forked from wrf-model/WRF
-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathmodule_convtrans_prep.F
204 lines (192 loc) · 7.66 KB
/
module_convtrans_prep.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
!------------------------------------------------------------------------
! original routinte by Georg Grell, adaptive timestepping by William Gustafson Jr. (PNNL), cloud fraction by Georg Grell (based on Stu's previous work wih so2so4 routine
!------------------------------------------------------------------------
Module module_convtrans_prep
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CONTAINS
subroutine convtrans_prep(gd_cloud,gd_cloud2,gd_cloud_a, &
gd_cloud_b,raincv,raincv_a,raincv_b, &
cldfr,moist,p_QV,p_QC,p_qi,T_PHY,P_PHY,num_moist, &
gd_cloud2_a,gd_cloud2_b,convtrans_avglen_m, &
adapt_step_flag,curr_secs, &
ktau,dt,cu_phys, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte,kts,kte )
REAL, PARAMETER :: coef_p = 0.25, coef_gamm = 0.49, coef_alph = 100.
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
p_QV,p_QC,p_qi,num_moist
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
OPTIONAL, &
INTENT(IN ) :: gd_cloud,gd_cloud2
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: t_phy,p_phy
REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
INTENT(IN ) :: moist
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
OPTIONAL, &
INTENT(INOUT ) :: gd_cloud_a,gd_cloud_b,gd_cloud2_a, &
cldfr,gd_cloud2_b
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: raincv
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(INOUT ) :: raincv_a,raincv_b
INTEGER, INTENT(IN) :: ktau,cu_phys
INTEGER :: stepave
INTEGER, SAVE :: stepave_count
REAL, INTENT(IN) :: curr_secs
REAL, INTENT(IN) :: convtrans_avglen_m, dt
LOGICAL, INTENT(IN) :: adapt_step_flag
LOGICAL :: avg_end_flag, first_period_flag
REAL :: satvp,rhgrid,h2oliq
real :: pmax,pmin
!
! Determine where we are in relation to the averaging period...
!
! convtrans_avglen_m = 30. !Averaging time for convective transport in min.
stepave=convtrans_avglen_m*60./dt
avg_end_flag = .false. !Initially assume we are not at the end
first_period_flag = .false. !Nor at the beginning
if( adapt_step_flag ) then
!If end of period...
if( curr_secs+real(dt,8)+0.01 >= &
( int( curr_secs/real(convtrans_avglen_m*60.,8) + 1_8, 8) &
*real(convtrans_avglen_m*60.,8) ) ) &
avg_end_flag = .true.
if( curr_secs <= real(convtrans_avglen_m*60.,8) ) first_period_flag = .true.
else
if( mod(ktau,stepave)==0 ) avg_end_flag = .true.
if( ktau <= stepave ) first_period_flag = .true.
end if
!
! Initialize the averaging arrays at the simulation start
!
if(ktau.le.1)then
stepave_count = 0
raincv_a(its:ite,jts:jte) = 0.
raincv_b(its:ite,jts:jte) = 0.
end if
if(present(gd_cloud2_a))then
if(ktau.le.1) gd_cloud2_a(its:ite,kts:kte,jts:jte)=0.
end if
if(present(cldfr))then
if(ktau.le.1) cldfr(its:ite,kts:kte,jts:jte)=0.
end if
!
! no time average available in first half hour
!
if( first_period_flag )then
do j=jts,jte
do i=its,ite
raincv_b(i,j)=raincv(i,j)
enddo
enddo
end if
!
! build time average, and stored in raincv_b to be used by convective transport routine
!
stepave_count = stepave_count+1
do j=jts,jte
do i=its,ite
raincv_a(i,j)=raincv_a(i,j)+raincv(i,j)
enddo
enddo
if( avg_end_flag )then
do j=jts,jte
do i=its,ite
raincv_b(i,j)=raincv_a(i,j)/real(stepave_count)
raincv_a(i,j)=0.
enddo
enddo
end if
!
! do the same for convective parameterization cloud water mix ratio,
! currently only for cu_physics=3,5, used by both photolysis and atmospheric radiation
!
if(cu_phys.eq.3.or.cu_phys.eq.5.or.cu_phys.eq.93)then
! if(config_flags%cu_physics == GDSCHEME .OR. &
! config_flags%cu_physics == GFSCHEME .OR. &
! config_flags%cu_physics == GFSCHEME ) THEN
! pmax=maxval(gd_cloud)
! pmin=maxval(gd_cloud2)
! print *,pmax,pmin
if( first_period_flag )then
do j=jts,jte
do k=kts,kte
do i=its,ite
gd_cloud_b(i,k,j)=gd_cloud(i,k,j)
gd_cloud2_b(i,k,j)=gd_cloud2(i,k,j)
enddo
enddo
enddo
end if ! stepave
!
!
!
do j=jts,jte
do k=kts,kte
do i=its,ite
gd_cloud_a(i,k,j)=gd_cloud_a(i,k,j)+gd_cloud(i,k,j)
gd_cloud2_a(i,k,j)=gd_cloud2_a(i,k,j)+gd_cloud2(i,k,j)
enddo
enddo
enddo
if( avg_end_flag )then
do j=jts,jte
do k=kts,kte
do i=its,ite
gd_cloud_b(i,k,j)=.1*gd_cloud_a(i,k,j)/real(stepave_count)
gd_cloud_a(i,k,j)=0.
gd_cloud2_b(i,k,j)=.1*gd_cloud2_a(i,k,j)/real(stepave_count)
gd_cloud2_a(i,k,j)=0.
enddo
enddo
enddo
! pmax=maxval(gd_cloud_b)
! pmin=maxval(gd_cloud2_b)
! print *,'avg_end_flag ',pmax,pmin
end if !stepave
end if ! cu_rad_feedback
!
! Clear the accumulator counter if we just finished an average...
!
if( avg_end_flag ) stepave_count = 0
! Siebesma et al., JAS, Vol. 60, no. 10, 1201-1219, 2003 (based on LES comparisons with trade-wind cumulus from BOMEX)
! SAM: Note units of liquid water and saturation vapor pressure must be in g/kg
! within the Siebesma et al. cloud fraction scheme
if( first_period_flag .or. avg_end_flag )then
do j=jts,jte
do k=kts,kte
do i=its,ite
cldfr(i,k,j)=0.
! if(gd_cloud_b(i,k,j).gt.0 .or. gd_cloud2_b(i,k,j).gt.0)then
if(p_qc.gt.1 .and. p_qi.le.1)then
satvp = 3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
(t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))
rhgrid = max(.1,MIN( .95, moist(i,k,j,p_qv) /satvp))
h2oliq=1000.*(gd_cloud_b(i,k,j) + moist(i,k,j,p_qc))
satvp=1000.*satvp
cldfr(i,k,j)=(1.-exp(-coef_alph*h2oliq/((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p)
cldfr(i,k,j)=max(0.,MIN(1.,cldfr(i,k,j)))
if(moist(i,k,j,p_qc).eq.0)cldfr(i,k,j)=cldfr(i,k,j)*.2
endif
if(p_qc.gt.1 .and. p_qi.gt.1)then
satvp = 3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
(t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))
rhgrid = max(.1,MIN( .95, moist(i,k,j,p_qv) /satvp))
h2oliq=1000.*(gd_cloud_b(i,k,j) + moist(i,k,j,p_qc) + &
gd_cloud2_b(i,k,j) + moist(i,k,j,p_qi))
satvp=1000.*satvp
cldfr(i,k,j)=(1.-exp(-coef_alph*h2oliq/((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p)
cldfr(i,k,j)=max(0.,MIN(1.,cldfr(i,k,j)))
if(moist(i,k,j,p_qc).eq.0 .and. moist(i,k,j,p_qi).eq.0)cldfr(i,k,j)=cldfr(i,k,j)*.2
endif
! endif
enddo
enddo
enddo
endif
END subroutine convtrans_prep
END MODULE MODULE_CONVTRANS_prep