| blob | 9300dcc30673ef4d9dd02f7263eba3bc37b06552 |
1 !#define DENG_SHCU1D
2 !#define DENG_SHCU1D_subsidence
3 MODULE module_shcu_deng
5 USE module_wrf_error
7 REAL , PARAMETER :: TO = 263.15
8 ! Lookup table variables:
9 INTEGER, PARAMETER :: KFNT=250,KFNP=220
10 REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
11 REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
12 REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
13 REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
15 CONTAINS
17 SUBROUTINE deng_shcu_driver( &
18 ids,ide, jds,jde, kds,kde &
19 ,ims,ime, jms,jme, kms,kme &
20 ,its,ite, jts,jte, kts,kte &
21 ,DT,KTAU,DX, xtime, gmt &
22 ,XLV, XLS,XLV0,XLV1,XLS0,XLS1,CP,R,G & ! Constants
23 ,SVP1,SVP2,SVP3,SVPT0 & ! Constants
24 ,ADAPT_STEP_FLAG, DSIGMA & ! in
25 ,XLONG, ht, PBLH & ! in 2D
26 ,U,V,w,TH,T,QV,QC,QR,dz8w,PCPS,rho,z_at_w,pi & ! in 3D
27 ,tke, kth, bbls &
28 ,ten_radl, ten_rads & ! in 3D
29 ,rainsh, rainshv, rainshvb, pblmax, capesave & ! inout 2D
30 ,xtime1, radsave, clddpthb, cldtopb, MAVAIL, PBLHAVG & ! inout 2D
31 ,ainckfsa, ltopb, kdcldtop, kdcldbas & ! inout 2D
32 ,W0AVG, TKEAVG, cldareaa, cldareab & ! inout 3D
33 ,cldliqa, cldliqb, cldfra_sh, ca_rad, cw_rad, wub & ! inout 3D
34 ,RUSHTEN, RVSHTEN, RTHSHTEN, RQVSHTEN, RQCSHTEN & ! inout 3D
35 ,RQRSHTEN, RDCASHTEN, RQCDCSHTEN & ! inout 3D
36 )
37 !-------------------------------------------------------------
38 IMPLICIT NONE
39 ! ------------ IN ---------------------------------------------
40 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
41 ims,ime, jms,jme, kms,kme, &
42 its,ite, jts,jte, kts,kte
43 REAL, INTENT(IN ) :: DT, DX, xtime, gmt
44 INTEGER, INTENT(IN ) :: KTAU
45 REAL, INTENT(IN ) :: XLV, XLS,XLV0,XLV1,XLS0,XLS1,CP,R,G, &
46 SVP1,SVP2,SVP3,SVPT0
47 LOGICAL , INTENT(IN ) :: adapt_step_flag
48 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: dsigma
50 REAL, DIMENSION( ims:ime , jms:jme ), &
51 INTENT(IN ) :: xlong, ht, pblh
52 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
53 INTENT(IN ) :: U, V, W, TH, T, QV, QR, dz8w, Pcps, &
54 rho, z_at_w, pi, tke, &
55 kth, bbls, ten_radl, ten_rads
56 ! ------------ INOUT ------------------------------------------
57 REAL, DIMENSION( ims:ime , jms:jme ), &
58 INTENT(INOUT) :: RAINSHV, RAINSH, pblmax, cldtopb, clddpthb, &
59 rainshvb, capesave, xtime1, radsave, &
60 MAVAIL, PBLHAVG
62 REAL, DIMENSION( ims:ime , 1:100, jms:jme ), &
63 INTENT(INOUT) :: ainckfsa
64 INTEGER, DIMENSION( ims:ime , jms:jme ), &
65 INTENT(INOUT) :: ltopb, kdcldtop, kdcldbas
66 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
67 INTENT(INOUT) :: W0AVG, TKEAVG, &
68 cldareaa, cldareab, cldliqa, cldliqb, wub
69 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
70 INTENT(IN ) :: QC
71 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
72 INTENT( OUT) :: ca_rad, cw_rad, cldfra_sh
73 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
74 INTENT(INOUT) :: RUSHTEN, RVSHTEN, RTHSHTEN, RQVSHTEN, RQCSHTEN, &
75 RQRSHTEN, RDCASHTEN, RQCDCSHTEN
76 !
77 ! Local Variables
78 !
79 REAL, DIMENSION( its:ite , jts:jte ) :: CAPEI, RADIUSC, AINCKFI
80 REAL, DIMENSION( kts:kte ) :: U1D, V1D, TH1D, T1D, DZ1D, QV1D, &
81 QC1D, RH1D, QR1D, P1D, z1d, RHO1D, W0AVG1D, &
82 kth1D, tke1d, bbls1D, &
83 #ifdef DENG_SHCU1D
84 ca_rad1d, &
85 #endif
86 ten_radl1d, ten_rads1d
87 REAL, DIMENSION( kts:kte+1 ) :: z_at_w1d
88 REAL, DIMENSION( kts:kte ) :: QVTEN, QCTEN, QRTEN, TTEN, DCATEN, QCDCTEN
89 REAL, DIMENSION( its:ite , jts:jte ) :: pblh_tile, pblhavg_tile
90 REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: w_at_half, w_at_half_mean
91 REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: tke_tile, tkeavg_tile
92 REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: cldareac, cldliqc, rh
94 REAL :: dt2, DXSQ
95 INTEGER :: i, j, k, ii, ntst, kk
96 REAL :: AMOIS, time_period_avg_min ! time window to perform a running mean value
97 INTEGER :: if_avg_w=1, if_avg_pbl=0, if_avg_tke=0
98 REAL :: GNUHF, OMUHF
99 REAL :: es, qes, alf1, alf2, alf3, c1, c2, cs, qctot
100 #ifdef DENG_SHCU1D_subsidence
101 REAL, DIMENSION(100) :: P_P, WSUBS_P
102 REAL, DIMENSION(kte) :: P_SIGMA, WSUBSID, WSUBSID_td, thten_sub
103 REAL, DIMENSION(kte+1) :: th_full
104 #endif
105 INTEGER :: kx, kxp1, kl
107 dxsq = dx*dx
108 dt2 = 2.0*dt
109 GNUHF = 0.1
110 OMUHF = 1.-2.*GNUHF
111 kx = kte
112 kl = kte
113 kxp1 = kx + 1
115 ! print*, ids, ide, jds, jde, kds, kde
116 ! print*, ims, ime, jms, jme, kms, kme
117 ! print*, its, ite, jts, jte, kts, kte
119 #ifdef DENG_SHCU1D
120 write(301,*) xtime/60.,pblh(1,1)
121 #ifdef DENG_SHCU1D_subsidence
122 OPEN(3,FILE='input03.dat',STATUS='old')
123 READ(3,*)
124 DO K=1,35
125 READ(3,*) P_P(K),WSUBS_P(K) ! SUBSIDENCE IN CM/SEC. (from Betts 1994)
126 ENDDO
127 CLOSE(3)
128 #endif
129 #endif
130 !
131 ! Time averaging w on the half levels
132 !
133 IF( if_avg_w == 1 ) THEN
134 time_period_avg_min = 10.0
135 DO J = jts,jte
136 DO I = its,ite
137 DO k = kts,kte
138 w_at_half(i,k,j) = 0.5* ( w(i,k,j) + w(i,k+1,j) )
139 w_at_half_mean(i,k,j) = w0avg(i,k,j)
140 ENDDO
141 ENDDO
142 ENDDO
143 CALL time_avg_3d( dt, ADAPT_STEP_FLAG, time_period_avg_min &
144 ,w_at_half, w_at_half_mean &
145 ,its,ite, jts,jte, kts,kte &
146 )
147 DO J = jts,jte
148 DO I = its,ite
149 DO k = kts,kte
150 w0avg(i,k,j) = w_at_half_mean(i,k,j)
151 ENDDO
152 ENDDO
153 ENDDO
154 ELSE
155 DO J = jts,jte
156 DO I = its,ite
157 DO k = kts,kte
158 w_at_half(i,k,j) = 0.5* ( w(i,k,j) + w(i,k+1,j) )
159 w0avg(i,k,j) = w_at_half(i,k,j)
160 ENDDO
161 ENDDO
162 ENDDO
163 ENDIF
164 !
165 ! Time averaging PBLH
166 !
167 IF( if_avg_pbl == 1 ) THEN
168 time_period_avg_min = 120.0
169 DO J = jts,jte
170 DO I = its,ite
171 pblh_tile(i,j) = pblh(i,j)
172 pblhavg_tile(i,j) = pblhavg(i,j)
173 ENDDO
174 ENDDO
175 CALL time_avg_2d( dt, ADAPT_STEP_FLAG, time_period_avg_min &
176 ,pblh_tile, pblhavg_tile &
177 ,its,ite, jts,jte &
178 )
179 DO J = jts,jte
180 DO I = its,ite
181 pblhavg(i,j) = pblhavg_tile(i,j)
182 ENDDO
183 ENDDO
184 ELSE
185 DO J = jts,jte
186 DO I = its,ite
187 pblh_tile(i,j) = pblh(i,j)
188 pblhavg(i,j) = pblh_tile(i,j)
189 ENDDO
190 ENDDO
191 ENDIF
192 #ifdef DENG_SHCU1D
193 write(302,*) xtime/60.,PBLHAVG(1,1)
194 #endif
195 !
196 ! Time averaging TKE
197 !
198 IF( if_avg_tke == 1 ) THEN
199 time_period_avg_min = 30.0
200 DO J = jts,jte
201 DO I = its,ite
202 DO k = kts,kte
203 tke_tile(i,k,j) = tke(i,k,j)
204 tkeavg_tile(i,k,j) = tkeavg(i,k,j)
205 ENDDO
206 ENDDO
207 ENDDO
208 CALL time_avg_3d( dt, ADAPT_STEP_FLAG, time_period_avg_min &
209 ,tke_tile, tkeavg_tile &
210 ,its,ite, jts,jte, kts,kte &
211 )
212 DO J = jts,jte
213 DO I = its,ite
214 DO k = kts,kte
215 tkeavg(i,k,j) = tkeavg_tile(i,k,j)
216 ENDDO
217 ENDDO
218 ENDDO
219 ELSE
220 DO J = jts,jte
221 DO I = its,ite
222 DO k = kts,kte
223 tke_tile(i,k,j) = tke(i,k,j)
224 tkeavg(i,k,j) = tke_tile(i,k,j)
225 ENDDO
226 ENDDO
227 ENDDO
228 ENDIF
229 !
230 ! Calculate the RH-based cloud fraction, CS, based on Xu and Randall (1996, JAS)
231 ! As shown in Deng et al. 2003a, the virtual cloud fraction CAREA is a wighted average
232 ! of Deng scheme-predicted cloud fraction, CLDLIQA and CS. This effective cloud
233 ! fraction is used in the radiation calculation.
234 !
235 DO J = jts,jte
236 DO I = its,ite
237 DO k = kts,kte
238 IF( T(i,k,j) .GE. TO ) THEN
239 ES=1.E3*SVP1*EXP(SVP2*(T(i,k,j)-SVPT0)/(T(i,k,j)-SVP3))
240 ELSE
241 ES=611.0*EXP(22.514-6.15E3/T(i,k,j))
242 ENDIF
243 QES = 0.622*ES/(PCPS(i,k,j)-ES)
244 RH(i,k,j)=AMAX1(QV(i,k,j)/QES,1.0E-25)
246 alf1 = 0.25
247 alf2 = 0.49
248 alf3 = 100.0
249 IF( RH(I,K,J) .LT. 0.999) THEN
250 C1=(1.0-RH(I,K,J))*QES**alf2
251 QCTOT=CLDLIQA(I,K,J)*CLDAREAA(I,K,J)+QC(I,K,J)
252 C2=-alf3*QCTOT/C1
253 CS=RH(I,K,J)**alf1*(1.0-EXP(C2))
254 ELSE
255 CS=1.0
256 ENDIF
258 CA_RAD(I,K,J)=(1.0-CLDAREAA(I,K,J))*CS+CLDAREAA(I,K,J)
259 CLDFRA_SH(I,K,J)= CA_RAD(I,K,J)
260 CW_RAD(I,K,J)=CLDLIQA(I,K,J)*CLDAREAA(I,K,J)
261 ENDDO
262 ENDDO
263 ENDDO
265 DO J = jts,jte
266 DO I = its,ite
267 DO k=kts,kte
268 QVTEN(k)=0.
269 QCTEN(k)=0.
270 QRTEN(k)=0.
271 TTEN(k)=0.
272 DCATEN(k)=0.
273 QCDCTEN(k)=0.
274 ENDDO
275 RAINSHV(I,J)=0.
276 RAINSH(I,J)=0.
278 DO K=kts,kte
279 U1D(K) =U(I,K,J)
280 V1D(K) =V(I,K,J)
281 TH1D(K) =TH(I,K,J)
282 T1D(K) =T(I,K,J)
283 W0AVG1D(k) = W0AVG(i,k,j)
284 tke1d(k) = tkeavg(i,k,j)
285 RHO1D(K) =rho(I,K,J)
286 QV1D(K)=MAX(QV(I,K,J),1.E-12)
287 RH1D(K)=RH(I,K,J)
288 #ifdef DENG_SHCU1D
289 ca_rad1d(k) = ca_rad(i,k,j)
290 #endif
291 QC1D(K)=QC(I,K,J)
292 QR1D(K)=QR(I,K,J)
293 P1D(K) =Pcps(I,K,J)
294 DZ1D(k)=dz8w(I,K,J)
295 kth1d(k) = kth(i,k,j)
296 bbls1d(k) = bbls(i,k,j)
297 ten_radl1d(k) = ten_radl(i,k,j)
298 ten_rads1d(k) = ten_rads(i,k,j)
299 ENDDO
301 DO K=kts,kte+1
302 z_at_w1d(k) = z_at_w(I,K,J)
303 ENDDO
305 #ifdef DENG_SHCU1D_subsidence
306 DO K=1,KL
307 kk = KL - K + 1
308 P_SIGMA(K) = p1d(kk)/100.0
309 ENDDO
310 CALL INTERP1D( P_SIGMA, P_P, WSUBS_P, WSUBSID_td, KX, 35 )
311 DO K = 1,KX
312 kk = KX - K + 1
313 WSUBSID(K) = WSUBSID_td(KK)/100.0 ! convert to m/s (READ IN INPUT2.F)
314 ENDDO
315 DO K = 2,KX
316 th_full(K) = 0.5*( th1d(k)+th1d(k-1) )
317 ENDDO
318 DO K = 2,KX-1
319 thten_sub(K) = - WSUBSID(K) * ( th_full(K+1)-th_full(K) ) / dz1d(k)
320 ENDDO
321 thten_sub(1) = 0.0
322 thten_sub(kx) = 0.0
323 #endif
324 CALL deng_shcu(I, J, &
325 ids,ide, jds,jde, kds,kde, &
326 ims,ime, jms,jme, kms,kme, &
327 its,ite, jts,jte, kts,kte, &
328 XLV,XLS,XLV0,XLV1,XLS0,XLS1,CP,R,G, &
329 SVP1,SVP2,SVP3,SVPT0, &
330 DT, ktau, dt2, DX, DXSQ, xtime, gmt, &
331 PBLHAVG(i,j), ht(i,j), xlong, capesave, radsave, ainckfsa, &
332 U1D,V1D,T1D,QV1D,rh1d,QC1D,QR1D,p1d,RHO1D,z_at_w1d,DZ1D,W0AVG1D, &
333 #ifdef DENG_SHCU1D
334 ca_rad1d, & ! for 1d printing only
335 #endif
336 tke1d, kth1d, bbls1d, ten_radl1d, ten_rads1d, dsigma, &
337 RAINSHV, RAINSH, pblmax, CLDTOPB, CLDDPTHB, &
338 ltopb, kdcldtop, kdcldbas, &
339 cldareab, cldliqb, wub, &
340 CAPEI, RADIUSC, AINCKFI, &
341 QVTEN,QCTEN,QRTEN,TTEN,DCATEN,QCDCTEN &
342 )
344 DO K=kts,kte
345 RUSHTEN(i,k,j) = 0.0
346 RVSHTEN(i,k,j) = 0.0
347 RTHSHTEN(i,k,j) = TTEN(k)/pi(I,K,J) * 1.0
348 #ifdef DENG_SHCU1D_subsidence
349 RTHSHTEN(i,k,j) = RTHSHTEN(i,k,j) + thten_sub(k)
350 #endif
351 RQVSHTEN(i,k,j) = QVTEN(k) * 1.0
352 RQCSHTEN(i,k,j) = QCTEN(k) * 1.0
353 RQRSHTEN(i,k,j) = QRTEN(k) * 1.0
354 RDCASHTEN(i,k,j) = DCATEN(K) * 1.0
355 RQCDCSHTEN(i,k,j) = QCDCTEN(K) * 1.0
356 ENDDO
357 !
358 ! IF THERE IS CONVECTIVE PRECIPITATION, UPDATE THE SURFACE MOISTURE
359 ! AVAILABILITY THAT THE SURFACE IS WET FOR THE FOLLOWING HOUR
360 !
361 IF(RAINSHV(I,J) .EQ. 0.0 .AND. RAINSHVB(I,J) .NE. 0.0) THEN
362 XTIME1(I,J)=XTIME+60.
363 ENDIF
365 AMOIS = MAVAIL(I,J)
366 IF(RAINSHV(I,J) .NE. 0.0 .OR. XTIME .LE. XTIME1(I,J)) THEN
367 MAVAIL(I,J)=AMOIS+0.67*(1.0-AMOIS)
368 ELSE
369 MAVAIL(I,J)=AMOIS
370 ENDIF
372 RAINSHVB(I,J)=RAINSHV(I,J)
373 DO II = 100,2,-1
374 AINCKFSA(i,ii,j) = AINCKFSA(i,ii-1,j)
375 ENDDO
376 RADSAVE(I,J)=RADIUSC(i,j)
377 AINCKFSA(i,1,j) = AINCKFI(i,j)
378 CAPESAVE(i,j) = CAPEI(i,j)
380 ENDDO ! i-loop
381 ENDDO ! j-loop
382 !
383 ! NBC prediction using the MM5 leapfog method
384 !
385 DO J = jts,jte
386 DO I = its,ite
387 DO K = kts,kte
388 CLDAREAC(I,K,J) = CLDAREAB(I,K,J) + RDCASHTEN(I,K,J) * DT
389 CLDLIQC(I,K,J) = CLDLIQB(I,K,J) + RQCDCSHTEN(I,K,J) * DT
390 CLDAREAC(I,K,J) = AMAX1(0.0,CLDAREAC(I,K,J))
391 CLDLIQC(I,K,J) = AMAX1(0.0,CLDLIQC(I,K,J))
392 IF( CLDAREAC(I,K,J) .GT. 0.99999 ) THEN
393 CLDLIQC(I,K,J) = CLDLIQC(I,K,J)*CLDAREAC(I,K,J)
394 CLDAREAC(I,K,J) = 1.0
395 ENDIF
396 IF( CLDAREAC(I,K,J) .LE. 1.0e-17 .OR. CLDLIQC(I,K,J) .LE. 1.0e-17 ) THEN
397 CLDAREAC(I,K,J) = 0.0
398 CLDLIQC(I,K,J) = 0.0
399 ENDIF
400 ENDDO
401 ENDDO
402 ENDDO
404 DO J = jts,jte
405 DO I = its,ite
406 DO K = kts,kte
407 CLDAREAB(I,K,J)=OMUHF*CLDAREAA(I,K,J)+GNUHF*(CLDAREAB(I,K,J)+CLDAREAC(I,K,J))
408 CLDLIQB(I,K,J)=OMUHF*CLDLIQA(I,K,J) +GNUHF*(CLDLIQB(I,K,J)+CLDLIQC(I,K,J))
409 CLDAREAA(I,K,J)=CLDAREAC(I,K,J)
410 CLDLIQA(I,K,J)= CLDLIQC(I,K,J)
412 IF( CLDAREAA(I,K,J) .LE. 1.0e-17 .OR. CLDLIQA(I,K,J) .LE. 1.0e-17 ) THEN
413 CLDAREAA(I,K,J) = 0.0
414 CLDLIQA(I,K,J) = 0.0
415 ENDIF
417 IF(CLDDPTHB(I,J) .LE. 0.0) THEN ! No active updraft
418 RQCSHTEN(I,K,J)=RQCSHTEN(I,K,J)+CLDAREAB(I,K,J)*CLDLIQB(I,K,J)/DT
419 CLDAREAB(I,K,J) = 0.0
420 CLDLIQB(I,K,J) = 0.0
421 ENDIF
422 ENDDO
423 ENDDO
424 ENDDO
426 END SUBROUTINE deng_shcu_driver
428 !================
429 SUBROUTINE deng_shcu(I,J, & ! in
430 ids,ide, jds,jde, kds,kde, & ! in
431 ims,ime, jms,jme, kms,kme, & ! in
432 its,ite, jts,jte, kts,kte, & ! in
433 XLV,XLS,XLV0,XLV1,XLS0,XLS1,CP,R,G, & ! in
434 SVP1,SVP2,SVP3,SVPT0, & ! in
435 DT, ktau, dt2, DX, DXSQ, xtime, gmt, zpbl, ht, & ! in
436 xlong, capesave, radsave, ainckfsa, & ! in 2D
437 U0,V0,T0,Q0,RH0, QC0_expl,QR0,P0, & ! in 3D
438 RHOE,z_at_w0, DZQ,W0AVG0, & ! in 1D
439 #ifdef DENG_SHCU1D
440 ca_rad0, & ! for 1d printing only
441 #endif
442 tke0, kth0, bbls0, ten_radl0, ten_rads0, dsigma, & ! in 1D
443 RAINSHV, RAINSH, pblmax, CLDTOPB, CLDDPTHB, & ! inout 2D
444 ltopb, kdcldtop, kdcldbas, & ! inout 2D
445 cldareab, cldliqb, wub, & ! inout 3D
446 CAPEI,RADIUSC, AINCKFI, & ! out 2D
447 QVTEN,QCTEN,QRTEN,TTEN,DCATEN,QCDCTEN & ! out 1D (output tendencie)
448 )
449 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450 IMPLICIT NONE
451 ! ----------- IN -------------------------------------------
452 INTEGER, INTENT(IN ) :: i, j, &
453 ids,ide, jds,jde, kds,kde, &
454 ims,ime, jms,jme, kms,kme, &
455 its,ite, jts,jte, kts,kte
456 REAL, INTENT(IN ) :: XLV,XLS,XLV0,XLV1,XLS0,XLS1, &
457 CP,R,G,SVP1,SVP2,SVP3,SVPT0
458 REAL, INTENT(IN ) :: DT, dt2, DX, DXSQ, xtime, gmt
459 INTEGER, INTENT(IN ) :: KTAU
460 REAL, INTENT(IN ) :: zpbl, ht
461 REAL, DIMENSION( ims:ime , jms:jme ), &
462 INTENT(IN ) :: xlong, capesave, radsave
463 REAL, DIMENSION( ims:ime , 1:100, jms:jme ), &
464 INTENT(IN ) :: ainckfsa
465 REAL, DIMENSION( kts:kte ), &
466 INTENT(IN ) :: U0, V0, T0, QR0, RH0, P0, rhoe, DZQ, W0AVG0, &
467 TKE0, KTH0, BBLS0, ten_radl0, ten_rads0
468 REAL, DIMENSION( kms:kme ), &
469 INTENT(IN ) :: dsigma
470 REAL, DIMENSION( kts:kte+1 ), &
471 INTENT(IN ) :: z_at_w0
472 ! ------------ INOUT ---------------------------------------------
473 REAL, DIMENSION( kts:kte ), &
474 INTENT(INOUT) :: Q0, QC0_expl
475 REAL, DIMENSION( ims:ime, jms:jme ), &
476 INTENT(INOUT) :: RAINSHV, RAINSH, pblmax, cldtopb, clddpthb
477 INTEGER, DIMENSION( ims:ime , jms:jme ), &
478 INTENT(INOUT) :: ltopb, kdcldtop, kdcldbas
479 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
480 INTENT(INOUT) :: cldareab, cldliqb, wub
481 REAL, DIMENSION( kts:kte ), &
482 INTENT(INOUT) :: QVTEN, QCTEN, QRTEN, TTEN, DCATEN, QCDCTEN
483 ! ----------- OUT ------------------------------------------------
484 REAL, DIMENSION( its:ite, jts:jte ), &
485 INTENT( OUT) :: CAPEI, RADIUSC, AINCKFI
486 !
487 ! Local variables
488 !
490 REAL, DIMENSION( kts:kte ) :: DTDT, DQDT, DQRDT, DQLDT, DDCADT, DQCDCDT
492 REAL , PARAMETER :: PIE = 3.141592654
493 REAL , PARAMETER :: TTFRZ = 268.16
494 REAL , PARAMETER :: TBFRZ = 248.16
495 REAL , PARAMETER :: RHBC = 0.90
496 REAL , PARAMETER :: P00 = 100000.0
497 REAL , PARAMETER :: T00 = 273.16
498 REAL , PARAMETER :: RLF = 3.339E5
499 REAL , PARAMETER :: AVTS = 11.72
500 REAL , PARAMETER :: BVTS = 0.41
501 REAL , PARAMETER :: N0R = 8.E6
502 ! REAL , PARAMETER :: TO = 263.15
503 REAL , PARAMETER :: XN0 = 1.E-2
504 REAL , PARAMETER :: XMMAX = 9.4E-10 ! from MM5 param.F for autoconvertion of qc to qs
505 REAL , PARAMETER :: N0S = 2.E7
506 REAL , PARAMETER :: ESI = 0.1
507 REAL , PARAMETER :: QCTH = 0.5E-3 ! from MM5 param.F for autoconvertion of qc to qr
508 REAL , PARAMETER :: QCK1 = 1.0E-3 ! from MM5 param.F for autoconvertion of qc to qr
509 REAL , PARAMETER :: AVT = 841.99667
510 REAL , PARAMETER :: BVT = 0.8
511 REAL , PARAMETER :: RHCRIT = 0.999
513 REAL, DIMENSION( kts:kte ) :: SCR3
515 REAL :: nupdraft, ZLCL,ZLFC, TRPPT, CLDTOP, CLDDPTH, & ! nupdraft is for sure a real number
516 TTLCL, PPLCL, ainc, ainc2, ainc3, ainc4, tenv, qenv, tven, tvbar, &
517 AINCM1, AINCM2, AINCMX, timeunv, timeloc, wbar, wtke, wpbltp, wnh
519 INTEGER :: LLC, KPBL, KKPBL,KLCL,KDCB,KDCT,LTOP, LTOP1, LTOPM1, LTOP2, &
520 kpblmax, kl, klm, kx, kxp1, kp1, LCL, LET, LFS, L, n
522 REAL, DIMENSION( kts:kte ) :: z0, prc, pra, tlong,tshort, &
523 QU,TU,TVU,UMF,UER,TZ,QD,TV0, &
524 WD,TVD,DMF,DER,TG,QQG,TVG, &
525 EMS,EMSD,THETEU,THETED,THETEE,THTAU, &
526 THTAD,THTA0,RLIQ,RICE,QLQOUT,QICOUT, &
527 PPTLIQ,PPTICE,DETLQ,DETIC,UMF2,DMF2, &
528 DETLQ2,DETIC2,UDR,UDR2,DDR,UER2, &
529 DER2,RATIO2,DOMGDP,EXN, &
530 DZA,TVQU,DP,DMS,EQFRC,WSPD,QDT,FXM, &
531 THTES,THTESG,THTAG,DDR2,THPA,QPA, &
532 THFXIN,THFXOUT,QFXIN,QFXOUT,DTFM, &
533 QC0,QCG,QCPA,QCFXIN,QCFXOUT, &
534 DCA,DCAFXIN,DCAFXOUT,DCA0, &
535 RGVC,FALOUTC,wu,qes
537 REAL, DIMENSION( kts:kte+1 ) :: OMG, KV,KR, LFLUX, LFLUX1, LFLUX2
539 REAL :: zagl_bot, zagl_top, &
540 tmix, thmix, qmix, zmix, pmix, dpthmx, &
541 qmix1, zmix1, dpthmx1, &
542 tmix2, thmix2, qmix2, pmix2, emix, emix2, &
543 c1, c2, c3, AREA, AREA1, AREA2, CLQ, CLQ1, CLQ2, &
544 emax, zval, TLCL, PLCL, TVLCL, ES, QS, TVAVG, QESE, WTW, RHOLCL, &
545 sum, val, val0, val1, val2, ff, h, b, rad, &
546 PPBL, TENVPBL, QENVPBL, TVENPBL, TPBL, TVPBL, RHOPBL, ROCPQ, &
547 RHO_INI, T_INI, TV_INI, TVEN_INI, Z_INI, P_INI, &
548 DPTH, PPLSDP, AU0, VMFPBLCL, VMFLCL, UPOLD, UPNEW, &
549 abe, THTUDL, TUDL, TTEMP, RATE, FRC, FRC1, QNEWLQ, QNEWIC, RL, &
550 QNWFRZ, EFFQ, BE, BOTERM, ENTERM, DZZ, WSQ, &
551 WABS, UDLBE, ee, ee1, ee2, ud1, ud2, rei, ttmp, f1, f2, THTTMP, &
552 RTMP, TMPLIQ, TMPICE, TU95, TU10, EQMAX, EQMIN, CTP1, CTP2, WUMAX, &
553 DPTT, DUMFDP, TSAT, THTA, P165, PPTMLT, PPTML2, CLVF, USR, VCONV, &
554 TIMEC, TADVEC, TCMIN, TCMAX, SHSIGN, VWS, PEF, PEFMX, PEFMN, &
555 CBH, RCBH, PEFCBH, CBHMX, PEFF, TDER, TDER2, THTMIN, DTMLTD, &
556 DPT, DPDD, a1, rdd, DQSSDT, DTMP, QSRH, PPTFLX, PPTFL2, CPR, CNDTNF, &
557 UPDINC, DMFMIN, DEVDMF, PPR, RCED, DPPTDF, DMFLFS, &
558 DDINC, FABE, STAB, TKEMAX, DSOURCE, TIMECS, &
559 DTT, DTT1, GD, CPM, GM, DTEMPDT, DTHETADT, &
560 TMA, TMB, TMM, BCOEFF, ACOEFF, TOPOMG, DTMLTE, TDP, &
561 DQ, QMIN, TLOG, TDPT, CPORQ, DLP, ZLCL1, ABEG, THATA, &
562 THTFC, DABE, AINCMIN, RF, UPDIN2, &
563 DR, UDFRC, UEFRC, DDFRC, DEFRC, TUC, TDC, QGS, RH_0, RHG, &
564 QINIT, QFNL, QCFINL, ERR2, RELERR, &
565 DIFFK, RLC, E2, EOVLC, RCP, RLOVCP, &
566 RHCR, AREADIFF, FX, FX1, FX2, &
567 EFOLD, EX, DEPTHMAX, HIN, HOUT, DELTH, QIN, QOUT, &
568 DELTQ, RKAPA, DELTA, RATIO, RATIOMIN, RATIOMAX, &
569 GAMA, BEITA, SIGMMA, WFUN, BVT3, BVTS3, PPI, G3PB, G3PBS, PRAC, PRACS, &
570 TOUT, SC2, SC3, SC7, RHOS, PPIS, XNC, &
571 RHO2, VT2C, FALTNDC, gdry, p300, factor1, &
572 z1, p1, t1, w1, z2, p2, t2, w2, r1, &
573 t1rh, dtime, ROVCP
575 INTEGER :: ii, k, ks, k1, k2, kk, k0, ndk, nj, nm, nk, nk1, &
576 kval, kval1, kval2, iflag, kflag, nic, &
577 KMIN, KSTART, ND, NT, ND1, LDB, LDT, LMAX, NCOUNT, ISTOP, &
578 k100, k1000, KAMAX, KLCMAX, dbg_level, i0, j0, KRADLMAX, &
579 NTC, NSTEP, ML, L5, LLFC, LM, LM1, LM2, LC, lvf, kradsmax
581 REAL :: ALIQ, BLIQ, CLIQ, DLIQ
583 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
585 ! REAL :: tpdd ! Function declaration
586 ! REAL :: GAMMA ! Function declaration
588 CHARACTER*1024 message
590 #ifdef DENG_SHCU1D
591 REAL, DIMENSION( its:ite , jts:jte ) :: &
592 TPBLG,PPBLG,TLCLG,PLCLG, & ! used to be 1-d mm5 graphics
593 CLDHGTG,TCLDTPG,PCLDTPG,TLFSG,PLFSG,TLDTG,PLDTG, &
594 TLDBG,PLDBG,TLLCG,PLLCG
596 REAL, DIMENSION( its:ite , jts:jte ) :: &
597 CLDBMFLX, &
598 nca_ctr1,nca_ctr2,nca_ctr3,nca_ctr4, &
599 dcldtop, dcldbase
601 REAL, DIMENSION( kts:kte ) :: tt4, tt6, tt7, wsubsub, wutemp, cldaten0, cldaten, &
602 cldlctn0, cldlctn1, cldlctn2, cldlctn3, cldlctn4, &
603 cldlctn5, thtetmp, thtatmp, cldua, tug, wsp_plot, &
604 wdr_plot, td_plot, ca_rad0
606 REAL :: PLET,TLET, aa=60.0*60.0, bb=60.0*60.0*1000.0, qqq, ees, &
607 PLOTFRQ = 30.0
609 INTEGER :: KTAU1HUR, KCBASE
610 #endif
611 !******************************************************************************
612 !
613 CALL get_wrf_debug_level( dbg_level )
615 ISTOP= 0
616 NT = 0
617 ALIQ = SVP1*1000.
618 BLIQ = SVP2
619 CLIQ = SVP2*SVPT0
620 DLIQ = SVP3
622 GDRY = -G/CP
623 ROVCP=R/CP
625 KL=kte
626 klm = kl - 1
627 KX=kte
628 kxp1 = kx + 1
630 i0 = (ite-its)/2+its
631 j0 = (jte-jts)/2+jts
633 #ifdef DENG_SHCU1D
634 kcbase=1
635 KTAU1HUR = 3600/INT(DT)
636 IF( i == 1 .AND. j == 1 ) write(31,*) xtime/60.,zpbl
637 do k=1,kl
638 tt4(k)=0.0
639 tt6(k)=0.0
640 tt7(k)=0.0
641 wutemp(k)=0.0
642 wsubsub(k)=0.0
643 cldaten0(k)=0.0
644 cldaten(k)=0.0
645 cldlctn0(k)=0.0
646 cldlctn1(k)=0.0
647 cldlctn2(k)=0.0
648 cldlctn3(k)=0.0
649 cldlctn4(k)=0.0
650 cldlctn5(k)=0.0
651 thtetmp(k)=0.0
652 thtatmp(k)=0.0
653 enddo
654 #endif
656 LTOP=1
657 KPBLMAX = 1
658 CLDDPTH=0.0
659 CLDTOP=0.0
660 DO K=1,KL
661 #ifdef DENG_SHCU1D
662 CLDUA(K)=0.0
663 TUG(K)=0.0
664 #endif
665 WU(K)=0.0
666 OMG(K)=0.0
667 UMF(K)=0.0
668 THTAD(K)=0.0
669 WD(K)=0.0
670 KV(K)=0.0
671 KR(K)=0.0
672 ENDDO
673 RAINSHV(I,J)=0.0
674 AINCKFI(I,J)=0.0
675 CAPEI(I,J)=0.0
676 #ifdef DENG_SHCU1D
677 DCLDTOP(I,J)=0.0
678 DCLDBASE(I,J)=0.0
679 CLDBMFLX(I,J)=0.0
680 #endif
681 RADIUSC(I,J)=0.0
682 !
683 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED FROM THE
684 !...BOTTOM-UP IN THE SHALLOW CONVECTION SCHEME...
685 !
687 P300=P0(1)-30000.
688 ML = 0
690 DO K=1,KX
691 QC0(K)=0.0
692 DCA0(K)=0.0
693 !
694 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
695 !
696 ! ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
697 IF( T0(K) .GT. TO ) THEN
698 ES=1.E3*SVP1*EXP(SVP2*(T0(K)-SVPT0)/(T0(K)-SVP3))
699 ELSE
700 ES=611.0*EXP(22.514-6.15E3/T0(K))
701 ENDIF
702 QES(K)=0.622*ES/(P0(K)-ES)
703 Q0(K)=AMIN1(QES(K),Q0(K))
704 IF(Q0(K).GT.QES(K))Q0(K)=QES(K)
706 TV0(K)=T0(K)*(1.+.608*Q0(K))
707 ! RHOE(K)=P0(K)/(R*TV0(K))
708 !
709 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVELS,
710 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
711 !
712 DP(K)=rhoe(k)*g*DZQ(k)
713 IF(P0(K).GE.500E2) L5=K
714 IF(P0(K).GE.P300) LLFC=K
715 IF(T0(K).GT.T00) ML=K
716 ENDDO
718 Z0(1)=.5*DZQ(1)
719 DO K=2,KL
720 Z0(K) = Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
721 DZA(K-1) = Z0(K)-Z0(K-1)
722 ENDDO
723 DZA(KL) = 0.
725 LC = 1
726 KFLAG=0
727 AINC=999.0
728 CTP2=99999.9
729 LTOP2=99999
730 ! ZPBL = PBL(I,J)
731 !
732 !...FIND CONVECTIVE SOURCE LAYER WHERE THE INITIAL PARCEL CHARACTERISTICS
733 ! ARE DEFINED. IT IS THE MAXIMUM OF 20% OF PBL AND 2 (IF PBL IS AT LEAST
734 ! 4 LATERS DEEP. DEFINE THE PARCEL PURTERBATION PROPERTITES FROM THE
735 ! SOURCE LAYER (THMIX,QMIX,ZMIX,PMIX)
736 !
737 LM1 = 1
738 FACTOR1 = 0.20 ! 20% of PBL
739 source: DO NK = 1,KL
740 IF( (Z0(NK)+0.5*DZQ(NK)) .LT. FACTOR1*ZPBL ) THEN
741 LM1=LM1+1
742 ELSE
743 EXIT source
744 ENDIF
745 ENDDO source
747 kpbl = 1
749 loop_ku: DO k = 1, kl
750 zagl_bot = z_at_w0(k)-ht
751 zagl_top = z_at_w0(k+1)-ht
752 IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
753 kpbl = k
754 EXIT loop_ku
755 ENDIF
756 ENDDO loop_ku
758 LM2=1
759 IF( KPBL .GE. 4 ) LM2=2
760 LM=MAX(LM1,LM2)
762 THMIX=0.
763 QMIX=0.
764 ZMIX=0.
765 PMIX=0.
766 DPTHMX=0.
767 DO NK = 1, LM
768 DPTHMX=DPTHMX+DP(NK)
769 ROCPQ=0.2854*(1.-0.28*Q0(NK))
770 THMIX=THMIX+DP(NK)*T0(NK)*(P00/P0(NK))**ROCPQ
771 QMIX=QMIX+DP(NK)*Q0(NK)
772 ZMIX=ZMIX+DP(NK)*Z0(NK)
773 PMIX=PMIX+DP(NK)*P0(NK)
774 ENDDO
775 THMIX=THMIX/DPTHMX
776 QMIX=QMIX/DPTHMX
777 ZMIX=ZMIX/DPTHMX
778 PMIX=PMIX/DPTHMX
779 !
780 ! FIND THE LCL (TLCL,TVLCL,PLCL,KLCL)
781 !
782 ROCPQ=0.2854*(1.-0.28*QMIX)
783 TMIX=THMIX*(PMIX/P00)**ROCPQ
784 EMIX=QMIX*PMIX/(0.622+QMIX)
785 ! TLOG=ALOG(EMIX/ALIQ)
786 ! TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
787 TLOG=ALOG(EMIX/SVP1/1000.)
788 TDPT=(SVP2*SVPT0-SVP3*TLOG)/(SVP2-TLOG)
789 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))* &
790 (TMIX-TDPT)
791 TLCL=AMIN1(TLCL,TMIX)
792 TVLCL=TLCL*(1.+0.608*QMIX)
793 CPORQ=1./ROCPQ
794 PLCL=P00*(TLCL/THMIX)**CPORQ
796 DO NK = 1,KL
797 KLCL=NK
798 IF( PLCL >= P0(NK) ) GO TO 35
799 ENDDO
800 GOTO 425
801 35 CONTINUE
802 K=KLCL-1
804 IF(KLCL .EQ. 1 ) THEN
805 #ifdef DENG_SHCU1D
806 write(*,*) 'KL layer is saturated and skip over active con.'
807 PLET = PLCL
808 TLET = TLCL
809 PLFSG(I,J) = PLCL
810 TLFSG(I,J) = TLCL
811 PLDTG(I,J) = PLCL
812 TLDTG(I,J) = TLCL
813 PLDBG(I,J) = PLCL
814 TLDBG(I,J) = TLCL
815 PLLCG(I,J) = PLCL
816 TLLCG(I,J) = TLCL
817 CLDHGTG(I,J)=0.0
818 PCLDTPG(I,J)=PLCL
819 TCLDTPG(I,J)=TLCL
820 #endif
822 IF(ZPBL .LE. Z0(1)) THEN
823 PPBL=P0(1)
824 ELSE
825 loop12: DO KK=1, KL-1
826 IF(ZPBL .GE. Z0(KK) .AND. ZPBL .LT. Z0(KK+1)) THEN
827 K0=KK
828 EXIT loop12
829 END IF
830 ENDDO loop12
832 C1=(ZPBL-Z0(K0))/(Z0(K0+1)-Z0(K0))
833 PPBL=P0(K0)*EXP(C1*(ALOG(P0(K0+1))-ALOG(P0(K0))))
834 END IF
835 TPBL=TLCL*(PPBL/PLCL)**ROCPQ
837 ZLCL=Z0(KLCL)
838 CLDTOP=ZLCL
839 ZLFC=Z0(KL)
840 #ifdef DENG_SHCU1D
841 nca_ctr1(i,j)=nca_ctr1(i,j)+1
842 #endif
843 GOTO 425
845 ! KLCL=KLCL+1
846 ! K=KLCL-1
847 ENDIF
849 DLP=LOG(PLCL/P0(K))/LOG(P0(KLCL)/P0(K))
850 !
851 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
852 !
853 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
854 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
855 TVEN=TENV*(1.+0.608*QENV)
856 TVBAR=0.5*(TV0(K)+TVEN)
857 ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
858 !
859 ! DETERMINE PURTERBATION VERTICAL VELOCITIES FOR INITIATING PARCEL,
860 ! WPBLTP=WBAR+WTKE+WNH, WHERE WBAR IS THE RESOLVABLE-SCALE MEAN VERTICAL
861 ! MOTION, WTKE IS THE CONTRIBUTION FROM THE PBL TKE AND WNH IS THE NON-HYDROSTAIC
862 ! PUMPING EFFECT
863 !
864 IF( zpbl .GT. PBLMAX(I,J) ) THEN
865 PBLMAX(I,J) = zpbl
866 KPBLMAX = KPBL
867 ENDIF
868 !
869 ! CALCULATE LOCAL TIME USING UTC (DATE) AND THE LONGITUDE.
870 !
871 TIMEUNV = gmt + XTIME/60.
872 IF( TIMEUNV .GE. 24.0 ) TIMEUNV = TIMEUNV - 24.0
873 TIMELOC = TIMEUNV + XLONG(I,J)/15.0
874 IF( TIMELOC .LT. 0.0 ) TIMELOC = TIMELOC + 24.0
875 IF( TIMELOC .GE.24.0 ) TIMELOC = TIMELOC - 24.0
877 IF( ABS(TIMELOC-6.0) .LT. 0.01 ) THEN
878 PBLMAX(I,J) = 0.0 ! SET PBLMAX=0 IF LOCAL TIME=6AM
879 KPBLMAX = 1
880 END IF
881 KVAL=MIN(KLCL,KPBLMAX)
883 EMAX = 0.0
884 loop154: DO nk = 1, KVAL
885 val = TKE0(nK)
886 IF( val.LT.EMAX ) CYCLE loop154
887 EMAX = val
888 ENDDO loop154
890 KVAL=MIN(KLCL-1,KPBL)
892 WBAR=W0AVG0(KVAL)
893 WTKE=SQRT(2.0*EMAX/3.0)
895 WPBLTP=WBAR+WTKE
897 WNH=0.0
898 IF( CLDDPTHB(I,J) .GT. 4.0E3 ) THEN
899 ZVAL = ZLCL + 2000.0
900 DO NK = 1,KL
901 IF( ABS(Z0(NK)-ZVAL) .LT. DZQ(NK)/2.0 ) KVAL = NK
902 ENDDO
903 IF( WUB(I,KVAL,J) .GT. WPBLTP ) THEN
904 WNH=0.25*(WUB(I,KVAL,J)-WPBLTP)
905 ELSE
906 WNH=0.0
907 ENDIF
908 ENDIF
910 WPBLTP=WPBLTP+WNH
911 #ifdef DENG_SHCU1D
912 IF( i == 1 .AND. j == 1 ) write(45,*) xtime/60.,WBAR,WTKE,WNH,WPBLTP
913 #endif
914 !
915 !...THE PARCEL IS BUOYANT; COMPUTE EQUIVALENT POTENTIAL TEMPERATURE (THETEU)
916 !...AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...USE BOLTONS (1980)
917 !...FORMULA FOR THETE...
918 !
919 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
920 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
921 ES=1.E3*SVP1*EXP(SVP2*(TENV-SVPT0)/(TENV-SVP3))
922 TVAVG=0.5*(TV0(KLCL)+TENV*(1.+0.608*QENV))
923 QESE=0.622*ES/(PLCL-ES)
924 THTES(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
925 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
926 WTW=WPBLTP*WPBLTP
927 TVLCL=TLCL*(1.+0.608*QMIX)
928 RHOLCL=PLCL/(R*TVLCL)
929 LCL=KLCL
930 LET=LCL
931 !
932 ! DEFINE THE UPDRAFT RADIUS, RAD, AS A FUNCTION OF DEPTHS OF THE PBL
933 ! AND UPDRAFT AT THE PREVIOUS TIME STEP. RAD IS IN THE RANGE OF 150
934 ! AND 1500 M
935 !
936 VAL=CLDDPTHB(I,J)
937 VAL1=ZPBL
938 IF(VAL/1000.0 .GT. 3.9) THEN
939 RAD=1500.0
940 GO TO 1001
941 END IF
942 FF=12.0*VAL/(4.0-VAL/1000.)/1000.
943 H=VAL1*FF
944 B=-0.5*(7.0+2.0*H/1000.)
945 RAD=0.5*(-B-SQRT(B*B-12.0*H/1000.0))
946 RAD=RAD*0.5*1000.
947 IF(RAD .LT. 150.0) RAD=150.0
948 IF(RAD .GT.1500.0) RAD=1500.0
949 1001 CONTINUE
950 ! write(*,'(i7, 2i3, 4f10.3)') ktau, i, j, xtime/60.0, rad, zpbl, CLDDPTHB(I,J)
952 ! Averaging the radius between current and the previous time value
953 !
954 IF(CAPESAVE(I,J) .LT. 100.0) RAD=(RAD+RADSAVE(I,J))/2.0
956 RADIUSC(i,j)=RAD
957 !
958 ! FIND PBL PROPERTIES
959 !
960 IF( ZPBL .LE. Z0(1) ) THEN
961 K = 1
962 PPBL = P0(K)
963 TENVPBL = T0(K)
964 QENVPBL = Q0(K)
965 TVENPBL = TV0(K)
966 ELSE
967 loop13: DO KK = 1, KL-1
968 IF( ZPBL .GE. Z0(KK) .AND. ZPBL .LT. Z0(KK+1) ) THEN
969 K = KK
970 EXIT loop13
971 END IF
972 ENDDO loop13
974 C1 = ( ZPBL-Z0(K) )/( Z0(K+1)-Z0(K) )
975 PPBL = P0(K)*EXP( C1 * ( ALOG(P0(K+1))-ALOG(P0(K)) ) )
976 TENVPBL = C1* (T0(K+1)-T0(K)) + T0(K)
977 QENVPBL = C1* (Q0(K+1)-Q0(K)) + Q0(K)
978 TVENPBL = C1* (TV0(K+1)-TV0(K)) + TV0(K)
979 END IF
980 KKPBL = K
981 TPBL=TLCL*(PPBL/PLCL)**ROCPQ
982 TVPBL = TPBL*(1.+0.608*QMIX)
983 RHOPBL=PPBL/(R*TVPBL)
984 !
985 ! THE LOWER OF LCL AND PBL IS DETERMINED AS THE PARCEL RELEASING POINT
986 ! WHERE K,RHO_INI, T_INI, TV_INI, TVEN_IN, Z_INI AND P_INI ARE DEFINED
987 !
988 IF( PPBL .GE. PLCL ) THEN
989 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
990 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
991 ES=1.E3*SVP1*EXP(SVP2*(TENVPBL-SVPT0)/(TENVPBL-SVP3))
992 QESE=0.622*ES/(PPBL-ES)
993 THTES(K)=TENVPBL*(1.E5/PPBL)**(0.2854*(1.-0.28*QESE))* &
994 EXP((3374.6525/TENVPBL-2.5403)*QESE*(1.+0.81*QESE))
995 RHO_INI = RHOPBL
996 T_INI = TPBL
997 TV_INI = TVPBL
998 TVEN_INI = TVENPBL
999 Z_INI = ZPBL
1000 P_INI = PPBL
1001 LET = KKPBL
1002 ELSE
1003 K = KLCL-1
1004 RHO_INI = RHOLCL
1005 T_INI = TLCL
1006 TV_INI = TVLCL
1007 TVEN_INI= TVEN
1008 Z_INI = ZLCL
1009 P_INI = PLCL
1010 END IF
1011 !
1012 ! DETERMIN THE DEPTH OF SOURCE LAYER WHERE THE MASS IS IS REMOVED AS A RESULTS OF
1013 ! CONVECTION. THE DEPTH IS DEFINED AS A FUNCTION OF CLOUD RADIUS AND IS IN THE
1014 ! RANGE OF 10 MB TO 60 MB
1015 !
1016 DPTHMX1 = 0.0
1017 DO NK = 1, K
1018 DPTHMX1 = DPTHMX1 + DP(NK)
1019 ENDDO
1020 IF( K .EQ. 1 ) THEN
1021 LLC=1
1022 GO TO 897
1023 ENDIF
1024 DPTH = 100.0*(RAD-150.0)/27.0 + 1000.0
1025 IF( DPTH .GE. DPTHMX1 ) THEN
1026 LLC = 1
1027 GO TO 897
1028 END IF
1029 PPLSDP = P0(K+1) + DP(K+1)/2.0 + DPTH
1030 LLC = K
1031 DO NK = K-1, 1, -1
1032 IF( ABS(PPLSDP-P0(NK)) .LE. ABS(PPLSDP-P0(NK+1)) ) LLC = NK
1033 ENDDO
1035 897 CONTINUE
1037 #ifdef DENG_SHCU1D
1038 KCBASE=KLCL
1039 IF( i == 1 .AND. j == 1 ) write(300,*) xtime/60.,K, RHO_INI, T_INI, TV_INI, TVEN_INI, Z_INI, P_INI
1040 #endif
1041 !
1042 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
1043 ! BEGINING OF UPDRAFT CACULATION
1044 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
1045 !
1046 ! COMPUTE INITIAL UPDRAFT MASS FLUX UMF(K)
1047 !
1048 WU(K)=WPBLTP
1049 AU0=PIE*RAD*RAD
1050 UMF(K)=RHO_INI*AU0*WPBLTP
1051 VMFPBLCL=UMF(K)
1052 VMFLCL=UMF(K)
1053 UPOLD=VMFPBLCL
1054 UPNEW=UPOLD
1055 RATIO2(K)=0.
1056 !
1057 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE BUOYANT ENERGY,
1058 ! TRPPT IS THE TOTAL RATE OF PRECIPITATION PRODUCTION...
1059 !
1060 UER(K)=0.
1061 ABE=0.
1062 TRPPT=0.
1063 TU(K)=T_INI
1064 TVU(K)=TV_INI
1065 QU(K)=QMIX
1066 EQFRC(K)=1.
1067 RLIQ(K)=0.
1068 RICE(K)=0.
1069 QLQOUT(K)=0.
1070 QICOUT(K)=0.
1071 DETLQ(K)=0.
1072 DETIC(K)=0.
1073 PPTLIQ(K)=0.
1074 PPTICE(K)=0.
1075 IFLAG = 0
1076 ! KFRZ=LC
1077 !
1078 !...THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH RESPECT TO
1079 !...UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILT PARCEL IS THTUDL, UNDILT
1080 ! TEMP IS GIVEN BY TUDL...
1081 !
1082 THTUDL=THETEU(K)
1083 TUDL=T_INI
1084 !
1085 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION PROCESS; IT IS
1086 ! INITIALLY SET TO THE TEMPERATURE AT WHICH FREEZING IS SPECIFIED TO BEGIN;
1087 ! WITHIN THE GLACIATION INTERVAL, IT IS SET EQUAL TO THE UPDR TEMP AT THE
1088 ! PREVIOUS MODEL LEVEL...
1089 !
1090 TTEMP=TTFRZ
1091 !
1092 ! FOR SHALLOW CONVECTION, SET RATE=0 SO NO LIQUID FALLS OUT
1093 !
1094 IF( CLDDPTHB(I,J) .LT. 4.E3) THEN
1095 RATE = 0.0 ! TURN OFF DOWNDRAFT ALSO
1096 ELSE
1097 RATE = 0.01
1098 ENDIF
1099 !
1100 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP, MIXING
1101 ! RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND MOISTURE, AND
1102 ! PRECIPITATION RATES AT EACH MODEL LEVEL...
1103 !
1104 EE1=1.
1105 UD1=0.
1106 REI = 0.
1107 loop60: DO NK=K,KLM
1108 NK1=NK+1
1109 RATIO2(NK1)=RATIO2(NK)
1110 !
1111 !...UPDATE UPDRAFT PROPERTIES AT THE NEXT MODEL LVL TO REFLECT ENTRAINMNT OF
1112 ! ENVIRONMENTAL AIR...
1113 !
1114 FRC1=0.
1115 TU(NK1)=T0(NK1)
1116 THETEU(NK1)=THETEU(NK)
1117 QU(NK1)=QU(NK)
1118 RLIQ(NK1)=RLIQ(NK)
1119 RICE(NK1)=RICE(NK)
1120 CALL TPMIX(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),RLIQ(NK1), &
1121 RICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1, &
1122 SVP1,SVP2,SVPT0,SVP3)
1123 ! CALL TPMIXBG(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),RLIQ(NK1), &
1124 ! RICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1, &
1125 ! SVP1,SVP2,SVPT0,SVP3)
1127 ! CALL TPMIX2(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),RLIQ(NK1), &
1128 ! RICE(NK1),QNEWLQ,QNEWIC,XLV1,XLV0,ktau,i,j,nk1, 1)
1129 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
1130 !
1131 ! CHECK TO SEE IF UPDRAFT TEMP IS WITHIN THE FREEZING INTERVAL, IF IT IS,
1132 ! CALCULATE THE FRACTIONAL CONVERSION TO GLACIATION AND ADJUST QNEWLQ TO
1133 ! REFLECT THE GRADUAL CHANGE IN THETAU SINCE THE LAST MODEL LEVEL...THE
1134 ! GLACIATION EFFECTS WILL BE DETERMINED AFTER THE AMOUNT OF CONDENSATE
1135 ! AVAILABLE AFTER PRECIP FALLOUT IS DETERMINED...TTFRZ IS THE TEMP AT WHICH
1136 ! GLACIATION BEGINS, TBFRZ THE TEMP AT WHICH IT ENDS...
1137 !
1138 IF(TU(NK1).LE.TTFRZ.AND.IFLAG.LT.1)THEN
1139 IF(TU(NK1).GT.TBFRZ)THEN
1140 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
1141 FRC1=(TTEMP-TU(NK1))/(TTFRZ-TBFRZ)
1142 R1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
1143 ELSE
1144 FRC1=(TTEMP-TBFRZ)/(TTFRZ-TBFRZ)
1145 R1=1.
1146 IFLAG=1
1147 ENDIF
1148 QNWFRZ=QNEWLQ
1149 QNEWIC=QNEWIC+QNEWLQ*R1*0.5
1150 QNEWLQ=QNEWLQ-QNEWLQ*R1*0.5
1151 IF( ABS(TTEMP-TBFRZ) < 1.E-3 ) THEN
1152 EFFQ= 1.0
1153 ELSE
1154 EFFQ=(TTFRZ-TBFRZ)/(TTEMP-TBFRZ)
1155 ENDIF
1156 TTEMP=TU(NK1)
1157 ENDIF
1158 !
1159 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
1160 !
1161 IF(NK.EQ.K)THEN
1162 BE=(TV_INI+TVU(NK1))/(TVEN_INI+TV0(NK1))-1.
1163 BOTERM=2.*(Z0(NK1)-Z_INI)*G*BE/1.5
1164 ENTERM=0.
1165 DZZ=Z0(NK1)-Z_INI
1166 ELSE
1167 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
1168 BOTERM=2.*DZA(NK)*G*BE/1.5
1169 ENTERM=2.*UER(NK)*WTW/UPOLD
1170 DZZ=DZA(NK)
1171 ENDIF
1172 WSQ=WTW
1173 CALL CONDLOAD(RLIQ(NK1),RICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
1174 RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1), G)
1175 IF(WTW < 1.E-3) THEN
1176 EXIT loop60
1177 ELSE
1178 WU(NK1)=SQRT(WTW)
1179 ENDIF
1180 !
1181 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND, IF CLOUD
1182 ! IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
1183 !
1184 ! IF(WU(NK1).LT.0.) EXIT loop60
1185 !
1186 ! UPDATE THE ABE FOR UNDILUTE ASCENT...
1187 !
1188 THTES(NK1)=T0(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QES(NK1)))* &
1189 EXP((3374.6525/T0(NK1)-2.5403)*QES(NK1)*(1.+0.81*QES(NK1)))
1190 UDLBE=((2.*THTUDL)/(THTES(NK)+THTES(NK1))-1.)*DZZ
1191 IF(UDLBE.GT.0.)ABE=ABE+UDLBE*G
1192 !
1193 ! DETERMINE THE EFFECTS OF CLOUD GLACIATION IF WITHIN THE SPECIFIED
1194 ! TEMP INTERVAL...
1195 !
1196 IF(FRC1.GT.1.E-6)THEN
1197 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),RLIQ(NK1), &
1198 RICE(NK1),RATIO2(NK1),QNWFRZ,RL,FRC1,EFFQ, &
1199 IFLAG,XLV0,XLV1,XLS0,XLS1, &
1200 SVP1,SVP2,SVPT0,SVP3)
1201 ENDIF
1202 !
1203 ! CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMP...
1204 ! WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED WITH RESPECT TO THE
1205 ! SAME DEGREE OF GLACIATION FOR ALL ENTRAINING AIR...
1206 !
1207 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1), &
1208 THETEE(NK1),RATIO2(NK1),RL, &
1209 SVP1,SVP2,SVPT0,SVP3,I,J)
1210 ! CALL ENVIRTHT2(P0(NK1),T0(NK1),Q0(NK1), &
1211 ! THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1212 !
1213 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
1214 !
1215 REI=VMFPBLCL*DP(NK1)*0.03/RAD
1216 !
1217 !...TVQU IS THE PARCEL TEMP INCLUDING WATER LOADING...
1218 !...IF IT BECOMES GREATER THAN TV0, THEN YOU ARE PAST LFC...
1219 !...
1220 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-RLIQ(NK1)-RICE(NK1))
1221 !
1222 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, NO
1223 ! ENTRAINMENT IS ALLOWED AT THIS LEVEL...
1224 !
1225 IF(TVQU(NK1).LE.TV0(NK1))THEN
1226 UER(NK1)=0.0
1227 UDR(NK1)=REI
1228 EE2=0.
1229 UD2=1.
1230 EQFRC(NK1)=0.
1231 GOTO 55
1232 ENDIF
1233 LET=NK1
1234 TTMP=TVQU(NK1)
1235 !
1236 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
1237 !
1238 F1=0.95
1239 F2=1.-F1
1240 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
1241 RTMP=F1*Q0(NK1)+F2*QU(NK1)
1242 TMPLIQ=F2*RLIQ(NK1)
1243 TMPICE=F2*RICE(NK1)
1244 CALL TPMIX(P0(NK1),THTTMP,TTMP,RTMP,TMPLIQ,TMPICE,QNEWLQ, &
1245 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1, &
1246 SVP1,SVP2,SVPT0,SVP3)
1247 ! CALL TPMIX2(P0(NK1),THTTMP,TTMP,RTMP,TMPLIQ,TMPICE,QNEWLQ, &
1248 ! QNEWIC,XLV1,XLV0,ktau,i,j,nk1,2)
1249 ! CALL TPMIXBG(P0(NK1),THTTMP,TTMP,RTMP,TMPLIQ,TMPICE,QNEWLQ, &
1250 ! QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1, &
1251 ! SVP1,SVP2,SVPT0,SVP3)
1252 TU95=TTMP*(1.+0.608*RTMP-TMPLIQ-TMPICE)
1253 IF(TU95.GT.TV0(NK1))THEN
1254 EE2=1.
1255 UD2=0.
1256 EQFRC(NK1)=1.0
1257 GOTO 50
1258 ENDIF
1259 F1=0.10
1260 F2=1.-F1
1261 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
1262 RTMP=F1*Q0(NK1)+F2*QU(NK1)
1263 TMPLIQ=F2*RLIQ(NK1)
1264 TMPICE=F2*RICE(NK1)
1265 CALL TPMIX(P0(NK1),THTTMP,TTMP,RTMP,TMPLIQ,TMPICE,QNEWLQ, &
1266 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1, &
1267 SVP1,SVP2,SVPT0,SVP3)
1268 ! CALL TPMIX2(P0(NK1),THTTMP,TTMP,RTMP,TMPLIQ,TMPICE,QNEWLQ, &
1269 ! QNEWIC,XLV1,XLV0,ktau,i,j,nk1,3)
1270 ! CALL TPMIXBG(P0(NK1),THTTMP,TTMP,RTMP,TMPLIQ,TMPICE,QNEWLQ, &
1271 ! QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1, &
1272 ! SVP1,SVP2,SVPT0,SVP3)
1273 TU10=TTMP*(1.+0.608*RTMP-TMPLIQ-TMPICE)
1274 IF( ABS(TU10-TVQU(NK1)) < 1.E-3) THEN
1275 EQFRC(NK1)= 1.0
1276 ELSE
1277 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1)+1.E-20)
1278 ENDIF
1279 EQMAX = 1.
1280 EQMIN = 0.
1281 EQFRC(NK1)=MAX(EQMIN,EQFRC(NK1))
1282 EQFRC(NK1)=MIN(EQMAX,EQFRC(NK1))
1283 IF(EQFRC(NK1).EQ.1)THEN
1284 EE2=1.
1285 UD2=0.
1286 GOTO 50
1287 ELSEIF(EQFRC(NK1).EQ.0.)THEN
1288 EE2=0.
1289 UD2=1.
1290 GOTO 50
1291 ELSE
1292 !
1293 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
1294 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
1295 !
1296 CALL PROF5(EQFRC(NK1),EE2,UD2)
1297 ENDIF
1299 50 CONTINUE
1301 ! IF( NK .EQ. K) THEN
1302 ! EE1=1.
1303 ! UD1=0.
1304 ! ENDIF
1305 !
1306 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
1307 ! VALUES IN THE LAYER...
1308 !
1309 UER(NK1)=0.5*REI*(EE1+EE2)
1310 UDR(NK1)=0.5*REI*(UD1+UD2)
1311 !
1312 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
1313 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
1314 !
1315 55 CONTINUE
1317 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
1318 !
1319 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
1320 ! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
1321 !
1322 IF(UDLBE.GT.0.)ABE=ABE-UDLBE*G
1323 LET=NK
1324 EXIT loop60
1325 ENDIF
1327 EE1=EE2
1328 UD1=UD2
1329 UPOLD=UMF(NK)-UDR(NK1)
1330 UPNEW=UPOLD+UER(NK1)
1331 UMF(NK1)=UPNEW
1332 !
1333 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND ICE IN THE
1334 ! DETRAINING UPDRAFT MASS...
1335 !
1336 DETLQ(NK1)=RLIQ(NK1)*UDR(NK1)
1337 DETIC(NK1)=RICE(NK1)*UDR(NK1)
1338 QDT(NK1)=QU(NK1)
1339 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
1340 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
1341 RLIQ(NK1)=RLIQ(NK1)*UPOLD/UPNEW
1342 RICE(NK1)=RICE(NK1)*UPOLD/UPNEW
1343 !
1344 !...KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE IS GENERATED...
1345 ! PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF LIQUID PRECIP AT A GIVEN
1346 ! MODEL LVL, PPTICE THE SAME FOR ICE, TRPPT IS THE TOTAL RATE OF PRODUCTION
1347 ! OF PRECIP UP TO THE CURRENT MODEL LEVEL...
1348 !
1349 ! IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
1350 PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
1351 PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
1352 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
1353 ENDDO loop60
1355 CAPEI(I,J) = ABE
1356 #ifdef DENG_SHCU1D
1357 IF( i == 1 .AND. j == 1 ) write(88,*) xtime/60.,ABE
1358 #endif
1359 !
1360 ! UPDRAFT TOP IS DEFINED AS THE LAGER VALUE OF LTOP1 (CTP1) AND LTOP1 (CTP2)
1361 ! WHERE LTOP1 IS THE KAIN-FRITSCH UPDRAFT TOP AND LTOP2 IS THE HEIGHT CALCULATED
1362 ! USING THE EQUATION: CTP2=CLDTOPB(I,J)+0.2*WUMAX*DT2/2.0
1363 !
1364 LTOP1=MIN(NK,KLM)
1365 CTP1=Z0(LTOP1)
1367 WUMAX=WU(1)
1368 DO NK=2,LTOPB(I,J)
1369 IF(WU(NK) .GT. WUMAX) WUMAX=WU(NK)
1370 ENDDO
1372 IF(LTOPB(I,J) .LE. 1) GOTO 666
1373 CTP2=CLDTOPB(I,J)+0.2*WUMAX*DT2/2.0
1374 LTOP2=1
1375 IF(CTP2 .GE. (Z0(KL)+DZQ(KL)/2.0)) THEN
1376 CTP2=Z0(KL)
1377 LTOP2=KL
1378 ELSE
1379 DO NK=1,KL
1380 IF(CTP2 .LT. (Z0(NK)+DZQ(NK)/2.0) .AND. &
1381 CTP2 .GE. (Z0(NK)-DZQ(NK)/2.0)) LTOP2 = NK
1382 ENDDO
1383 CTP2=MAX(CTP2,ZLCL)
1384 LTOP2=MAX(LTOP2,KLCL)
1385 ENDIF
1387 CLDTOP=MIN(CTP1,CTP2)
1388 LTOP=MIN(LTOP1,LTOP2)
1389 GOTO 667
1390 666 CONTINUE
1391 CLDTOP=CTP1
1392 LTOP=LTOP1
1393 667 CONTINUE
1394 #ifdef DENG_SHCU1D
1395 IF( i == 1 .AND. j == 1 ) write(89,'(8f10.3,i3)') xtime/60.,CLDTOPB(I,J),CTP1,CTP2, &
1396 CLDTOP,ZLCL,WUMAX,ABE, LTOPB(I,J)
1397 #endif
1398 !
1399 ! MODIFY TRPPT CONSIDERING ACTUAL CLOUD TOP WHICH IS
1400 ! LOWER IN MOST CASES
1401 !
1402 IF(LTOP .EQ. LTOP2 .AND. LTOP1 .GT. LTOP2) THEN
1403 DO NK=LTOP2+1,LTOP1
1404 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1405 ENDDO
1406 ENDIF
1407 CLDDPTH=CLDTOP-ZLCL
1408 IF(CLDDPTH .LE. 0.0) CLDDPTH=0.0
1409 #ifdef DENG_SHCU1D
1410 CLDHGTG(I,J)=CLDDPTH ! Graphics
1411 PCLDTPG(I,J)=P0(LTOP) ! Graphics
1412 TCLDTPG(I,J)=TU(LTOP) ! Graphics
1413 #endif
1414 !
1415 ! FIND LEVEL OF FREE CONVECTION. WHEN THE PARCCEL IS WARMER THAN THE ENVIRONMENT
1416 ! AT ITS LCL, IT NEED AT LEAST 400 M TO SAY THAT LFC IF AT LCL
1417 !
1418 ZLFC=Z0(KL)
1419 IF(CLDDPTH .LE. 0.0) THEN
1420 ZLFC=Z0(KL)
1421 GOTO 67
1422 ELSE IF(TVQU(KLCL) .GE. TV0(KLCL)) THEN
1423 SUM=0.0
1424 KVAL=MAX(KLCL,KKPBL)
1425 loop61: DO KS=KVAL,LTOP1
1426 IF(TVQU(KS).GT.TV0(KS)) THEN
1427 SUM=SUM+DZQ(KS)
1428 ELSE
1429 EXIT loop61
1430 ENDIF
1431 ENDDO loop61
1433 IF(SUM .GE. 300.) THEN
1434 ZLFC=ZLCL
1435 GO TO 67
1436 ENDIF
1437 ENDIF
1439 KVAL=KLCL
1440 loop66: DO KK =KVAL,LTOP1-1
1441 IF(TVQU(KK) .LT. TV0(KK) .AND. TVQU(KK+1).GT.TV0(KK+1)) THEN
1442 K1=KK
1443 K2=KK+1
1444 ! VAL0=TVQU(K2)-TV0(K2)+TV0(K1)-TVQU(K1)
1445 ! IF( ABS(VAL0) < 0.001 ) then
1446 ! if( dbg_level > 100 ) then
1447 ! write(message,*) kval, k1, k2, TVQU(K2), TV0(K2), TVQU(K1), TV0(K1), TVQU(K2)-TV0(K2)+TV0(K1)-TVQU(K1)
1448 ! call wrf_message( message )
1449 ! endif
1450 ! STOP 999 ! aj
1451 ! ENDIF
1452 ! VAL1=(TVQU(K2)-TV0(K2))/VAL0
1453 ! VAL2=(TV0(K1)-TVQU(K1))/VAL0
1454 ! VAL=VAL1*Z0(K1)+VAL2*Z0(K2)
1455 VAL=Z0(K2)
1457 SUM=0.0
1458 loop62: DO KS=K2,LTOP1
1459 IF(TVQU(KS).GT.TV0(KS)) THEN
1460 SUM=SUM+DZQ(KS)
1461 ELSE
1462 EXIT loop62
1463 ENDIF
1464 ENDDO loop62
1466 IF(SUM .GE. 300.) THEN
1467 ZLFC=VAL
1468 EXIT loop66
1469 ELSE
1470 CYCLE loop66
1471 ENDIF
1472 ENDIF
1474 ENDDO loop66
1476 67 CONTINUE
1478 ZLFC=MAX(ZLFC,ZPBL)
1479 ! ZLFC=ZLCL
1480 ! ZLFC=Z0(KL)
1481 !=================================================================
1482 #ifdef DENG_SHCU1D
1483 DO NK=LLC,K
1484 TUG(NK)=TLCL*(P0(NK)/PLCL)**ROCPQ
1485 ENDDO
1486 DO NK=K+1,LTOP
1487 TUG(NK)=TU(NK)
1488 ENDDO
1489 #endif
1491 IF( CLDTOP .LE. ZLCL .OR. LTOP .LT. KLCL) THEN
1492 #ifdef DENG_SHCU1D
1493 PLET = P0(1)
1494 TLET = TLCL*(PLET/PLCL)**ROVCP
1495 PLFSG(I,J) = P0(1)
1496 TLFSG(I,J) = TLCL*(PLFSG(I,J)/PLCL)**ROVCP
1497 PLDTG(I,J) = P0(1)
1498 TLDTG(I,J) = TLCL*(PLDTG(I,J)/PLCL)**ROVCP
1499 PLDBG(I,J) = P0(1)
1500 TLDBG(I,J) = TLCL*(PLDBG(I,J)/PLCL)**ROVCP
1501 #endif
1502 CLDTOP = Z0(1)
1503 LTOP = 1
1504 GO TO 425
1505 END IF
1506 IF( LET .GT. LTOP ) LET = LTOP
1507 !
1508 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FLUX AT
1509 ! THIS LEVEL...
1510 !
1511 IF(LET.EQ.LTOP)THEN
1512 UPOLD=UMF(LTOP-1)-UDR(LTOP)
1513 UPNEW=UPOLD+UER(LTOP)
1514 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
1515 IF( ABS(UPOLD) < 0.001 ) STOP 907
1516 DETLQ(LTOP)=RLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1517 DETIC(LTOP)=RICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1518 UER(LTOP)=0.
1519 UMF(LTOP)=0.
1520 GOTO 85
1521 ENDIF
1522 !
1523 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1524 !
1525 DPTT=0.
1526 DO NJ=LET+1,LTOP
1527 DPTT=DPTT+DP(NJ)
1528 ENDDO
1529 DUMFDP=UMF(LET)/DPTT
1530 !
1531 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALLOUT
1532 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND LTOP
1533 !
1534 DO NK=LET+1,LTOP
1535 UDR(NK)=DP(NK)*DUMFDP
1536 UMF(NK)=UMF(NK-1)-UDR(NK)
1537 DETLQ(NK)=RLIQ(NK)*UDR(NK)
1538 DETIC(NK)=RICE(NK)*UDR(NK)
1539 IF(NK.GE.LET+2)THEN
1540 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1541 PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
1542 PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
1543 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1544 ENDIF
1545 ENDDO
1547 85 CONTINUE
1548 !
1549 !...EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR THE
1550 ! UPDRAFT AIR...ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL...
1551 !
1552 QMIX1 = 0.0
1553 ZMIX1 = 0.0
1554 DPTHMX1 = 0.0
1555 DO NK = LLC, K
1556 DPTHMX1 = DPTHMX1 + DP(NK)
1557 QMIX1 = QMIX1 + DP(NK)*Q0(NK)
1558 ZMIX1 = ZMIX1 + DP(NK)*Z0(NK)
1559 ENDDO
1560 QMIX1 = QMIX1 / DPTHMX1
1561 ZMIX1 = ZMIX1 / DPTHMX1
1563 DO NK=1,K
1564 IF(NK.GE.LLC)THEN
1565 IF(NK.EQ.LLC)THEN
1566 UMF(NK)=VMFPBLCL*DP(NK)/DPTHMX1
1567 UER(NK)=VMFPBLCL*DP(NK)/DPTHMX1
1568 ELSEIF(NK .LE. K)THEN
1569 UER(NK)=VMFPBLCL*DP(NK)/DPTHMX1
1570 UMF(NK)=UMF(NK-1)+UER(NK)
1571 ELSE
1572 UMF(NK)=VMFPBLCL
1573 UER(NK)=0.
1574 ENDIF
1575 QU(NK)=QMIX1
1576 TU(NK)=TLCL+(Z0(NK)-ZLCL)*GDRY
1577 WU(NK)=WPBLTP
1578 ELSE
1579 TU(NK)=0.
1580 QU(NK)=0.
1581 UMF(NK)=0.
1582 WU(NK)=0.
1583 UER(NK)=0.
1584 ENDIF
1585 UDR(NK)=0.
1586 QDT(NK)=0.
1587 RLIQ(NK)=0.
1588 RICE(NK)=0.
1589 QLQOUT(NK)=0.
1590 QICOUT(NK)=0.
1591 PPTLIQ(NK)=0.
1592 PPTICE(NK)=0.
1593 DETLQ(NK)=0.
1594 DETIC(NK)=0.
1595 RATIO2(NK)=0.
1596 EE=Q0(NK)*P0(NK)/(0.622+Q0(NK))
1597 TLOG=ALOG(EE/SVP1/1000.)
1598 TDPT=(SVP2*SVPT0-SVP3*TLOG)/(SVP2-TLOG)
1599 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T0(NK)-T00))*(T0(NK)-TDPT)
1600 THTA=T0(NK)*(1.E5/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1601 THETEE(NK)=THTA*EXP((3374.6525/TSAT-2.5403)*Q0(NK)*(1.+0.81*Q0(NK)))
1602 THTES(NK)=THTA*EXP((3374.6525/T0(NK)-2.5403)*QES(NK)*(1.+0.81*QES(NK)))
1603 EQFRC(NK)=1.0
1604 ENDDO
1606 LTOP1=LTOP+1
1607 LTOPM1=LTOP-1
1608 !
1609 !...DEFINE VARIABLES ABOVE CLOUD TOP...
1610 !
1611 DO NK = LTOP1, KX
1612 UMF(NK) = 0.
1613 UDR(NK)=0.
1614 UER(NK)=0.
1615 QDT(NK)=0.
1616 RLIQ(NK)=0.
1617 RICE(NK)=0.
1618 QLQOUT(NK)=0.
1619 QICOUT(NK)=0.
1620 DETLQ(NK)=0.
1621 DETIC(NK)=0.
1622 PPTLIQ(NK)=0.
1623 PPTICE(NK)=0.
1624 IF(NK.GT.LTOP1)THEN
1625 TU(NK)=0.
1626 QU(NK)=0.
1627 WU(NK)=0.
1628 ENDIF
1629 THTA0(NK)=0.
1630 THTAU(NK)=0.
1631 EMS(NK)=DP(NK)*DXSQ/G
1632 EMSD(NK)=1./EMS(NK)
1633 TG(NK)=T0(NK)
1634 QQG(NK)=Q0(NK)
1635 QCG(NK)=0.
1636 DCA(NK)=0.
1637 DMS(NK)=0.
1638 DTFM(NK)=0.
1639 FXM(NK)=0.0
1640 OMG(NK)=0.
1641 ENDDO
1642 OMG(KXP1)=0.
1644 #ifdef DENG_SHCU1D
1645 !
1646 ! 1-D GRAPHICS
1647 !
1648 DO NK=1,KL
1649 wutemp(nk)=wu(NK)
1650 ENDDO
1651 #endif
1653 P165=P0(KLCL)-1.65E4
1654 PPTMLT=0.
1655 DO NK=1,LTOP
1656 EMS(NK)=DP(NK)*DXSQ/G
1657 EMSD(NK)=1./EMS(NK)
1658 DTFM(NK)=0.
1659 !
1660 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCHEME.
1661 !
1662 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1663 THTAU(NK)=TU(NK)*EXN(NK)
1664 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1665 THTA0(NK)=T0(NK)*EXN(NK)
1667 IF(P0(NK).GT.P165) LVF=NK
1668 PPTMLT=PPTMLT+PPTICE(NK)
1670 OMG(NK)=0.
1671 ENDDO
1673 CLVF=0.
1674 DO NK=KLCL,LVF+1
1675 CLVF=CLVF+PPTLIQ(NK)+PPTICE(NK)
1676 ENDDO
1677 USR=UMF(LVF+1)*QU(LVF+1)+CLVF
1678 USR=MIN(USR,TRPPT)
1679 ! WRITE(28,1025)KLCL,ZLCL,LTOP,P0(LTOP),IFLAG, &
1680 ! TMIX-T00,PMIX,QMIX,ABE
1681 ! WRITE(28,1030)P0(LET)/100.,P0(LTOP)/100.,VMFPBLCL,XTIME/60. &
1682 ! ,PLCL/100.,WPBLTP,CLDDPTH
1683 !1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
1684 ! ' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
1685 ! I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
1686 ! ' CAPE=',0PF7.1)
1687 !1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
1688 ! E12.3,' XTIME =',F10.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
1689 ! F8.1)
1690 !
1691 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1692 !...AND MIDTROPOSPHERE IS USED.
1693 !
1694 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1695 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1696 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1697 VCONV = .5*(WSPD(KLCL)+WSPD(L5))
1698 TIMEC = DX/VCONV
1699 TADVEC = TIMEC
1700 TCMIN = 1800.
1701 TCMAX = 3600.
1702 TIMEC = MAX(TCMIN,TIMEC)
1703 TIMEC = MIN(TCMAX,TIMEC)
1704 NIC = NINT(TIMEC/(.5*DT2))
1705 TIMEC = FLOAT(NIC)*.5*DT2
1706 !
1707 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1708 !
1709 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1710 SHSIGN=1.
1711 ELSE
1712 SHSIGN=-1.
1713 ENDIF
1714 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))*(V0(LTOP)-V0(KLCL))
1715 IF(LTOP .LE. LCL) GO TO 425
1716 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1717 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1718 PEFMX=0.9
1719 PEFMN=0.2
1720 PEF=MAX(PEF,PEFMN)
1721 PEF=MIN(PEF,PEFMX)
1722 !
1723 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1724 !
1725 CBH=(ZLCL-Z0(1))*3.281E-3
1726 IF(CBH.LT.3.) THEN
1727 RCBH=.02
1728 ELSE
1729 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
1730 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1731 ENDIF
1732 IF (CBH.GT.25) RCBH=2.4
1733 PEFCBH=1./(1.+RCBH)
1734 CBHMX=0.9
1735 PEFCBH=MIN(PEFCBH,CBHMX)
1736 !
1737 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1738 ! DOWNDRAFT IS TURNED OFF BY RATE=0.0
1739 !
1740 PEFF=.5*(PEF+PEFCBH)
1741 IF(CLDDPTHB(I,J) .LT. 4.E3) THEN
1742 USR=TRPPT
1743 PEFF=0.99999
1744 ENDIF
1745 ! IF ( wrf_dm_on_monitor()) THEN
1746 ! CALL get_wrf_debug_level( dbg_level )
1747 ! IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
1748 ! WRITE(message,1035)PEF,PEFCBH,LC,LET,VWS
1749 !1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'VWS=',F5.2)
1750 ! CALL wrf_message( message )
1751 ! ENDIF
1752 ! ENDIF
1753 !
1754 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
1755 ! BEGINING OF DOWNDRAFT CACULATION
1756 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
1757 !
1758 !...LET DOWNDRAFT ORIGINATE AT THE LEVEL OF MINIMUM SATURATION EQUIVALENT
1759 !...POTENTIAL TEMPERATURE (SEQT) IN THE CLOUD LAYER, EXTEND DOWNWARD TO THE
1760 !...SURFACE, OR TO THE LAYER BELOW CLOUD BASE AT WHICH ENVIR SEQT IS LESS
1761 !...THAN MIN SEQT IN THE CLOUD LAYER...LET DOWNDRAFT DETRAIN OVER A LAYER
1762 !...OF SPECIFIED PRESSURE-DEPTH (DPDD)...
1763 !
1764 TDER=0.
1765 KSTART=MAX0(KKPBL,KLCL)
1766 THTMIN=THTES(KSTART+1)
1767 KMIN=KSTART+1
1768 DO NK=KSTART+2,LTOP-1
1769 THTMIN=MIN(THTMIN,THTES(NK))
1770 IF(THTMIN.EQ.THTES(NK))KMIN=NK
1771 ENDDO
1772 LFS=KMIN
1773 IF(LFS .GE. LTOP) LFS=LTOP-1
1774 IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT(P0(LFS),T0(LFS),Q0(LFS), &
1775 THETEE(LFS),0.,RL, &
1776 SVP1,SVP2,SVPT0,SVP3,I,J)
1777 ! IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT2(P0(LFS),T0(LFS),Q0(LFS), &
1778 ! THETEE(LFS),ALIQ,BLIQ,CLIQ,DLIQ)
1780 IF( ABS(THETEE(LFS)-THETEU(LFS)) < 1.E-3 ) THEN ! DENG 20140110 fix to avoid dividing by 0
1781 EQFRC(LFS) = 1.0
1782 ELSE
1783 EQFRC(LFS)=(THTES(LFS)-THETEU(LFS))/(THETEE(LFS)-THETEU(LFS))
1784 ENDIF
1785 EQFRC(LFS)=MAX(EQFRC(LFS),0.)
1786 EQFRC(LFS)=MIN(EQFRC(LFS),1.)
1787 THETED(LFS)=THTES(LFS)
1788 !
1789 !...ESTIMATE THE EFFECT OF MELTING ON THE DOWNDRAFT...
1790 !
1791 IF(ML.GT.0)THEN
1792 DTMLTD=0.5*(QU(KLCL)-QU(LTOP))*RLF/CP
1793 ELSE
1794 DTMLTD=0.
1795 ENDIF
1796 TZ(LFS)=T0(LFS)-DTMLTD
1797 ES=1.E3*SVP1*EXP(SVP2*(TZ(LFS)-SVPT0)/(TZ(LFS)-SVP3))
1798 QS = 0.622*ES/(P0(LFS)-ES)
1799 QD(LFS)=EQFRC(LFS)*Q0(LFS)+(1.-EQFRC(LFS))*QU(LFS)
1800 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QD(LFS)))
1801 IF( QD(LFS) .GE. QS ) THEN
1802 THETED(LFS)=THTAD(LFS)* &
1803 EXP((3374.6525/TZ(LFS)-2.5403)*QS*(1.+0.81*QS))
1804 ELSE
1805 CALL ENVIRTHT(P0(LFS),TZ(LFS),QD(LFS), &
1806 THETED(LFS),0.,RL, &
1807 SVP1,SVP2,SVPT0,SVP3,I,J)
1808 ! CALL ENVIRTHT2(P0(LFS),TZ(LFS),QD(LFS), &
1809 ! THETED(LFS),ALIQ,BLIQ,CLIQ,DLIQ)
1810 ENDIF
1811 !
1812 !...DETERMINE THE LOWEST LEVEL TO WHICH DOWNDRAFT CAN PENETRATE (LDB)...
1813 !
1814 loop107: DO NK=1,LFS
1815 ND=LFS-NK
1816 IF(THETED(LFS).GT.THTES(ND) .OR. ND.EQ.1)THEN
1817 LDB=ND
1818 EXIT loop107
1819 ENDIF
1820 ENDDO loop107
1821 !
1822 !...ASSUME ALL OF DOWNDRAFT DETRAINMENT OCCURS AT ITS BASE, SO THE TOP OF THE
1823 !...DETRAINMENT LAYER (LDT) IS THE SAME AS LDB...
1824 !
1825 DPDD= 20.0E2
1826 DPT=0.
1827 loop115: DO NK=LDB,LFS
1828 DPT=DPT+DP(NK)
1829 IF(DPT.GT.DPDD)THEN
1830 LDT=NK
1831 FRC=(DPDD+DP(NK)-DPT)/DP(NK)
1832 EXIT loop115
1833 ENDIF
1834 IF(NK.EQ.LFS-1)THEN
1835 LDT=NK
1836 FRC=1.
1837 DPDD=DPT
1838 EXIT loop115
1839 ENDIF
1840 ENDDO loop115
1842 IF(LDB.EQ.LFS) THEN
1843 TDER=0.
1844 GOTO 141
1845 ENDIF
1846 !
1847 !...MAKE A FIRST GUESS AT INITIAL DOWNDRAFT MASS FLUX, DETERMINE ITS VERTICAL
1848 !...PROFILE...
1849 !
1850 TVD(LFS)=T0(LFS)*(1.+0.608*QES(LFS))
1851 RDD=P0(LFS)/(R*TVD(LFS))
1852 A1=(1.-PEFF)*AU0
1853 DMF(LFS)=-A1*RDD
1854 DER(LFS)=EQFRC(LFS)*DMF(LFS)
1855 DDR(LFS)=0.
1856 DO ND=LFS-1,LDB,-1
1857 ND1=ND+1
1858 IF(ND.LE.LDT)THEN
1859 DER(ND)=0.
1860 DDR(ND)=-DMF(LDT+1)*DP(ND)*FRC/DPDD
1861 DMF(ND)=DMF(ND1)+DDR(ND)
1862 FRC=1.
1863 THETED(ND)=THETED(ND1)
1864 QD(ND)=QD(ND1)
1865 ELSE
1866 DER(ND)=DMF(LFS)*0.03*DP(ND)/RAD
1867 DDR(ND)=0.
1868 DMF(ND)=DMF(ND1)+DER(ND)
1869 IF(RATIO2(ND).GT.0.)CALL ENVIRTHT(P0(ND),T0(ND),Q0(ND), &
1870 THETEE(ND),0.,RL, &
1871 SVP1,SVP2,SVPT0,SVP3,I,J)
1872 ! IF(RATIO2(ND).GT.0.)CALL ENVIRTHT2(P0(ND),T0(ND),Q0(ND), &
1873 ! THETEE(ND),ALIQ,BLIQ,CLIQ,DLIQ)
1874 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1875 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
1876 ENDIF
1877 ENDDO
1878 !
1879 !...DETERMINE TOTAL DOWNDRAFT EVAPORATION RATE (TDER)
1880 !
1881 TDER=0.
1882 DO ND=LDB,LDT
1883 TZ(ND)=TPDDBG(P0(ND),THETED(LDT),T0(ND),QS,QD(ND),1.0, &
1884 XLV0,XLV1, &
1885 SVP1,SVP2,SVPT0,SVP3)
1886 ! CALL TPMIX2DD(p0(ND),theted(LDT),tz(ND),qs,ktau,i,j,nd)
1888 ES=1.E3*SVP1*EXP(SVP2*(TZ(ND)-SVPT0)/(TZ(ND)-SVP3))
1889 QS=0.622*ES/(P0(ND)-ES)
1890 DQSSDT=SVP2*(SVPT0-SVP3)/((TZ(ND)-SVP3)*(TZ(ND)-SVP3))
1891 RL=XLV0-XLV1*TZ(ND)
1892 DTMP=RL*QS*(1.-RHBC)/(CP+RL*RHBC*QS*DQSSDT)
1893 T1RH=TZ(ND)+DTMP
1894 ES=RHBC*1.E3*SVP1*EXP(SVP2*(T1RH-SVPT0)/(T1RH-SVP3))
1895 QSRH=0.622*ES/(P0(ND)-ES)
1896 !
1897 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1898 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1899 !
1900 IF(QSRH.LT.QD(ND))THEN
1901 QSRH=QD(ND)
1902 T1RH=TZ(ND)
1903 ENDIF
1904 TZ(ND)=T1RH
1905 QS=QSRH
1906 TDER=TDER + (QS-QD(ND))*DDR(ND)
1907 QD(ND)=QS
1908 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1909 ENDDO
1910 !
1911 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1912 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1913 !
1914 141 CONTINUE
1915 IF(TDER.LT.1.)THEN
1916 PPTFLX=TRPPT
1917 CPR=TRPPT
1918 TDER=0.
1919 CNDTNF=0.
1920 UPDINC=1.
1921 LDB=LFS
1923 DO NDK=1,KX
1924 DMS(NDK)=0.
1925 DMF(NDK)=0.
1926 DER(NDK)=0.
1927 DDR(NDK)=0.
1928 THTAD(NDK)=0.
1929 WD(NDK)=0.
1930 TZ(NDK)=0.
1931 QD(NDK)=0.
1932 ENDDO
1934 AINCM2=100.
1935 DMFMIN=0.
1936 GOTO 165
1937 ENDIF
1938 !
1939 !...ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT IS
1940 !...CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP...
1941 !
1942 DEVDMF=TDER/DMF(LFS)
1943 PPR=0.
1944 PPTFLX=PEFF*USR
1945 RCED=TRPPT-PPTFLX
1946 !
1947 !...PPR IS THE TOTAL AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
1948 !...FROM CLOUD BASE TO THE LFS...UPDRAFT MASS FLUX WILL BE INCREASED UP TO
1949 !...THE LFS TO ACCOUNT FOR UPDRAFT AIR MIXING WITH ENVIRONMENTAL AIR TO MAKE
1950 !...THE UPDRAFT, SO PPR WILL INCREASE PROPORTIONATELY...
1951 !
1952 DO NM=KLCL,LFS
1953 PPR=PPR+PPTLIQ(NM)+PPTICE(NM)
1954 ENDDO
1955 IF(LFS.GE.KLCL)THEN
1956 DPPTDF=(1.-PEFF)*PPR*(1.-EQFRC(LFS))/UMF(LFS)
1957 ELSE
1958 DPPTDF=0.
1959 ENDIF
1960 !
1961 !...CNDTNF IS THE AMOUNT OF CONDENSATE TRANSFERRED ALONG WITH UPDRAFT MASS TO
1962 !...THE DOWNDRAFT AT THE LFS...
1963 !
1964 CNDTNF=(RLIQ(LFS)+RICE(LFS))*(1.-EQFRC(LFS))
1965 DMFLFS=RCED/(DEVDMF+DPPTDF+CNDTNF)
1966 IF(DMFLFS.GT.0.)THEN
1967 TDER=0.
1968 GOTO 141
1969 ENDIF
1970 !...DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT MASS FLUX
1971 !...TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP, UPDINC IS THE FACTOR BY WHICH
1972 !...TO INCREASE THE UPDRAFT MASS FLUX BELOW THE LFS TO ACCOUNT FOR THE TRANSFER
1973 !...OF MASS FROM UPDRAFT TO DOWNDRAFT...
1974 !
1975 DDINC=DMFLFS/DMF(LFS)
1976 IF(LFS.GE.KLCL)THEN
1977 UPDINC=(UMF(LFS)-(1.-EQFRC(LFS))*DMFLFS)/UMF(LFS)
1978 ELSE
1979 UPDINC=1.
1980 ENDIF
1981 DO NK=LDB,LFS
1982 DMF(NK)=DMF(NK)*DDINC
1983 DER(NK)=DER(NK)*DDINC
1984 DDR(NK)=DDR(NK)*DDINC
1985 ENDDO
1986 CPR=TRPPT+PPR*(UPDINC-1.)
1987 PPTFLX = PPTFLX+PEFF*PPR*(UPDINC-1.)
1988 TDER=TDER*DDINC
1989 !
1990 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AND ICE
1991 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATED MASS
1992 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1993 !
1994 DO NK=LC,LFS
1995 UMF(NK)=UMF(NK)*UPDINC
1996 UDR(NK)=UDR(NK)*UPDINC
1997 UER(NK)=UER(NK)*UPDINC
1998 PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1999 PPTICE(NK)=PPTICE(NK)*UPDINC
2000 DETLQ(NK)=DETLQ(NK)*UPDINC
2001 DETIC(NK)=DETIC(NK)*UPDINC
2002 ENDDO
2003 IF(LDB.GT.1)THEN
2004 DO NK=1,LDB-1
2005 DMF(NK)=0.
2006 DER(NK)=0.
2007 DDR(NK)=0.
2008 WD(NK)=0.
2009 TZ(NK)=0.
2010 QD(NK)=0.
2011 THTAD(NK)=0.
2012 ENDDO
2013 ENDIF
2015 DO NK=LFS+1,KX
2016 DMF(NK)=0.
2017 DER(NK)=0.
2018 DDR(NK)=0.
2019 WD(NK)=0.
2020 TZ(NK)=0.
2021 QD(NK)=0.
2022 THTAD(NK)=0.
2023 ENDDO
2025 DO NK=LDT+1,LFS-1
2026 TZ(NK)=0.
2027 QD(NK)=0.
2028 ENDDO
2029 !
2030 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFLOW
2031 ! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILABLE
2032 ! IN THAT LAYER INITIALLY...
2033 !
2034 165 CONTINUE
2035 AINCMX=1000.
2036 LMAX=MAX0(KLCL,LFS)
2037 DO NK=LLC,LMAX
2038 IF((UER(NK)-DER(NK)).GT.0.) &
2039 AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
2040 AINCMX=MIN(AINCMX,AINCM1)
2041 ENDDO
2043 AINC=1.
2044 IF(AINCMX.LT.AINC)AINC=AINCMX
2045 !
2046 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT...THEY WILL BE
2047 ! ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION CLOSURE...
2048 !
2049 NCOUNT=0
2050 PPTMLT=0.
2051 TDER2=TDER
2052 PPTFL2=PPTFLX
2053 DO NK=1,LTOP
2054 DETLQ2(NK)=DETLQ(NK)
2055 DETIC2(NK)=DETIC(NK)
2056 UDR2(NK)=UDR(NK)
2057 UER2(NK)=UER(NK)
2058 DDR2(NK)=DDR(NK)
2059 DER2(NK)=DER(NK)
2060 UMF2(NK)=UMF(NK)
2061 DMF2(NK)=DMF(NK)
2062 DMS(NK)=0.
2063 IF(NK.GT.ML .AND. NK.LT.LTOP)THEN
2064 PPTMLT=PPTMLT+PPTICE(NK+1)
2065 ENDIF
2066 ENDDO
2068 PPTML2=PPTMLT
2069 FABE=1.
2070 STAB=0.95
2071 AINC2=AINC
2072 IF(AINC/AINCMX.GT.0.999)GOTO 255
2073 ISTOP=0
2075 KKPBL=K
2076 !
2077 !...FIND THE MAXIMUM TKE VALUE BETWEEN LLC AND KLCL... AND USE IT TO
2078 ! CALCULATE AINC2
2079 !
2080 TKEMAX=0.
2081 DO KK=LLC,KLCL
2082 TKEMAX=AMAX1(TKEMAX,TKE0(KK))
2083 ENDDO
2084 TKEMAX=AMIN1(TKEMAX,10.)
2085 TKEMAX=AMAX1(TKEMAX,1.0)
2086 DSOURCE=Z0(K)-Z0(LLC)+0.5*(DZQ(K)+DZQ(LLC))
2087 TIMECS=DSOURCE/WPBLTP*45.0
2088 AINC2=TKEMAX*DPTHMX1*DXSQ/(VMFLCL*G*TIMECS)
2089 #ifdef DENG_SHCU1D
2090 IF( i == 1 .AND. j == 1 ) write(808,'(i7,f10.3, 2i4)') ktau, xtime/60.,LLC, KLCL
2091 IF( i == 1 .AND. j == 1 ) write(809,'(i7,4f10.3)') ktau, xtime/60.,TKEMAX, &
2092 DSOURCE, AINC2
2093 #endif
2094 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
2095 ! BEGINING OF CLOSURE ASSUMPTION SELECTIONS
2096 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
2097 !
2098 IF(CAPEI(I,J) .LT. 100.0) THEN ! averaging KF number of clouds
2099 NT = NINT(15.0*60.0/(0.5*DT))-1 ! for 15 min.
2100 NT = MIN(NT, 100)
2101 ELSE
2102 NT = 0
2103 ENDIF
2105 175 NCOUNT=NCOUNT+1
2106 !
2107 !
2108 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
2109 !...SATISFY MASS CONTINUITY...
2110 !
2111 DTT=TIMEC
2112 DO NK=1,LTOP
2113 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
2114 IF(NK.GT.1)THEN
2115 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
2116 DTT1 = 0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
2117 DTT=MIN(DTT,DTT1)
2118 ENDIF
2119 ENDDO
2120 !
2121 ! NOTE:
2122 ! KFLAG=2: CLOUD TOP LOWER THAN LFC
2123 ! KFLAG=3: CLOUD TOP HIGHER THAN LFC BUT LOWER THAN 4 KM
2124 ! KFLAG=4: CLOUD TOP HIGHER THAN 4 KM
2125 !
2126 IF(KFLAG .EQ. 2) GO TO 176
2128 IF((CLDDPTH+ZLCL) .LE. ZLFC) THEN
2129 AINC=AINC2
2130 KFLAG=2
2131 GO TO 255
2132 ENDIF
2134 176 CONTINUE
2136 DO NK=1,LTOP
2137 THPA(NK)=THTA0(NK)
2138 QPA(NK)=Q0(NK)
2139 QCPA(NK)= QC0(NK)
2140 DCA(NK) = DCA0(NK)
2141 NSTEP=NINT(TIMEC/DTT+1)
2142 DTIME=TIMEC/FLOAT(NSTEP)
2143 FXM(NK)=OMG(NK)*DXSQ/G
2144 #ifdef DENG_SHCU1D
2145 wsubsub(nk)=-omg(nk)/rhoe(nk)/g
2146 #endif
2147 ENDDO
2149 loop495: DO NTC=1,NSTEP
2150 !
2151 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED ON THE
2152 !...SIGN OF OMEGA...
2153 !
2154 DO NK=1,LTOP
2155 THFXIN(NK)=0.
2156 THFXOUT(NK)=0.
2157 QFXIN(NK)=0.
2158 QFXOUT(NK)=0.
2159 QCFXIN(NK)=0.
2160 QCFXOUT(NK)=0.
2161 DCAFXIN(NK)=0.0
2162 DCAFXOUT(NK)=0.0
2163 ENDDO
2165 DO NK=2,LTOP
2166 !
2167 ! THIS PART IS ALSO MODIFIED TO INCLUDE THE SUBSIDANCE EFFECT
2168 ! ON DEAD CLOUD. THE FLUXES INTO THE LAYER ARE MODIFIED
2169 ! TO CONSIDER THE EFFECT ON QC,THETA,QV, GD ID ADIABATIC LAPSE
2170 ! RATE AND GM IS MOIST ADIABATIC LAPSE RATE
2171 !
2172 IF(OMG(NK).LE.0.)THEN
2173 THFXIN(NK)=-FXM(NK)*THPA(NK-1)
2174 QFXIN(NK)=-FXM(NK)*QPA(NK-1)
2175 QCFXIN(NK)=-FXM(NK)*QCPA(NK-1)
2176 DCAFXIN(NK) = -FXM(NK)*DCA(NK-1)
2177 THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
2178 QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
2179 QCFXOUT(NK-1)=QCFXOUT(NK-1)+QCFXIN(NK)
2180 DCAFXOUT(NK-1)=DCAFXOUT(NK-1)+DCAFXIN(NK)
2181 ELSE
2182 GD=G/CP
2183 CPM=CP*(1.+0.887*Q0(NK))
2184 IF(T0(NK) .GT. TO) THEN
2185 RL=XLV0-XLV1*T0(NK)
2186 DQSSDT=QES(NK)*SVP2*(SVPT0-SVP3) &
2187 /((T0(NK)-SVP3)*(T0(NK)-SVP3))
2188 ELSE
2189 RL=XLS
2190 DQSSDT=QES(NK)*6.15E3/(T0(NK)*T0(NK))
2191 ENDIF
2192 GM=GD/(1.0+RL/CPM*DQSSDT)
2193 DTEMPDT=(GD-GM)*OMG(NK)/RHOE(NK)/G
2194 IF((QCPA(NK)-DTEMPDT*CP/RL*DTIME) .LE. 0.0) &
2195 DTEMPDT=QCPA(NK)/(CP/RL*DTIME)
2196 DTHETADT=DTEMPDT*(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
2197 THFXOUT(NK)=FXM(NK)*(THPA(NK)-DTHETADT*DTIME)
2198 QFXOUT(NK)=FXM(NK)*(QPA(NK) +DTEMPDT*CP/RL*DTIME)
2199 QCFXOUT(NK)=FXM(NK)*(QCPA(NK)-DTEMPDT*CP/RL*DTIME)
2200 DCAFXOUT(NK)=FXM(NK)*DCA(NK)
2201 THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
2202 QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
2203 QCFXIN(NK-1)=QCFXIN(NK-1)+QCFXOUT(NK)
2204 DCAFXIN(NK-1)=DCAFXIN(NK-1)+DCAFXOUT(NK)
2205 ENDIF
2206 ENDDO
2207 !
2208 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
2209 !
2210 DO NK=LTOP,1,-1
2211 IF( NK .GT. KKPBL .OR. NK .LT. LLC) THEN
2212 THPA(NK)=THPA(NK)+ &
2213 (THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)*THTAD(NK)- &
2214 THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))*DTIME*EMSD(NK)
2215 QPA(NK)=QPA(NK)+ &
2216 (QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK) &
2217 -QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
2218 ELSE
2219 THPA(NK)=THPA(NK)+ &
2220 (THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)*THTAD(NK)- &
2221 THFXOUT(NK)-UER(NK)*THMIX+DER(NK)*THTA0(NK))* &
2222 DTIME*EMSD(NK)
2223 QPA(NK)=QPA(NK)+ &
2224 (QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK) &
2225 -QFXOUT(NK)-UER(NK)*QMIX+DER(NK)*Q0(NK))* &
2226 DTIME*EMSD(NK)
2227 END IF
2228 QCPA(NK)=QCPA(NK)+ &
2229 (QCFXIN(NK)+DETLQ(NK)+DETIC(NK)-QCFXOUT(NK))*DTIME*EMSD(NK)
2230 !
2231 ! ADDING DEAD CLOUD AREA DCA WHICH WILL BE USED IN LOOP 320 TO CALCULATED
2232 ! THE TENDENCY DDCADT IN LOOP 320, THE TENDENCY (DQCDCDT) FOR DEAD CLOUD
2233 ! LIQUID WATER CONTENT QCDC IS SOLVE AS A RESIDUAL OF DQLDT-DDCADT
2234 !
2235 DCA(NK)=DCA(NK)+(DCAFXIN(NK)-DCAFXOUT(NK)+UDR(NK))*DTIME*EMSD(NK)
2236 ENDDO
2238 ENDDO loop495
2240 DO NK=LTOP,1,-1
2241 THTAG(NK)=THPA(NK)
2242 QQG(NK)=QPA(NK)
2243 QCG(NK)=QCPA(NK)
2244 ENDDO
2245 !
2246 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORROW
2247 !...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO; WITH
2248 !...INDEX NK-1, MINIMUM VALUE FOR NK=2...
2249 !
2250 DO NK = 2,LTOP
2251 IF(QQG(NK).LT.0.)THEN
2252 TMA = QQG(NK+1)*EMS(NK+1)
2253 TMB = QQG(NK-1)*EMS(NK-1)
2254 TMM = (QQG(NK)-1.E-9)*EMS(NK )
2255 IF(TMA .NE. 0.0 .AND. TMB .NE. 0.0)THEN
2256 BCOEFF = -TMM/((TMA*TMA)/TMB+TMB)
2257 ACOEFF = BCOEFF*TMA/TMB
2258 ELSE IF(TMA .EQ. 0.0 .AND. TMB .EQ. 0.0) THEN
2259 BCOEFF = 0.0
2260 ACOEFF = 0.0
2261 ELSE IF(TMA .EQ. 0.0) THEN
2262 BCOEFF = -TMM/TMB
2263 ACOEFF = 0.0
2264 ELSE IF(TMB .EQ. 0.0) THEN
2265 BCOEFF = 0.0
2266 ACOEFF = -TMM/TMA
2267 ENDIF
2268 TMB = TMB*(1.-BCOEFF)
2269 TMA = TMA*(1.-ACOEFF)
2270 QQG(NK) = 1.E-9
2271 QQG(NK+1) = TMA*EMSD(NK+1)
2272 QQG(NK-1) = TMB*EMSD(NK-1)
2273 ENDIF
2274 ENDDO
2275 TOPOMG = (UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
2276 IF(ABS(TOPOMG-OMG(LTOP)).GT. 1.E-3)THEN
2277 IF ( wrf_dm_on_monitor()) THEN
2278 CALL get_wrf_debug_level( dbg_level )
2279 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2280 WRITE(message,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;', &
2281 ' TOPOMG, OMG =',TOPOMG,OMG(LTOP)
2282 CALL wrf_message( message )
2283 ENDIF
2284 ENDIF
2285 ISTOP=1
2286 GOTO 265
2287 ENDIF
2288 !
2289 !...CONVERT THETA TO T, FREEZE ALL SUPERCOOLED DETRAINED LIQUID WATER BECAUSE
2290 !...SUPERCOOLED WATER NOT ALLOWED IN EXMOIS...FREEZE SUPERCOOLED PRECIP IF
2291 !...FEEDING IT BACK TO EXMOIS...
2292 !
2293 DO NK=1,LTOP
2294 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QQG(NK)))
2295 TG(NK)=THTAG(NK)/EXN(NK)
2296 IF(NK.GT.ML)DTFM(NK)=DETLQ(NK)*RLF*EMSD(NK)/CP
2297 TG(NK)=TG(NK)+DTFM(NK)*timec
2298 TVG(NK)=TG(NK)*(1.+0.608*QQG(NK))
2299 ENDDO
2300 !
2301 !...ALLOW FROZEN PRECIP TO MELT OVER A LAYER 200mb DEEP IF PRECIP IS NOT FED
2302 !...BACK TO EXMOIS...
2303 !
2304 DTMLTE=0.
2305 TDP=0.
2306 loop231: DO K = 1,ML
2307 NK = ML-K+1
2308 TDP=TDP+DP(NK)
2309 IF(TDP.LT.2.E4)THEN
2310 DTFM(NK)=DTMLTE
2311 TG(NK)=TG(NK)+DTMLTE*TIMEC
2312 ELSE
2313 DTFM(NK)=DTMLTE*(2.E4+DP(NK)-TDP)/DP(NK)
2314 TG(NK)=TG(NK)+DTMLTE*TIMEC*(2.E4+DP(NK)-TDP)/DP(NK)
2315 EXIT loop231
2316 ENDIF
2317 ENDDO loop231
2319 IF(KFLAG .EQ. 2 .OR. KFLAG .EQ. 3) GO TO 265
2320 !
2321 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY.
2322 !
2323 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
2324 !
2325 THMIX2=0.
2326 QMIX2=0.
2327 PMIX2=0.
2328 DO NK = 1,LM
2329 ROCPQ=0.2854*(1.-0.28*QQG(NK))
2330 THMIX2=THMIX2+DP(NK)*TG(NK)*(P00/P0(NK))**ROCPQ
2331 QMIX2=QMIX2+DP(NK)*QQG(NK)
2332 PMIX2=PMIX2+DP(NK)*P0(NK)
2333 ENDDO
2334 THMIX2=THMIX2/DPTHMX
2335 QMIX2=QMIX2/DPTHMX
2336 PMIX2=PMIX2/DPTHMX
2338 ROCPQ=0.2854*(1.-0.28*QMIX2)
2339 TMIX2=THMIX2*(PMIX2/P00)**ROCPQ
2340 ES=1.E3*SVP1*EXP(SVP2*(TMIX2-SVPT0)/(TMIX2-SVP3))
2341 QS=0.622*ES/(PMIX2-ES)
2342 IF(QMIX2.GT.QS)THEN
2343 RL=XLV0-XLV1*TMIX2
2344 CPM=CP*(1.+0.887*QMIX2)
2345 DQSSDT=QS*SVP2*(SVPT0-SVP3)/((TMIX2-SVP3)*(TMIX2-SVP3))
2346 DQ=(QMIX2-QS)/(1.+RL*DQSSDT/CPM)
2347 TMIX2=TMIX2+RL/CP*DQ
2348 QMIX2=QMIX2-DQ
2349 ROCPQ=0.2854*(1.-0.28*QMIX2)
2350 THMIX2=TMIX2*(P00/PMIX2)**ROCPQ
2351 TTLCL=TMIX2
2352 PPLCL=PMIX2
2353 ELSE
2354 QMIN=0.
2355 QMIX2=MAX(QMIX2,QMIN)
2356 EMIX2=QMIX2*PMIX2/(0.622+QMIX2)
2357 TLOG=LOG(EMIX2/SVP1/1000.)
2358 TDPT=(SVP2*SVPT0-SVP3*TLOG)/(SVP2-TLOG)
2359 TTLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX2-T00))* &
2360 (TMIX2-TDPT)
2361 TTLCL=MIN(TTLCL,TMIX2)
2362 CPORQ=1./ROCPQ
2363 PPLCL=P00*(TTLCL/THMIX2)**CPORQ
2364 ENDIF
2365 TVLCL=TTLCL*(1.+0.608*QMIX2)
2367 loop235: DO NK=LC,KL
2368 KLCL=NK
2369 IF(PPLCL.GE.P0(NK)) EXIT loop235
2370 ENDDO loop235
2372 K=KLCL-1
2374 DLP=LOG(PPLCL/P0(K))/LOG(P0(KLCL)/P0(K))
2375 !
2376 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
2377 !
2378 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
2379 QENV=QQG(K)+(QQG(KLCL)-QQG(K))*DLP
2380 TVEN=TENV*(1.+0.608*QENV)
2381 TVBAR=0.5*(TVG(K)+TVEN)
2382 ZLCL1=Z0(K)+R*TVBAR*LOG(P0(K)/PPLCL)/G
2383 TVAVG=0.5*(TVEN+TG(KLCL)*(1.+0.608*QQG(KLCL)))
2384 PPLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL1))
2385 THETEU(K)=TMIX2*(1.E5/PMIX2)**(0.2854*(1.-0.28*QMIX2))* &
2386 EXP((3374.6525/TTLCL-2.5403)*QMIX2* &
2387 (1.+0.81*QMIX2))
2388 ES=1.E3*SVP1*EXP(SVP2*(TENV-SVPT0)/(TENV-SVP3))
2389 QESE=0.622*ES/(PPLCL-ES)
2390 THTESG(K)=TENV*(1.E5/PPLCL)**(0.2854*(1.-0.28*QESE))* &
2391 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
2392 !
2393 !...COMPUTE ADJUSTED ABE(ABEG).
2394 !
2395 ABEG=0.
2396 THTUDL=THETEU(K)
2397 THATA = TMIX2*(1.E5/PMIX2)**0.286
2398 THTFC = THATA*EXP((XLV0-XLV1*TTLCL)*QMIX2/(CP*TTLCL))
2399 DO NK=K,LTOPM1
2400 NK1=NK+1
2401 ES=1.E3*SVP1*EXP(SVP2*(TG(NK1)-SVPT0)/(TG(NK1)-SVP3))
2402 QESE=0.622*ES/(P0(NK1)-ES)
2403 THTESG(NK1)=TG(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QESE))* &
2404 EXP((3374.6525/TG(NK1)-2.5403)*QESE*(1.+0.81*QESE))
2405 IF(NK.EQ.K)THEN
2406 DZZ=Z0(KLCL)-ZLCL1
2407 ELSE
2408 DZZ=DZA(NK)
2409 ENDIF
2410 BE=((2.*THTUDL)/(THTESG(NK1)+THTESG(NK))-1.)*DZZ
2411 IF(BE.GT.0.)ABEG=ABEG+BE*G
2412 ENDDO
2413 !
2414 !...FIND OUT HOW MUCH CAPE HAS BEEN REMOVED...
2415 !
2416 C1=999.9
2417 C2=999.9
2418 C3=999.9
2419 AINC3=999.0
2420 IF(ABEG.EQ.0)THEN
2421 AINC=AINC*0.5
2422 IF(NCOUNT .LT. 10) THEN
2423 GOTO 255
2424 ELSE
2425 AINC4=AINC
2426 IF(CLDDPTH .GE. 4.E3) THEN
2427 KFLAG=4
2428 GOTO 265
2429 ELSE IF((CLDDPTH) .LT. 4.E3 .AND. &
2430 CLDDPTH+ZLCL .GT. ZLFC) THEN
2431 C1=CLDDPTH+ZLCL-ZLFC
2432 C2=4000.0-ZLFC+ZLCL
2433 IF( ABS(C2) < 0.001 ) STOP 901
2434 C3=C1/C2
2435 VAL = 0.0
2436 DO II = 1, NT
2437 VAL = VAL + AINCKFSA(I,II,J)
2438 ENDDO
2439 AINC=(AINC+VAL)/FLOAT(NT+1)
2440 AINCKFI(i,j)=AINC
2441 AINC3=C3*AINC+(1-C3)*AINC2
2442 AINC=AINC3
2443 KFLAG=3
2444 GO TO 255
2445 ENDIF
2446 ENDIF
2447 ENDIF
2448 DABE=MAX(ABE-ABEG,0.1*ABE)
2449 FABE=ABEG/(ABE+1.E-8)
2450 IF(AINC/AINCMX.GT.0.999 .AND. FABE.GT.1.05-STAB) THEN
2451 AINC4=AINC
2452 IF(CLDDPTH .GE. 4.E3) THEN
2453 KFLAG=4
2454 GOTO 265
2455 ELSE IF((CLDDPTH) .LT. 4.E3 .AND. &
2456 CLDDPTH+ZLCL .GT. ZLFC) THEN
2457 C1=CLDDPTH+ZLCL-ZLFC
2458 C2=4000.0-ZLFC+ZLCL
2459 IF( ABS(C2) < 0.001 ) STOP 902
2460 C3=C1/C2
2461 VAL = 0.0
2462 DO II = 1, NT
2463 VAL = VAL + AINCKFSA(I,II,J)
2464 ENDDO
2465 AINC=(AINC+VAL)/FLOAT(NT+1)
2466 AINCKFI(i,j)=AINC
2467 AINC3=C3*AINC+(1-C3)*AINC2
2468 AINC=AINC3
2469 KFLAG=3
2470 GO TO 255
2471 ENDIF
2472 ENDIF
2474 IF(NCOUNT.GT.10)THEN
2475 AINC4=AINC
2476 IF(CLDDPTH .GE. 4.E3) THEN
2477 KFLAG=4
2478 GOTO 265
2479 ELSE IF((CLDDPTH) .LT. 4.E3 .AND. &
2480 CLDDPTH+ZLCL .GT. ZLFC) THEN
2481 C1=CLDDPTH+ZLCL-ZLFC
2482 C2=4000.0-ZLFC+ZLCL
2483 IF( ABS(C2) < 0.001 ) STOP 903
2484 C3=C1/C2
2485 VAL = 0.0
2486 DO II = 1, NT
2487 VAL = VAL + AINCKFSA(I,II,J)
2488 ENDDO
2489 AINC=(AINC+VAL)/FLOAT(NT+1)
2490 AINCKFI(i,j)=AINC
2491 AINC3=C3*AINC+(1-C3)*AINC2
2492 AINC=AINC3
2493 KFLAG=3
2494 GO TO 255
2495 ENDIF
2496 ENDIF
2498 IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) THEN
2499 AINC4=AINC
2500 IF(CLDDPTH .GE. 4.E3) THEN
2501 KFLAG=4
2502 GOTO 265
2503 ELSE IF((CLDDPTH) .LT. 4.E3 .AND. &
2504 CLDDPTH+ZLCL .GT. ZLFC) THEN
2505 C1=CLDDPTH+ZLCL-ZLFC
2506 C2=4000.0-ZLFC+ZLCL
2507 IF( ABS(C2) < 0.001 ) STOP 904
2508 C3=C1/C2
2509 VAL = 0.0
2510 DO II = 1, NT
2511 VAL = VAL + AINCKFSA(I,II,J)
2512 ENDDO
2513 AINC=(AINC+VAL)/FLOAT(NT+1)
2514 AINCKFI(i,j)=AINC
2515 AINC3=C3*AINC+(1-C3)*AINC2
2516 AINC=AINC3
2517 KFLAG=3
2518 GO TO 255
2519 ENDIF
2520 ENDIF
2521 !
2522 !...ADJUST MASS FLUX BY THE FACTOR AINC TO CONVERGE TO SPECIFIED DEGREE OF STAB-
2523 !...ILIZATION...
2524 !
2525 IF ( wrf_dm_on_monitor()) THEN
2526 CALL get_wrf_debug_level( dbg_level )
2527 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2528 WRITE(message,1080) LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
2529 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, NSTEP=',F5.0,I3, &
2530 'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2531 CALL wrf_message(message)
2532 ENDIF
2533 ENDIF
2535 AINC=AINC*STAB*ABE/(DABE+1.E-8)
2537 255 CONTINUE
2539 AINC=MIN(AINCMX,AINC)
2540 AINCMIN = 0.
2541 AINC=MAX(AINC,AINCMIN)
2542 PPTMLT=PPTML2*AINC
2543 TDER=TDER2*AINC
2544 PPTFLX=PPTFL2*AINC
2545 DO NK=1,LTOP
2546 UMF(NK)=UMF2(NK)*AINC
2547 DMF(NK)=DMF2(NK)*AINC
2548 DETLQ(NK)=DETLQ2(NK)*AINC
2549 DETIC(NK)=DETIC2(NK)*AINC
2550 UDR(NK)=UDR2(NK)*AINC
2551 UER(NK)=UER2(NK)*AINC
2552 DER(NK)=DER2(NK)*AINC
2553 DDR(NK)=DDR2(NK)*AINC
2554 DTFM(NK)=0.
2555 DMS(NK)=0.
2556 ENDDO
2557 !
2558 !...IF THE DOWNDRAFT OVERDRAWS THE INITIAL MASS AVAILABLE AT THE LFS, ALLOW
2559 !...PEFF TO CHANGE SO THAT DOWNDRAFT MASS FLUX IS NOT THE LIMITING FACTOR
2560 !...IN THE TOTAL CONVECTIVE MASS FLUX. THIS IS USUALLY NECESSARY ONLY WHEN
2561 !...THE DOWNDRAFT IS VERY SHALLOW...
2562 !
2563 DMFMIN=UER(LFS)-EMS(LFS)/TIMEC
2564 IF(DER(LFS).LT.DMFMIN .AND. DMFMIN.LT.0.)THEN
2565 RF=DMFMIN/DER(LFS)
2566 DO NK=LDB,LFS
2567 DER2(NK)=DER2(NK)*RF
2568 DDR2(NK)=DDR2(NK)*RF
2569 DMF2(NK)=DMF2(NK)*RF
2570 DER(NK)=DER2(NK)*AINC
2571 DDR(NK)=DDR2(NK)*AINC
2572 DMF(NK)=DMF2(NK)*AINC
2573 ENDDO
2574 TDER2=TDER2*RF
2575 TDER=TDER2*AINC
2576 IF(LFS.GE.KLCL)THEN
2577 UPDIN2=1.-(1.-EQFRC(LFS))*DMF(LFS)*UPDINC/UMF(LFS)
2578 ELSE
2579 UPDIN2=1.
2580 ENDIF
2581 CPR=TRPPT+PPR*(UPDIN2-1.)
2582 PEFF=1.-(TDER+(1.-EQFRC(LFS))*DMF(LFS)*(RLIQ(LFS)+RICE(LFS)))/ &
2583 (CPR*AINC)
2584 PPTFL2=PEFF*CPR
2585 PPTFLX=PPTFL2*AINC
2586 F1=UPDIN2/(AINC*UPDINC)
2587 DO NK=LC,LFS
2588 UMF2(NK)=UMF(NK)*F1
2589 UMF(NK)=UMF2(NK)*AINC
2590 UDR2(NK)=UDR(NK)*F1
2591 UDR(NK)=UDR2(NK)*AINC
2592 UER2(NK)=UER(NK)*F1
2593 UER(NK)=UER2(NK)*AINC
2594 DETLQ2(NK)=DETLQ(NK)*F1
2595 DETLQ(NK)=DETLQ2(NK)*AINC
2596 DETIC2(NK)=DETIC(NK)*F1
2597 DETIC(NK)=DETIC2(NK)*AINC
2598 PPTML2=PPTML2-PPTICE(NK)*(1.-UPDIN2/UPDINC)
2599 PPTLIQ(NK)=PPTLIQ(NK)*F1*AINC
2600 PPTICE(NK)=PPTICE(NK)*F1*AINC
2601 ENDDO
2602 PPTMLT=PPTML2*AINC
2603 UPDINC=UPDIN2
2604 ENDIF
2606 GO TO 175
2608 265 CONTINUE
2610 #ifdef DENG_SHCU1D
2611 DO K=KLCL,LTOP
2612 CLDUA(K)= AINC*AU0/DXSQ ! updraft fraction
2613 ENDDO
2614 #endif
2616 ! IF ( wrf_dm_on_monitor()) THEN
2617 ! CALL get_wrf_debug_level( dbg_level )
2618 ! IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2619 ! WRITE(message,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
2620 ! CALL wrf_message(message)
2621 ! ENDIF
2622 ! ENDIF
2623 !
2624 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
2625 !
2626 ! IF ( wrf_dm_on_monitor()) THEN
2627 ! CALL get_wrf_debug_level( dbg_level )
2628 ! IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2629 ! WRITE(message,3071)EQFRC(LFS),UPDINC
2630 ! CALL wrf_message(message)
2631 !3071 FORMAT('EQFRC(LFS) =',F7.4,'UPDINC =',F8.4)
2632 ! WRITE(message,969)
2633 ! CALL wrf_message(message)
2634 !969 FORMAT(1x,' P ',' DP ',' DT K/D ', &
2635 ! ' DR K/D ',' OMG ', &
2636 ! ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
2637 ! ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC')
2638 ! ENDIF
2639 ! ENDIF
2641 DO NK=1,LTOP
2642 K=LTOP-NK+1
2643 DTT=(TG(K)-T0(K))*86400./TIMEC
2644 RL=XLV0-XLV1*TG(K)
2645 DR=-(QQG(K)-Q0(K))*RL*86400./(TIMEC*CP)
2646 UDFRC=UDR(K)*TIMEC*EMSD(K)
2647 UEFRC=UER(K)*TIMEC*EMSD(K)
2648 DDFRC=DDR(K)*TIMEC*EMSD(K)
2649 DEFRC=-DER(K)*TIMEC*EMSD(K)
2651 ! IF ( wrf_dm_on_monitor()) THEN
2652 ! CALL get_wrf_debug_level( dbg_level )
2653 ! IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2654 ! WRITE(message,'(F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F7.3)') &
2655 ! P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4, &
2656 ! UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
2657 ! W0AVG0(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)* &
2658 ! TIMEC*EMSD(K)*1.E3
2659 ! CALL wrf_message(message)
2660 ! ENDIF
2661 ! ENDIF
2663 ENDDO
2665 #ifdef DENG_SHCU1D
2666 DO k = 1, KL
2667 tt6(k)=UMF(k)*QU(k)
2668 tt7(k)=OMG(k)/G*DXSQ*Q0(k)
2669 ENDDO
2670 CLDBMFLX(I,J)=UMF(KCBASE)/1.0E+6
2671 #endif
2673 ! IF ( wrf_dm_on_monitor()) THEN
2674 ! CALL get_wrf_debug_level( dbg_level )
2675 ! IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2676 ! WRITE(message,'(A3,15A7,2A8)') &
2677 ! 'K','P','Z','T0','TG','DT','TU','TD','Q0','QQG', &
2678 ! 'DQ','QU','QD','QCG','WU','WD','RH0','RHG'
2679 ! CALL wrf_message(message)
2680 ! ENDIF
2681 ! ENDIF
2683 DO NK=1,KX
2684 K=KX-NK+1
2685 DTT=TG(K)-T0(K)
2686 TUC=TU(K)-T00
2687 IF(K.LT.LC .OR. K.GT.LTOP)TUC=0.
2688 TDC=TZ(K)-T00
2689 IF((K.LT.LDB .OR. K.GT.LDT) .AND. K.NE.LFS)TDC=0.
2690 ES=1.E3*SVP1*EXP(SVP2*(TG(K)-SVPT0)/(TG(K)-SVP3))
2691 QGS=ES*0.622/(P0(K)-ES)
2692 RH_0=Q0(K)/QES(K)
2693 RHG=QQG(K)/QGS
2695 ! IF ( wrf_dm_on_monitor()) THEN
2696 ! CALL get_wrf_debug_level( dbg_level )
2697 ! IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2698 ! WRITE(message,'(I3,F7.2,F7.0,10F7.2,F7.3,2F7.2,2F8.3)') &
2699 ! K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC, &
2700 ! TDC,Q0(K)*1000.,QQG(K)*1000.,(QQG(K)-Q0(K))*1000.,QU(K)* &
2701 ! 1000.,QD(K)*1000.,QCG(K)*1000.,WU(K),WD(K),RH_0,RHG
2702 ! CALL wrf_message(message)
2703 ! ENDIF
2704 ! ENDIF
2706 ENDDO
2707 !
2708 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A SOUNDING
2709 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
2710 !
2711 IF(ISTOP.EQ.1)THEN
2712 IF ( wrf_dm_on_monitor()) THEN
2713 CALL get_wrf_debug_level( dbg_level )
2714 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2715 DO K = 1,KX
2716 WRITE(message,'(2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)') &
2717 Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2718 U0(K),V0(K),DP(K)/100.,W0AVG0(K)
2719 CALL wrf_message(message)
2720 ENDDO
2721 ENDIF
2722 ENDIF
2723 ENDIF
2725 CNDTNF=(1.-EQFRC(LFS))*(RLIQ(LFS)+RICE(LFS))*DMF(LFS)
2727 ! IF ( wrf_dm_on_monitor()) THEN
2728 ! CALL get_wrf_debug_level( dbg_level )
2729 ! IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2730 ! WRITE(message,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
2731 !1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
2732 ! CALL wrf_message(message)
2733 ! ENDIF
2734 ! ENDIF
2735 !
2736 ! EVALUATE MOISTURE BUDGET...
2737 !
2738 QINIT=0.
2739 QFNL=0.
2740 QCFINL=0.
2741 DO NK=1,LTOP
2742 QINIT=QINIT+Q0(NK)*EMS(NK)
2743 QFNL=QFNL+(QQG(NK)+QCG(NK))*DP(NK)*DXSQ/G
2744 QCFINL=QCFINL+QCG(NK)*DP(NK)*DXSQ/G
2745 ENDDO
2746 QFNL=QFNL+PPTFLX*TIMEC
2747 ERR2=(QFNL-QINIT)*100./QINIT
2748 RELERR=ERR2*QINIT/(PPTFLX*TIMEC+1.E-10)
2750 IF ( wrf_dm_on_monitor()) THEN
2751 CALL get_wrf_debug_level( dbg_level )
2752 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2753 WRITE(message,*)'QFNL-QINIT, QCFINL =',QFNL-QINIT,QCFINL
2754 CALL wrf_message(message)
2755 WRITE(message,1110)QINIT,QFNL,ERR2
2756 CALL wrf_message(message)
2757 WRITE(message,1200)RELERR
2758 CALL wrf_message(message)
2759 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2760 ' TOTAL WATER CHANGE =',F8.2,'%')
2761 1200 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
2762 ! WRITE(message,*)'TDER, CPR, USR, TRPPT =',TDER,CPR*AINC,USR*AINC,TRPPT*AINC
2763 ! CALL wrf_message(message)
2764 ENDIF
2765 ENDIF
2766 !
2767 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM FOR
2768 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
2769 !
2770 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/(0.5*DT2))
2771 !
2772 !...CALCULATE CONVECTIVE TENDENCIES...
2773 !
2774 loop320: DO K = KX,1,-1
2775 DTDT(K) = (TG(K)-T0(K))/TIMEC
2776 DQDT(K) = (QQG(K)-Q0(K))/TIMEC
2777 DQLDT(K)= (QCG(K)-QC0(K))/TIMEC
2778 DQRDT(K)= 0.0 ! set this to zero for now, unlike kfpara (don't remember why yet)
2779 DDCADT(K)=(DCA(K)-DCA0(K))/TIMEC
2780 IF(DQLDT(K) .LE. 1.0E-17 &
2781 .OR. RH0(K) .GT. RHCRIT) DDCADT(K)=0.0
2782 AREA=AMAX1((CLDAREAB(I,K,J)+DDCADT(K)*DT),0.0)
2783 IF(AREA .LE. 1.0E-17 &
2784 .OR. (DQLDT(K) .LE. 1.0E-17) &
2785 .OR. RH0(K) .GT. RHCRIT) THEN
2786 DQCDCDT(K)=0.0
2787 ELSE
2788 DQCDCDT(K)=(DQLDT(K)-CLDLIQB(I,K,J)*DDCADT(K))/AREA
2789 ENDIF
2790 #ifdef DENG_SHCU1D
2792 CLDATEN0(K)=DDCADT(K)
2793 CLDLCTN0(K)=DQCDCDT(K)
2794 tt4(k) = RLIQ(K)+RICE(K)
2795 #endif
2796 ENDDO loop320
2798 !
2799 ! CALCULATE THE CONVECTIVE RAINFALL
2800 !
2801 RAINSHV(I,J)=.5*DT2*PPTFLX/DXSQ
2803 IF ( wrf_dm_on_monitor()) THEN
2804 CALL get_wrf_debug_level( dbg_level )
2805 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2806 WRITE(message,909) RAINSHV(I,J)*NIC
2807 909 FORMAT(' CONVECTIVE RAINFALL =',F8.4,' CM')
2808 CALL wrf_message(message)
2809 ENDIF
2810 ENDIF
2811 !
2812 !...FEED BACK CONVECTIVE TENDENCIES TO RESOLVABLE SCALE EXCEPT THAT QC IS NOT
2813 ! FED BACK TO THE GRID SCALE WHEN THE GRID IS SUB-SATURATED. INSTEAD, IT
2814 ! IS DETRAINED TO FORM NBC. DCATEN(I,K): NBC AREA TENDENCY.
2815 ! QCDCTEN(I,K): NBC LIQUID/ICE WATER CONTENT.
2816 !
2817 DO K=1,KX
2818 TTEN(K)=TTEN(K)+ DTDT(K)
2819 QVTEN(K)=QVTEN(K)+ DQDT(K)
2820 QRTEN(K)=QRTEN(K)+DQRDT(K)
2821 IF(RH0(K) .GT. RHCRIT) THEN
2822 QCTEN(K)=QCTEN(K)+DQLDT(K)
2823 ELSE
2824 DCATEN(K)=DCATEN(K)+DDCADT(K)
2825 QCDCTEN(K)=QCDCTEN(K)+DQCDCDT(K)
2826 ENDIF
2827 ENDDO
2829 !
2830 ! raincn(i,j) is in cm. We can convert it to KG/m^2 in following way:
2831 !
2832 ! 1cm=1cm*1g/cm^3=1g/cm^2=10^-3Kg/10^-4m^2=10Kg/m^2
2833 !
2834 RAINSH(I,J) = RAINSH(I,J)+RAINSHV(I,J)
2836 #ifdef DENG_SHCU1D
2838 PLET = P0(LET)
2839 TLET = TU(LET)
2840 PLFSG(I,J) = P0(LFS)
2841 TLFSG(I,J) = TU(LFS)
2842 PLDTG(I,J) = P0(LDT)
2843 TLDTG(I,J) = TLCL*(PLDTG(I,J)/PLCL)**ROVCP
2844 PLDBG(I,J) = P0(LDB)
2845 TLDBG(I,J) = TLCL*(PLDBG(I,J)/PLCL)**ROVCP
2846 PLLCG(I,J) = P0(LLC)
2847 TLLCG(I,J) = TU(LLC)
2849 nca_ctr2(i,j)=nca_ctr2(i,j)+1
2851 #endif
2853 425 CONTINUE
2855 #ifdef DENG_SHCU1D
2856 nca_ctr4(i,j)=nca_ctr4(i,j)+1
2857 PPBLG(I,J)=PPBL
2858 TPBLG(I,J)=TPBL
2859 PLCLG(I,J)=PLCL
2860 TLCLG(I,J)=TLCL
2862 IF( i == 1 .AND. j == 1 ) THEN
2863 write(32,*) xtime/60.,ZLCL
2864 write(33,*) xtime/60.,CLDDPTH+zlcl
2865 write(34,*) xtime/60.0,zlfc
2866 ! CALL MAXIM(wutemp,KX,KLCL+1,KX,KVAL1)
2867 ! CALL MINIM(wsubsub,KX,KLCL+1,KX,KVAL2)
2868 ! write(82,*) xtime/60.0,wutemp(kval1),wsubsub(kval2)
2869 write(63,*) xtime/60.,CLDBMFLX(I,J)
2870 write(98,*) xtime/60.,tt6(KCBASE),tt7(KCBASE)
2871 ENDIF
2872 #endif
2873 !
2874 ! FIND NBC TOP (KDCT), NBC BASE (KDCB), LOCATION OF MAXIMUM NBC AREA
2875 ! (KAMAX) AND MAXIMUM CLD LIQUID WATER CONTENT (KLCMAX), AND THE
2876 ! LOCATION 1000M AWAY FROM KAMAX (K1000)
2877 !
2878 KDCB=1
2879 KDCT=1
2880 K1000=KDCT
2881 loop429: DO K=KL-1,2,-1
2882 IF(CLDAREAB(I,K,J) .GT. 0.001 .OR. RH0(K) .GT. RHCRIT) THEN
2883 KDCT=K
2884 EXIT loop429
2885 ENDIF
2886 ENDDO loop429
2888 loop441: DO K=2,KL-1 ! UPWARD
2889 IF(CLDAREAB(I,K,J) .GT. 0.001 .OR. RH0(K) .GT. RHCRIT ) THEN
2890 KDCB=K
2891 EXIT loop441
2892 ENDIF
2893 ENDDO loop441
2895 #ifdef DENG_SHCU1D
2896 DCLDTOP(I,J)=Z0(KDCT)
2897 DCLDBASE(I,J)=Z0(KDCB)
2898 #endif
2899 KDCLDTOP(I,J)=KDCT
2900 KDCLDBAS(I,J)=KDCB
2902 IF(AINC .NE. 999.0) THEN
2903 NUPDRAFT=AINC
2904 ELSE
2905 NUPDRAFT=0
2906 ENDIF
2908 #ifdef DENG_SHCU1D
2909 IF( i == 1 .AND. j == 1 ) THEN
2910 write(93,*) xtime/60.,DCLDTOP(I,J)
2911 write(94,*) xtime/60.,DCLDBASE(I,J)
2912 IF(AINC .NE. 999.0) THEN
2913 write(19,*) xtime/60.,RAD
2914 write(25,*) xtime/60.,NUPDRAFT
2915 ELSE
2916 VAL=0.0
2917 write(19,*) xtime/60.,VAL
2918 write(25,*) xtime/60.,VAL
2919 ENDIF
2920 ENDIF
2921 #endif
2923 IF(KDCT .LE. 1) THEN
2924 #ifdef DENG_SHCU1D
2925 nca_ctr3(i,j)=nca_ctr3(i,j)+1
2926 #endif
2927 GOTO 345 ! END OF DISSIPATION CALCULATION
2928 ENDIF
2930 KAMAX=KDCB
2931 KLCMAX=KDCB
2932 DO K=KDCB,KDCT
2933 IF( CLDAREAB(I,K,J) .GT. CLDAREAB(I,KAMAX,J) ) KAMAX = K
2934 IF( CLDLIQB(I,K,J) .GT. CLDLIQB(I,KLCMAX,J) ) KLCMAX = K
2935 ENDDO
2937 loop430: DO K=KAMAX,1,-1
2938 IF(Z0(KAMAX)-Z0(K) .GE. 1000.0) THEN
2939 K1000 = K
2940 EXIT loop430
2941 ENDIF
2942 ENDDO loop430
2944 K1000=MAX(K1000,KDCB)
2945 !
2946 ! WHEN THE GRID IS SATURATED, CONVERT NBC TO RESOLVABLE-SCALE CLOUD
2947 !
2948 DO K=1,KL
2949 IF(RH0(K) .GT. RHCRIT) THEN
2950 CLQ=AMAX1(CLDLIQB(I,K,J)+QCDCTEN(K)*DT,0.0)
2951 AREA=AMAX1(CLDAREAB(I,K,J)+DCATEN(K)*DT,0.0)
2952 QCDCTEN(K)=QCDCTEN(K)-CLQ/DT
2953 DCATEN(K)=DCATEN(K)-AREA/DT
2954 QCTEN(K)=QCTEN(K)+AREA*CLQ/DT
2955 ENDIF
2956 ENDDO
2957 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
2958 ! BEGINING OF DISSIPATION CACULATION
2959 ! goto 345
2960 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
2961 !
2962 ! EVAPORATION OF NBC AT THE CLOUD EDGE *******************************
2963 !
2964 #ifdef DENG_SHCU1D
2965 IF ( wrf_dm_on_monitor()) THEN
2966 CALL get_wrf_debug_level( dbg_level )
2967 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
2968 if( mod(ktau,KTAU1HUR) .eq. 0 ) then
2969 write(*, *) ' ------- a evaporation'
2970 endif
2971 ENDIF
2972 ENDIF
2973 #endif
2975 DIFFK=1.0E-5
2976 loop332: DO K=1,KX
2977 RLC=AMAX1(CLDLIQB(I,K,J)+QCDCTEN(K)*DT,0.0)
2978 VAL=QES(K)-Q0(K)
2979 VAL=AMAX1(0.0,VAL)
2980 E2=(0.01+CLDAREAB(I,K,J))*DIFFK*VAL*SQRT(NUPDRAFT)
2981 IF(RLC .LE. 1.0E-17) CYCLE loop332
2982 EOVLC=AMIN1(E2/RLC,CLDAREAB(I,K,J)/DT)
2983 DCATEN(K)=DCATEN(K)-EOVLC
2984 #ifdef DENG_SHCU1D
2985 cldaten(k)=-EOVLC
2986 #endif
2987 QVTEN(K)=QVTEN(K)+EOVLC*RLC
2988 IF(T0(K) .LE. TO) THEN
2989 RL=XLS
2990 ELSE
2991 RL=XLV0-XLV1*T0(K)
2992 ENDIF
2993 RCP=CP*(1.+0.887*Q0(K))
2994 RLOVCP=RL/RCP
2995 TTEN(K)=TTEN(K)-EOVLC*RLC*RLOVCP
2996 ENDDO loop332
2997 !
2998 ! VERTICAL DIFFUSION OF NBC LIQUID/ICE WATER CONTENT **********************
2999 !
3000 ! CALCULAT KH (VERTICAL DIFF COEF ) AND KR ( RADIATIVE INDUCED VERTICAL
3001 ! DIFF COEF). KR IS CALCULATED AT KAMAX+1 (FULL LAYER AT CLOUD TOP), THEN
3002 ! CONSTANT UNTIL CLOUD TOP. IT LINEARLY DROPS TO ZERO AT K1000, WHICH IS
3003 ! THE LARGER ONE OF KDCB AND K1000.
3004 !
3005 DO K=KDCB,KDCT+1 ! UPWARD
3006 KV(K)=KTH0(K)
3007 ENDDO
3009 DO K=1,KL
3010 tlong(k)=TEN_RADL0(K)
3011 tshort(k)=TEN_RADS0(K)
3012 ENDDO
3013 K=KAMAX
3014 KP1=K+1
3015 CALL MINIM(tlong, KX,KDCB,KDCT,KRADLMAX)
3016 CALL MAXIM(tshort,KX,KDCB,KDCT,KRADSMAX)
3018 KR(KP1)=BBLS0(K)*BBLS0(K)/(T0(K)*(P00/P0(K))**ROCPQ) &
3019 *(ABS(TEN_RADL0(KRADLMAX))*DZQ(K)/15.0 &
3020 +ABS(TEN_RADS0(KRADSMAX))*DZQ(K)/50.0)
3022 DO K=K1000,KAMAX+1
3023 IF(KDCT+1 .NE. K1000) KR(K)=KR(KAMAX+1)*(FLOAT(K-K1000) &
3024 /FLOAT(KAMAX+1-K1000))
3025 ENDDO
3026 DO K=KAMAX+1,KDCT+1
3027 KR(K)=KR(KAMAX+1)
3028 ENDDO
3029 !
3030 ! CALCULATE IN-CLOUD DIFF FLUXES ASSOCIATED WITH KV and KR. MAKE SURE
3031 ! THAT UPWARD FLUXES AT CLOUD TOP(KDCT+1) VANISH, SO THAT THERE IS NO
3032 ! LIQUID DIFFUSE UPWARD
3033 !
3034 DO K=1,KDCB-1
3035 LFLUX(K)=0.0
3036 LFLUX1(K)=0.0
3037 LFLUX2(K)=0.0
3038 ENDDO
3040 DO K=KDCT+1,KXP1
3041 LFLUX(K)=0.0
3042 LFLUX1(K)=0.0
3043 LFLUX2(K)=0.0
3044 ENDDO
3046 DO K=KDCT,KDCB,-1
3047 IF(K .EQ. KDCB) THEN
3048 C1=15.0
3049 ELSE
3050 C1=DZA(K-1)
3051 ENDIF
3052 AREA1=AMAX1(CLDAREAB(I,K,J)+DCATEN(K)*DT,0.0)
3053 AREA2=AMAX1(CLDAREAB(I,K-1,J)+DCATEN(K-1)*DT,0.0)
3054 CLQ1=AMAX1(CLDLIQB(I,K,J)+QCDCTEN(K)*DT,0.0)
3055 CLQ2=AMAX1(CLDLIQB(I,K-1,J)+QCDCTEN(K-1)*DT,0.0)
3056 LFLUX1(K)=-RHOE(K)*KV(K)*(AREA1*CLQ1-AREA2*CLQ2)/C1
3057 LFLUX2(K)=RHOE(K)*KR(K)*(AREA1*CLQ1-AREA2*CLQ2)/C1
3058 LFLUX1(K)=AMIN1(LFLUX1(K),0.0)
3059 LFLUX2(K)=AMIN1(LFLUX2(K),0.0)
3060 LFLUX(K)=LFLUX1(K)+LFLUX2(K)
3061 ! 0.05 IS USED HERE TO
3062 ! MAKE SURE THAT NOT MORE THAN 5% OF THE TOTAL LIQUID MASS AT THIS LAYER
3063 ! IS REMOVED BY IN-CLOUD MIXING IN A SINGLE TIME STEP.
3064 C1=LFLUX(K+1)-CLQ1*DP(K)/(G*DT)*0.05
3065 LFLUX(K)=AMAX1(LFLUX(K),C1)
3066 ENDDO
3067 DO K=KDCT+1,KDCB,-1
3068 IF(RH0(K) .GT. RHCRIT .OR. &
3069 (CLDAREAB(I,K,J)+DCATEN(K)*DT) .LE. 1.0e-6) then
3070 LFLUX(K)=0.0
3071 LFLUX(K+1)=0.0
3072 ENDIF
3073 ENDDO
3074 !
3075 ! CALCULATE DIFF TENDENCY ASSOCIATED WITH IN-CLOUD FLUXES CALCULATED ABOVE.
3076 ! FLUXES AT CLD BASE (KDCB) CARRY LIQUID TO THE LAYER BELOW THE CLD BASE
3077 ! AND EVAPORATE IN THAT LAYER
3078 !
3079 #ifdef DENG_SHCU1D
3080 IF ( wrf_dm_on_monitor()) THEN
3081 CALL get_wrf_debug_level( dbg_level )
3082 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
3083 if( mod(ktau,KTAU1HUR) .eq. 0 ) then
3084 write(*, *) ' --- lc in-cloud mixing'
3085 endif
3086 ENDIF
3087 ENDIF
3088 #endif
3089 !
3090 ! THIS IS TO ENSURE THAT THE LAYER BELOW THE CLOUD BASE (KDCB-1) DOES
3091 ! NOT REACH SATURATION ONLY BECAUSE OF MIXING (USE 95% HERE). ADJUST
3092 ! LFLUX(KDCB) SO THAT THE AMOUNT OF MASS TRANPORTED TO THIS LAYER IS
3093 ! LESS THAN THE AMOUNT OF MASS NEEDED TO RAISE THE RH TO 95%.
3094 !
3095 RHCR=0.95
3096 IF(Q0(KDCB-1)/QES(KDCB-1) .GT. RHCR) THEN
3097 LFLUX(KDCB)=0.0
3098 ELSE
3099 VAL1=-G*(LFLUX(KDCB)-LFLUX(KDCB-1))/DP(KDCB-1)
3100 VAL2=(RHCR*QES(KDCB-1)-Q0(KDCB-1))/DT
3101 VAL=AMIN1(VAL1,VAL2)
3102 LFLUX(KDCB)=-VAL*DP(KDCB-1)/G
3103 ENDIF
3105 loop436: DO K = KDCT, KDCB, -1
3106 AREA=AMAX1(CLDAREAB(I,K,J)+DCATEN(K)*DT,0.0)
3107 IF( RH0(K) .GT. RHCRIT .OR. AREA .LE. 1.0e-6 ) CYCLE loop436
3108 QCDCTEN(K)=QCDCTEN(K)-G*(LFLUX(K+1)-LFLUX(K))/DP(K)/AREA
3109 #ifdef DENG_SHCU1D
3110 cldlctn1(k)=-G*(LFLUX(K+1)-LFLUX(K))/DP(K)/AREA
3111 #endif
3112 ENDDO loop436
3114 K=KDCB-1
3115 IF(T0(K) .LE. TO) THEN
3116 RL=XLS
3117 ELSE
3118 RL=XLV0-XLV1*T0(K)
3119 ENDIF
3120 RCP=CP*(1.+0.887*Q0(K))
3121 RLOVCP=RL/RCP
3122 QVTEN(K)=QVTEN(K)-G*(LFLUX(K+1)-LFLUX(K))/DP(K)
3123 TTEN(K)=TTEN(K)+G*(LFLUX(K+1)-LFLUX(K))/DP(K)*RLOVCP
3124 !
3125 ! CALCULATE DIFF TENDENCY ASSOCIATED WITH DOWNWARD FLUXES AT
3126 ! EXPOSED CLD BASE IF CLD AREA INCREASE WITH HEIGHT. THIS FLUXE WILL
3127 ! CAUSE LIQUID DECREASE IN THE CLD LAYER AND WATER VAPOR INCREASE IN THE
3128 ! LAYER BELOW THE EXPOSED CLD BASE DUE TO EVAPORATION
3129 !
3130 #ifdef DENG_SHCU1D
3131 IF ( wrf_dm_on_monitor()) THEN
3132 CALL get_wrf_debug_level( dbg_level )
3133 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
3134 if(mod(ktau,KTAU1HUR) .eq. 0) then
3135 write(*, *) ' ----lc exp-cld-base mixing'
3136 endif
3137 ENDIF
3138 ENDIF
3139 #endif
3141 loop440: DO K=KDCT,KDCB+1,-1
3142 AREADIFF=CLDAREAB(I,K,J)-CLDAREAB(I,K-1,J)
3143 AREADIFF=AMAX1(AREADIFF,0.0)
3144 CLQ=AMAX1(CLDLIQB(I,K,J)+QCDCTEN(K)*DT,0.0)
3145 C1=-CLQ*DP(K)/(G*DT)*0.05 ! limit of flux that remove 5% of
3146 ! all the liquid in the layer above
3147 FX1=-RHOE(K)*KV(K)*CLDLIQB(I,K,J)/15.*AREADIFF
3148 FX2=-RHOE(K)*KR(K)*CLDLIQB(I,K,J)/15.*AREADIFF
3149 FX=FX1+FX2
3150 FX=AMAX1(FX,C1)
3152 IF((CLDAREAB(I,K,J)+DCATEN(K)*DT) .LE. 1.0E-6) CYCLE loop440
3154 IF(RH0(K) .GT. RHCRIT) CYCLE loop440
3155 !
3156 ! THIS IS TO ENSURE THAT THE LAYER BELOW THE EXPOSED CLOUD BASE (K-1) DOES
3157 ! NOT REACH SATURATION ONLY BECAUSE OF MIXING (USE 95% HERE). ADJUST
3158 ! LFLUX(KDCB) SO THAT THE AMOUNT OF MASS TRANPORTED TO THIS LAYER IS
3159 ! LESS THAN THE AMOUNT OF MASS NEEDED TO RAISE THE RH TO 95%.
3160 !
3161 RHCR=0.95
3162 IF(Q0(K-1)/QES(K-1) .GT. RHCR) THEN
3163 FX=0.0
3164 ELSE
3165 VAL1=-G*FX/DP(K-1)
3166 VAL2=(RHCR*QES(K-1)-Q0(K-1))/DT
3167 VAL=AMIN1(VAL1,VAL2)
3168 FX=-VAL*DP(K-1)/G
3169 ENDIF
3171 AREA=AMAX1(CLDAREAB(I,K,J)+DCATEN(K)*DT,0.0)
3172 IF(AREA .LE. 1.0E-6) CYCLE loop440 ! this line was not in MM5
3173 QCDCTEN(K)=QCDCTEN(K)+G*FX/DP(K)/AREA
3174 #ifdef DENG_SHCU1D
3175 cldlctn2(k)=G*FX/DP(K)/AREA
3176 #endif
3177 IF(T0(K-1) .LE. TO) THEN
3178 RL=XLS
3179 ELSE
3180 RL=XLV0-XLV1*T0(K-1)
3181 ENDIF
3182 RCP=CP*(1.+0.887*Q0(K-1))
3183 RLOVCP=RL/RCP
3184 QVTEN(K-1)=QVTEN(K-1)-G*FX/DP(K-1)
3185 TTEN(K-1)=TTEN(K-1)+G*FX/DP(K-1)*RLOVCP
3186 ENDDO loop440
3187 !
3188 ! CLOUD TOP ENTRAINMENT INSTABILITY (CTEI) ********************************
3189 ! RANDALL (JAS 1980), DEARDOFF (JAS 1980) AND DEL GENIO (J. CLIMATE 1996)
3190 !
3191 RATIO=0.0
3192 RATIOMIN=0.0
3193 RATIOMAX=0.0
3194 IF(RH0(KDCT) .GT. RHCRIT) GOTO 340
3195 EFOLD=1.0E-4
3196 EX=1.0
3197 DEPTHMAX=100.0
3198 IF(KDCT .LE. 0) GOTO 340
3199 IF(T0(KDCT) .GE. TO) THEN
3200 HIN=CP*T0(KDCT)+G*Z0(KDCT)+XLV*QES(KDCT)
3201 ELSE
3202 HIN=CP*T0(KDCT)+G*Z0(KDCT)+XLS*QES(KDCT)
3203 ENDIF
3204 IF(T0(KDCT+1) .GE. TO) THEN
3205 HOUT=CP*T0(KDCT+1)+G*Z0(KDCT+1)+XLV*Q0(KDCT+1)
3206 ELSE
3207 HOUT=CP*T0(KDCT+1)+G*Z0(KDCT+1)+XLS*Q0(KDCT+1)
3208 ENDIF
3209 DELTH=HIN-HOUT
3210 QIN=QES(KDCT)+CLDLIQB(I,KDCT,J)
3211 QOUT=Q0(KDCT+1)
3212 DELTQ=QIN-QOUT
3213 RATIO=DELTH/(XLV*DELTQ+1.E-8)
3214 RKAPA=CP*0.5*(T0(KDCT)+T0(KDCT+1))/XLV
3215 DELTA=0.608
3216 IF(T0(KDCT) .GT. TO) THEN
3217 DQSSDT=QES(KDCT)*SVP2*(SVPT0-SVP3) &
3218 /((T0(KDCT)-SVP3)*(T0(KDCT)-SVP3))
3219 RL=XLV0-XLV1*T0(KDCT)
3220 ELSE
3221 DQSSDT=QES(KDCT)*6.15E3/(T0(KDCT)*T0(KDCT))
3222 RL=XLS
3223 ENDIF
3224 RCP=CP*(1.+0.887*Q0(KDCT))
3225 RLOVCP=RL/RCP
3227 GAMA=RLOVCP*DQSSDT
3228 BEITA=(1+(1+DELTA)*GAMA*RKAPA)/(1+GAMA)
3229 RATIOMIN=RKAPA/BEITA
3230 RATIOMAX=(1+GAMA)*(1+(1-DELTA)*RKAPA) &
3231 /(2+(1+(1+DELTA)*RKAPA)*GAMA)
3233 IF(RATIO .LE. RATIOMIN) THEN
3234 SIGMMA=0.0
3235 GOTO 340
3236 ELSE IF(RATIO .GE. RATIOMAX) THEN
3237 SIGMMA=EFOLD
3238 ELSE
3239 SIGMMA=EFOLD*((RATIO-RATIOMIN)/(RATIOMAX-RATIOMIN))**EX
3240 ENDIF
3242 #ifdef DENG_SHCU1D
3243 IF ( wrf_dm_on_monitor()) THEN
3244 CALL get_wrf_debug_level( dbg_level )
3245 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
3246 if( mod(ktau,KTAU1HUR) .eq. 0 ) then
3247 write(* ,*) ' -------- lc CTEI '
3248 endif
3249 ENDIF
3250 ENDIF
3251 #endif
3252 !
3253 ! LOCATE K WHERE ITS DISTANCE FROM KDCT IS 100M
3254 !
3255 K100=KDCT
3256 loop335: DO L=KDCT-1,KDCB,-1
3257 IF((Z0(KDCT)+0.5*DZQ(KDCT))-(Z0(L)-0.5*DZQ(L)) .GE. DEPTHMAX) THEN
3258 K100=K100-1
3259 EXIT loop335
3260 ENDIF
3261 K100=L
3262 ENDDO loop335
3264 RHCR=0.95
3265 IF(RH0(KDCT+1) .GT. RHCR) GOTO 340
3267 loop344: DO L=KDCT,K100,-1
3268 IF(RH0(L) .GT. RHCRIT) CYCLE loop344
3269 WFUN=FLOAT(L-K100+1)/FLOAT(KDCT-K100+1)
3270 AREA=AMAX1(CLDAREAB(I,L,J)+DCATEN(L)*DT,0.0)
3271 IF(T0(L) .LE. TO) THEN
3272 RL=XLS
3273 ELSE
3274 RL=XLV0-XLV1*T0(L)
3275 ENDIF
3276 RCP=CP*(1.+0.887*Q0(L))
3277 RLOVCP=RL/RCP
3279 VAL1=SIGMMA*CLDLIQB(I,L,J)*WFUN
3280 QCDCTEN(L)=QCDCTEN(L)-VAL1
3281 #ifdef DENG_SHCU1D
3282 cldlctn5(L)=-VAL1
3283 #endif
3284 QVTEN(KDCT+1)=QVTEN(KDCT+1) &
3285 +VAL1*AREA*DP(L)/DP(KDCT+1)
3286 TTEN(L)=TTEN(L)-VAL1*AREA*RLOVCP
3287 VAL=Q0(KDCT+1) + QVTEN(KDCT+1)*DT
3288 !
3289 ! ENSURING THAT THE LAYER ABOVE CLOUD TOP DOES NOT REACH SATURATION (RH=95%)
3290 ! ONLY BECAUSE OF CTEI. IF IT HAPPENS DUE TO CTEI THEN REMOVE THIS EFFECT
3291 !
3292 IF( VAL/QES(KDCT+1) .GT. RHCR ) THEN
3293 QVTEN(KDCT+1)=QVTEN(KDCT+1) &
3294 -VAL1*DP(L)/DP(KDCT+1)
3295 TTEN(L)=TTEN(L)+VAL1*AREA*RLOVCP
3296 ENDIF
3297 ENDDO loop344
3299 340 CONTINUE
3300 !
3301 ! PRECIPITATION (DRIZZLE) PROCESS ****************************************
3302 !
3303 #ifdef DENG_SHCU1D
3304 IF ( wrf_dm_on_monitor()) THEN
3305 CALL get_wrf_debug_level( dbg_level )
3306 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
3307 if( mod(ktau,KTAU1HUR) .eq. 0 ) then
3308 write( *,*) ' ------ lc drizzle'
3309 endif
3310 ENDIF
3311 ENDIF
3312 #endif
3314 DO K=1,KL
3315 PRC(K)=0.
3316 PRA(K)=0.
3317 ENDDO
3319 BVTS3=3.+BVTS
3320 PPI=1./(PIE*N0R)
3321 G3PBS=GAMMA(3.+BVTS)
3322 PRACS=PIE*N0S*AVTS*G3PBS*.25*ESI ! from MM5
3324 DO K=1,KX
3325 TOUT=T0(K)
3326 SC2=AMAX1(0.0,CLDLIQB(I,K,J))
3327 SC3=AMAX1(0.0,QR0(K))
3328 IF(TOUT.GT.TO)THEN
3329 SC7=(RHOE(K)*SC3*PPI/1000.)**0.25
3330 ELSE
3331 RHOS=0.1
3332 PPIS=1./(PIE*N0S*RHOS)
3333 SC7=(RHOE(K)*SC3*PPIS/1000.)**0.25
3334 ENDIF
3335 IF(TOUT .LT. TO)THEN
3336 XNC=XN0*EXP(0.6*(TO-TOUT))/RHOE(K) ! XN0 = 1.E-2 in MM5 domain/initial/param.F
3337 XNC=AMAX1(XNC,10000./RHOE(K))
3338 !... AUTOCONVERSION OF CLOUD ICE TO SNOW:
3339 PRC(K)=AMAX1(0.,(SC2-XMMAX*XNC)/DT)
3340 !... ACCRETION OF CLOUD ICE BY SNOW:
3341 IF(SC3 .LT. 1.0E-17)THEN
3342 PRA(K)=0.
3343 ELSE
3344 PRA(K)=PRACS*SC7**BVTS3*SC2
3345 ENDIF
3346 ELSE
3347 !--- AUTOCONVERSION OF CLOUD WATER TO RAINWATER:
3348 PRC(K)=AMAX1(0.,QCK1*(SC2-QCTH))
3349 !--- ACCRETION OF CLOUD WATER BY RAINWATER:
3350 BVT3=3.+BVT
3351 IF (SC3 .LT. 1.0E-17) THEN
3352 PRA(K)=0.
3353 ELSE
3354 G3PB=GAMMA(3.+BVT)
3355 PRAC=PIE*N0R*AVT*G3PB*0.25
3356 PRA(K)=PRAC*SC7**BVT3*SC2
3357 ENDIF
3358 ENDIF
3360 VAL=PRC(K)+PRA(K)
3361 VAL=MIN(VAL,CLDLIQB(I,K,J)/DT)
3362 AREA=AMAX1(CLDAREAB(I,K,J)+DCATEN(K)*DT,0.0)
3363 IF(AREA .LT. 1.0E-17) VAL=0.0
3364 QRTEN(K)=QRTEN(K)+VAL*AREA
3365 QCDCTEN(K)=QCDCTEN(K)-VAL
3366 #ifdef DENG_SHCU1D
3367 cldlctn3(k)=-VAL
3368 #endif
3370 ENDDO
3371 !
3372 ! ICE SEDIMENTATION AS A SINK TERM TO CLDLIQ ******************************
3373 !
3374 #ifdef DENG_SHCU1D
3375 IF ( wrf_dm_on_monitor()) THEN
3376 CALL get_wrf_debug_level( dbg_level )
3377 IF( dbg_level >= 10 .AND. i == i0 .AND. j == j0 ) THEN
3378 if( mod(ktau,KTAU1HUR) .eq. 0 ) then
3379 write( *,*) ' ----- lc ice settling'
3380 write(18,*) ' ----- lc ice settling'
3381 endif
3382 ENDIF
3383 ENDIF
3384 #endif
3386 DO K=1,KL
3387 CLQ=AMAX1(CLDLIQB(I,K,J)+QCDCTEN(K)*DT,0.0)
3388 AREA=AMAX1(CLDAREAB(I,K,J)+DCATEN(K)*DT,0.0)
3389 SCR3(K)=CLQ*AREA
3390 ENDDO
3392 NSTEP=1
3393 DO K=1,KL
3394 !
3395 ! NOTE: RHO2=RHO/(PSB*1000) WHERE RHO IS THE AIR DENSITY IN KG/M^3
3396 ! NSTEP IS THE # OF TIME STEPS NEEDED FOR A DROPLET TO MOVE FROM THE LAYER-TOP
3397 ! TO THE LAYER-BOTTOM
3398 !
3399 RHO2=P0(K)/(R*( T0(K)+DT*TTEN(K)))
3400 IF( (T0(K)+DT*TTEN(K)) .GT.TO)THEN
3401 VT2C=0.
3402 ELSE
3403 !
3404 ! NOTE: SEDIMENTATION FORMULA OF HEYMSFIELD AND DONNER (1990, JAS)
3405 !
3406 VT2C=3.29*(RHO2* SCR3(K) )**0.16 ! m/s
3407 ENDIF
3409 RGVC(K)=G*RHO2*VT2C ! rho*g*vt2c/(1000*ps)
3410 NSTEP=MAX0(IFIX(RGVC(K)*DT/DSIGMA(K)+1.),NSTEP) ! dvt2c*dt/DZ
3411 ENDDO
3413 DO N=1,NSTEP
3414 IF(N.GT.1000)STOP &
3415 'IN SUB. SHALLOW (ICE SETTLING), NSTEP TOO LARGE, PROBABLY NAN'
3417 DO K=1,KL
3418 FALOUTC(K)=RGVC(K)*SCR3(K)
3419 ENDDO
3421 loop341: DO K=KL-1,2,-1
3422 AREA=AMAX1(CLDAREAB(I,K,J)+DCATEN(K)*DT,0.0)
3423 CLQ=AMAX1(CLDLIQB(I,K,J)+QCDCTEN(K)*DT,0.0)
3424 IF(RH0(K) .GT. RHCRIT) CYCLE loop341
3425 IF(AREA .LE. 1.E-4) CYCLE loop341
3426 FALTNDC=(FALOUTC(K)-FALOUTC(K+1))/DSIGMA(K)
3427 C1=AMIN1(FALTNDC,AREA*CLQ/DT*NSTEP)
3428 C1=AMAX1(0.0,C1)
3429 QCDCTEN(K)=QCDCTEN(K)-C1/AREA/NSTEP
3430 #ifdef DENG_SHCU1D
3431 cldlctn4(k) = cldlctn4(k)-C1/AREA/NSTEP ! 1-D graphics
3432 #endif
3433 SCR3(K)=SCR3(K)-C1*DT/NSTEP/AREA
3434 QCTEN(K-1)=QCTEN(K-1)+C1*DSIGMA(K)/DSIGMA(K-1)/NSTEP
3435 ! if(i.eq.82.and.j.eq.74) then
3436 ! endif
3437 RGVC(K)=AMAX1(RGVC(K)/DSIGMA(K),RGVC(K+1)/DSIGMA(K+1))*DSIGMA(K)
3438 ENDDO loop341
3440 ENDDO
3441 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
3442 345 CONTINUE ! END OF DISSIPATION CACULATION
3443 !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
3444 !
3445 ! PUT INTO TAU-1 ARRAYS
3446 !
3447 CLDDPTHB(I,J)=CLDDPTH
3448 CLDTOPB(I,J)=CLDTOP
3449 LTOPB(I,J)=LTOP
3450 DO K=1,KL
3451 WUB(I,K,J)=WU(K)
3452 ENDDO
3454 #ifdef DENG_SHCU1D
3456 IF( i == 1 .AND. j == 1 ) THEN
3458 DO NK = 1,KX
3459 EE=Q0(NK)*P0(NK)/(0.622+Q0(NK))
3460 TLOG=LOG(EE/SVP1/1000.)
3461 TDPT=(SVP2*SVPT0-SVP3*TLOG)/(SVP2-TLOG)
3462 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T0(NK)-T00))* &
3463 (T0(NK)-TDPT)
3464 THTA=T0(NK)*(1.E5/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
3465 THETEE(NK)=THTA* &
3466 EXP((3374.6525/TSAT-2.5403)*Q0(NK)*(1.+0.81*Q0(NK)))
3467 thtetmp(nk)=THETEE(NK)
3468 thtatmp(nk)=THTA
3469 ENDDO
3471 if(mod(ktau,KTAU1HUR) .eq. 0) then
3472 write(*,878) xtime/60.0, zpbl, Zlcl,Zlfc,CLDDPTHB(I,J)
3473 write(18,878) Zlcl,Zlfc,CLDDPTH
3474 878 format(5x,'xtime=',f6.1,' Zpbl=',f5.1,' Zlcl=',f6.1,1x,'Zlfc=',f8.1, &
3475 1x,'CLD_Depth=',f7.1)
3476 endif
3478 KK=INT(PLOTFRQ*(KTAU1HUR/60.))
3479 IF( MOD(KTAU-1,KK) .EQ. 0) THEN
3480 write(17,*) xtime/60.0,ktau,ZPBL,ZLCL,LLC,KPBL,KLCL, &
3481 CLDTOP,KX-LTOP+1,KX-KDCB+1,KX-KDCT+1,RHCRIT
3482 ! write(43,*) xtime/60.0,ktau,ZPBL,ZLCL,LLC,KPBL,KLCL, &
3483 ! CLDTOP,LTOP,KDCB,KDCT
3484 DO k = KL, 1, -1
3485 write(17,*) kl-k+1,p0(k)/100.0,z0(k),CLDUA(k)*100.0, &
3486 CLDAREAB(I,K,J)*100.0, &
3487 CLDLIQB(I,K,J)*1000.0, &
3488 CLDAREAB(I,K,J)*CLDLIQB(I,K,J)*1000.0, &
3489 tt4(k)*1000.0, &
3490 RH0(K)*100.0, &
3491 ca_rad0(K)*100.0, &
3492 QC0_expl(K)*1000.
3493 ENDDO
3495 ! DO k = 1, KL
3496 ! write(43,*) k,p0(k),z0(k), &
3497 ! cldaten0(k)*aa, &
3498 ! cldaten(k)*aa, &
3499 ! cldlctn0(k)*bb, &
3500 ! cldlctn1(k)*bb, &
3501 ! cldlctn2(k)*bb, &
3502 ! cldlctn3(k)*bb, &
3503 ! cldlctn4(k)*bb, &
3504 ! cldlctn5(k)*bb ! g/hour
3505 ! ENDDO
3506 ENDIF
3509 KK=INT(PLOTFRQ*(KTAU1HUR/60.))
3510 IF( MOD(KTAU-1,KK) .EQ. 0) THEN
3512 write(10,*) xtime/60.0, ktau-1
3514 DO k = kl, 1, -1
3515 wsp_plot(k) = sqrt(u0(k)**2+v0(k)**2)
3516 if( ABS(wsp_plot(k)) < 0.001 ) then
3517 wdr_plot(k)=0.
3518 else if(v0(k) .gt. 0.) then
3519 wdr_plot(k)=atan(u0(k)/v0(k))
3520 wdr_plot(k)=wdr_plot(k)*180./pie+180.
3521 else if (v0(k) .lt. 0.) then
3522 if( u0(k) .gt. 0.0 ) then
3523 wdr_plot(k)=atan(u0(k)/v0(k))
3524 wdr_plot(k)=wdr_plot(k)*180./PIE+360.
3525 else if( u0(k) .lt. 0.0 ) then
3526 wdr_plot(k)=atan(u0(k)/v0(k))
3527 wdr_plot(k)=wdr_plot(k)*180./PIE
3528 end if
3529 else if (u0(k) .gt.0.) then
3530 wdr_plot(k)=270.
3531 else
3532 wdr_plot(k)=90.
3533 end if
3535 qqq= Q0(K)
3536 qqq = AMAX1(qqq,1.0E-25)
3537 EES=qqq*p0(k)/100.0/(0.622+qqq)
3538 IF( t0(k) .GT. TO ) THEN
3539 TLOG=LOG(EES/SVP1/10.)
3540 td_plot(k) = (SVP2*SVPT0-SVP3*TLOG)/(SVP2-TLOG)-SVPT0
3541 ELSE
3542 TLOG=LOG(EES/0.611/10.)
3543 td_plot(k) = 6150./(22.514-TLOG)-SVPT0
3544 ENDIF
3546 write(10,*) p0(k)/100.0, t0(k)-SVPT0, td_plot(k), wsp_plot(k), wdr_plot(k)
3547 ENDDO
3549 write(10,*) (TUG(K),K=1,KL),TPBLG(I,J),PPBLG(I,J),TLCLG(I,J), &
3550 PLCLG(I,J), PCLDTPG(I,J), TCLDTPG(I,J), &
3551 CLDHGTG(I,J),PLET,TLET,PLFSG(I,J),TLFSG(I,J),PLDTG(I,J), &
3552 TLDTG(I,J),PLDBG(I,J),TLDBG(I,J),PLLCG(I,J),TLLCG(I,J), &
3553 TRPPT
3554 ENDIF
3556 KK=INT(PLOTFRQ*(KTAU1HUR/60.))
3557 IF( MOD(KTAU-1,KK) .EQ. 0) THEN
3559 write(7,*) xtime/60.0,ktau,ZPBL,ZLCL,LLC,KKPBL,KLCL, &
3560 CLDTOP,kx-LTOP+1,kx-KDCB+1,kx-KDCT+1
3561 do k = kx, 1, -1
3562 write(7,*) kx-k+1, p0(k)/100.0, z_at_w0(k)-ht, tke0(k)
3563 enddo
3565 write(67,*) (TKE0(K),K=KXP1-KL/3+1,1,-1), (z_at_w0(k)-ht,k=KXP1-KL/3+1,1,-1)
3566 ENDIF
3568 ENDIF
3569 #endif
3571 END SUBROUTINE deng_shcu
3573 !====================================================================
3574 SUBROUTINE deng_shcu_init(RTHSHTEN,RQVSHTEN,RQCSHTEN,RQRSHTEN, &
3575 RUSHTEN,RVSHTEN,RDCASHTEN,RQCDCSHTEN,W0AVG, &
3576 PBLHAVG, TKEAVG, cldareaa, cldareab, &
3577 cldliqa, cldliqb, ca_rad, cw_rad, &
3578 wub, pblmax, ltopb, clddpthb, cldtopb, &
3579 capesave, ainckfsa, radsave, &
3580 rainsh, rainshvb, kdcldtop, kdcldbas, &
3581 xtime1, &
3582 restart, &
3583 SVP1,SVP2,SVP3,SVPT0, &
3584 ids, ide, jds, jde, kds, kde, &
3585 ims, ime, jms, jme, kms, kme, &
3586 its, ite, jts, jte, kts, kte )
3587 !--------------------------------------------------------------------
3588 IMPLICIT NONE
3589 !--------------------------------------------------------------------
3590 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
3591 ims, ime, jms, jme, kms, kme, &
3592 its, ite, jts, jte, kts, kte
3594 LOGICAL , INTENT(IN) :: restart
3595 REAL , INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
3597 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
3598 RUSHTEN, &
3599 RVSHTEN, &
3600 RTHSHTEN, &
3601 RQVSHTEN, &
3602 RQCSHTEN, &
3603 RQRSHTEN, &
3604 RDCASHTEN, &
3605 RQCDCSHTEN
3607 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG, &
3608 TKEAVG
3610 REAL , DIMENSION( ims:ime , jms:jme), INTENT(OUT) :: pblmax, &
3611 rainsh, rainshvb, clddpthb, cldtopb, capesave, radsave, xtime1, &
3612 pblhavg
3614 REAL , DIMENSION( ims:ime , 1:100, jms:jme ), INTENT(OUT) :: ainckfsa
3616 INTEGER , DIMENSION( ims:ime , jms:jme), INTENT(OUT) :: ltopb &
3617 ,kdcldtop, kdcldbas
3619 REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(OUT) :: wub, &
3620 cldareaa, cldareab, cldliqa, cldliqb, ca_rad, cw_rad
3623 INTEGER :: i, j, k, itf, jtf, ktf
3625 jtf=min0(jte,jde-1)
3626 ktf=min0(kte,kde-1)
3627 itf=min0(ite,ide-1)
3629 IF(.NOT.RESTART)THEN
3630 DO j=jts,jtf
3631 DO i=its,itf
3632 DO k=kts,ktf
3633 RUSHTEN(i,k,j) = 0.
3634 RVSHTEN(i,k,j) = 0.
3635 RTHSHTEN(i,k,j) = 0.
3636 RQVSHTEN(i,k,j) = 0.
3637 RQCSHTEN(i,k,j) = 0.
3638 RQRSHTEN(i,k,j) = 0.
3639 RDCASHTEN(i,k,j) = 0.
3640 RQCDCSHTEN(i,k,j) = 0.
3641 W0AVG(i,k,j) = 0.
3642 TKEAVG(i,k,j) = 0.
3643 wub(i,k,j) = 0.
3644 cldareaa(i,k,j) = 0.0
3645 cldareab(i,k,j) = 0.0
3646 cldliqa(i,k,j) = 0.0
3647 cldliqb(i,k,j) = 0.0
3648 ca_rad(i,k,j) = 0.0
3649 cw_rad(i,k,j) = 0.0
3650 ENDDO
3651 PBLHAVG(i,j) = 0.0
3652 pblmax(i,j) = 0.0
3653 clddpthb(I,J) = 0.0
3654 cldtopb(I,J) = 0.0
3655 rainsh(i,j) = 0.0
3656 rainshvb(i,j) = 0.0
3657 capesave(i,j) = 0.0
3658 radsave(i,j) = 0.0
3659 xtime1(i,j) = 0.0
3660 ltopb(i,j) = 1
3661 kdcldtop(i,j) = 1
3662 kdcldbas(i,j) = 1
3663 ENDDO
3664 ENDDO
3666 DO j=jts,jtf
3667 DO i=its,itf
3668 DO k=1,100
3669 ainckfsa(i,k,j) = 0.0
3670 ENDDO
3671 ENDDO
3672 ENDDO
3674 ENDIF
3676 CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
3677 ! CALL KF_LUTABBG(SVP1,SVP2,SVP3,SVPT0)
3679 END SUBROUTINE deng_shcu_init
3681 !==========================================================
3682 SUBROUTINE time_avg_2d( dt, ADAPT_STEP_FLAG, time_period_avg_min &
3683 ,array2d, array2d_avg &
3684 ,its,ite, jts,jte &
3685 )
3686 IMPLICIT NONE
3687 INTEGER, INTENT(IN ) :: its,ite, jts,jte
3689 REAL, INTENT(IN ) :: time_period_avg_min, dt ! time window to perform
3690 !a running mean value
3691 LOGICAL ,INTENT(IN) :: adapt_step_flag
3692 REAL, DIMENSION( its:ite , jts:jte ), INTENT(IN ) :: array2d
3693 REAL, DIMENSION( its:ite , jts:jte ), INTENT(INOUT) :: array2d_avg
3695 INTEGER :: i, j, ntst
3696 REAL :: tst, AVGfctr, fctr, den
3698 NTST = NINT( time_period_avg_min * 60.0 / dt ) ! number of time steps in time_period_avg_min
3699 TST = FLOAT(NTST)
3700 IF( TST <= 0 ) STOP 'Denominator has to be greater than zero - deng_shcu_driver'
3702 if (ADAPT_STEP_FLAG) then
3703 AVGfctr = MAX( time_period_avg_min*60,dt ) - dt
3704 fctr = dt
3705 den = MAX( time_period_avg_min*60,dt )
3706 else
3707 AVGfctr = (TST-1.)
3708 fctr = 1.
3709 den = TST
3710 endif
3712 DO J = jts,jte
3713 DO I = its,ite
3714 array2d_avg(i,j) = ( array2d_avg(i,j) * AVGfctr + array2d(i,j) * fctr ) / den
3715 ENDDO
3716 ENDDO
3718 ! write(*,'(5f10.4)') array2d(1,1), array2d_avg(1,1), AVGfctr, fctr, den
3720 END SUBROUTINE time_avg_2d
3722 !==========================================================
3723 SUBROUTINE time_avg_3d( dt, ADAPT_STEP_FLAG, time_period_avg_min &
3724 ,array3d, array3d_avg &
3725 ,its,ite, jts,jte, kts,kte &
3726 )
3727 IMPLICIT NONE
3728 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
3730 REAL, INTENT(IN ) :: time_period_avg_min, dt ! time window to perform
3731 !a running mean value
3732 LOGICAL ,INTENT(IN) :: adapt_step_flag
3733 REAL, DIMENSION( its:ite , kts:kte , jts:jte ), INTENT(IN ) :: array3d
3734 REAL, DIMENSION( its:ite , kts:kte , jts:jte ), INTENT(INOUT) :: array3d_avg
3736 INTEGER :: i, j, k, ntst
3737 REAL :: tst, AVGfctr, fctr, den
3739 NTST = NINT( time_period_avg_min * 60.0 / dt ) ! number of time steps in time_period_avg_min
3740 TST = FLOAT(NTST)
3741 IF( TST <= 0 ) STOP 'Denominator has to be greater than zero - deng_shcu_driver'
3743 if (ADAPT_STEP_FLAG) then
3744 AVGfctr = MAX( time_period_avg_min*60,dt ) - dt
3745 fctr = dt
3746 den = MAX( time_period_avg_min*60,dt )
3747 else
3748 AVGfctr = (TST-1.)
3749 fctr = 1.
3750 den = TST
3751 endif
3753 DO J = jts,jte
3754 DO I = its,ite
3755 DO k = kts,kte
3756 array3d_avg(i,k,j) = ( array3d_avg(i,k,j) * AVGfctr + array3d(i,k,j) * fctr ) / den
3757 ENDDO
3758 ENDDO
3759 ENDDO
3761 ! do k = 1,kte
3762 ! write(*,'(i3,5f10.4)') k, array3d(1,k,1), array3d_avg(1,k,1), AVGfctr, fctr, den
3763 ! enddo
3765 END SUBROUTINE time_avg_3d
3767 !====================================================================
3768 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
3769 QNEWIC,QLQOUT,QICOUT,G)
3770 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
3771 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
3772 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
3773 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
3774 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
3775 !----------------------
3776 IMPLICIT NONE
3777 !----------------------
3778 REAL, INTENT(IN ) :: G
3779 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
3780 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
3781 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
3783 IF(ABS(WTW).LT.1.E-4) WTW=1.E-4
3784 QTOT=QLIQ+QICE
3785 QNEW=QNEWLQ+QNEWIC
3786 !
3787 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
3788 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
3789 ! LEVELS...
3790 !
3791 QEST=0.5*(QTOT+QNEW)
3792 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
3793 IF(G1.LT.0.0)G1=0.
3794 WAVG=0.5*(SQRT(WTW)+SQRT(G1))
3795 CONV=RATE*DZ/WAVG
3796 !
3797 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
3798 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
3799 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
3800 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
3801 !
3802 RATIO3=QNEWLQ/(QNEW+1.E-8)
3803 ! OLDQ=QTOT
3804 QTOT=QTOT+0.6*QNEW
3805 OLDQ=QTOT
3806 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
3807 QTOT=QTOT*EXP(-CONV)
3808 !
3809 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
3810 ! PARCEL AT THIS LEVEL...
3811 !
3812 DQ=OLDQ-QTOT
3813 QLQOUT=RATIO4*DQ
3814 QICOUT=(1.-RATIO4)*DQ
3815 !
3816 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
3817 ! LATE VERTICAL VELOCITY
3818 !
3819 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
3820 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
3821 IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
3822 !
3823 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
3824 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
3825 !
3826 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
3827 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
3828 QNEWLQ=0.
3829 QNEWIC=0.
3831 END SUBROUTINE CONDLOAD
3832 !====================================================================
3833 SUBROUTINE DTFRZNEW(TU,P,THTEU,QVAP,QLIQ,QICE,RATIO2, &
3834 QNWFRZ,RL,FRC1,EFFQ,IFLAG,XLV0,XLV1,XLS0,XLS1, &
3835 ! ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE, &
3836 SVP1,SVP2,SVPT0,SVP3)
3837 !-----------------------------------------------------------------------
3838 IMPLICIT NONE
3839 !-----------------------------------------------------------------------
3840 REAL, INTENT(IN ) :: P, SVP1,SVP2,SVPT0,SVP3, &
3841 XLV0,XLV1,XLS0,XLS1, EFFQ
3843 REAL, INTENT(INOUT) :: TU, THTEU, QVAP, QLIQ, QICE, RATIO2, &
3844 QNWFRZ, FRC1, RL
3845 INTEGER, INTENT(INOUT) :: IFLAG
3846 REAL :: A, B, C, C5, RLC, RLS, RLF, QVAP1, &
3847 RV, CS, QLQFRZ, QNEW, ESLIQ, ESICE, &
3848 CP, DQVAP, DTFRZ, TU1, ES, pi
3849 !
3850 !...ALLOW GLACIATION OF THE UPDRAFT TO OCCUR AS AN APPROXIMATELY LINEAR
3851 ! FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE TTFRZ TO TBFRZ...
3852 !
3853 RV=461.5
3854 C5=1.0723E-3
3855 !
3856 !...ADJUST THE LIQUID WATER CONCENTRATIONS FROM FRESH CONDENSATE AND THA
3857 ! BROUGHT UP FROM LOWER LEVELS TO AN AMOUNT THAT WOULD BE PRESENT IF N
3858 ! LIQUID WATER HAD FROZEN THUS FAR...THIS IS NECESSARY BECAUSE THE
3859 ! EXPRESSION FOR TEMP CHANGE IS MULTIPLIED BY THE FRACTION EQUAL TO TH
3860 ! PARCEL TEMP DECREASE SINCE THE LAST MODEL LEVEL DIVIDED BY THE TOTAL
3861 ! GLACIATION INTERVAL, SO THAT EFFECTIVELY THIS APPROXIMATELY ALLOWS A
3862 ! AMOUNT OF LIQUID WATER TO FREEZE WHICH IS EQUAL TO THIS SAME FRACTIO
3863 ! OF THE LIQUID WATER THAT WAS PRESENT BEFORE THE GLACIATION PROCESS W
3864 ! INITIATED...ALSO, TO ALLOW THETAU TO CONVERT APPROXIMATELY LINEARLY
3865 ! ITS VALUE WITH RESPECT TO ICE, WE NEED TO ALLOW A PORTION OF THE FRE
3866 ! CONDENSATE TO CONTRIBUTE TO THE GLACIATION PROCESS; THE FRACTIONAL
3867 ! AMOUNT THAT APPLIES TO THIS PORTION IS 1/2 OF THE FRACTIONAL AMOUNT
3868 ! FROZEN OF THE "OLD" CONDENSATE BECAUSE THIS FRESH CONDENSATE IS ONLY
3869 ! PRODUCED GRADUALLY OVER THE LAYER...NOTE THAT IN TERMS OF THE DYNAMI
3870 ! OF THE PRECIPITATION PROCESS, IE. PRECIPITATION FALLOUT, THIS FRACTI
3871 ! AMNT OF FRESH CONDENSATE HAS ALREADY BEEN INCLUDED IN THE ICE CATEGO
3872 !
3873 QLQFRZ=QLIQ*EFFQ
3874 QNEW=QNWFRZ*EFFQ*0.5
3875 ! ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
3876 ESLIQ=1.E3*SVP1*EXP(SVP2*(TU-SVPT0)/(TU-SVP3))
3877 ! ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
3878 ESICE=611.0*EXP(22.514-6.15E3/TU)
3879 RLC=2.5E6-2369.276*(TU-273.16)
3880 RLS=2833922.-259.532*(TU-273.16)
3881 RLF=RLS-RLC
3882 CP=1005.7*(1.+0.89*QVAP)
3883 !
3884 ! A = D(ES)/DT IS THAT CALCULATED FROM BUCK`S (1981) EMPIRICAL FORMULAS
3885 ! FOR SATURATION VAPOR PRESSURE...
3886 !
3887 ! A=(CICE-BICE*DICE)/((TU-DICE)*(TU-DICE))
3888 A=6.15E+3/(TU*TU)
3889 B=RLS*0.622/P
3890 C=A*B*ESICE/CP
3891 DQVAP=B*(ESLIQ-ESICE)/(RLS+RLS*C)-RLF*(QLQFRZ+QNEW)/(RLS+RLS/C)
3892 DTFRZ=(RLF*(QLQFRZ+QNEW)+B*(ESLIQ-ESICE))/(CP+A*B*ESICE)
3893 TU1=TU
3894 QVAP1=QVAP
3895 TU=TU+FRC1*DTFRZ
3896 QVAP=QVAP-FRC1*DQVAP
3897 ES=QVAP*P/(0.622+QVAP)
3898 ! ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
3899 ESLIQ=1.E3*SVP1*EXP(SVP2*(TU-SVPT0)/(TU-SVP3))
3900 ! ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
3901 ESICE=611.0*EXP(22.514-6.15E3/TU)
3902 RATIO2=(ESLIQ-ES)/(ESLIQ-ESICE)
3903 !
3904 ! TYPICALLY, RATIO2 IS VERY CLOSE TO (TTFRZ-TU)/(TTFRZ-TBFRZ), USUALLY
3905 ! WITHIN 1% (USING TU BEFORE GALCIATION EFFECTS ARE APPLIED); IF THE
3906 ! INITIAL UPDRAFT TEMP IS BELOW TBFRZ AND RATIO2 IS STILL LESS THAN 1,
3907 ! AN ADJUSTMENT TO FRC1 AND RATIO2 IS INTRODUCED SO THAT GLACIATION
3908 ! EFFECTS ARE NOT UNDERESTIMATED; CONVERSELY, IF RATIO2 IS GREATER THAN
3909 ! FRC1 IS ADJUSTED SO THAT GLACIATION EFFECTS ARE NOT OVERESTIMATED...
3910 !
3911 IF(IFLAG.GT.0.AND.RATIO2.LT.1)THEN
3912 FRC1=FRC1+(1.-RATIO2)
3913 TU=TU1+FRC1*DTFRZ
3914 QVAP=QVAP1-FRC1*DQVAP
3915 RATIO2=1.
3916 IFLAG=1
3917 GOTO 20
3918 ENDIF
3919 IF(RATIO2.GT.1.)THEN
3920 FRC1=FRC1-(RATIO2-1.)
3921 FRC1=AMAX1(0.0,FRC1)
3922 TU=TU1+FRC1*DTFRZ
3923 QVAP=QVAP1-FRC1*DQVAP
3924 RATIO2=1.
3925 IFLAG=1
3926 ENDIF
3927 !
3928 ! CALCULATE A HYBRID VALUE OF THETAU, ASSUMING THAT THE LATENT HEAT OF
3929 ! VAPORIZATION/SUBLIMATION CAN BE ESTIMATED USING THE SAME WEIGHTING
3930 ! FUNCTION AS THAT USED TO CALCULATE SATURATION VAPOR PRESSURE, CALCU-
3931 ! LATE NEW LIQUID WATER AND ICE CONCENTRATIONS...
3932 !
3933 20 RLC=XLV0-XLV1*TU
3934 RLS=XLS0-XLS1*TU
3935 RL=RATIO2*RLS+(1.-RATIO2)*RLC
3936 PI=(1.E5/P)**(0.2854*(1.-0.28*QVAP))
3937 THTEU=TU*PI*EXP(RL*QVAP*C5/TU*(1.+0.81*QVAP))
3938 IF(IFLAG.EQ.1)THEN
3939 QICE=QICE+FRC1*DQVAP+QLIQ
3940 QLIQ=0.
3941 ELSE
3942 QICE=QICE+FRC1*(DQVAP+QLQFRZ)
3943 QLIQ=QLIQ-FRC1*QLQFRZ
3944 ENDIF
3945 QNWFRZ=0.
3947 END SUBROUTINE DTFRZNEW
3948 ! ***********************************************************************
3949 !====================================================================
3950 !====================================================================
3951 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,R1,RL, &
3952 ! ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE, &
3953 SVP1,SVP2,SVPT0,SVP3,I,J)
3954 !-----------------------------------------------------------------------
3955 IMPLICIT NONE
3956 !-----------------------------------------------------------------------
3957 REAL, INTENT(IN ) :: P1,T1,Q1,R1,RL,SVP1,SVP2,SVPT0,SVP3
3958 INTEGER, INTENT(IN ) :: I,J
3959 REAL, INTENT( OUT) :: THT1
3960 REAL :: EE,TLOG,TDPT,TSAT,THT,TFPT
3961 REAL :: T00,P00,C1,C2,C3,C4,C5,tlogic, tsatlq, tsatic
3962 !-----------------------------------------------------------------------
3963 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
3964 0.278296,1.0723E-3/
3965 !
3966 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
3967 !
3968 IF(R1.LT.1.E-6)THEN
3969 ! EE=Q1*P1/(0.622+Q1)
3970 EE=max(Q1*P1/(0.622+Q1),1.e-9)
3971 ! TLOG=ALOG(EE/ALIQ)
3972 ! TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
3973 TLOG=ALOG(EE/SVP1/1000.)
3974 TDPT=(SVP2*SVPT0-SVP3*TLOG)/(SVP2-TLOG)
3975 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
3976 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
3977 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
3978 ELSEIF(ABS(R1-1.).LT.1.E-6)THEN
3979 ! EE=Q1*P1/(0.622+Q1)
3980 EE=max(Q1*P1/(0.622+Q1),1.e-9)
3981 ! TLOG=ALOG(EE/AICE)
3982 ! TFPT=(CICE-DICE*TLOG)/(BICE-TLOG)
3983 TLOG=ALOG(EE/611.0)
3984 TFPT=6150./(22.514-TLOG)
3985 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
3986 TSAT=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
3987 THT1=THT*EXP((C3/TSAT-C4)*Q1*(1.+0.81*Q1))
3988 ELSE
3989 ! EE=Q1*P1/(0.622+Q1)
3990 EE=max(Q1*P1/(0.622+Q1),1.e-9)
3991 ! TLOG=ALOG(EE/ALIQ)
3992 ! TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
3993 ! TLOGIC=ALOG(EE/AICE)
3994 ! TFPT=(CICE-DICE*TLOGIC)/(BICE-TLOGIC)
3995 TLOG=ALOG(EE/SVP1/1000.)
3996 TDPT=(SVP2*SVPT0-SVP3*TLOG)/(SVP2-TLOG)
3997 TLOGIC=ALOG(EE/611.0)
3998 TFPT=6150./(22.514-TLOGIC)
3999 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
4000 TSATLQ=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
4001 TSATIC=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
4002 TSAT=R1*TSATIC+(1.-R1)*TSATLQ
4003 THT1=THT*EXP(RL*Q1*C5/TSAT*(1.+0.81*Q1))
4004 ! NaN 326. NaN 30.88963 NaN 0.0 NaN -Infinity -Infinity 0
4005 ENDIF
4007 END SUBROUTINE ENVIRTHT
4008 !====================================================================
4009 SUBROUTINE ENVIRTHT2(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
4010 !
4011 !-----------------------------------------------------------------------
4012 IMPLICIT NONE
4013 !-----------------------------------------------------------------------
4014 REAL, INTENT(IN ) :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
4015 REAL, INTENT(INOUT) :: THT1
4016 REAL :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT, &
4017 T00,P00,C1,C2,C3,C4,C5
4018 INTEGER :: INDLU
4019 !-----------------------------------------------------------------------
4020 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
4021 0.278296,1.0723E-3/
4022 !
4023 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
4024 !
4025 ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
4026 ! For example, KF90 Eq. 10 no longer used
4027 !
4028 EE=Q1*P1/(0.622+Q1)
4029 ! TLOG=ALOG(EE/ALIQ)
4030 ! ...calculate LOG term using lookup table...
4031 !
4032 astrt=1.e-3
4033 ainc=0.075
4034 a1=ee/aliq
4035 tp=(a1-astrt)/ainc
4036 indlu=int(tp)+1
4037 value=(indlu-1)*ainc+astrt
4038 aintrp=(a1-value)/ainc
4039 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
4040 !
4041 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
4042 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
4043 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
4044 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
4045 !
4046 END SUBROUTINE ENVIRTHT2
4047 ! ***********************************************************************
4048 !====================================================================
4049 SUBROUTINE PROF5(EQ,EE,UD)
4050 !-----------------------------------------------------------------------
4051 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
4052 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
4053 ! `HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICAL
4054 ! TABLES` ED. BY ABRAMOWITZ AND STEGUN, NAT`L BUREAU OF STANDARDS APPLIED
4055 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
4056 ! JACK KAIN
4057 ! 7/6/89
4058 !-----------------------------------------------------------------------
4059 IMPLICIT NONE
4060 !-----------------------------------------------------------------------
4061 REAL, INTENT(IN ) :: EQ
4062 REAL, INTENT( OUT) :: EE,UD
4063 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
4064 !-----------------------------------------------------------------------
4065 !*********************************************************************
4066 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
4067 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
4068 0.9372980,0.33267,0.166666667,0.202765151/
4069 X=(EQ-0.5)/SIGMA
4070 Y=6.*EQ-3.
4071 EY=EXP(Y*Y/(-2))
4072 E45=EXP(-4.5)
4073 T2=1./(1.+P*ABS(Y))
4074 T1=0.500498
4075 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
4076 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
4077 IF(Y.GE.0.)THEN
4078 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
4079 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.-EQ)
4080 ELSE
4081 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
4082 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
4083 EQ/2.-EQ)
4084 ENDIF
4085 EE=EE/FE
4086 UD=UD/FE
4088 END SUBROUTINE PROF5
4089 ! ***********************************************************************
4090 !====================================================================
4091 FUNCTION TPDD(P,THTED,TGS,RS,RD,RH,XLV0,XLV1, &
4092 ! ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE, &
4093 SVP1,SVP2,SVPT0,SVP3)
4094 !************************* TPDD.FOR ************************************
4095 ! THIS SUBROUTINE ITERATIVELY EXTRACTS TEMPERATURE FROM EQUIVALENT *
4096 ! POTENTIAL TEMP. IT IS DESIGNED FOR USE WITH DOWNDRAFT CALCULATIONS.
4097 ! IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, PARCEL *
4098 ! TEMP, SPECIFIC HUMIDITY, AND LIQUID WATER CONTENT ARE ITERATIVELY *
4099 ! CALCULATED. *
4100 !-----------------------------------------------------------------------
4101 IMPLICIT NONE
4102 !-----------------------------------------------------------------------
4103 REAL, INTENT(IN ) :: P,THTED,TGS,RD,RH,XLV0,XLV1,SVP1,SVP2,SVPT0,SVP3
4104 ! REAL, INTENT( OUT ) :: tpdd
4105 REAL :: tpdd
4106 REAL, INTENT(INOUT) :: RS
4107 REAL :: ES, PI, THTGS, F0, T1, T0, CP, F1, DT, DSSDT, T1RH, RSRH, RL
4108 INTEGER :: ITCNT
4109 !-----------------------------------------------------------------------
4110 ! ES=ALIQ*EXP((BLIQ*TGS-CLIQ)/(TGS-DLIQ))
4111 ES=1.E3*SVP1*EXP(SVP2*(TGS-SVPT0)/(TGS-SVP3))
4112 RS=0.622*ES/(P-ES)
4113 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
4114 THTGS=TGS*PI*EXP((3374.6525/TGS-2.5403)*RS*(1.+0.81*RS))
4115 F0=THTGS-THTED
4116 T1=TGS-0.5*F0
4117 T0=TGS
4118 CP=1005.7
4119 !
4120 !...ITERATE TO FIND WET-BULB TEMPERATURE...
4121 !
4122 ITCNT=0
4123 ! 90 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
4124 90 ES=1.E3*SVP1*EXP(SVP2*(T1-SVPT0)/(T1-SVP3))
4125 RS=0.622*ES/(P-ES)
4126 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
4127 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*RS*(1.+0.81*RS))
4128 F1=THTGS-THTED
4129 IF(ABS(F1).LT.0.05)GOTO 50
4130 ITCNT=ITCNT+1
4131 IF(ITCNT.GT.10)GOTO 50
4132 DT=F1*(T1-T0)/(F1-F0)
4133 T0=T1
4134 F0=F1
4135 T1=T1-DT
4136 GOTO 90
4137 50 RL=XLV0-XLV1*T1
4138 !
4139 !...IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, ESTIMATE THE
4140 ! TEMPERATURE AND MIXING RATIO WHICH WILL YIELD THE APPROPRIATE VALUE.
4141 !
4142 IF(RH.EQ.1.)GOTO 110
4143 DSSDT=SVP2*(SVPT0-SVP3)/((T1-SVP3)*(T1-SVP3))
4144 DT=RL*RS*(1.-RH)/(CP+RL*RH*RS*DSSDT)
4145 T1RH=T1+DT
4146 ! ES=RH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
4147 ES=RH*1.E3*SVP1*EXP(SVP2*(T1RH-SVPT0)/(T1RH-SVP3))
4148 RSRH=0.622*ES/(P-ES)
4149 !
4150 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
4151 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
4152 !
4153 IF(RSRH.LT.RD)THEN
4154 RSRH=RD
4155 T1RH=T1+(RS-RSRH)*RL/CP
4156 ENDIF
4157 T1=T1RH
4158 RS=RSRH
4159 110 TPDD=T1
4160 ! IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPDD = ',
4161 ! + ITCNT
4162 RETURN
4163 END FUNCTION TPDD
4164 !====================================================================
4165 FUNCTION TPDDBG(P,THTED,TGS,RS,RD,RH,XLV0,XLV1, &
4166 ! ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE, &
4167 SVP1,SVP2,SVPT0,SVP3)
4168 !************************* TPDD.FOR ************************************
4169 ! THIS SUBROUTINE ITERATIVELY EXTRACTS TEMPERATURE FROM EQUIVALENT *
4170 ! POTENTIAL TEMP. IT IS DESIGNED FOR USE WITH DOWNDRAFT CALCULATIONS.
4171 ! IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, PARCEL *
4172 ! TEMP, SPECIFIC HUMIDITY, AND LIQUID WATER CONTENT ARE ITERATIVELY *
4173 ! CALCULATED. *
4174 !-----------------------------------------------------------------------
4175 !
4176 ! Modified BJG 1/17/2014
4177 !
4178 IMPLICIT NONE
4179 !-----------------------------------------------------------------------
4180 REAL, INTENT(IN ) :: P,THTED,TGS,RD,RH,XLV0,XLV1,SVP1,SVP2,SVPT0,SVP3
4181 ! REAL, INTENT( OUT ) :: tpddbg
4182 REAL :: tpddbg
4183 REAL, INTENT(INOUT) :: RS
4184 REAL :: ES, PI, THTGS, F0, T1, T0, CP, F1, DT, DSSDT, T1RH, RSRH, RL
4185 real :: f2, t2
4186 INTEGER :: ITCNT
4187 !-----------------------------------------------------------------------
4188 ! ES=ALIQ*EXP((BLIQ*TGS-CLIQ)/(TGS-DLIQ))
4189 ES=1.E3*SVP1*EXP(SVP2*(TGS-SVPT0)/(TGS-SVP3))
4190 RS=0.622*ES/(P-ES)
4191 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
4192 THTGS=TGS*PI*EXP((3374.6525/TGS-2.5403)*RS*(1.+0.81*RS))
4193 F0=THTGS-THTED
4194 !c T1=TGS-0.5*F0
4195 T2=TGS-0.5*F0
4196 T0=TGS
4198 t1 = t0
4199 f1 = f0
4201 CP=1005.7
4202 !
4203 !...ITERATE TO FIND WET-BULB TEMPERATURE...
4204 !
4205 ITCNT=0
4206 ! 90 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
4207 90 ES=1.E3*SVP1*EXP(SVP2*(t2-SVPT0)/(t2-SVP3))
4208 RS=0.622*ES/(P-ES)
4209 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
4210 THTGS=t2*PI*EXP((3374.6525/t2-2.5403)*RS*(1.+0.81*RS))
4212 f2=THTGS-THTED
4213 if((f1 * f2).lt.0.0) then
4214 f0 = f1
4215 t0 = t1
4216 else
4217 f0 = f0
4218 t0 = t0
4219 endif
4220 f1 = f2
4221 t1 = t2
4223 IF(ABS(F1).LT.0.05)GOTO 50
4224 ITCNT=ITCNT+1
4225 IF(ITCNT.GT.10)GOTO 50
4226 DT=F1*(T1-T0)/(F1-F0)
4228 if(abs(dt).lt.abs( (t1-t0)/2.0 )) then
4229 dt = (t1 - t0) / 2.0
4230 endif
4232 !c T0=T1
4233 !c F0=F1
4234 t2=T1-DT
4236 GOTO 90
4237 50 RL=XLV0-XLV1*T1
4238 !
4239 !...IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, ESTIMATE THE
4240 ! TEMPERATURE AND MIXING RATIO WHICH WILL YIELD THE APPROPRIATE VALUE.
4241 !
4242 IF(RH.EQ.1.)GOTO 110
4243 DSSDT=SVP2*(SVPT0-SVP3)/((T1-SVP3)*(T1-SVP3))
4244 DT=RL*RS*(1.-RH)/(CP+RL*RH*RS*DSSDT)
4245 T1RH=T1+DT
4246 ! ES=RH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
4247 ES=RH*1.E3*SVP1*EXP(SVP2*(T1RH-SVPT0)/(T1RH-SVP3))
4248 RSRH=0.622*ES/(P-ES)
4249 !
4250 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
4251 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
4252 !
4253 IF(RSRH.LT.RD)THEN
4254 RSRH=RD
4255 T1RH=T1+(RS-RSRH)*RL/CP
4256 ENDIF
4257 T1=T1RH
4258 RS=RSRH
4259 110 TPDDBG=T1
4260 ! IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPDD = ',
4261 ! + ITCNT
4262 RETURN
4263 END FUNCTION TPDDBG
4264 !====================================================================
4265 SUBROUTINE TPMIX(P,THTU,TU,QU,QLIQ,QICE,QNEWLQ,QNEWIC,RATIO2,RL, &
4266 XLV0,XLV1,XLS0,XLS1, &
4267 ! ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE, &
4268 SVP1,SVP2,SVPT0,SVP3)
4269 !-----------------------------------------------------------------------
4270 IMPLICIT NONE
4271 !-----------------------------------------------------------------------
4272 REAL, INTENT(IN ) :: P,THTU,RATIO2,RL, &
4273 XLV0,XLV1,XLS0,XLS1,SVP1,SVP2,SVPT0,SVP3
4274 REAL, INTENT( OUT) :: QNEWLQ,QNEWIC
4275 REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
4276 REAL :: C5, RV, ES, QS, PI, THTGS, ESLIQ, ESICE, F0, T1, T0, &
4277 F1, DT, QNEW, QTOT, DQICE, DQLIQ, RLL, CP, dq
4278 INTEGER :: ITCNT
4279 !-----------------------------------------------------------------------
4280 !
4281 !...THIS SUBROUTINE ITERATIVELY EXTRACTS WET-BULB TEMPERATURE FROM EQUIV
4282 ! POTENTIAL TEMPERATURE, THEN CHECKS TO SEE IF SUFFICIENT MOISTURE IS
4283 ! AVAILABLE TO ACHIEVE SATURATION...IF NOT, TEMPERATURE IS ADJUSTED
4284 ! ACCORDINGLY, IF SO, THE RESIDUAL LIQUID WATER/ICE CONCENTRATION IS
4285 ! DETERMINED...
4287 C5=1.0723E-3
4288 RV=461.5
4289 !
4290 ! ITERATE TO FIND WET BULB TEMPERATURE AS A FUNCTION OF EQUIVALENT POT
4291 ! TEMP AND PRS, ASSUMING SATURATION VAPOR PRESSURE...RATIO2 IS THE DEG
4292 ! OF GLACIATION...
4293 !
4294 IF(RATIO2.LT.1.E-6)THEN
4295 ! ES=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
4296 ES=1.E3*SVP1*EXP(SVP2*(TU-SVPT0)/(TU-SVP3))
4297 QS=0.622*ES/(P-ES)
4298 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4299 THTGS=TU*PI*EXP((3374.6525/TU-2.5403)*QS*(1.+0.81*QS))
4300 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
4301 ! ES=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
4302 ES=611.0*EXP(22.514-6.15E3/TU)
4303 QS=0.622*ES/(P-ES)
4304 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4305 THTGS=TU*PI*EXP((3114.834/TU-0.278296)*QS*(1.+0.81*QS))
4306 ELSE
4307 ! ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
4308 ! ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
4309 ESLIQ=1.E3*SVP1*EXP(SVP2*(TU-SVPT0)/(TU-SVP3))
4310 ESICE=611.0*EXP(22.514-6.15E3/TU)
4311 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
4312 QS=0.622*ES/(P-ES)
4313 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4314 THTGS=TU*PI*EXP(RL*QS*C5/TU*(1.+0.81*QS))
4315 ENDIF
4316 F0=THTGS-THTU
4317 T1=TU-0.5*F0
4318 T0=TU
4319 ITCNT=0
4320 90 IF(RATIO2.LT.1.E-6)THEN
4321 ! ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
4322 ES=1.E3*SVP1*EXP(SVP2*(T1-SVPT0)/(T1-SVP3))
4323 QS=0.622*ES/(P-ES)
4324 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4325 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*QS*(1.+0.81*QS))
4326 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
4327 ! ES=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
4328 ES=611.0*EXP(22.514-6.15E3/T1)
4329 QS=0.622*ES/(P-ES)
4330 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4331 THTGS=T1*PI*EXP((3114.834/T1-0.278296)*QS*(1.+0.81*QS))
4332 ELSE
4333 ! ESLIQ=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
4334 ! ESICE=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
4335 ESLIQ=1.E3*SVP1*EXP(SVP2*(T1-SVPT0)/(T1-SVP3))
4336 ESICE=611.0*EXP(22.514-6.15E3/T1)
4337 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
4338 QS=0.622*ES/(P-ES)
4339 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4340 THTGS=T1*PI*EXP(RL*QS*C5/T1*(1.+0.81*QS))
4341 ENDIF
4342 F1=THTGS-THTU
4343 IF(ABS(F1).LT.0.01)GOTO 50
4344 ITCNT=ITCNT+1
4345 IF(ITCNT.GT.10)GOTO 50
4346 DT=F1*(T1-T0)/(F1-F0)
4347 T0=T1
4348 F0=F1
4349 T1=T1-DT
4350 GOTO 90
4351 !
4352 ! IF THE PARCEL IS SUPERSATURATED, CALCULATE CONCENTRATION OF FRESH
4353 ! CONDENSATE...
4354 !
4355 50 IF(QS.LE.QU)THEN
4356 QNEW=QU-QS
4357 QU=QS
4358 GOTO 96
4359 ENDIF
4360 !
4361 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
4362 ! ADJUSTED...IF LIQUID WATER OR ICE IS PRESENT, IT IS ALLOWED TO EVAPO
4363 ! SUBLIMATE.
4364 !
4365 QNEW=0.
4366 DQ=QS-QU
4367 QTOT=QLIQ+QICE
4368 !
4369 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS
4370 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MI
4371 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURA
4372 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPR
4373 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE
4374 !
4375 !...NOTE THAT THE LIQ AND ICE MAY BE PRESENT IN PROPORTIONS SLIGHTLY DIF
4376 ! THAN SUGGESTED BY THE VALUE OF RATIO2...CHECK TO MAKE SURE THAT LIQ
4377 ! ICE CONCENTRATIONS ARE NOT REDUCED TO BELOW ZERO WHEN EVAPORATION/
4378 ! SUBLIMATION OCCURS...
4379 !
4380 IF(QTOT.GE.DQ)THEN
4381 DQICE=0.0
4382 DQLIQ=0.0
4383 QLIQ=QLIQ-(1.-RATIO2)*DQ
4384 IF(QLIQ.LT.0.)THEN
4385 DQICE=0.0-QLIQ
4386 QLIQ=0.0
4387 ENDIF
4388 QICE=QICE-RATIO2*DQ+DQICE
4389 IF(QICE.LT.0.)THEN
4390 DQLIQ=0.0-QICE
4391 QICE=0.0
4392 ENDIF
4393 QLIQ=QLIQ+DQLIQ
4394 QU=QS
4395 GOTO 96
4396 ELSE
4397 IF(RATIO2.LT.1.E-6)THEN
4398 RLL=XLV0-XLV1*T1
4399 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
4400 RLL=XLS0-XLS1*T1
4401 ELSE
4402 RLL=RL
4403 ENDIF
4404 CP=1005.7*(1.+0.89*QU)
4405 IF(QTOT.LT.1.E-10)THEN
4406 !
4407 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
4408 T1=T1+RLL*(DQ/(1.+DQ))/CP
4409 GOTO 96
4410 ELSE
4411 !
4412 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURA
4413 ! THE TEMPERATURE IS GIVEN BY:
4414 T1=T1+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CP
4415 QU=QU+QTOT
4416 QTOT=0.
4417 ENDIF
4418 QLIQ=0
4419 QICE=0.
4420 ENDIF
4421 96 TU=T1
4422 QNEWLQ=(1.-RATIO2)*QNEW
4423 QNEWIC=RATIO2*QNEW
4424 ! IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPMIX =',
4425 ! + ITCNT
4427 END SUBROUTINE TPMIX
4429 !===============================================================================
4431 SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0,ktau,i,j,nk,tracker)
4432 !
4433 ! Lookup table variables:
4434 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
4435 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
4436 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
4437 ! REAL, SAVE, DIMENSION(1:200) :: ALU
4438 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
4439 ! End of Lookup table variables:
4440 !-----------------------------------------------------------------------
4441 IMPLICIT NONE
4442 !-----------------------------------------------------------------------
4443 REAL, INTENT(IN ) :: P,THES,XLV1,XLV0
4444 REAL, INTENT(OUT ) :: QNEWLQ,QNEWIC
4445 REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
4446 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11, &
4447 TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
4448 INTEGER :: IPTB,ITHTB
4449 INTEGER, INTENT(IN ) :: ktau,i,j,nk,tracker ! aj
4450 !-----------------------------------------------------------------------
4452 !c******** LOOKUP TABLE VARIABLES... ****************************
4453 ! parameter(kfnt=250,kfnp=220)
4454 !c
4455 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
4456 ! * alu(200),rdpr,rdthk,plutop
4457 !C***************************************************************
4458 !c
4459 !c***********************************************************************
4460 !c scaling pressure and tt table index
4461 !c***********************************************************************
4462 !c
4463 tp=(p-plutop)*rdpr
4464 qq=tp-aint(tp)
4465 iptb=int(tp)+1
4467 !
4468 !***********************************************************************
4469 ! base and scaling factor for the
4470 !***********************************************************************
4471 !
4472 ! scaling the and tt table index
4473 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
4474 tth=(thes-bth)*rdthk
4475 pp =tth-aint(tth)
4476 ithtb=int(tth)+1
4477 IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
4478 write(98,'(a,6i5,2f9.2,i3)')'*** OUT OF BOUNDS ***', ktau,i,j,nk,IPTB, ITHTB, &
4479 p/100.0, thes, tracker
4480 flush(98)
4481 ENDIF
4482 !
4483 t00=ttab(ithtb ,iptb )
4484 t10=ttab(ithtb+1,iptb )
4485 t01=ttab(ithtb ,iptb+1)
4486 t11=ttab(ithtb+1,iptb+1)
4487 !
4488 q00=qstab(ithtb ,iptb )
4489 q10=qstab(ithtb+1,iptb )
4490 q01=qstab(ithtb ,iptb+1)
4491 q11=qstab(ithtb+1,iptb+1)
4492 !
4493 !***********************************************************************
4494 ! parcel temperature
4495 !***********************************************************************
4496 !
4497 temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
4498 !
4499 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
4500 !
4501 DQ=QS-QU
4502 IF(DQ.LE.0.)THEN
4503 QNEW=QU-QS
4504 QU=QS
4505 ELSE
4506 !
4507 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
4508 ! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
4509 !
4510 QNEW=0.
4511 QTOT=QLIQ+QICE
4512 !
4513 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
4514 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
4515 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
4516 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
4517 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
4518 !
4519 !...subsaturated values only occur in calculations involving various mixtures of
4520 !...updraft and environmental air for estimation of entrainment and detrainment.
4521 !...For these purposes, assume that reasonable estimates can be given using
4522 !...liquid water saturation calculations only - i.e., ignore the effect of the
4523 !...ice phase in this process only...will not affect conservative properties...
4524 !
4525 IF(QTOT.GE.DQ)THEN
4526 qliq=qliq-dq*qliq/(qtot+1.e-10)
4527 qice=qice-dq*qice/(qtot+1.e-10)
4528 QU=QS
4529 ELSE
4530 RLL=XLV0-XLV1*TEMP
4531 CPP=1004.5*(1.+0.89*QU)
4532 IF(QTOT.LT.1.E-10)THEN
4533 !
4534 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
4535 TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
4536 ELSE
4537 !
4538 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
4539 ! THE TEMPERATURE IS GIVEN BY:
4540 !
4541 TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
4542 QU=QU+QTOT
4543 QTOT=0.
4544 QLIQ=0.
4545 QICE=0.
4546 ENDIF
4547 ENDIF
4548 ENDIF
4549 TU=TEMP
4550 qnewlq=qnew
4551 qnewic=0.
4552 !
4553 END SUBROUTINE TPMIX2
4555 !===============================================================================
4556 SUBROUTINE TPMIXBG(P,THTU,TU,QU,QLIQ,QICE,QNEWLQ,QNEWIC,RATIO2,RL, &
4557 XLV0,XLV1,XLS0,XLS1, &
4558 ! ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE, &
4559 SVP1,SVP2,SVPT0,SVP3)
4560 !-----------------------------------------------------------------------
4561 IMPLICIT NONE
4562 !-----------------------------------------------------------------------
4563 REAL, INTENT(IN ) :: P,THTU,RATIO2,RL, &
4564 XLV0,XLV1,XLS0,XLS1,SVP1,SVP2,SVPT0,SVP3
4565 REAL, INTENT( OUT) :: QNEWLQ,QNEWIC
4566 REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
4567 REAL :: C5, RV, ES, QS, PI, THTGS, ESLIQ, ESICE, F0, T1, T0, &
4568 F1, DT, QNEW, QTOT, DQICE, DQLIQ, RLL, CP, dq
4569 real :: f2, t2
4570 INTEGER :: ITCNT
4571 !-----------------------------------------------------------------------
4572 !
4573 !...THIS SUBROUTINE ITERATIVELY EXTRACTS WET-BULB TEMPERATURE FROM EQUIV
4574 ! POTENTIAL TEMPERATURE, THEN CHECKS TO SEE IF SUFFICIENT MOISTURE IS
4575 ! AVAILABLE TO ACHIEVE SATURATION...IF NOT, TEMPERATURE IS ADJUSTED
4576 ! ACCORDINGLY, IF SO, THE RESIDUAL LIQUID WATER/ICE CONCENTRATION IS
4577 ! DETERMINED...
4578 !
4579 ! Modified BJG 1/17/2014
4580 !
4582 C5=1.0723E-3
4583 RV=461.5
4585 tu = thtu / (1.e5/p) ** 0.2854
4588 !
4589 ! ITERATE TO FIND WET BULB TEMPERATURE AS A FUNCTION OF EQUIVALENT POT
4590 ! TEMP AND PRS, ASSUMING SATURATION VAPOR PRESSURE...RATIO2 IS THE DEG
4591 ! OF GLACIATION...
4592 !
4593 IF(RATIO2.LT.1.E-6)THEN
4594 ! ES=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
4595 ES=1.E3*SVP1*EXP(SVP2*(TU-SVPT0)/(TU-SVP3))
4596 QS=0.622*ES/(P-ES)
4597 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4598 THTGS=TU*PI*EXP((3374.6525/TU-2.5403)*QS*(1.+0.81*QS))
4599 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
4600 ! ES=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
4601 ES=611.0*EXP(22.514-6.15E3/TU)
4602 QS=0.622*ES/(P-ES)
4603 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4604 THTGS=TU*PI*EXP((3114.834/TU-0.278296)*QS*(1.+0.81*QS))
4605 ELSE
4606 ! ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
4607 ! ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
4608 ESLIQ=1.E3*SVP1*EXP(SVP2*(TU-SVPT0)/(TU-SVP3))
4609 ESICE=611.0*EXP(22.514-6.15E3/TU)
4610 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
4611 QS=0.622*ES/(P-ES)
4612 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4613 THTGS=TU*PI*EXP(RL*QS*C5/TU*(1.+0.81*QS))
4614 ENDIF
4615 F0=THTGS-THTU
4616 !c T1=TU-0.5*F0
4617 T0=TU
4618 t1 = t0
4619 f1 = f0
4620 t2 = 50.0
4624 ITCNT=0
4625 90 IF(RATIO2.LT.1.E-6)THEN
4626 ! ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
4627 ES=1.E3*SVP1*EXP(SVP2*(t2-SVPT0)/(t2-SVP3))
4628 QS=0.622*ES/(P-ES)
4629 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4630 THTGS=t2*PI*EXP((3374.6525/t2-2.5403)*QS*(1.+0.81*QS))
4631 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
4632 ! ES=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
4633 ES=611.0*EXP(22.514-6.15E3/t2)
4634 QS=0.622*ES/(P-ES)
4635 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4636 THTGS=t2*PI*EXP((3114.834/t2-0.278296)*QS*(1.+0.81*QS))
4637 ELSE
4638 ! ESLIQ=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
4639 ! ESICE=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
4640 ESLIQ=1.E3*SVP1*EXP(SVP2*(t2-SVPT0)/(t2-SVP3))
4641 ESICE=611.0*EXP(22.514-6.15E3/t2)
4642 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
4643 QS=0.622*ES/(P-ES)
4644 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
4645 THTGS=t2*PI*EXP(RL*QS*C5/t2*(1.+0.81*QS))
4646 ENDIF
4648 f2=THTGS-THTU
4649 if((f1 * f2).lt.0.0) then
4650 f0 = f1
4651 t0 = t1
4652 else
4653 f0 = f0
4654 t0 = t0
4655 endif
4656 f1 = f2
4657 t1 = t2
4659 IF(ABS(F1).LT.0.01)GOTO 50
4660 ITCNT=ITCNT+1
4661 IF(ITCNT.GT.10)GOTO 50
4662 DT=F1*(T1-T0)/(F1-F0)
4664 if(abs(dt).lt.abs( (t1-t0)/2.0 )) then
4665 dt = (t1 - t0) / 2.0
4666 endif
4668 !c T0=T1
4669 !c F0=F1
4670 t2=T1-DT
4671 GOTO 90
4672 !
4673 ! IF THE PARCEL IS SUPERSATURATED, CALCULATE CONCENTRATION OF FRESH
4674 ! CONDENSATE...
4675 !
4676 50 IF(QS.LE.QU)THEN
4677 QNEW=QU-QS
4678 QU=QS
4679 GOTO 96
4680 ENDIF
4681 !
4682 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
4683 ! ADJUSTED...IF LIQUID WATER OR ICE IS PRESENT, IT IS ALLOWED TO EVAPO
4684 ! SUBLIMATE.
4685 !
4686 QNEW=0.
4687 DQ=QS-QU
4688 QTOT=QLIQ+QICE
4689 !
4690 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS
4691 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MI
4692 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURA
4693 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPR
4694 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE
4695 !
4696 !...NOTE THAT THE LIQ AND ICE MAY BE PRESENT IN PROPORTIONS SLIGHTLY DIF
4697 ! THAN SUGGESTED BY THE VALUE OF RATIO2...CHECK TO MAKE SURE THAT LIQ
4698 ! ICE CONCENTRATIONS ARE NOT REDUCED TO BELOW ZERO WHEN EVAPORATION/
4699 ! SUBLIMATION OCCURS...
4700 !
4701 IF(QTOT.GE.DQ)THEN
4702 DQICE=0.0
4703 DQLIQ=0.0
4704 QLIQ=QLIQ-(1.-RATIO2)*DQ
4705 IF(QLIQ.LT.0.)THEN
4706 DQICE=0.0-QLIQ
4707 QLIQ=0.0
4708 ENDIF
4709 QICE=QICE-RATIO2*DQ+DQICE
4710 IF(QICE.LT.0.)THEN
4711 DQLIQ=0.0-QICE
4712 QICE=0.0
4713 ENDIF
4714 QLIQ=QLIQ+DQLIQ
4715 QU=QS
4716 GOTO 96
4717 ELSE
4718 IF(RATIO2.LT.1.E-6)THEN
4719 RLL=XLV0-XLV1*T1
4720 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
4721 RLL=XLS0-XLS1*T1
4722 ELSE
4723 RLL=RL
4724 ENDIF
4725 CP=1005.7*(1.+0.89*QU)
4726 IF(QTOT.LT.1.E-10)THEN
4727 !
4728 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
4729 T1=T1+RLL*(DQ/(1.+DQ))/CP
4730 GOTO 96
4731 ELSE
4732 !
4733 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURA
4734 ! THE TEMPERATURE IS GIVEN BY:
4735 T1=T1+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CP
4736 QU=QU+QTOT
4737 QTOT=0.
4738 ENDIF
4739 QLIQ=0
4740 QICE=0.
4741 ENDIF
4742 96 TU=T1
4743 QNEWLQ=(1.-RATIO2)*QNEW
4744 QNEWIC=RATIO2*QNEW
4745 ! IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPMIX =',
4746 ! + ITCNT
4748 END SUBROUTINE TPMIXBG
4750 !===============================================================================
4751 SUBROUTINE TPMIX2DD(p,thes,ts,qs,ktau,i,j,nk)
4752 !
4753 ! Lookup table variables:
4754 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
4755 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
4756 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
4757 ! REAL, SAVE, DIMENSION(1:200) :: ALU
4758 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
4759 ! End of Lookup table variables:
4760 !-----------------------------------------------------------------------
4761 IMPLICIT NONE
4762 !-----------------------------------------------------------------------
4763 REAL, INTENT(IN ) :: P,THES
4764 REAL, INTENT(INOUT) :: TS,QS
4765 INTEGER, INTENT(IN ) :: ktau,i,j,nk ! avail for debugging aj
4766 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
4767 INTEGER :: IPTB,ITHTB
4768 CHARACTER*256 :: MESS
4769 !-----------------------------------------------------------------------
4771 !
4772 !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
4773 ! parameter(kfnt=250,kfnp=220)
4774 !
4775 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
4776 ! alu(200),rdpr,rdthk,plutop
4777 !***************************************************************
4778 !
4779 !***********************************************************************
4780 ! scaling pressure and tt table index
4781 !***********************************************************************
4782 !
4783 tp=(p-plutop)*rdpr
4784 qq=tp-aint(tp)
4785 iptb=int(tp)+1
4786 !
4787 !***********************************************************************
4788 ! base and scaling factor for the
4789 !***********************************************************************
4790 !
4791 ! scaling the and tt table index
4792 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
4793 tth=(thes-bth)*rdthk
4794 pp =tth-aint(tth)
4795 ithtb=int(tth)+1
4796 IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
4797 write(97,'(a,6i5,2f9.2)')'*** OUT OF BOUNDS ***', ktau,i,j,nk,IPTB, ITHTB, p/100., thes
4798 flush(97)
4799 ENDIF
4800 !
4801 t00=ttab(ithtb ,iptb )
4802 t10=ttab(ithtb+1,iptb )
4803 t01=ttab(ithtb ,iptb+1)
4804 t11=ttab(ithtb+1,iptb+1)
4805 !
4806 q00=qstab(ithtb ,iptb )
4807 q10=qstab(ithtb+1,iptb )
4808 q01=qstab(ithtb ,iptb+1)
4809 q11=qstab(ithtb+1,iptb+1)
4810 !
4811 !***********************************************************************
4812 ! parcel temperature and saturation mixing ratio
4813 !***********************************************************************
4814 !
4815 ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
4816 !
4817 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
4818 !
4819 END SUBROUTINE TPMIX2DD
4821 !===============================================================================
4823 REAL FUNCTION GAMMA(X)
4824 !D DOUBLE PRECISION FUNCTION DGAMMA(X)
4825 !-----------------------------------------------------------------------
4826 IMPLICIT NONE
4827 !-----------------------------------------------------------------------
4828 REAL, INTENT(IN ) :: X
4829 ! REAL, INTENT( OUT) :: GAMMA
4830 !----------------------------------------------------------------------
4831 !
4832 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
4833 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
4834 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
4835 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
4836 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
4837 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
4838 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
4839 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
4840 ! MACHINE-DEPENDENT CONSTANTS.
4841 !
4842 !
4843 !*******************************************************************
4844 !*******************************************************************
4845 !
4846 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
4847 !
4848 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
4849 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
4850 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
4851 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
4852 ! GAMMA(XBIG) = BETA**MAXEXP
4853 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
4854 ! APPROXIMATELY BETA**MAXEXP
4855 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
4856 ! 1.0+EPS .GT. 1.0
4857 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
4858 ! 1/XMININ IS MACHINE REPRESENTABLE
4859 !
4860 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
4861 !
4862 ! BETA MAXEXP XBIG
4863 !
4864 ! CRAY-1 (S.P.) 2 8191 966.961
4865 ! CYBER 180/855
4866 ! UNDER NOS (S.P.) 2 1070 177.803
4867 ! IEEE (IBM/XT,
4868 ! SUN, ETC.) (S.P.) 2 128 35.040
4869 ! IEEE (IBM/XT,
4870 ! SUN, ETC.) (D.P.) 2 1024 171.624
4871 ! IBM 3033 (D.P.) 16 63 57.574
4872 ! VAX D-FORMAT (D.P.) 2 127 34.844
4873 ! VAX G-FORMAT (D.P.) 2 1023 171.489
4874 !
4875 ! XINF EPS XMININ
4876 !
4877 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
4878 ! CYBER 180/855
4879 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
4880 ! IEEE (IBM/XT,
4881 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
4882 ! IEEE (IBM/XT,
4883 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
4884 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
4885 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
4886 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
4887 !
4888 !**************************************************************
4889 !*************************************************************
4890 !
4891 ! ERROR RETURNS
4892 !
4893 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
4894 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
4895 ! TO BE FREE OF UNDERFLOW AND OVERFLOW.
4896 !
4897 !
4898 ! INTRINSIC FUNCTIONS REQUIRED ARE:
4899 !
4900 ! INT, DBLE, EXP, LOG, REAL, SIN
4901 !
4902 !
4903 ! REFERENCES: `AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
4904 ! FUNCTIONS`, W. J. CODY, LECTURE NOTES IN MATHEMATICS,
4905 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
4906 ! (ED.), SPRINGER VERLAG, BERLIN, 1976.
4907 !
4908 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
4909 ! SONS, NEW YORK, 1968.
4910 !
4911 ! LATEST MODIFICATION: OCTOBER 12, 1989
4912 !
4913 ! AUTHORS: W. J. CODY AND L. STOLTZ
4914 ! APPLIED MATHEMATICS DIVISION
4915 ! ARGONNE NATIONAL LABORATORY
4916 ! ARGONNE, IL 60439
4917 !
4918 !------------------------------------------------------
4919 INTEGER I,N
4920 LOGICAL PARITY
4921 REAL :: &
4922 !D DOUBLE PRECISION
4923 C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, &
4924 TWO,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
4925 DIMENSION C(7),P(8),Q(8)
4926 !-------------------------------------------------------------
4927 ! MATHEMATICAL CONSTANTS
4928 !-----------------------------------------------------------
4929 DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/, &
4930 SQRTPI/0.9189385332046727417803297E0/, &
4931 PI/3.1415926535897932384626434E0/
4932 !D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
4933 !D 1 SQRTPI/0.9189385332046727417803297D0/,
4934 !D 2 PI/3.1415926535897932384626434D0/
4935 !----------------------------------------------------------------
4936 ! MACHINE DEPENDENT PARAMETERS
4937 !--------------------------------------------------------------
4938 DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/, &
4939 XINF/3.4E38/
4940 !D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
4941 !D 1 XINF/1.79D308/
4942 !---------------------------------------------------------
4943 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
4944 ! APPROXIMATION OVER (1,2).
4945 !-----------------------------------------------------------
4946 DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
4947 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
4948 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
4949 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
4950 DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
4951 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
4952 2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
4953 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
4954 !D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
4955 !D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
4956 !D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
4957 !D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
4958 !D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
4959 !D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
4960 !D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3,
4961 !D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
4962 !----------------------------------------------------------------------
4963 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
4964 !--------------------------------------------------------------------
4965 DATA C/-1.910444077728E-03,8.4171387781295E-04, &
4966 -5.952379913043012E-04,7.93650793500350248E-04, &
4967 -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
4968 5.7083835261E-03/
4969 !D DATA C/-1.910444077728D-03,8.4171387781295D-04,
4970 !D 1 -5.952379913043012D-04,7.93650793500350248D-04,
4971 !D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02,
4972 !D 3 5.7083835261D-03/
4973 !--------------------------------------------------------------------
4974 ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
4975 !------------------------------------------------------------------
4976 CONV(I) = REAL(I)
4977 !D CONV(I) = DBLE(I)
4978 PARITY=.FALSE.
4979 FACT=ONE
4980 N=0
4981 Y=X
4982 IF(Y.LE.ZERO)THEN
4983 !----------------------------------------------------------
4984 ! ARGUMENT IS NEGATIVE
4985 !--------------------------------------------------------
4986 Y=-X
4987 Y1=AINT(Y)
4988 RES=Y-Y1
4989 IF(RES.NE.ZERO)THEN
4990 IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
4991 FACT=-PI/SIN(PI*RES)
4992 Y=Y+ONE
4993 ELSE
4994 RES=XINF
4995 GOTO 900
4996 ENDIF
4997 ENDIF
4998 !-------------------------------------------
4999 ! ARGUMENT IS POSITIVE
5000 !-----------------------------------------
5001 IF(Y.LT.EPS)THEN
5002 !---------------------------------------
5003 ! ARGUMENT .LT. EPS
5004 !-------------------------------------
5005 IF(Y.GE.XMININ)THEN
5006 RES=ONE/Y
5007 ELSE
5008 RES=XINF
5009 GOTO 900
5010 ENDIF
5011 ELSEIF(Y.LT.TWELVE)THEN
5012 Y1=Y
5013 IF(Y.LT.ONE)THEN
5014 !---------------------------
5015 ! 0.0 .LT. ARGUMENT .LT. 1.0
5016 !--------------------------
5017 Z=Y
5018 Y=Y+ONE
5019 ELSE
5020 !--------------------------
5021 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
5022 !-----------------------------------------------------------
5023 N=INT(Y)-1
5024 Y=Y-CONV(N)
5025 Z=Y-ONE
5026 ENDIF
5027 !------------------------------------------------------
5028 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
5029 !-----------------------------------------------------
5030 XNUM=ZERO
5031 XDEN=ONE
5032 DO 260 I=1,8
5033 XNUM=(XNUM+P(I))*Z
5034 XDEN=XDEN*Z+Q(I)
5035 260 CONTINUE
5036 RES=XNUM/XDEN+ONE
5037 IF(Y1.LT.Y)THEN
5038 !-----------------------------------------------
5039 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
5040 !----------------------------------------------
5041 RES=RES/Y1
5042 ELSEIF(Y1.GT.Y)THEN
5043 !-------------------------------------------------
5044 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
5045 !------------------------------------------------
5046 DO 290 I=1,N
5047 RES=RES*Y
5048 Y=Y+ONE
5049 290 CONTINUE
5050 ENDIF
5051 ELSE
5052 !----------------------------------------------
5053 ! EVALUATE FOR ARGUMENT .GE. 12.0,
5054 !--------------------------------------------
5055 IF(Y.LE.XBIG)THEN
5056 YSQ=Y*Y
5057 SUM=C(7)
5058 DO 350 I=1,6
5059 SUM=SUM/YSQ+C(I)
5060 350 CONTINUE
5061 SUM=SUM/Y-Y+SQRTPI
5062 SUM=SUM+(Y-HALF)*LOG(Y)
5063 RES=EXP(SUM)
5064 ELSE
5065 RES=XINF
5066 GOTO 900
5067 ENDIF
5068 ENDIF
5069 !-----------------------------
5070 ! FINAL ADJUSTMENTS AND RETURN
5071 !------------------------------
5072 IF(PARITY)RES=-RES
5073 IF(FACT.NE.ONE)RES=FACT/RES
5074 900 GAMMA=RES
5075 !D900 DGAMMA = RES
5076 RETURN
5077 ! ---------- LAST LINE OF GAMMA -
5078 END FUNCTION GAMMA
5079 !====================================================================
5080 SUBROUTINE MINIM(ARRAY,KDIM,KS,KEND,KT)
5081 DIMENSION ARRAY(KDIM)
5082 KT=KS
5083 X=ARRAY(KS)
5085 DO K=KS+1,KEND
5086 IF(ARRAY(K).LT.X)THEN
5087 X=ARRAY(K)
5088 KT=K
5089 ENDIF
5090 ENDDO
5092 END SUBROUTINE MINIM
5093 !====================================================================
5094 SUBROUTINE MAXIM(ARRAY,KDIM,KS,KE,MAX)
5095 DIMENSION ARRAY(KDIM)
5096 MAX=KS
5097 X=ARRAY(KS)
5099 DO K=KS,KE
5100 XAR=ARRAY(K)
5101 IF(XAR.GE.X)THEN
5102 X=XAR
5103 MAX=K
5104 ENDIF
5105 ENDDO
5107 END SUBROUTINE MAXIM
5108 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5109 !
5110 ! Interpolating from p surface to sigma surface
5111 subroutine interp1d(pp,p,datain,dataout,kl,km)
5112 !
5113 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5114 IMPLICIT NONE
5115 !---------------------------------
5116 INTEGER, INTENT(IN ) :: kl, km
5117 REAL, dimension(km), INTENT(IN ) :: p, datain
5118 REAL, dimension(kl), INTENT(IN ) :: pp
5119 REAL, dimension(kl), INTENT( OUT) :: dataout
5121 INTEGER :: i, j, k1, k2, kmin
5122 REAL :: pmin, pdif, pmindif
5124 DO i = 1, kl
5125 pmin=2000.
5126 DO j = 1, km
5127 pdif = pp(i) - p(j)
5128 if( abs( pdif ) .le. pmin ) then
5129 pmin = abs( pdif )
5130 pmindif = pdif
5131 kmin = j
5132 endif
5133 ENDDO
5135 if( pmindif .le. 0.0 ) then
5136 k1 = kmin - 1
5137 k2 = kmin
5138 else
5139 k1 = kmin
5140 k2 = kmin + 1
5141 endif
5143 dataout(i) = (datain(k2)-datain(k1))/LOG(p(k2)/p(k1)) &
5144 *LOG(pp(i)/p(k1))+datain(k1)
5146 ENDDO
5148 end subroutine interp1d
5149 !====================================================================
5152 subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
5153 !
5154 ! This subroutine is a lookup table.
5155 ! Given a series of series of saturation equivalent potential
5156 ! temperatures, the temperature is calculated.
5157 !
5158 !--------------------------------------------------------------------
5159 IMPLICIT NONE
5160 !--------------------------------------------------------------------
5161 ! Lookup table variables
5162 ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
5163 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
5164 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
5165 ! REAL, SAVE, DIMENSION(1:200) :: ALU
5166 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
5167 ! End of Lookup table variables
5169 INTEGER :: KP,IT,ITCNT,I
5170 REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
5171 TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
5172 ASTRT,AINC,A1,THTGS
5173 ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
5174 REAL :: ALIQ,BLIQ,CLIQ,DLIQ
5175 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
5176 !
5177 ! equivalent potential temperature increment
5178 data dth/1./
5179 ! minimum starting temp
5180 data tmin/150./
5181 ! tolerance for accuracy of temperature
5182 data toler/0.001/
5183 ! top pressure (pascals)
5184 plutop=5000.0
5185 ! bottom pressure (pascals)
5186 pbot=110000.0
5188 ALIQ = SVP1*1000.
5189 BLIQ = SVP2
5190 CLIQ = SVP2*SVPT0
5191 DLIQ = SVP3
5193 !
5194 ! compute parameters
5195 !
5196 ! 1._over_(sat. equiv. theta increment)
5197 rdthk=1./dth
5198 ! pressure increment
5199 !
5200 DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
5201 ! dpr=(pbot-plutop)/REAL(kfnp-1)
5202 ! 1._over_(pressure increment)
5203 rdpr=1./dpr
5204 ! compute the spread of thes
5205 ! thespd=dth*(kfnt-1)
5206 !
5207 ! calculate the starting sat. equiv. theta
5208 !
5209 temp=tmin
5210 p=plutop-dpr
5211 do kp=1,kfnp
5212 p=p+dpr
5213 es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
5214 qs=0.622*es/(p-es)
5215 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
5216 the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs* &
5217 (1.+0.81*qs))
5218 enddo
5219 !
5220 ! compute temperatures for each sat. equiv. potential temp.
5221 !
5222 p=plutop-dpr
5223 do kp=1,kfnp
5224 thes=the0k(kp)-dth
5225 p=p+dpr
5226 do it=1,kfnt
5227 ! define sat. equiv. pot. temp.
5228 thes=thes+dth
5229 ! iterate to find temperature
5230 ! find initial guess
5231 if(it.eq.1) then
5232 tgues=tmin
5233 else
5234 tgues=ttab(it-1,kp)
5235 endif
5236 es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
5237 qs=0.622*es/(p-es)
5238 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
5239 thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs* &
5240 (1.+0.81*qs))
5241 f0=thgues-thes
5242 t1=tgues-0.5*f0
5243 t0=tgues
5244 itcnt=0
5245 ! iteration loop
5246 do itcnt=1,11
5247 es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
5248 qs=0.622*es/(p-es)
5249 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
5250 thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
5251 f1=thtgs-thes
5252 if(abs(f1).lt.toler)then
5253 exit
5254 endif
5255 ! itcnt=itcnt+1
5256 dt=f1*(t1-t0)/(f1-f0)
5257 t0=t1
5258 f0=f1
5259 t1=t1-dt
5260 enddo
5261 ttab(it,kp)=t1
5262 qstab(it,kp)=qs
5263 enddo
5264 enddo
5265 !
5266 ! lookup table for tlog(emix/aliq)
5267 !
5268 ! set up intial values for lookup tables
5269 !
5270 astrt=1.e-3
5271 ainc=0.075
5272 !
5273 a1=astrt-ainc
5274 do i=1,200
5275 a1=a1+ainc
5276 alu(i)=alog(a1)
5277 enddo
5278 !
5279 END SUBROUTINE KF_LUTAB
5281 !====================================================================
5282 subroutine kf_lutabbg(SVP1,SVP2,SVP3,SVPT0)
5283 !
5284 ! This subroutine is a lookup table.
5285 ! Given a series of series of saturation equivalent potential
5286 ! temperatures, the temperature is calculated.
5287 !
5288 !--------------------------------------------------------------------
5289 IMPLICIT NONE
5290 !--------------------------------------------------------------------
5291 ! Lookup table variables
5292 ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
5293 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
5294 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
5295 ! REAL, SAVE, DIMENSION(1:200) :: ALU
5296 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
5297 ! End of Lookup table variables
5299 INTEGER :: KP,IT,ITCNT,I
5300 REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
5301 TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
5302 ASTRT,AINC,A1,THTGS
5303 real :: t2, f2
5304 ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
5305 REAL :: ALIQ,BLIQ,CLIQ,DLIQ
5306 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
5307 !
5308 ! equivalent potential temperature increment
5309 data dth/1./
5310 ! minimum starting temp
5311 data tmin/150./
5312 ! tolerance for accuracy of temperature
5313 data toler/0.001/
5314 ! top pressure (pascals)
5315 plutop=5000.0
5316 ! bottom pressure (pascals)
5317 pbot=110000.0
5319 ALIQ = SVP1*1000.
5320 BLIQ = SVP2
5321 CLIQ = SVP2*SVPT0
5322 DLIQ = SVP3
5324 !
5325 ! compute parameters
5326 !
5327 ! 1._over_(sat. equiv. theta increment)
5328 rdthk=1./dth
5329 ! pressure increment
5330 !
5331 DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
5332 ! dpr=(pbot-plutop)/REAL(kfnp-1)
5333 ! 1._over_(pressure increment)
5334 rdpr=1./dpr
5335 ! compute the spread of thes
5336 ! thespd=dth*(kfnt-1)
5337 !
5338 ! calculate the starting sat. equiv. theta
5339 !
5340 temp=tmin
5341 p=plutop-dpr
5342 do kp=1,kfnp
5343 p=p+dpr
5344 es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
5345 qs=0.622*es/(p-es)
5346 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
5347 the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs* &
5348 (1.+0.81*qs))
5349 enddo
5350 !
5351 ! compute temperatures for each sat. equiv. potential temp.
5352 !
5353 p=plutop-dpr
5354 do kp=1,kfnp
5355 thes=the0k(kp)-dth
5356 p=p+dpr
5357 do it=1,kfnt
5358 ! define sat. equiv. pot. temp.
5359 thes=thes+dth
5360 ! iterate to find temperature
5361 ! find initial guess
5362 if(it.eq.1) then
5363 tgues=tmin
5364 else
5365 tgues=ttab(it-1,kp)
5366 endif
5367 es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
5368 qs=0.622*es/(p-es)
5369 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
5370 thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs* &
5371 (1.+0.81*qs))
5372 f0=thgues-thes
5373 ! t1=tgues-0.5*f0
5374 t2=tgues-0.5*f0
5375 t0=tgues
5376 t1 = t0
5377 f1 = f0
5380 itcnt=0
5381 ! iteration loop
5382 do itcnt=1,11
5383 ! es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
5384 es=aliq*exp((bliq*t2-cliq)/(t2-dliq))
5385 qs=0.622*es/(p-es)
5386 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
5387 ! thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
5388 thtgs=t2*pi*exp((3374.6525/t2-2.5403)*qs*(1.+0.81*qs))
5390 ! f1=thtgs-thes
5391 f2=thtgs-thes
5393 if((f1 * f2).lt.0.0) then
5394 f0 = f1
5395 t0 = t1
5396 else
5397 f0 = f0
5398 t0 = t0
5399 endif
5400 f1 = f2
5401 t1 = t2
5404 if(abs(f1).lt.toler)then
5405 exit
5406 endif
5407 ! itcnt=itcnt+1
5408 dt=f1*(t1-t0)/(f1-f0)
5410 if(abs(dt).lt.abs( (t1-t0)/2.0 )) then
5411 dt = (t1 - t0) / 2.0
5412 endif
5414 ! t0=t1
5415 ! f0=f1
5416 ! t1=t1-dt
5417 t2=t1-dt
5418 enddo
5419 ttab(it,kp)=t1
5420 qstab(it,kp)=qs
5421 enddo
5422 enddo
5423 !
5424 ! lookup table for tlog(emix/aliq)
5425 !
5426 ! set up intial values for lookup tables
5427 !
5428 astrt=1.e-3
5429 ainc=0.075
5430 !
5431 a1=astrt-ainc
5432 do i=1,200
5433 a1=a1+ainc
5434 alu(i)=alog(a1)
5435 enddo
5436 !
5437 END SUBROUTINE KF_LUTABBG
5438 !====================================================================
5440 END MODULE module_shcu_deng