module_mp_kessler.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3
4 MODULE module_mp_kessler
5
6 CONTAINS
7 !----------------------------------------------------------------
8 SUBROUTINE kessler( t, qv, qc, qr, rho, pii &
9 ,dt_in, z, xlv, cp &
10 ,EP2,SVP1,SVP2,SVP3,SVPT0,rhowater &
11 ,dz8w &
12 ,RAINNC, RAINNCV &
13 ,ids,ide, jds,jde, kds,kde & ! domain dims
14 ,ims,ime, jms,jme, kms,kme & ! memory dims
15 ,its,ite, jts,jte, kts,kte & ! tile dims
16 )
17 !----------------------------------------------------------------
18 IMPLICIT NONE
19 !----------------------------------------------------------------
20 ! taken from the COMMAS code - WCS 10 May 1999.
21 ! converted from FORTRAN 77 to 90, tiled, WCS 10 May 1999.
22 !----------------------------------------------------------------
23 REAL , PARAMETER :: c1 = .001
24 REAL , PARAMETER :: c2 = .001
25 REAL , PARAMETER :: c3 = 2.2
26 REAL , PARAMETER :: c4 = .875
27 REAL , PARAMETER :: fudge = 1.0
28 REAL , PARAMETER :: mxfall = 10.0
29 !----------------------------------------------------------------
30 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
31 ims,ime, jms,jme, kms,kme, &
32 its,ite, jts,jte, kts,kte
33 REAL , INTENT(IN ) :: xlv, cp
34 REAL , INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
35 REAL , INTENT(IN ) :: rhowater
36
37 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
38 INTENT(INOUT) :: &
39 t , &
40 qv, &
41 qc, &
42 qr
43
44 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
45 INTENT(IN ) :: &
46 rho, &
47 pii, &
48 dz8w
49
50 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
51 INTENT(IN ) :: z
52
53 REAL, INTENT(IN ) :: dt_in
54
55 REAL, DIMENSION( ims:ime , jms:jme ), &
56 INTENT(INOUT) :: RAINNC, &
57 RAINNCV
58
59
60
61 ! local variables
62
63 REAL :: qrprod, ern, gam, rcgs, rcgsi
64 REAL, DIMENSION( its:ite , kts:kte, jts:jte ) :: prod
65 REAL, DIMENSION(kts:kte) :: vt, prodk, vtden,rdzk,rhok,factor,rdzw
66 INTEGER :: i,j,k
67 INTEGER :: nfall, n, nfall_new
68 REAL :: qrr, pressure, temp, es, qvs, dz, dt
69 REAL :: f5, dtfall, rdz, product
70 REAL :: max_heating, max_condense, max_rain, maxqrp
71 REAL :: vtmax, ernmax, crmax, factorn, time_sediment
72 REAL :: qcr, factorr, ppt
73 REAL, PARAMETER :: max_cr_sedimentation = 0.75
74 !----------------------------------------------------------------
75
76 INTEGER :: imax, kmax
77
78 dt = dt_in
79
80 ! f5 = 237.3 * 17.27 * 2.5e6 / cp
81 f5 = svp2*(svpt0-svp3)*xlv/cp
82 ernmax = 0.
83 maxqrp = -100.
84
85 !------------------------------------------------------------------------------
86 ! parameters for the time split terminal advection
87 !------------------------------------------------------------------------------
88
89 max_heating = 0.
90 max_condense = 0.
91 max_rain = 0.
92
93 !-----------------------------------------------------------------------------
94 ! outer J loop for entire microphysics, outer i loop for sedimentation
95 !-----------------------------------------------------------------------------
96
97 microphysics_outer_j_loop: DO j = jts, jte
98
99 sedimentation_outer_i_loop: DO i = its,ite
100
101 ! vtmax = 0.
102 crmax = 0.
103
104
105 !------------------------------------------------------------------------------
106 ! Terminal velocity calculation and advection, set up coefficients and
107 ! compute stable timestep
108 !------------------------------------------------------------------------------
109
110 DO k = 1, kte
111 prodk(k) = qr(i,k,j)
112 rhok(k) = rho(i,k,j)
113 qrr = amax1(0.,qr(i,k,j)*0.001*rhok(k))
114 vtden(k) = sqrt(rhok(1)/rhok(k))
115 vt(k) = 36.34*(qrr**0.1364) * vtden(k)
116 ! vtmax = amax1(vt(k), vtmax)
117 rdzw(k) = 1./dz8w(i,k,j)
118 crmax = amax1(vt(k)*dt*rdzw(k),crmax)
119 ENDDO
120 DO k = 1, kte-1
121 rdzk(k) = 1./(z(i,k+1,j) - z(i,k,j))
122 ENDDO
123 rdzk(kte) = 1./(z(i,kte,j) - z(i,kte-1,j))
124
125 nfall = max(1,nint(0.5+crmax/max_cr_sedimentation)) ! courant number for big timestep.
126 dtfall = dt / float(nfall) ! splitting so courant number for sedimentation
127 time_sediment = dt ! is stable
128
129 !------------------------------------------------------------------------------
130 ! Terminal velocity calculation and advection
131 ! Do a time split loop on this for stability.
132 !------------------------------------------------------------------------------
133
134 column_sedimentation: DO WHILE ( nfall > 0 )
135
136 time_sediment = time_sediment - dtfall
137 DO k = 1, kte-1
138 factor(k) = dtfall*rdzk(k)/rhok(k)
139 ENDDO
140 factor(kte) = dtfall*rdzk(kte)
141
142 ppt=0.
143
144 k = 1
145 ppt=rhok(k)*prodk(k)*vt(k)*dtfall/rhowater
146 RAINNCV(i,j)=ppt*1000.
147 RAINNC(i,j)=RAINNC(i,j)+ppt*1000. ! unit = mm
148
149 !------------------------------------------------------------------------------
150 ! Time split loop, Fallout done with flux upstream
151 !------------------------------------------------------------------------------
152
153 DO k = kts, kte-1
154 prodk(k) = prodk(k) - factor(k) &
155 * (rhok(k)*prodk(k)*vt(k) &
156 -rhok(k+1)*prodk(k+1)*vt(k+1))
157 ENDDO
158
159 k = kte
160 prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
161
162 !------------------------------------------------------------------------------
163 ! compute new sedimentation velocity, and check/recompute new
164 ! sedimentation timestep if this isn't the last split step.
165 !------------------------------------------------------------------------------
166
167 IF( nfall > 1 ) THEN ! this wasn't the last split sedimentation timestep
168
169 nfall = nfall - 1
170 crmax = 0.
171 DO k = kts, kte
172 qrr = amax1(0.,prodk(k)*0.001*rhok(k))
173 vt(k) = 36.34*(qrr**0.1364) * vtden(k)
174 ! vtmax = amax1(vt(k), vtmax)
175 crmax = amax1(vt(k)*time_sediment*rdzw(k),crmax)
176 ENDDO
177
178 nfall_new = max(1,nint(0.5+crmax/max_cr_sedimentation))
179 if (nfall_new /= nfall ) then
180 nfall = nfall_new
181 dtfall = time_sediment/nfall
182 end if
183
184 ELSE ! this was the last timestep
185
186 DO k=kts,kte
187 prod(i,k,j) = prodk(k)
188 ENDDO
189 nfall = 0 ! exit condition for sedimentation loop
190
191 END IF
192
193 ENDDO column_sedimentation
194
195 ENDDO sedimentation_outer_i_loop
196
197 !------------------------------------------------------------------------------
198 ! Production of rain and deletion of qc
199 ! Production of qc from supersaturation
200 ! Evaporation of QR
201 !------------------------------------------------------------------------------
202
203 DO k = kts, kte
204 DO i = its, ite
205 factorn = 1.0 / (1.+c3*dt*amax1(0.,qr(i,k,j))**c4)
206 qrprod = qc(i,k,j) * (1.0 - factorn) &
207 + factorn*c1*dt*amax1(qc(i,k,j)-c2,0.)
208 rcgs = 0.001*rho(i,k,j)
209
210 qc(i,k,j) = amax1(qc(i,k,j) - qrprod,0.)
211 qr(i,k,j) = (qr(i,k,j) + prod(i,k,j)-qr(i,k,j))
212 qr(i,k,j) = amax1(qr(i,k,j) + qrprod,0.)
213
214 temp = pii(i,k,j)*t(i,k,j)
215 pressure = 1.000e+05 * (pii(i,k,j)**(1004./287.))
216 gam = 2.5e+06/(1004.*pii(i,k,j))
217 ! qvs = 380.*exp(17.27*(temp-273.)/(temp- 36.))/pressure
218 es = 1000.*svp1*exp(svp2*(temp-svpt0)/(temp-svp3))
219 qvs = ep2*es/(pressure-es)
220 ! prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+qvs*f5/(temp-36.)**2)
221 prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+pressure/(pressure-es)*qvs*f5/(temp-svp3)**2)
222 ern = amin1(dt*(((1.6+124.9*(rcgs*qr(i,k,j))**.2046) &
223 *(rcgs*qr(i,k,j))**.525)/(2.55e8/(pressure*qvs) &
224 +5.4e5))*(dim(qvs,qv(i,k,j))/(rcgs*qvs)), &
225 amax1(-prod(i,k,j)-qc(i,k,j),0.),qr(i,k,j))
226
227 ! Update all variables
228
229 product = amax1(prod(i,k,j),-qc(i,k,j))
230 t (i,k,j) = t(i,k,j) + gam*(product - ern)
231 qv(i,k,j) = amax1(qv(i,k,j) - product + ern,0.)
232 qc(i,k,j) = qc(i,k,j) + product
233 qr(i,k,j) = qr(i,k,j) - ern
234
235 ENDDO
236 ENDDO
237
238 ENDDO microphysics_outer_j_loop
239
240 RETURN
241
242 END SUBROUTINE kessler
243
244 END MODULE module_mp_kessler