| blob | 59707e297971a5309a86101491e51aa66cf67804 |
1 MODULE module_cu_kfeta
3 USE module_wrf_error
5 !
6 ! V3.3: A new trigger function is added based Ma and Tan (2009):
7 ! Ma, L.-M. and Z.-M. Tan, 2009: Improving the behavior of
8 ! the cumulus parameterization for tropical cyclone prediction:
9 ! Convection trigger. Atmospheric Research, 92, 190 - 211.
10 !
11 !ckay
12 ! WRF v3.5 with diagnosed deep and shallow KF cloud fraction using
13 ! CAM3-CAM5 methodology, along with captured liquid and ice condensates.
14 ! JAH & KA (U.S. EPA) -- May 2013
15 !
16 !--------------------------------------------------------------------
17 ! Lookup table variables:
18 INTEGER, PARAMETER, PRIVATE :: KFNT=250,KFNP=220!wig, 24-Aug-2006: added private to prevent CuP conflicts !BSINGH - For WRFCuP scheme
19 REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
20 REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
21 REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
22 REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
23 ! Note: KF Lookup table is used by subroutines KF_eta_PARA, TPMIX2,
24 ! TPMIX2DD, ENVIRTHT
25 ! End of Lookup table variables:
27 CONTAINS
29 SUBROUTINE KF_eta_CPS( &
30 ids,ide, jds,jde, kds,kde &
31 ,ims,ime, jms,jme, kms,kme &
32 ,its,ite, jts,jte, kts,kte &
33 ,trigger &
34 ,DT,KTAU,DX,CUDT,ADAPT_STEP_FLAG &
35 ,rho,RAINCV,PRATEC,NCA &
36 ,U,V,TH,T,W,dz8w,Pcps,pi &
37 ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
38 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
39 ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
40 ,QV &
41 ,shall & !BSINGH - For WRFCuP scheme added "shall"
42 ! optionals
43 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
44 ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
45 ,RQICUTEN,RQSCUTEN, RQVFTEN &
46 !ckay
47 ,cldfra_dp_KF,cldfra_sh_KF &
48 ,qc_KF,qi_KF &
49 !kf_edrates
50 ,UDR_KF,DDR_KF &
51 ,UER_KF,DER_KF &
52 ,TIMEC_KF,KF_EDRATES &
53 )
54 !
55 !-------------------------------------------------------------
56 IMPLICIT NONE
57 !-------------------------------------------------------------
58 INTEGER, INTENT(IN ) :: &
59 ids,ide, jds,jde, kds,kde, &
60 ims,ime, jms,jme, kms,kme, &
61 its,ite, jts,jte, kts,kte
63 INTEGER, INTENT(IN ) :: trigger
64 INTEGER, INTENT(IN ) :: STEPCU
65 LOGICAL, INTENT(IN ) :: warm_rain
67 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
68 REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
69 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
71 INTEGER, INTENT(IN ) :: KTAU
73 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
74 INTENT(IN ) :: &
75 U, &
76 V, &
77 W, &
78 TH, &
79 T, &
80 QV, &
81 dz8w, &
82 Pcps, &
83 rho, &
84 pi
85 !
86 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
87 INTENT(INOUT) :: &
88 W0AVG
90 REAL, INTENT(IN ) :: DT, DX
91 REAL, INTENT(IN ) :: CUDT
92 LOGICAL,OPTIONAL,INTENT(IN ) :: ADAPT_STEP_FLAG
93 !
94 REAL, DIMENSION( ims:ime , jms:jme ), &
95 INTENT(INOUT) :: RAINCV
97 REAL, DIMENSION( ims:ime , jms:jme ), &
98 INTENT(INOUT) :: PRATEC
100 REAL, DIMENSION( ims:ime , jms:jme ), &
101 INTENT(INOUT) :: NCA
103 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
104 INTENT(INOUT) :: SHALL !BSINGH - For WRFCuP scheme
106 REAL, DIMENSION( ims:ime , jms:jme ), &
107 INTENT(OUT) :: CUBOT, &
108 CUTOP
110 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
111 INTENT(INOUT) :: CU_ACT_FLAG
113 !
114 ! Optional arguments
115 !
117 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
118 OPTIONAL, &
119 INTENT(INOUT) :: &
120 RTHCUTEN, &
121 RQVCUTEN, &
122 RQCCUTEN, &
123 RQRCUTEN, &
124 RQICUTEN, &
125 RQSCUTEN, &
126 RQVFTEN
127 !
128 ! Flags relating to the optional tendency arrays declared above
129 ! Models that carry the optional tendencies will provdide the
130 ! optional arguments at compile time; these flags all the model
131 ! to determine at run-time whether a particular tracer is in
132 ! use or not.
133 !
134 LOGICAL, OPTIONAL :: &
135 F_QV &
136 ,F_QC &
137 ,F_QR &
138 ,F_QI &
139 ,F_QS
141 !ckay
142 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
143 INTENT(INOUT) :: &
144 cldfra_dp_KF, &
145 cldfra_sh_KF, &
146 qc_KF, &
147 qi_KF
149 !kf_edrates
150 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
151 INTENT(INOUT) :: &
152 UDR_KF, &
153 DDR_KF, &
154 UER_KF, &
155 DER_KF
157 REAL, DIMENSION( ims:ime , jms:jme ) , &
158 INTENT(INOUT) :: &
159 TIMEC_KF
161 INTEGER, INTENT(IN) :: KF_EDRATES
163 ! LOCAL VARS
165 LOGICAL :: flag_qr, flag_qi, flag_qs
167 REAL, DIMENSION( kts:kte ) :: &
168 U1D, &
169 V1D, &
170 T1D, &
171 DZ1D, &
172 QV1D, &
173 P1D, &
174 RHO1D, &
175 tpart_v1D, &
176 tpart_h1D, &
177 W0AVG1D
179 REAL, DIMENSION( kts:kte ):: &
180 DQDT, &
181 DQIDT, &
182 DQCDT, &
183 DQRDT, &
184 DQSDT, &
185 DTDT
187 REAL, DIMENSION (its-1:ite+1,kts:kte,jts-1:jte+1) :: aveh_t, aveh_q
188 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: aveh_qmax, aveh_qmin
189 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_t, avev_q
190 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_qmax, avev_qmin
191 REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: coef_v, coef_h, tpart_h, tpart_v
192 INTEGER :: ii,jj,kk
194 REAL :: ttop
195 REAL, DIMENSION (kts:kte) :: z0
197 REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
198 integer :: ibegh,iendh,jbegh,jendh
199 integer :: istart,iend,jstart,jend
200 INTEGER :: i,j,k,NTST
201 INTEGER :: ishall !BSINGH - For WRFCuP Scheme
202 REAL :: lastdt = -1.0
203 REAL :: W0AVGfctr, W0fctr, W0den
205 !
206 DXSQ=DX*DX
208 !----------------------
209 NTST=STEPCU
210 TST=float(NTST*2)
211 flag_qr = .FALSE.
212 flag_qi = .FALSE.
213 flag_qs = .FALSE.
214 IF ( PRESENT(F_QR) ) flag_qr = F_QR
215 IF ( PRESENT(F_QI) ) flag_qi = F_QI
216 IF ( PRESENT(F_QS) ) flag_qs = F_QS
217 !
218 if (lastdt < 0) then
219 lastdt = dt
220 endif
222 if (ADAPT_STEP_FLAG) then
223 W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
224 W0fctr = dt
225 W0den = 2 * MAX(CUDT*60,dt)
226 else
227 W0AVGfctr = (TST-1.)
228 W0fctr = 1.
229 W0den = TST
230 endif
232 DO J = jts,jte
233 DO K=kts,kte
234 DO I= its,ite
235 ! SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
236 ! TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
237 ! RHOE=Pcps(I,K,J)/(R*TV)
238 ! W0=-101.9368*SCR1/RHOE
239 W0=0.5*(w(I,K,J)+w(I,K+1,J))
241 ! Old:
242 !
243 ! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
244 !
245 ! New, to support adaptive time step:
246 !
247 W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
248 ENDDO
249 ENDDO
250 ENDDO
251 lastdt = dt
253 ! New trigger function
254 IF (trigger.eq.2) THEN
255 !
256 ! calculate 9-point average of moisture advection and temperature using halo (Horizontal)
257 !
258 aveh_t=-999 ! horizontal 9-point ave
259 aveh_q=-999
260 avev_t=0 ! vertical 3-level ave
261 avev_q=0
262 avev_qmax=0
263 avev_qmin=0
264 aveh_qmax=0
265 aveh_qmin=0
266 tpart_h=0
267 tpart_v=0
268 coef_h=0
269 coef_v=0
270 ibegh=max(its-1, ids+1) ! start from 2
271 jbegh=max(jts-1, jds+1)
272 iendh=min(ite+1, ide-2) ! end at ide-2
273 jendh=min(jte+1, jde-2)
274 DO J = jbegh,jendh
275 DO K = kts,kte
276 DO I = ibegh,iendh
277 aveh_t(i,k,j)=(T(i-1,k,j-1)+T(i-1,k,j) +T(i-1,k,j+1)+ &
278 T(i,k,j-1) +T(i,k,j) +T(i,k,j+1)+ &
279 T(i+1,k,j-1) +T(i+1,k,j) +T(i+1,k,j+1))/9.
280 aveh_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) +rqvften(i-1,k,j+1)+ &
281 rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
282 rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
283 ENDDO
284 ENDDO
285 ENDDO
286 ! boundary value ( all processors will do the following? Or just those processsors handling sub-area including boundary)
287 DO K = kts,kte
288 DO J = jts-1,jte+1
289 DO I = its-1,ite+1
291 if(i.eq.ids) then
292 aveh_t(i,k,j)=aveh_t(i+1,k,j)
293 aveh_q(i,k,j)=aveh_q(i+1,k,j)
294 elseif(i.eq.ide-1) then
295 aveh_t(i,k,j)=aveh_t(i-1,k,j)
296 aveh_q(i,k,j)=aveh_q(i-1,k,j)
297 endif
299 if(j.eq.jds) then
300 aveh_t(i,k,j)=aveh_t(i,k,j+1)
301 aveh_q(i,k,j)=aveh_q(i,k,j+1)
302 elseif(j.eq.jde-1) then
303 aveh_t(i,k,j)=aveh_t(i,k,j-1)
304 aveh_q(i,k,j)=aveh_q(i,k,j-1)
305 endif
307 if(j.eq.jds.and.i.eq.ids) then
308 aveh_q(i,k,j)=aveh_q(i+1,k,j+1)
309 aveh_t(i,k,j)=aveh_t(i+1,k,j+1)
310 endif
312 if(j.eq.jde-1.and.i.eq.ids) then
313 aveh_q(i,k,j)=aveh_q(i+1,k,j-1)
314 aveh_t(i,k,j)=aveh_t(i+1,k,j-1)
315 endif
317 if(j.eq.jde-1.and.i.eq.ide-1) then
318 aveh_q(i,k,j)=aveh_q(i-1,k,j-1)
319 aveh_t(i,k,j)=aveh_t(i-1,k,j-1)
320 endif
322 if(j.eq.jds.and.i.eq.ide-1) then
323 aveh_q(i,k,j)=aveh_q(i-1,k,j+1)
324 aveh_t(i,k,j)=aveh_t(i-1,k,j+1)
325 endif
327 ENDDO
328 ENDDO
329 ENDDO
330 ! search for max/min moisture advection in 9-point range, calculate horizontal T-perturbation (tpart_h)
331 istart=max(its, ids+1) ! start from 2
332 jstart=max(jts, jds+1)
333 iend=min(ite, ide-2) ! end at ide-2
334 jend=min(jte, jde-2)
335 DO K = kts,kte
336 DO J = jstart,jend
337 DO I = istart,iend
338 aveh_qmax(i,k,j)=aveh_q(i,k,j)
339 aveh_qmin(i,k,j)=aveh_q(i,k,j)
340 DO ii=-1, 1
341 DO jj=-1,1
342 if(aveh_q(i+II,k,j+JJ).gt.aveh_qmax(i,k,j)) aveh_qmax(i,k,j)=aveh_q(i+II,k,j+JJ)
343 if(aveh_q(i+II,k,j+JJ).lt.aveh_qmin(i,k,j)) aveh_qmin(i,k,j)=aveh_q(i+II,k,j+JJ)
344 ENDDO
345 ENDDO
346 if(aveh_qmax(i,k,j).gt.aveh_qmin(i,k,j))then
347 coef_h(i,k,j)=(aveh_q(i,k,j)-aveh_qmin(i,k,j))/(aveh_qmax(i,k,j)-aveh_qmin(i,k,j))
348 else
349 coef_h(i,k,j)=0.
350 endif
351 coef_h(i,k,j)=amin1(coef_h(i,k,j),1.0)
352 coef_h(i,k,j)=amax1(coef_h(i,k,j),0.0)
353 tpart_h(i,k,j)=coef_h(i,k,j)*(T(i,k,j)-aveh_t(i,k,j))
354 ENDDO
355 ENDDO
356 ENDDO
357 89 continue
358 ! vertical 3-layer calculation
359 DO J = jts, jte
360 DO I = its, ite
361 z0(1) = 0.5 * dz8w(i,1,j)
362 DO K = 2, kte
363 Z0(K) = Z0(K-1) + .5 * (DZ8W(i,K,j) + DZ8W(i,K-1,j))
364 ENDDO
365 DO K = kts+1,kte-1
366 ttop = t(i,k,j) + ((t(i,k,j) - t(i,k+1,j)) / (z0(k) - z0(k+1))) * (z0(k)-z0(k-1))
367 avev_t(i,k,j)=(T(i,k-1,j) + T(i,k,j) + ttop)/3.
368 ! avev_t(i,k,j)=(T(i,k-1,j)+T(i,k,j) + T(i,k+1,j))/3.
369 avev_q(i,k,j)=(rqvften(i,k-1,j)+rqvften(i,k,j) + rqvften(i,k+1,j))/3.
370 ENDDO
371 avev_t(i,kts,j)=avev_t(i,kts+1,j) ! lowest level value, is it the same as avev_t(i,kds,j)=avev_t(i,kds+1,j)?
372 avev_q(i,kts,j)=avev_q(i,kts+1,j)
373 avev_t(i,kte,j)=avev_t(i,kte-1,j) ! highest level value
374 avev_q(i,kte,j)=avev_q(i,kte-1,j)
375 ENDDO
376 ENDDO
377 ! max /min value
378 DO J = jts, jte
379 DO I = its, ite
380 DO K = kts+1,kte-1
381 avev_qmax(i,k,j)=avev_q(i,k,j)
382 avev_qmin(i,k,j)=avev_q(i,k,j)
383 DO kk=-1,1
384 if(avev_q(i,k+kk,j).gt.avev_qmax(i,k,j)) avev_qmax(i,k,j)=avev_q(i,k+kk,j)
385 if(avev_q(i,k+kk,j).lt.avev_qmin(i,k,j)) avev_qmin(i,k,j)=avev_q(i,k+kk,j)
386 ENDDO
387 if(avev_qmax(i,k,j).gt.avev_qmin(i,k,j)) then
388 coef_v(i,k,j)=(avev_q(i,k,j)-avev_qmin(i,k,j))/(avev_qmax(i,k,j)-avev_qmin(i,k,j))
389 else
390 coef_v(i,k,j)=0
391 endif
392 tpart_v(i,k,j)=coef_v(i,k,j)*(T(i,k,j)-avev_t(i,k,j))
393 ENDDO
394 tpart_v(i,kts,j)= tpart_v(i,kts+1,j) ! lowest level
395 tpart_v(i,kte,j)= tpart_v(i,kte-1,j) ! highest level
396 ENDDO
397 ENDDO
398 ENDIF ! endif (trigger.eq.2)
399 !
400 DO J = jts,jte
401 DO I= its,ite
402 CU_ACT_FLAG(i,j) = .true.
403 ENDDO
404 ENDDO
406 DO J = jts,jte
407 DO I=its,ite
410 IF ( NCA(I,J) .ge. 0.5*DT ) then
411 CU_ACT_FLAG(i,j) = .false.
412 ELSE
414 DO k=kts,kte
415 DQDT(k)=0.
416 DQIDT(k)=0.
417 DQCDT(k)=0.
418 DQRDT(k)=0.
419 DQSDT(k)=0.
420 DTDT(k)=0.
421 !ckay
422 cldfra_dp_KF(I,k,J)=0.
423 cldfra_sh_KF(I,k,J)=0.
424 qc_KF(I,k,J)=0.
425 qi_KF(I,k,J)=0.
426 ENDDO
427 !kf_edrates
428 IF (KF_EDRATES == 1) THEN
429 DO k=kts,kte
430 UDR_KF(I,k,J)=0.
431 DDR_KF(I,k,J)=0.
432 UER_KF(I,k,J)=0.
433 DER_KF(I,k,J)=0.
434 ENDDO
435 TIMEC_KF(I,J)=0.
436 ENDIF
437 RAINCV(I,J)=0.
438 CUTOP(I,J)=KTS
439 CUBOT(I,J)=KTE+1
440 PRATEC(I,J)=0.
441 !
442 ! assign vars from 3D to 1D
444 DO K=kts,kte
445 U1D(K) =U(I,K,J)
446 V1D(K) =V(I,K,J)
447 T1D(K) =T(I,K,J)
448 RHO1D(K) =rho(I,K,J)
449 QV1D(K)=QV(I,K,J)
450 P1D(K) =Pcps(I,K,J)
451 W0AVG1D(K) =W0AVG(I,K,J)
452 DZ1D(k)=dz8w(I,K,J)
453 IF (trigger.eq.2) THEN
454 tpart_h1D(K) =tpart_h(I,K,J)
455 tpart_v1D(K) =tpart_v(I,K,J)
456 ELSE
457 tpart_h1D(K) = 0.
458 tpart_v1D(K) = 0.
459 ENDIF
460 ENDDO
461 CALL KF_eta_PARA(I, J, &
462 U1D,V1D,T1D,QV1D,P1D,DZ1D,W0AVG1D, &
463 tpart_h1D,tpart_v1D, &
464 trigger, &
465 DT,DX,DXSQ,RHO1D, &
466 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
467 EP2,SVP1,SVP2,SVP3,SVPT0, &
468 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
469 RAINCV,PRATEC,NCA, &
470 flag_QI,flag_QS,warm_rain, &
471 CUTOP,CUBOT,ishall,CUDT, & !BSINGH - Added ishall as arg for WRFCuP
472 ids,ide, jds,jde, kds,kde, &
473 ims,ime, jms,jme, kms,kme, &
474 its,ite, jts,jte, kts,kte, &
475 !ckay
476 cldfra_dp_KF,cldfra_sh_KF, &
477 qc_KF,qi_KF, &
478 !kf_edrates
479 UDR_KF,DDR_KF, &
480 UER_KF,DER_KF, &
481 TIMEC_KF,KF_EDRATES )
482 if(present(shall))shall(i,j) = ishall !,BSINGH - For WRFCuP scheme
485 IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
486 DO K=kts,kte
487 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
488 RQVCUTEN(I,K,J)=DQDT(K)
489 ENDDO
490 ENDIF
492 IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
493 IF( F_QR )THEN
494 DO K=kts,kte
495 RQRCUTEN(I,K,J)=DQRDT(K)
496 RQCCUTEN(I,K,J)=DQCDT(K)
497 ENDDO
498 ELSE
499 ! This is the case for Eta microphysics without 3d rain field
500 DO K=kts,kte
501 RQRCUTEN(I,K,J)=0.
502 RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
503 ENDDO
504 ENDIF
505 ENDIF
507 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
510 IF(PRESENT( rqicuten )) THEN
511 IF ( F_QI ) THEN
512 DO K=kts,kte
513 RQICUTEN(I,K,J)=DQIDT(K)
514 ENDDO
515 ENDIF
516 ENDIF
518 IF(PRESENT( rqscuten )) THEN
519 IF ( F_QS ) THEN
520 DO K=kts,kte
521 RQSCUTEN(I,K,J)=DQSDT(K)
522 ENDDO
523 ENDIF
524 ENDIF
525 !
526 ENDIF
527 ENDDO ! i-loop
528 ENDDO ! j-loop
529 !
530 END SUBROUTINE KF_eta_CPS
531 ! ****************************************************************************
532 !-----------------------------------------------------------
533 SUBROUTINE KF_eta_PARA (I, J, &
534 U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
535 TPART_H0,TPART_V0, &
536 trigger, &
537 DT,DX,DXSQ,rhoe, &
538 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
539 EP2,SVP1,SVP2,SVP3,SVPT0, &
540 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
541 RAINCV,PRATEC,NCA, &
542 F_QI,F_QS,warm_rain, &
543 CUTOP,CUBOT,ishall,CUDT, & !BSINGH - For WRFCuP scheme
544 ids,ide, jds,jde, kds,kde, &
545 ims,ime, jms,jme, kms,kme, &
546 its,ite, jts,jte, kts,kte, &
547 !ckay
548 cldfra_dp_KF,cldfra_sh_KF, &
549 qc_KF,qi_KF, &
550 !kf_edrates
551 UDR_KF,DDR_KF, &
552 UER_KF,DER_KF, &
553 TIMEC_KF,KF_EDRATES )
554 !-----------------------------------------------------------
555 !***** The KF scheme that is currently used in experimental runs of EMCs
556 !***** Eta model....jsk 8/00
557 !
558 IMPLICIT NONE
559 !-----------------------------------------------------------
560 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
561 ims,ime, jms,jme, kms,kme, &
562 its,ite, jts,jte, kts,kte, &
563 I,J
564 ! ,P_QI,P_QS,P_FIRST_SCALAR
565 INTEGER, INTENT(IN ) :: trigger
567 LOGICAL, INTENT(IN ) :: F_QI, F_QS
569 LOGICAL, INTENT(IN ) :: warm_rain
570 !
571 REAL, DIMENSION( kts:kte ), &
572 INTENT(IN ) :: U0, &
573 V0, &
574 TPART_H0, &
575 TPART_V0, &
576 T0, &
577 QV0, &
578 P0, &
579 rhoe, &
580 DZQ, &
581 W0AVG1D
582 !
583 REAL, INTENT(IN ) :: DT,DX,DXSQ
584 !
586 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
587 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
589 !
590 INTEGER, INTENT(INOUT) :: ISHALL
591 REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
592 DQDT, &
593 DQIDT, &
594 DQCDT, &
595 DQRDT, &
596 DQSDT, &
597 DTDT
599 REAL, DIMENSION( ims:ime , jms:jme ), &
600 INTENT(INOUT) :: NCA
602 !ckay
603 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
604 INTENT(INOUT) :: cldfra_dp_KF, &
605 cldfra_sh_KF, &
606 qc_KF, &
607 qi_KF
608 !kf_edrates
609 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
610 INTENT(INOUT) :: UDR_KF, &
611 DDR_KF, &
612 UER_KF, &
613 DER_KF
615 REAL, DIMENSION( ims:ime , jms:jme ), &
616 INTENT(INOUT) :: TIMEC_KF
618 INTEGER, INTENT(IN) :: KF_EDRATES
620 REAL, DIMENSION( ims:ime , jms:jme ), &
621 INTENT(INOUT) :: RAINCV
623 REAL, DIMENSION( ims:ime , jms:jme ), &
624 INTENT(INOUT) :: PRATEC
626 REAL, DIMENSION( ims:ime , jms:jme ), &
627 INTENT(OUT) :: CUBOT, &
628 CUTOP
629 REAL, INTENT(IN ) :: CUDT
630 !
631 !...DEFINE LOCAL VARIABLES...
632 !
633 REAL, DIMENSION( kts:kte ) :: &
634 Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
635 QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
636 UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
637 UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
638 THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
639 QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
640 DETLQ2,DETIC2,RATIO,RATIO2
643 REAL, DIMENSION( kts:kte ) :: &
644 DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD, &
645 QDT,FXM,THTAG,THPA,THFXOUT, &
646 THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN, &
647 QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
648 QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
649 QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
652 REAL, DIMENSION( kts:kte+1 ) :: OMG
653 REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
654 REAL, DIMENSION( kts:kte ) :: &
655 CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
657 ! LOCAL VARS
659 REAL :: P00,T00,RLF,RHIC,RHBC,PIE, &
660 TTFRZ,TBFRZ,C5,RATE
661 REAL :: GDRY,ROCP,ALIQ,BLIQ, &
662 CLIQ,DLIQ
663 REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
664 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
665 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
666 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
667 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
668 UPNEW,ABE,WKLCL,TTEMP,FRC1, &
669 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
670 DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
671 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
672 UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
673 THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
674 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
675 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
676 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
677 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
678 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
679 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
680 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
681 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
682 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
683 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
684 DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
685 REAL :: ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
686 QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
687 !
688 INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
689 REAL :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
690 REAL :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
691 !ckay
692 REAL :: xcldfra,UMF_new
694 INTEGER :: KX,K,KL
695 !
696 INTEGER :: NCHECK
697 INTEGER, DIMENSION (kts:kte) :: KCHECK
699 INTEGER :: ISTOP,ML,L5,KMIX,LOW, &
700 LC,MXLAYR,LLFC,NLAYRS,NK, &
701 KPBL,KLCL,LCL,LET,IFLAG, &
702 NK1,LTOP,NJ,LTOP1, &
703 LTOPM1,LVF,KSTART,KMIN,LFS, &
704 ND,NIC,LDB,LDT,ND1,NDK, &
705 NM,LMAX,NCOUNT,NOITR, &
706 NSTEP,NTC,NCHM,NSHALL
707 LOGICAL :: IPRNT
708 REAL :: u00,qslcl,rhlcl,dqssdt !jfb
709 CHARACTER*1024 message
710 !
711 DATA P00,T00/1.E5,273.16/
712 DATA RLF/3.339E5/
713 DATA RHIC,RHBC/1.,0.90/
714 DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
715 DATA RATE/0.03/ ! wrf default
716 ! DATA RATE/0.01/ ! value used in NRCM
717 ! DATA RATE/0.001/ ! effectively turn off autoconversion
718 !-----------------------------------------------------------
719 IPRNT=.FALSE.
720 GDRY=-G/CP
721 ROCP=R/CP
722 NSHALL = 0
723 KL=kte
724 KX=kte
725 !
726 ! ALIQ = 613.3
727 ! BLIQ = 17.502
728 ! CLIQ = 4780.8
729 ! DLIQ = 32.19
730 ALIQ = SVP1*1000.
731 BLIQ = SVP2
732 CLIQ = SVP2*SVPT0
733 DLIQ = SVP3
734 !
735 !
736 !****************************************************************************
737 ! ! PPT FB MODS
738 !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER ! PPT FB MODS
739 !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL) ! PPT FB MODS
740 !...FIELD. "FBFRC" IS THE FRACTION OF AVAILABLE ! PPT FB MODS
741 !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)... ! PPT FB MODS
742 FBFRC=0.0 ! PPT FB MODS
743 !...mods to allow shallow convection...
744 NCHM = 0
745 ISHALL = 0
746 DPMIN = 5.E3
747 !...
748 P300=P0(1)-30000.
749 !
750 !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
751 !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
752 !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
753 !
754 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
755 !...FROM BOTTOM-UP IN THE KF SCHEME...
756 !
757 ML=0
758 !SUE tmprpsb=1./PSB(I,J)
759 !SUE CELL=PTOP*tmprpsb
760 !
761 DO K=1,KX
762 !
763 ! Saturation vapor pressure (ES) is calculated following Buck (1981)
764 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
765 !
766 ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
767 QES(K)=0.622*ES/(P0(K)-ES)
768 Q0(K)=AMIN1(QES(K),QV0(K))
769 Q0(K)=AMAX1(0.000001,Q0(K))
770 QL0(K)=0.
771 QI0(K)=0.
772 QR0(K)=0.
773 QS0(K)=0.
774 RH(K) = Q0(K)/QES(K)
775 DILFRC(K) = 1.
776 TV0(K)=T0(K)*(1.+0.608*Q0(K))
777 ! RHOE(K)=P0(K)/(R*TV0(K))
778 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
779 DP(K)=rhoe(k)*g*DZQ(k)
780 ! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
781 ! use it for shallow convection...For now, assume it is not available....
782 ! TKE(K) = Q2(I,J,NK)
783 TKE(K) = 0.
784 CLDHGT(K) = 0.
785 ! IF(P0(K).GE.500E2)L5=K
786 IF(P0(K).GE.0.5*P0(1))L5=K
787 IF(P0(K).GE.P300)LLFC=K
788 ENDDO
789 !
790 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
791 Z0(1)=.5*DZQ(1)
792 !cdir novector
793 DO K=2,KL
794 Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
795 DZA(K-1)=Z0(K)-Z0(K-1)
796 ENDDO
797 DZA(KL)=0.
798 !
799 !
800 ! To save time, specify a pressure interval to move up in sequential
801 ! check of different ~50 mb deep groups of adjacent model layers in
802 ! the process of identifying updraft source layer (USL). Note that
803 ! this search is terminated as soon as a buoyant parcel is found and
804 ! this parcel can produce a cloud greater than specifed minimum depth
805 ! (CHMIN)...For now, set interval at 15 mb...
806 !
807 NCHECK = 1
808 KCHECK(NCHECK)=1
809 PM15 = P0(1)-15.E2
810 DO K=2,LLFC
811 IF(P0(K).LT.PM15)THEN
812 NCHECK = NCHECK+1
813 KCHECK(NCHECK) = K
814 PM15 = PM15-15.E2
815 ENDIF
816 ENDDO
817 !
818 NU=0
819 NUCHM=0
820 usl: DO
821 NU = NU+1
822 IF(NU.GT.NCHECK)THEN
823 IF(ISHALL.EQ.1)THEN
824 CHMAX = 0.
825 NCHM = 0
826 DO NK = 1,NCHECK
827 NNN=KCHECK(NK)
828 IF(CLDHGT(NNN).GT.CHMAX)THEN
829 NCHM = NNN
830 NUCHM = NK
831 CHMAX = CLDHGT(NNN)
832 ENDIF
833 ENDDO
834 NU = NUCHM-1
835 FBFRC=1.
836 CYCLE usl
837 ELSE
838 !BSINGH - For WRFCuP scheme
839 ! wig, 29-Aug-2006: I think this is where no convecion occurs. So, force
840 ! ishall to a flag value to indicate this for accounting purposes.
841 ishall = 2
842 !BSINGH -ENDS
843 RETURN
844 ENDIF
845 ENDIF
846 KMIX = KCHECK(NU)
847 LOW=KMIX
848 !...
849 LC = LOW
850 !
851 !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
852 !...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
853 !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
854 !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
855 !
856 NLAYRS=0
857 DPTHMX=0.
858 NK=LC-1
859 IF ( NK+1 .LT. KTS ) THEN
860 WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
861 CALL wrf_message (TRIM(message))
862 ELSE
863 DO
864 NK=NK+1
865 IF ( NK .GT. KTE ) THEN
866 WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
867 CALL wrf_message (TRIM(message))
868 EXIT
869 ENDIF
870 DPTHMX=DPTHMX+DP(NK)
871 NLAYRS=NLAYRS+1
872 IF(DPTHMX.GT.DPMIN)THEN
873 EXIT
874 ENDIF
875 END DO
876 ENDIF
877 IF(DPTHMX.LT.DPMIN)THEN
878 !BSINGH - For WRFCuP scheme
879 ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
880 ishall = 2
881 !BSINGH -ENDS
882 RETURN
883 ENDIF
884 KPBL=LC+NLAYRS-1
885 !
886 !...********************************************************
887 !...for computational simplicity without much loss in accuracy,
888 !...mix temperature instead of theta for evaluating convective
889 !...initiation (triggering) potential...
890 ! THMIX=0.
891 TMIX=0.
892 QMIX=0.
893 ZMIX=0.
894 PMIX=0.
895 !
896 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
897 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
898 !...LAYERS...
899 !
900 !cdir novector
901 DO NK=LC,KPBL
902 TMIX=TMIX+DP(NK)*T0(NK)
903 QMIX=QMIX+DP(NK)*Q0(NK)
904 ZMIX=ZMIX+DP(NK)*Z0(NK)
905 PMIX=PMIX+DP(NK)*P0(NK)
906 ENDDO
907 ! THMIX=THMIX/DPTHMX
908 TMIX=TMIX/DPTHMX
909 QMIX=QMIX/DPTHMX
910 ZMIX=ZMIX/DPTHMX
911 PMIX=PMIX/DPTHMX
912 EMIX=QMIX*PMIX/(0.622+QMIX)
913 !
914 !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
915 !
916 ! TLOG=ALOG(EMIX/ALIQ)
917 ! ...calculate dewpoint using lookup table...
918 !
919 astrt=1.e-3
920 ainc=0.075
921 a1=emix/aliq
922 tp=(a1-astrt)/ainc
923 indlu=int(tp)+1
924 value=(indlu-1)*ainc+astrt
925 aintrp=(a1-value)/ainc
926 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
927 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
928 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
929 TLCL=AMIN1(TLCL,TMIX)
930 TVLCL=TLCL*(1.+0.608*QMIX)
931 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
932 ! NK = LC-1
933 ! DO
934 ! NK = NK+1
935 ! KLCL=NK
936 ! IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
937 ! EXIT
938 ! ENDIF
939 ! ENDDO
940 ! IF(NK.GT.KL)THEN
941 ! RETURN
942 ! ENDIF
944 DO NK = LC, KL
945 KLCL = NK
946 IF ( ZLCL.LE.Z0(NK) ) EXIT
947 END DO
948 IF ( ZLCL.GT.Z0(KL) ) then
949 !BSINGH - For WRFCuP scheme
950 !BSINGH - Improvised based on the following reason
951 !Before every "RETURN" statement assign ishall=2
952 ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
953 ishall = 2
954 !BSINGH -ENDS
955 RETURN
956 endif
958 K=KLCL-1
959 ! calculate DLP using Z instead of log(P)
960 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
961 !
962 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
963 !
964 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
965 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
966 TVEN=TENV*(1.+0.608*QENV)
967 !
968 !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
969 !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
970 !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
971 !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
972 !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
973 !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
974 !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
975 !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
976 !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
977 !
978 IF(ZLCL.LT.2.E3)THEN ! Kain (2004) Eq. 2
979 WKLCL=0.02*ZLCL/2.E3
980 ELSE
981 WKLCL=0.02 ! units of m/s
982 ENDIF
983 WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3-WKLCL
984 IF(WKL.LT.0.0001)THEN
985 DTLCL=0.
986 ELSE
987 DTLCL=4.64*WKL**0.33 ! Kain (2004) Eq. 1
988 ENDIF
990 ! New trigger function
991 IF(trigger.eq.2) then
992 DTLCL=amax1(TPART_H0(KLCL)+TPART_V0(KLCL),0.0)
993 ENDIF
994 ! end for trigger function
995 !
996 dtrh = 0.
997 if (trigger .eq. 3) then
998 !...for ETA model, give parcel an extra temperature perturbation based
999 !...the threshold RH for condensation (U00)...
1000 ! as described in Narita and Ohmori (2007, 12th Mesoscale Conf.
1001 !
1002 !...for now, just assume U00=0.75...
1003 !...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
1004 U00 = 0.75
1005 IF(U00.lt.1.)THEN
1006 QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
1007 RHLCL = QENV/QSLCL
1008 DQSSDT = QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
1009 IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
1010 DTRH = 0.25*(RHLCL-0.75)*QMIX/DQSSDT
1011 ELSEIF(RHLCL.GT.0.95)THEN
1012 DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
1013 ELSE
1014 DTRH = 0.
1015 ENDIF
1016 ENDIF
1017 endif ! trigger 3
1018 ! IF(ISHALL.EQ.1)IPRNT=.TRUE.
1019 ! IPRNT=.TRUE.
1020 ! IF(TLCL+DTLCL.GT.TENV)GOTO 45
1022 trigger2: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
1023 !
1024 ! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
1025 !
1026 CYCLE usl
1027 !
1028 ELSE ! Parcel is buoyant, determine updraft
1029 !
1030 !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
1031 !...EQUIVALENT POTENTIAL TEMPERATURE
1032 !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
1033 !
1034 CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
1035 !
1036 !...modify calculation of initial parcel vertical velocity...jsk 11/26/97
1037 !
1038 DTTOT = DTLCL+DTRH
1039 IF(DTTOT.GT.1.E-4)THEN
1040 GDT=2.*G*DTTOT*500./TVEN ! Kain (2004) Eq. 3 (sort of)
1041 WLCL=1.+0.5*SQRT(GDT)
1042 WLCL = AMIN1(WLCL,3.)
1043 ELSE
1044 WLCL=1.
1045 ENDIF
1046 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
1047 WTW=WLCL*WLCL
1048 !
1049 TVLCL=TLCL*(1.+0.608*QMIX)
1050 RHOLCL=PLCL/(R*TVLCL)
1051 !
1052 LCL=KLCL
1053 LET=LCL
1054 ! make RAD a function of background vertical velocity... (Kain (2004) Eq. 6)
1055 IF(WKL.LT.0.)THEN
1056 RAD = 1000.
1057 ELSEIF(WKL.GT.0.1)THEN
1058 RAD = 2000.
1059 ELSE
1060 RAD = 1000.+1000*WKL/0.1
1061 ENDIF
1062 !
1063 !*******************************************************************
1064 ! *
1065 ! COMPUTE UPDRAFT PROPERTIES *
1066 ! *
1067 !*******************************************************************
1068 !
1069 !
1070 !...
1071 !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
1072 !
1073 WU(K)=WLCL
1074 AU0=0.01*DXSQ
1075 UMF(K)=RHOLCL*AU0
1076 VMFLCL=UMF(K)
1077 UPOLD=VMFLCL
1078 UPNEW=UPOLD
1079 !
1080 !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
1081 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
1082 !...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
1083 !...PRODUCTION...
1084 !
1085 RATIO2(K)=0.
1086 UER(K)=0.
1087 ABE=0.
1088 TRPPT=0.
1089 TU(K)=TLCL
1090 TVU(K)=TVLCL
1091 QU(K)=QMIX
1092 EQFRC(K)=1.
1093 QLIQ(K)=0.
1094 QICE(K)=0.
1095 QLQOUT(K)=0.
1096 QICOUT(K)=0.
1097 DETLQ(K)=0.
1098 DETIC(K)=0.
1099 PPTLIQ(K)=0.
1100 PPTICE(K)=0.
1101 IFLAG=0
1102 !
1103 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
1104 !...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
1105 !...FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
1106 !...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
1107 !...PREVIOUS MODEL LEVEL...
1108 !
1109 TTEMP=TTFRZ
1110 !
1111 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
1112 !...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
1113 !...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
1114 !
1115 ! **1 variables indicate the bottom of a model layer
1116 ! **2 variables indicate the top of a model layer
1117 !
1118 EE1=1.
1119 UD1=0.
1120 REI = 0.
1121 DILBE = 0.
1122 updraft: DO NK=K,KL-1
1123 NK1=NK+1
1124 RATIO2(NK1)=RATIO2(NK)
1125 FRC1=0.
1126 TU(NK1)=T0(NK1)
1127 THETEU(NK1)=THETEU(NK)
1128 QU(NK1)=QU(NK)
1129 QLIQ(NK1)=QLIQ(NK)
1130 QICE(NK1)=QICE(NK)
1131 call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
1132 qice(nk1),qnewlq,qnewic,XLV1,XLV0)
1133 !
1134 !
1135 !...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
1136 !...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
1137 !...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
1138 !...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
1139 !...LIQUID WATER IS FROZEN AT EACH LEVEL...
1140 !
1141 IF(TU(NK1).LE.TTFRZ)THEN
1142 IF(TU(NK1).GT.TBFRZ)THEN
1143 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
1144 FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
1145 ELSE
1146 FRC1=1.
1147 IFLAG=1
1148 ENDIF
1149 TTEMP=TU(NK1)
1150 !
1151 ! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
1152 !...IS BELOW TTFRZ...
1153 !
1154 QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
1155 QNEWIC=QNEWIC+QNEWLQ*FRC1
1156 QNEWLQ=QNEWLQ-QNEWLQ*FRC1
1157 QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
1158 QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
1159 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
1160 QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1161 ENDIF
1162 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
1163 !
1164 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
1165 !
1166 IF(NK.EQ.K)THEN
1167 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
1168 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
1169 DZZ=Z0(NK1)-ZLCL
1170 ELSE
1171 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
1172 BOTERM=2.*DZA(NK)*G*BE/1.5
1173 DZZ=DZA(NK)
1174 ENDIF
1175 ENTERM=2.*REI*WTW/UPOLD
1177 CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
1178 RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
1179 !
1180 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
1181 !...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
1182 !
1183 IF(WTW.LT.1.E-3)THEN
1184 EXIT
1185 ELSE
1186 WU(NK1)=SQRT(WTW)
1187 ENDIF
1188 !...Calculate value of THETA-E in environment to entrain into updraft...
1189 !
1190 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1191 !
1192 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
1193 !
1194 REI=VMFLCL*DP(NK1)*0.03/RAD ! KF (1990) Eq. 1; Kain (2004) Eq. 5
1195 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
1196 IF(NK.EQ.K)THEN
1197 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
1198 ELSE
1199 DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
1200 ENDIF
1201 IF(DILBE.GT.0.)ABE=ABE+DILBE*G
1202 !
1203 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
1204 !...ENTRAINMENT (0.5*REI) IS IMPOSED...
1205 !
1206 IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
1207 EE2=0.5 ! Kain (2004) Eq. 4
1208 UD2=1.
1209 EQFRC(NK1)=0.
1210 ELSE
1211 LET=NK1
1212 TTMP=TVQU(NK1)
1213 !
1214 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
1215 !
1216 F1=0.95
1217 F2=1.-F1
1218 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
1219 QTMP=F1*Q0(NK1)+F2*QU(NK1)
1220 TMPLIQ=F2*QLIQ(NK1)
1221 TMPICE=F2*QICE(NK1)
1222 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
1223 qnewlq,qnewic,XLV1,XLV0)
1224 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
1225 IF(TU95.GT.TV0(NK1))THEN
1226 EE2=1.
1227 UD2=0.
1228 EQFRC(NK1)=1.0
1229 ELSE
1230 F1=0.10
1231 F2=1.-F1
1232 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
1233 QTMP=F1*Q0(NK1)+F2*QU(NK1)
1234 TMPLIQ=F2*QLIQ(NK1)
1235 TMPICE=F2*QICE(NK1)
1236 call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
1237 qnewlq,qnewic,XLV1,XLV0)
1238 TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
1239 TVDIFF = ABS(TU10-TVQU(NK1))
1240 IF(TVDIFF.LT.1.e-3)THEN
1241 EE2=1.
1242 UD2=0.
1243 EQFRC(NK1)=1.0
1244 ELSE
1245 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
1246 EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
1247 EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
1248 IF(EQFRC(NK1).EQ.1)THEN
1249 EE2=1.
1250 UD2=0.
1251 ELSEIF(EQFRC(NK1).EQ.0.)THEN
1252 EE2=0.
1253 UD2=1.
1254 ELSE
1255 !
1256 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
1257 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
1258 !
1259 CALL PROF5(EQFRC(NK1),EE2,UD2)
1260 ENDIF
1261 ENDIF
1262 ENDIF
1263 ENDIF ! End of Entrain/Detrain IF BLOCK
1264 !
1265 !
1266 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
1267 ! VALUES IN THE LAYER...
1268 !
1269 EE2 = AMAX1(EE2,0.5)
1270 UD2 = 1.5*UD2
1271 UER(NK1)=0.5*REI*(EE1+EE2)
1272 UDR(NK1)=0.5*REI*(UD1+UD2)
1273 !
1274 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
1275 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
1276 !
1277 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
1278 !
1279 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
1280 ! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
1281 ! First, correct ABE calculation if needed...
1282 !
1283 IF(DILBE.GT.0.)THEN
1284 ABE=ABE-DILBE*G
1285 ENDIF
1286 LET=NK
1287 ! WRITE(98,1015)P0(NK1)/100.
1288 EXIT
1289 ELSE
1290 EE1=EE2
1291 UD1=UD2
1292 UPOLD=UMF(NK)-UDR(NK1)
1293 UPNEW=UPOLD+UER(NK1)
1294 UMF(NK1)=UPNEW
1295 DILFRC(NK1) = UPNEW/UPOLD
1296 !
1297 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
1298 !...ICE IN THE DETRAINING UPDRAFT MASS...
1299 !
1300 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
1301 DETIC(NK1)=QICE(NK1)*UDR(NK1)
1302 QDT(NK1)=QU(NK1)
1303 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
1304 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
1305 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
1306 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
1307 !
1308 !...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
1309 !...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
1310 !...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
1311 !...CURRENT MODEL LEVEL...
1312 !
1313 PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
1314 PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
1315 !
1316 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
1317 IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
1318 ENDIF
1319 !
1320 END DO updraft
1321 !
1322 !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIUM
1323 ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
1324 ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BETWEEN
1325 ! THE LET AND CLOUD TOP...
1326 !
1327 !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOCITY
1328 ! FIRST BECOMES NEGATIVE...
1329 !
1330 LTOP=NK
1331 CLDHGT(LC)=Z0(LTOP)-ZLCL
1332 !
1333 !...Instead of using the same minimum cloud height (for deep convection)
1334 !...everywhere, try specifying minimum cloud depth as a function of TLCL...
1335 !
1336 ! Kain (2004) Eq. 7
1337 !
1338 IF(TLCL.GT.293.)THEN
1339 CHMIN = 4.E3
1340 ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
1341 CHMIN = 2.E3 + 100.*(TLCL-273.)
1342 ELSEIF(TLCL.LT.273.)THEN
1343 CHMIN = 2.E3
1344 ENDIF
1345 !ckay
1346 DO NK=K,LTOP
1347 qc_KF(I,NK,J)=QLIQ(NK)
1348 qi_KF(I,NK,J)=QICE(NK)
1349 END DO
1351 !
1352 !...If cloud top height is less than the specified minimum for deep
1353 !...convection, save value to consider this level as source for
1354 !...shallow convection, go back up to check next level...
1355 !
1356 !...Try specifying minimum cloud depth as a function of TLCL...
1357 !
1358 !
1359 !...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
1360 !
1361 !... 1.) if there is no CAPE, or
1362 !... 2.) cloud top is at model level just above LCL, or
1363 !... 3.) cloud top is within updraft source layer, or
1364 !... 4.) cloud-top detrainment layer begins within
1365 !... updraft source layer.
1366 !
1367 IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN ! No Convection Allowed
1368 CLDHGT(LC)=0.
1369 DO NK=K,LTOP
1370 UMF(NK)=0.
1371 UDR(NK)=0.
1372 UER(NK)=0.
1373 DETLQ(NK)=0.
1374 DETIC(NK)=0.
1375 PPTLIQ(NK)=0.
1376 PPTICE(NK)=0.
1377 !ckay
1378 cldfra_dp_KF(I,NK,J)=0.
1379 cldfra_sh_KF(I,NK,J)=0.
1380 qc_KF(I,NK,J)=0.
1381 qi_KF(I,NK,J)=0.
1382 ENDDO
1383 !
1384 ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
1385 ISHALL=0
1386 !ckay
1387 DO NK=K,LTOP
1388 cldfra_sh_KF(I,NK,J)=0.
1389 ENDDO
1390 EXIT usl
1391 ELSE
1392 !
1393 !...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
1394 ISHALL = 1
1395 !ckay
1396 DO NK=K,LTOP
1397 cldfra_dp_KF(I,NK,J)=0.
1398 ENDDO
1399 IF(NU.EQ.NUCHM)THEN
1400 EXIT usl ! Shallow Convection from this layer
1401 ELSE
1402 ! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
1403 DO NK=K,LTOP
1404 UMF(NK)=0.
1405 UDR(NK)=0.
1406 UER(NK)=0.
1407 DETLQ(NK)=0.
1408 DETIC(NK)=0.
1409 PPTLIQ(NK)=0.
1410 PPTICE(NK)=0.
1411 !ckay
1412 cldfra_dp_KF(I,NK,J)=0.
1413 cldfra_sh_KF(I,NK,J)=0.
1414 qc_KF(I,NK,J)=0.
1415 qi_KF(I,NK,J)=0.
1416 ENDDO
1417 ENDIF
1418 ENDIF
1419 ENDIF trigger2
1420 END DO usl
1421 IF(ISHALL.EQ.1)THEN
1422 KSTART=MAX0(KPBL,KLCL)
1423 LET=KSTART
1424 endif
1425 !
1426 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
1427 ! THIS LEVEL...
1428 !
1429 IF(LET.EQ.LTOP)THEN
1430 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
1431 DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1432 DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
1433 UER(LTOP)=0.
1434 UMF(LTOP)=0.
1435 ELSE
1436 !
1437 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1438 !
1439 DPTT=0.
1440 DO NJ=LET+1,LTOP
1441 DPTT=DPTT+DP(NJ)
1442 ENDDO
1443 DUMFDP=UMF(LET)/DPTT
1444 !
1445 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
1446 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
1447 !
1448 DO NK=LET+1,LTOP
1449 !
1450 !...entrainment is allowed at every level except for LTOP, so disallow
1451 !...entrainment at LTOP and adjust entrainment rates between LET and LTOP
1452 !...so the the dilution factor due to entrainment is not changed but
1453 !...the actual entrainment rate will change due due forced total
1454 !...detrainment in this layer...
1455 !
1456 IF(NK.EQ.LTOP)THEN
1457 UDR(NK) = UMF(NK-1)
1458 UER(NK) = 0.
1459 DETLQ(NK) = UDR(NK)*QLIQ(NK)*DILFRC(NK)
1460 DETIC(NK) = UDR(NK)*QICE(NK)*DILFRC(NK)
1461 ELSE
1462 UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
1463 UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
1464 UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
1465 DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
1466 DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
1467 ENDIF
1468 IF(NK.GE.LET+2)THEN
1469 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1470 PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
1471 PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
1472 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1473 ENDIF
1474 ENDDO
1475 ENDIF
1476 !
1477 ! Initialize some arrays below cloud base and above cloud top...
1478 !
1479 DO NK=1,LTOP
1480 IF(T0(NK).GT.T00)ML=NK
1481 ENDDO
1482 DO NK=1,K
1483 IF(NK.GE.LC)THEN
1484 IF(NK.EQ.LC)THEN
1485 UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1486 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1487 ELSEIF(NK.LE.KPBL)THEN
1488 UER(NK)=VMFLCL*DP(NK)/DPTHMX
1489 UMF(NK)=UMF(NK-1)+UER(NK)
1490 ELSE
1491 UMF(NK)=VMFLCL
1492 UER(NK)=0.
1493 ENDIF
1494 TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1495 QU(NK)=QMIX
1496 WU(NK)=WLCL
1497 ELSE
1498 TU(NK)=0.
1499 QU(NK)=0.
1500 UMF(NK)=0.
1501 WU(NK)=0.
1502 UER(NK)=0.
1503 !ckay
1504 cldfra_dp_KF(I,NK,J)=0.
1505 cldfra_sh_KF(I,NK,J)=0.
1506 qc_KF(I,NK,J)=0.
1507 qi_KF(I,NK,J)=0.
1508 ENDIF
1509 UDR(NK)=0.
1510 QDT(NK)=0.
1511 QLIQ(NK)=0.
1512 QICE(NK)=0.
1513 QLQOUT(NK)=0.
1514 QICOUT(NK)=0.
1515 PPTLIQ(NK)=0.
1516 PPTICE(NK)=0.
1517 DETLQ(NK)=0.
1518 DETIC(NK)=0.
1519 RATIO2(NK)=0.
1520 CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
1521 EQFRC(NK)=1.0
1522 ENDDO
1523 !
1524 LTOP1=LTOP+1
1525 LTOPM1=LTOP-1
1526 !
1527 !...DEFINE VARIABLES ABOVE CLOUD TOP...
1528 !
1529 DO NK=LTOP1,KX
1530 UMF(NK)=0.
1531 UDR(NK)=0.
1532 UER(NK)=0.
1533 QDT(NK)=0.
1534 QLIQ(NK)=0.
1535 QICE(NK)=0.
1536 QLQOUT(NK)=0.
1537 QICOUT(NK)=0.
1538 DETLQ(NK)=0.
1539 DETIC(NK)=0.
1540 PPTLIQ(NK)=0.
1541 PPTICE(NK)=0.
1542 IF(NK.GT.LTOP1)THEN
1543 TU(NK)=0.
1544 QU(NK)=0.
1545 WU(NK)=0.
1546 !ckay
1547 cldfra_dp_KF(I,NK,J)=0.
1548 cldfra_sh_KF(I,NK,J)=0.
1549 qc_KF(I,NK,J)=0.
1550 qi_KF(I,NK,J)=0.
1551 ENDIF
1552 THTA0(NK)=0.
1553 THTAU(NK)=0.
1554 EMS(NK)=0.
1555 EMSD(NK)=0.
1556 TG(NK)=T0(NK)
1557 QG(NK)=Q0(NK)
1558 QLG(NK)=0.
1559 QIG(NK)=0.
1560 QRG(NK)=0.
1561 QSG(NK)=0.
1562 OMG(NK)=0.
1563 ENDDO
1564 OMG(KX+1)=0.
1565 DO NK=1,LTOP
1566 EMS(NK)=DP(NK)*DXSQ/G
1567 EMSD(NK)=1./EMS(NK)
1568 !
1569 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCHEME
1570 !
1571 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1572 THTAU(NK)=TU(NK)*EXN(NK)
1573 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1574 THTA0(NK)=T0(NK)*EXN(NK)
1575 DDILFRC(NK) = 1./DILFRC(NK)
1576 OMG(NK)=0.
1577 ENDDO
1578 ! IF (XTIME.LT.10.)THEN
1579 ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1580 ! * TMIX-T00,PMIX,QMIX,ABE
1581 ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1582 ! * WLCL,CLDHGT
1583 ! ENDIF
1584 !
1585 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1586 !...AND MIDTROPOSPHERE IS USED.
1587 !
1588 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1589 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1590 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1591 VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1592 !...for ETA model, DX is a function of location...
1593 ! TIMEC=DX(I,J)/VCONV
1594 TIMEC=DX/VCONV
1595 TADVEC=TIMEC
1596 TIMEC=AMAX1(1800.,TIMEC) ! 30 minutes >= TIMEC <= 60 minutes
1597 TIMEC=AMIN1(3600.,TIMEC)
1598 IF(ISHALL.EQ.1)TIMEC=2400. ! shallow convection TIMEC = 40 minutes
1599 NIC=NINT(TIMEC/DT)
1600 TIMEC=FLOAT(NIC)*DT
1601 !
1602 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1603 !
1604 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1605 SHSIGN=1.
1606 ELSE
1607 SHSIGN=-1.
1608 ENDIF
1609 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
1610 (V0(LTOP)-V0(KLCL))
1611 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1612 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1613 PEF=AMAX1(PEF,.2)
1614 PEF=AMIN1(PEF,.9)
1615 !
1616 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1617 !
1618 CBH=(ZLCL-Z0(1))*3.281E-3
1619 IF(CBH.LT.3.)THEN
1620 RCBH=.02
1621 ELSE
1622 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
1623 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1624 ENDIF
1625 IF(CBH.GT.25)RCBH=2.4
1626 PEFCBH=1./(1.+RCBH)
1627 PEFCBH=AMIN1(PEFCBH,.9)
1628 !
1629 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1630 !
1631 PEFF=.5*(PEF+PEFCBH)
1632 PEFF2 = PEFF ! JSK MODS
1633 IF(IPRNT)THEN
1634 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1635 WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1636 CALL wrf_message( message )
1637 ! flush(98)
1638 endif
1639 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1640 !*****************************************************************
1641 ! *
1642 ! COMPUTE DOWNDRAFT PROPERTIES *
1643 ! *
1644 !*****************************************************************
1645 !
1646 !
1647 TDER=0.
1648 devap:IF(ISHALL.EQ.1)THEN
1649 LFS = 1
1650 ELSE
1651 !
1652 !...start downdraft about 150 mb above cloud base...
1653 !
1654 ! KSTART=MAX0(KPBL,KLCL)
1655 ! KSTART=KPBL ! Changed 7/23/99
1656 KSTART=KPBL+1 ! Changed 7/23/99
1657 KLFS = LET-1
1658 DO NK = KSTART+1,KL
1659 DPPP = P0(KSTART)-P0(NK)
1660 ! IF(DPPP.GT.200.E2)THEN
1661 IF(DPPP.GT.150.E2)THEN
1662 KLFS = NK
1663 EXIT
1664 ENDIF
1665 ENDDO
1666 KLFS = MIN0(KLFS,LET-1)
1667 LFS = KLFS
1668 !
1669 !...if LFS is not at least 50 mb above cloud base (implying that the
1670 !...level of equil temp, LET, is just above cloud base) do not allow a
1671 !...downdraft...
1672 !
1673 IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
1674 THETED(LFS) = THETEE(LFS)
1675 QD(LFS) = Q0(LFS)
1676 !
1677 !...call tpmix2dd to find wet-bulb temp, qv...
1678 !
1679 call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
1680 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
1681 !
1682 !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
1683 !
1684 TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
1685 RDD=P0(LFS)/(R*TVD(LFS))
1686 A1=(1.-PEFF)*AU0
1687 DMF(LFS)=-A1*RDD
1688 DER(LFS)=DMF(LFS)
1689 DDR(LFS)=0.
1690 RHBAR = RH(LFS)*DP(LFS)
1691 DPTT = DP(LFS)
1692 DO ND = LFS-1,KSTART,-1
1693 ND1 = ND+1
1694 DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
1695 DDR(ND)=0.
1696 DMF(ND)=DMF(ND1)+DER(ND)
1697 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1698 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
1699 DPTT = DPTT+DP(ND)
1700 RHBAR = RHBAR+RH(ND)*DP(ND)
1701 ENDDO
1702 RHBAR = RHBAR/DPTT
1703 DMFFRC = 2.*(1.-RHBAR) ! Kain (2004) eq. 11
1704 DPDD = 0.
1705 !...Calculate melting effect
1706 !... first, compute total frozen precipitation generated...
1707 !
1708 pptmlt = 0.
1709 DO NK = KLCL,LTOP
1710 PPTMLT = PPTMLT+PPTICE(NK)
1711 ENDDO
1712 if(lc.lt.ml)then
1713 !...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
1714 !...if DMFFRC=1. Otherwise, for small DMFFRC, DTMELT gets too large!
1715 !...12/14/98 jsk...
1716 DTMELT = RLF*PPTMLT/(CP*UMF(KLCL))
1717 else
1718 DTMELT = 0.
1719 endif
1720 LDT = MIN0(LFS-1,KSTART-1)
1721 !
1722 call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
1723 !
1724 tz(kstart) = tz(kstart)-dtmelt
1725 ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
1726 QSS=0.622*ES/(P0(KSTART)-ES)
1727 THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**(0.2854*(1.-0.28*QSS))* &
1728 EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
1729 !....
1730 LDT = MIN0(LFS-1,KSTART-1)
1731 DO ND = LDT,1,-1
1732 DPDD = DPDD+DP(ND)
1733 THETED(ND) = THETED(KSTART)
1734 QD(ND) = QD(KSTART)
1735 !
1736 !...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
1737 !
1738 call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
1739 qsd(nd) = qss
1740 !
1741 !...specify RH decrease of 20%/km in downdraft...
1742 !
1743 RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
1744 !
1745 !...adjust downdraft TEMP, Q to specified RH:
1746 !
1747 IF(RHH.LT.1.)THEN
1748 DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
1749 RL=XLV0-XLV1*TZ(ND)
1750 DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
1751 T1RH=TZ(ND)+DTMP
1752 ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1753 QSRH=0.622*ES/(P0(ND)-ES)
1754 !
1755 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1756 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1757 !
1758 IF(QSRH.LT.QD(ND))THEN
1759 QSRH=QD(ND)
1760 T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
1761 ENDIF
1762 TZ(ND)=T1RH
1763 QSS=QSRH
1764 QSD(ND) = QSS
1765 ENDIF
1766 TVD(nd) = tz(nd)*(1.+0.608*qsd(nd))
1767 IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
1768 LDB=ND
1769 EXIT
1770 ENDIF
1771 ENDDO
1772 IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN ! minimum Downdraft depth!
1773 DO ND=LDT,LDB,-1
1774 ND1 = ND+1
1775 DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
1776 DER(ND) = 0.
1777 DMF(ND) = DMF(ND1)+DDR(ND)
1778 TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
1779 QD(ND)=QSD(nd)
1780 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1781 ENDDO
1782 ENDIF
1783 ENDIF
1784 ENDIF devap
1785 !
1786 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1787 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1788 !
1789 d_mf: IF(TDER.LT.1.)THEN
1790 ! WRITE(98,3004)I,J
1791 !3004 FORMAT(' ','No Downdraft!; I=',I3,2X,'J=',I3,'ISHALL =',I2)
1792 PPTFLX=TRPPT
1793 CPR=TRPPT
1794 TDER=0.
1795 CNDTNF=0.
1796 UPDINC=1.
1797 LDB=LFS
1798 DO NDK=1,LTOP
1799 DMF(NDK)=0.
1800 DER(NDK)=0.
1801 DDR(NDK)=0.
1802 THTAD(NDK)=0.
1803 WD(NDK)=0.
1804 TZ(NDK)=0.
1805 QD(NDK)=0.
1806 ENDDO
1807 AINCM2=100.
1808 ELSE
1809 DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
1810 UPDINC=1.
1811 IF(TDER*DDINC.GT.TRPPT)THEN
1812 DDINC = TRPPT/TDER
1813 ENDIF
1814 TDER = TDER*DDINC
1815 DO NK=LDB,LFS
1816 DMF(NK)=DMF(NK)*DDINC
1817 DER(NK)=DER(NK)*DDINC
1818 DDR(NK)=DDR(NK)*DDINC
1819 ENDDO
1820 CPR=TRPPT
1821 PPTFLX = TRPPT-TDER
1822 PEFF=PPTFLX/TRPPT
1823 IF(IPRNT)THEN
1824 ! write(98,*)'PRECIP EFFICIENCY =',PEFF
1825 write(message,*)'PRECIP EFFICIENCY =',PEFF
1826 CALL wrf_message(message)
1827 ! flush(98)
1828 ENDIF
1829 !
1830 !
1831 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1832 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1833 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1834 !
1835 ! DO NK=LC,LFS
1836 ! UMF(NK)=UMF(NK)*UPDINC
1837 ! UDR(NK)=UDR(NK)*UPDINC
1838 ! UER(NK)=UER(NK)*UPDINC
1839 ! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1840 ! PPTICE(NK)=PPTICE(NK)*UPDINC
1841 ! DETLQ(NK)=DETLQ(NK)*UPDINC
1842 ! DETIC(NK)=DETIC(NK)*UPDINC
1843 ! ENDDO
1844 !
1845 !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1846 !...DOWNDRAFT...
1847 !
1848 IF(LDB.GT.1)THEN
1849 DO NK=1,LDB-1
1850 DMF(NK)=0.
1851 DER(NK)=0.
1852 DDR(NK)=0.
1853 WD(NK)=0.
1854 TZ(NK)=0.
1855 QD(NK)=0.
1856 THTAD(NK)=0.
1857 ENDDO
1858 ENDIF
1859 DO NK=LFS+1,KX
1860 DMF(NK)=0.
1861 DER(NK)=0.
1862 DDR(NK)=0.
1863 WD(NK)=0.
1864 TZ(NK)=0.
1865 QD(NK)=0.
1866 THTAD(NK)=0.
1867 ENDDO
1868 DO NK=LDT+1,LFS-1
1869 TZ(NK)=0.
1870 QD(NK)=0.
1871 THTAD(NK)=0.
1872 ENDDO
1873 ENDIF d_mf
1874 !
1875 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFLOW
1876 ! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILABLE
1877 ! IN THAT LAYER INITIALLY...
1878 !
1879 AINCMX=1000.
1880 LMAX=MAX0(KLCL,LFS)
1881 DO NK=LC,LMAX
1882 IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
1883 AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
1884 AINCMX=AMIN1(AINCMX,AINCM1)
1885 ENDIF
1886 ENDDO
1887 AINC=1.
1888 IF(AINCMX.LT.AINC)AINC=AINCMX
1889 !
1890 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL
1891 !...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
1892 !...CLOSURE...
1893 !
1894 TDER2=TDER
1895 PPTFL2=PPTFLX
1896 DO NK=1,LTOP
1897 DETLQ2(NK)=DETLQ(NK)
1898 DETIC2(NK)=DETIC(NK)
1899 UDR2(NK)=UDR(NK)
1900 UER2(NK)=UER(NK)
1901 DDR2(NK)=DDR(NK)
1902 DER2(NK)=DER(NK)
1903 UMF2(NK)=UMF(NK)
1904 DMF2(NK)=DMF(NK)
1905 ENDDO
1906 FABE=1.
1907 STAB=0.95
1908 NOITR=0
1909 ISTOP=0
1910 !
1911 IF(ISHALL.EQ.1)THEN ! First for shallow convection
1912 !
1913 ! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
1914 ! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
1915 ! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
1916 !
1917 !...find the maximum TKE value between LC and KLCL...
1918 ! TKEMAX = 0.
1919 TKEMAX = 5.
1920 ! DO 173 K = LC,KLCL
1921 ! NK = KX-K+1
1922 ! TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
1923 ! 173 CONTINUE
1924 ! TKEMAX = AMIN1(TKEMAX,10.)
1925 ! TKEMAX = AMAX1(TKEMAX,5.)
1926 !c TKEMAX = 10.
1927 !c...3_24_99...DPMIN was changed for shallow convection so that it is the
1928 !c... the same as for deep convection (5.E3). Since this doubles
1929 !c... (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
1930 !c... lation of EVAC...
1931 !c EVAC = TKEMAX*0.1
1932 EVAC = 0.5*TKEMAX*0.1
1933 ! AINC = 0.1*DPTHMX*DXIJ*DXIJ/(VMFLCL*G*TIMEC)
1934 ! AINC = EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
1935 AINC = EVAC*DPTHMX*DXSQ/(VMFLCL*G*TIMEC)
1936 TDER=TDER2*AINC
1937 PPTFLX=PPTFL2*AINC
1938 DO NK=1,LTOP
1939 UMF(NK)=UMF2(NK)*AINC
1940 DMF(NK)=DMF2(NK)*AINC
1941 DETLQ(NK)=DETLQ2(NK)*AINC
1942 DETIC(NK)=DETIC2(NK)*AINC
1943 UDR(NK)=UDR2(NK)*AINC
1944 UER(NK)=UER2(NK)*AINC
1945 DER(NK)=DER2(NK)*AINC
1946 DDR(NK)=DDR2(NK)*AINC
1947 ENDDO
1948 ENDIF ! Otherwise for deep convection
1949 ! use iterative procedure to find mass fluxes...
1950 iter: DO NCOUNT=1,10
1951 !
1952 !*****************************************************************
1953 ! *
1954 ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1955 ! *
1956 !*****************************************************************
1957 !
1958 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1959 !...SATISFY MASS CONTINUITY...
1960 !
1961 DTT=TIMEC
1962 DO NK=1,LTOP
1963 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1964 IF(NK.GT.1)THEN
1965 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1966 ABSOMG = ABS(OMG(NK))
1967 ABSOMGTC = ABSOMG*TIMEC
1968 FRDP = 0.75*DP(NK-1)
1969 IF(ABSOMGTC.GT.FRDP)THEN
1970 DTT1 = FRDP/ABSOMG
1971 DTT=AMIN1(DTT,DTT1)
1972 ENDIF
1973 ENDIF
1974 ENDDO
1975 DO NK=1,LTOP
1976 THPA(NK)=THTA0(NK)
1977 QPA(NK)=Q0(NK)
1978 NSTEP=NINT(TIMEC/DTT+1)
1979 DTIME=TIMEC/FLOAT(NSTEP)
1980 FXM(NK)=OMG(NK)*DXSQ/G
1981 ENDDO
1982 !
1983 !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1984 !
1985 DO NTC=1,NSTEP
1986 !
1987 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED ON
1988 !...SIGN OF OMEGA...
1989 !
1990 DO NK=1,LTOP
1991 THFXIN(NK)=0.
1992 THFXOUT(NK)=0.
1993 QFXIN(NK)=0.
1994 QFXOUT(NK)=0.
1995 ENDDO
1996 DO NK=2,LTOP
1997 IF(OMG(NK).LE.0.)THEN
1998 THFXIN(NK)=-FXM(NK)*THPA(NK-1)
1999 QFXIN(NK)=-FXM(NK)*QPA(NK-1)
2000 THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
2001 QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
2002 ELSE
2003 THFXOUT(NK)=FXM(NK)*THPA(NK)
2004 QFXOUT(NK)=FXM(NK)*QPA(NK)
2005 THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
2006 QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
2007 ENDIF
2008 ENDDO
2009 !
2010 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
2011 !
2012 DO NK=1,LTOP
2013 THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
2014 THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
2015 DTIME*EMSD(NK)
2016 QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)- &
2017 QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
2018 ENDDO
2019 ENDDO
2020 DO NK=1,LTOP
2021 THTAG(NK)=THPA(NK)
2022 QG(NK)=QPA(NK)
2023 ENDDO
2024 !
2025 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORROW
2026 !...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
2027 !
2028 DO NK=1,LTOP
2029 IF(QG(NK).LT.0.)THEN
2030 IF(NK.EQ.1)THEN ! JSK MODS
2031 ! PRINT *,' PROBLEM WITH KF SCHEME: ' ! JSK MODS
2032 ! PRINT *,'QG = 0 AT THE SURFACE!!!!!!!' ! JSK MODS
2033 CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
2034 ENDIF ! JSK MODS
2035 NK1=NK+1
2036 IF(NK.EQ.LTOP)THEN
2037 NK1=KLCL
2038 ENDIF
2039 TMA=QG(NK1)*EMS(NK1)
2040 TMB=QG(NK-1)*EMS(NK-1)
2041 TMM=(QG(NK)-1.E-9)*EMS(NK )
2042 BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
2043 ACOEFF=BCOEFF*TMA/TMB
2044 TMB=TMB*(1.-BCOEFF)
2045 TMA=TMA*(1.-ACOEFF)
2046 IF(NK.EQ.LTOP)THEN
2047 QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
2048 ! IF(ABS(QVDIFF).GT.1.)THEN
2049 ! PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ', &
2050 ! QVDIFF, &
2051 ! '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ', &
2052 ! 'VALUES IN KAIN-FRITSCH'
2053 ! ENDIF
2054 ENDIF
2055 QG(NK)=1.E-9
2056 QG(NK1)=TMA*EMSD(NK1)
2057 QG(NK-1)=TMB*EMSD(NK-1)
2058 ENDIF
2059 ENDDO
2060 TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
2061 IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
2062 ! WRITE(99,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME; &
2063 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
2064 ! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
2065 ISTOP=1
2066 IPRNT=.TRUE.
2067 EXIT iter
2068 ENDIF
2069 !
2070 !...CONVERT THETA TO T...
2071 !
2072 DO NK=1,LTOP
2073 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
2074 TG(NK)=THTAG(NK)/EXN(NK)
2075 TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
2076 ENDDO
2077 IF(ISHALL.EQ.1)THEN
2078 EXIT iter
2079 ENDIF
2080 !
2081 !*******************************************************************
2082 ! *
2083 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
2084 ! *
2085 !*******************************************************************
2086 !
2087 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
2088 !
2089 ! THMIX=0.
2090 TMIX=0.
2091 QMIX=0.
2092 !
2093 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
2094 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
2095 !...LAYERS...
2096 !
2097 DO NK=LC,KPBL
2098 TMIX=TMIX+DP(NK)*TG(NK)
2099 QMIX=QMIX+DP(NK)*QG(NK)
2100 ENDDO
2101 TMIX=TMIX/DPTHMX
2102 QMIX=QMIX/DPTHMX
2103 ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
2104 QSS=0.622*ES/(PMIX-ES)
2105 !
2106 !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
2107 !
2108 IF(QMIX.GT.QSS)THEN
2109 RL=XLV0-XLV1*TMIX
2110 CPM=CP*(1.+0.887*QMIX)
2111 DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
2112 DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
2113 TMIX=TMIX+RL/CP*DQ
2114 QMIX=QMIX-DQ
2115 TLCL=TMIX
2116 ELSE
2117 QMIX=AMAX1(QMIX,0.)
2118 EMIX=QMIX*PMIX/(0.622+QMIX)
2119 astrt=1.e-3
2120 binc=0.075
2121 a1=emix/aliq
2122 tp=(a1-astrt)/binc
2123 indlu=int(tp)+1
2124 value=(indlu-1)*binc+astrt
2125 aintrp=(a1-value)/binc
2126 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
2127 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2128 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-TDPT)
2129 TLCL=AMIN1(TLCL,TMIX)
2130 ENDIF
2131 TVLCL=TLCL*(1.+0.608*QMIX)
2132 ZLCL = ZMIX+(TLCL-TMIX)/GDRY
2133 DO NK = LC,KL
2134 KLCL=NK
2135 IF(ZLCL.LE.Z0(NK))THEN
2136 EXIT
2137 ENDIF
2138 ENDDO
2139 K=KLCL-1
2140 DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
2141 !
2142 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
2143 !
2144 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
2145 QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
2146 TVEN=TENV*(1.+0.608*QENV)
2147 PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
2148 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
2149 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
2150 !
2151 !...COMPUTE ADJUSTED ABE(ABEG).
2152 !
2153 ABEG=0.
2154 DO NK=K,LTOPM1
2155 NK1=NK+1
2156 THETEU(NK1) = THETEU(NK)
2157 !
2158 call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
2159 !
2160 TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
2161 IF(NK.EQ.K)THEN
2162 DZZ=Z0(KLCL)-ZLCL
2163 DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
2164 ELSE
2165 DZZ=DZA(NK)
2166 DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
2167 ENDIF
2168 IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
2169 !
2170 !...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
2171 !
2172 CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
2173 THETEU(NK1)=THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
2174 ENDDO
2175 !
2176 !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
2177 !...THE PERIOD TIMEC...
2178 !
2179 IF(NOITR.EQ.1)THEN
2180 ! write(98,*)' '
2181 ! write(98,*)'TAU, I, J, =',NTSD,I,J
2182 ! WRITE(98,1060)FABE
2183 ! GOTO 265
2184 EXIT iter
2185 ENDIF
2186 DABE=AMAX1(ABE-ABEG,0.1*ABE)
2187 FABE=ABEG/ABE
2188 IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
2189 ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
2190 ! *GRID POINT; NO CONVECTION ALLOWED!'
2191 !BSINGH - For WRFCuP scheme
2192 !BSINGH - Improvised based on the following reason
2193 !Before every "RETURN" statement assign ishall=2
2194 ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
2195 ishall = 2
2196 !BSINGH -ENDS
2197 RETURN
2198 ENDIF
2199 IF(NCOUNT.NE.1)THEN
2200 IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
2201 NOITR=1
2202 AINC=AINCOLD
2203 CYCLE iter
2204 ENDIF
2205 DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
2206 IF(DFDA.GT.0.)THEN
2207 NOITR=1
2208 AINC=AINCOLD
2209 CYCLE iter
2210 ENDIF
2211 ENDIF
2212 AINCOLD=AINC
2213 FABEOLD=FABE
2214 IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
2215 ! write(98,*)' '
2216 ! write(98,*)'TAU, I, J, =',NTSD,I,J
2217 ! WRITE(98,1055)FABE
2218 ! GOTO 265
2219 EXIT
2220 ENDIF
2221 IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
2222 EXIT iter
2223 ELSE
2224 IF(NCOUNT.GT.10)THEN
2225 ! write(98,*)' '
2226 ! write(98,*)'TAU, I, J, =',NTSD,I,J
2227 ! WRITE(98,1060)FABE
2228 ! GOTO 265
2229 EXIT
2230 ENDIF
2231 !
2232 !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTIVE
2233 !...MASS FLUX BY THE FACTOR AINC:
2234 !
2235 IF(FABE.EQ.0.)THEN
2236 AINC=AINC*0.5
2237 ELSE
2238 IF(DABE.LT.1.e-4)THEN
2239 NOITR=1
2240 AINC=AINCOLD
2241 CYCLE iter
2242 ELSE
2243 AINC=AINC*STAB*ABE/DABE
2244 ENDIF
2245 ENDIF
2246 ! AINC=AMIN1(AINCMX,AINC)
2247 AINC=AMIN1(AINCMX,AINC)
2248 !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
2249 !...WILL BE MINIMAL SO JUST IGNORE IT... ! JSK MODS
2250 IF(AINC.LT.0.05)then
2251 !BSINGH - For WRFCuP scheme
2252 !BSINGH - Improvised based on the following reason
2253 !Before every "RETURN" statement assign ishall=2
2254 ! wig, 29-Aug-2006: Indicate no convection occurred in ishall.
2255 ishall = 2
2256 !BSINGH -ENDS
2258 RETURN ! JSK MODS
2259 ENDIF
2260 ! AINC=AMAX1(AINC,0.05) ! JSK MODS
2261 TDER=TDER2*AINC
2262 PPTFLX=PPTFL2*AINC
2263 ! IF (XTIME.LT.10.)THEN
2264 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
2265 ! * FABEOLD,AINCOLD
2266 ! ENDIF
2267 DO NK=1,LTOP
2268 UMF(NK)=UMF2(NK)*AINC
2269 DMF(NK)=DMF2(NK)*AINC
2270 DETLQ(NK)=DETLQ2(NK)*AINC
2271 DETIC(NK)=DETIC2(NK)*AINC
2272 UDR(NK)=UDR2(NK)*AINC
2273 UER(NK)=UER2(NK)*AINC
2274 DER(NK)=DER2(NK)*AINC
2275 DDR(NK)=DDR2(NK)*AINC
2276 ENDDO
2277 !
2278 !...GO BACK UP FOR ANOTHER ITERATION...
2279 !
2280 ENDIF
2281 ENDDO iter
2282 !ckay
2283 ! get the cloud fraction for layer NK+1=NK1
2284 IF(ISHALL.EQ.1) THEN
2285 DO NK=KLCL-1, LTOP1
2286 UMF_new = UMF(NK)/DXSQ
2287 xcldfra = 0.07*alog(1.+(500.*UMF_new))
2288 xcldfra = amax1(0.01,xcldfra)
2289 cldfra_sh_KF(I,NK,J) = amin1(0.2,xcldfra)
2290 ENDDO
2291 ELSE
2292 DO NK=KLCL-1, LTOP1
2293 UMF_new = UMF(NK)/DXSQ
2294 xcldfra = 0.14*alog(1.+(500.*UMF_new))
2295 xcldfra = amax1(0.01,xcldfra)
2296 cldfra_dp_KF(I,NK,J) = amin1(0.6,xcldfra)
2297 ENDDO
2298 ENDIF
2300 !kf_edrates
2301 !Save up/down entrainment/detrainment rates as 3D variables
2302 IF (KF_EDRATES == 1) THEN
2303 DO NK=1,LTOP
2304 UDR_KF(I,NK,J)=UDR(NK)
2305 DDR_KF(I,NK,J)=DDR(NK)
2306 UER_KF(I,NK,J)=UER(NK)
2307 DER_KF(I,NK,J)=DER(NK)
2308 ENDDO
2309 ENDIF
2310 !
2311 !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
2312 !
2313 !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE ! PPT FB MODS
2314 !...GENERATED THAT GOES INTO PRECIPITIATION ! PPT FB MODS
2315 !
2316 ! Redistribute hydormeteors according to the final mass-flux values:
2317 !
2318 IF(CPR.GT.0.)THEN
2319 FRC2=PPTFLX/(CPR*AINC) ! PPT FB MODS
2320 ELSE
2321 FRC2=0.
2322 ENDIF
2323 DO NK=1,LTOP
2324 QLPA(NK)=QL0(NK)
2325 QIPA(NK)=QI0(NK)
2326 QRPA(NK)=QR0(NK)
2327 QSPA(NK)=QS0(NK)
2328 RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
2329 SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2 ! PPT FB MODS
2330 ENDDO
2331 DO NTC=1,NSTEP
2332 !
2333 !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAYER
2334 !...BASED ON THE SIGN OF OMEGA...
2335 !
2336 DO NK=1,LTOP
2337 QLFXIN(NK)=0.
2338 QLFXOUT(NK)=0.
2339 QIFXIN(NK)=0.
2340 QIFXOUT(NK)=0.
2341 QRFXIN(NK)=0.
2342 QRFXOUT(NK)=0.
2343 QSFXIN(NK)=0.
2344 QSFXOUT(NK)=0.
2345 ENDDO
2346 DO NK=2,LTOP
2347 IF(OMG(NK).LE.0.)THEN
2348 QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
2349 QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
2350 QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
2351 QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
2352 QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
2353 QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
2354 QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
2355 QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
2356 ELSE
2357 QLFXOUT(NK)=FXM(NK)*QLPA(NK)
2358 QIFXOUT(NK)=FXM(NK)*QIPA(NK)
2359 QRFXOUT(NK)=FXM(NK)*QRPA(NK)
2360 QSFXOUT(NK)=FXM(NK)*QSPA(NK)
2361 QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
2362 QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
2363 QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
2364 QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
2365 ENDIF
2366 ENDDO
2367 !
2368 !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
2369 !
2370 DO NK=1,LTOP
2371 QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*EMSD(NK)
2372 QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*EMSD(NK)
2373 QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
2374 QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))*DTIME*EMSD(NK) ! PPT FB MODS
2375 ENDDO
2376 ENDDO
2377 DO NK=1,LTOP
2378 QLG(NK)=QLPA(NK)
2379 QIG(NK)=QIPA(NK)
2380 QRG(NK)=QRPA(NK)
2381 QSG(NK)=QSPA(NK)
2382 ENDDO
2384 !kf_edrates
2385 !Save convective timescale (TIMEC) as 2D variable
2386 IF (KF_EDRATES == 1) THEN
2387 TIMEC_KF(I,J)=TIMEC
2388 ENDIF
2389 !
2390 !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
2391 !...GRID POINT...
2392 !
2393 ! IF (XTIME.LT.10.)THEN
2394 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2395 ! ENDIF
2396 IF(IPRNT)THEN
2397 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2398 WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2399 CALL wrf_message(message)
2400 ! flush(98)
2401 endif
2402 !
2403 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
2404 !
2405 !297 IF(IPRNT)then
2406 IF(IPRNT)then
2407 ! if(I.eq.16 .and. J.eq.41)then
2408 ! IF(ISTOP.EQ.1)THEN
2409 write(98,*)
2410 ! write(98,*)'At t(h), I, J =',float(NTSD)*72./3600.,I,J
2411 write(message,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100., &
2412 TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
2413 call wrf_message(message)
2414 write(message,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL, &
2415 DTRH,TENV
2416 call wrf_message(message)
2417 WRITE(message,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG, &
2418 TMIX-T00,PMIX,QMIX,ABE
2419 call wrf_message(message)
2420 WRITE(message,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100., &
2421 WLCL,CLDHGT(LC)
2422 call wrf_message(message)
2423 WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
2424 call wrf_message(message)
2425 write(message,*)'PRECIP EFFICIENCY =',PEFF
2426 call wrf_message(message)
2427 WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2428 call wrf_message(message)
2429 ! ENDIF
2430 !!!!! HERE !!!!!!!
2431 WRITE(message,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
2432 ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
2433 ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
2434 call wrf_message(message)
2435 write(message,*)'just before DO 300...'
2436 call wrf_message(message)
2437 ! flush(98)
2438 DO NK=1,LTOP
2439 K=LTOP-NK+1
2440 DTT=(TG(K)-T0(K))*86400./TIMEC
2441 RL=XLV0-XLV1*TG(K)
2442 DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
2443 UDFRC=UDR(K)*TIMEC*EMSD(K)
2444 UEFRC=UER(K)*TIMEC*EMSD(K)
2445 DDFRC=DDR(K)*TIMEC*EMSD(K)
2446 DEFRC=-DER(K)*TIMEC*EMSD(K)
2447 WRITE(message,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4, &
2448 UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
2449 W0AVG1D(K)*1.E2,DETLQ(K)*TIMEC*EMSD(K)*1.E3,DETIC(K)* &
2450 TIMEC*EMSD(K)*1.E3
2451 call wrf_message(message)
2452 ENDDO
2453 WRITE(message,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
2454 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
2455 call wrf_message(message)
2456 DO NK=1,KL
2457 K=KX-NK+1
2458 DTT=TG(K)-T0(K)
2459 TUC=TU(K)-T00
2460 IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
2461 TDC=TZ(K)-T00
2462 IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
2463 IF(T0(K).LT.T00)THEN
2464 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2465 ELSE
2466 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
2467 ENDIF
2468 QGS=ES*0.622/(P0(K)-ES)
2469 RH0=Q0(K)/QES(K)
2470 RHG=QG(K)/QGS
2471 WRITE(message,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC, &
2472 TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))*1000.,QU(K)* &
2473 1000.,QD(K)*1000.,QLG(K)*1000.,QIG(K)*1000.,QRG(K)*1000., &
2474 QSG(K)*1000.,RH0,RHG
2475 call wrf_message(message)
2476 ENDDO
2477 !
2478 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
2479 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
2480 !
2481 ! IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
2483 ! IF(ISHALL.NE.1)THEN
2484 ! write(98,4421)i,j,iyr,imo,idy,ihr,imn
2485 ! write(98)i,j,iyr,imo,idy,ihr,imn,kl
2486 ! 4421 format(7i4)
2487 ! write(98,4422)kl
2488 ! 4422 format(i6)
2489 DO 310 NK = 1,KL
2490 k = kl - nk + 1
2491 write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2492 u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
2493 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2494 ! WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
2495 ! * U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
2496 310 CONTINUE
2497 IF(ISTOP.EQ.1)THEN
2498 CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
2499 ENDIF
2500 ! ENDIF
2501 4455 format(8f11.3)
2502 ENDIF
2503 CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
2504 PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2505 RAINCV(I,J)=DT*PRATEC(I,J) ! PPT FB MODS
2506 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
2507 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
2508 RNC=RAINCV(I,J)*NIC
2509 IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
2511 ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
2512 !
2513 ! EVALUATE MOISTURE BUDGET...
2514 !
2516 QINIT=0.
2517 QFNL=0.
2518 DPT=0.
2519 DO 315 NK=1,LTOP
2520 DPT=DPT+DP(NK)
2521 QINIT=QINIT+Q0(NK)*EMS(NK)
2522 QFNL=QFNL+QG(NK)*EMS(NK)
2523 QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
2524 315 CONTINUE
2525 QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC) ! PPT FB MODS
2526 ! QFNL=QFNL+PPTFLX*TIMEC ! PPT FB MODS
2527 ERR2=(QFNL-QINIT)*100./QINIT
2528 IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
2529 IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN
2530 ! write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
2531 ! WRITE(99,1110)QINIT,QFNL,ERR2
2532 IPRNT=.TRUE.
2533 ISTOP=1
2534 write(98,4422)kl
2535 4422 format(i6)
2536 DO 311 NK = 1,KL
2537 k = kl - nk + 1
2538 ! write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2539 ! u0(k),v0(k),W0AVG1D(K),dp(k)
2540 ! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2541 ! WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2542 ! U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2543 WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2544 U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2545 311 CONTINUE
2546 ! flush(98)
2548 ! GOTO 297
2549 ! STOP 'QVERR'
2550 ENDIF
2551 1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2552 4456 format(8f12.3)
2553 IF(PPTFLX.GT.0.)THEN
2554 RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
2555 ELSE
2556 RELERR=0.
2557 ENDIF
2558 IF(IPRNT)THEN
2559 WRITE(98,1120)RELERR
2560 WRITE(98,*)'TDER, CPR, TRPPT =', &
2561 TDER,CPR*AINC,TRPPT*AINC
2562 ENDIF
2563 !
2564 !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
2565 !
2566 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
2567 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
2568 !
2569 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
2570 NCA(I,J) = REAL(NIC)*DT
2571 IF(ISHALL.EQ.1)THEN
2572 TIMEC = 2400.
2573 NCA(I,J) = CUDT*60.
2574 NSHALL = NSHALL+1
2575 ENDIF
2577 DO K=1,KX
2578 ! IF(IMOIST(INEST).NE.2)THEN
2579 !
2580 !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMATED
2581 !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
2582 !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
2583 !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
2584 !...OF QG...
2585 !
2586 ! RLC=XLV0-XLV1*TG(K)
2587 ! RLS=XLS0-XLS1*TG(K)
2588 ! CPM=CP*(1.+0.887*QG(K))
2589 ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
2590 ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
2591 ! DQLDT(I,J,NK)=0.
2592 ! DQIDT(I,J,NK)=0.
2593 ! DQRDT(I,J,NK)=0.
2594 ! DQSDT(I,J,NK)=0.
2595 ! ELSE
2596 !
2597 !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
2598 !
2599 IF(warm_rain)THEN
2601 CPM=CP*(1.+0.887*QG(K))
2602 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2603 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2604 DQIDT(K)=0.
2605 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2606 DQSDT(K)=0.
2607 ELSEIF(.NOT. F_QS)THEN
2608 !
2609 !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROMETEORS
2610 !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
2611 !
2612 CPM=CP*(1.+0.887*QG(K))
2613 IF(K.LE.ML)THEN
2614 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2615 ELSEIF(K.GT.ML)THEN
2616 TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
2617 ENDIF
2618 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2619 DQIDT(K)=0.
2620 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2621 DQSDT(K)=0.
2622 ELSEIF(F_QS) THEN
2623 !
2624 !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
2625 !...OF HYDROMETEORS DIRECTLY...
2626 !
2627 DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
2628 DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
2629 DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
2630 IF (F_QI) THEN
2631 DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
2632 ELSE
2633 DQSDT(K)=DQSDT(K)+(QIG(K)-QI0(K))/TIMEC
2634 ENDIF
2635 ELSE
2636 ! PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
2637 CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS MICROPHYSICS CHOICE IS NOT ALLOWED' )
2638 ENDIF
2639 DTDT(K)=(TG(K)-T0(K))/TIMEC
2640 DQDT(K)=(QG(K)-Q0(K))/TIMEC
2641 ENDDO
2642 PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2643 RAINCV(I,J)=DT*PRATEC(I,J)
2644 ! RAINCV(I,J)=.1*.5*DT*PPTFLX/DXSQ ! PPT FB MODS
2645 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
2646 RNC=RAINCV(I,J)*NIC
2647 909 FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
2648 ! write (98,909)I,J,RNC
2649 ! write (6,909)I,J,RNC
2650 ! WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
2651 ! * NCCNT
2652 ! flush(98)
2653 1000 FORMAT(' ',10A8)
2654 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
2655 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
2656 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
2657 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
2658 ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
2659 I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
2660 ' CAPE=',0PF7.1)
2661 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
2662 E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
2663 F8.1)
2664 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
2665 ,F6.3,'VWS=',F5.2)
2666 !1055 FORMAT('*** DEGREE OF STABILIZATION =',F5.3, &
2667 ! ', NO MORE MASS FLUX IS ALLOWED!')
2668 !1060 FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED &
2669 ! &DEGREE OF STABILIZATION! FABE= ',F6.4)
2670 1070 FORMAT (16A8)
2671 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2672 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=', &
2673 2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2674 1085 FORMAT (A3,16A7,2A8)
2675 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2676 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
2677 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
2678 E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
2679 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2680 ' TOTAL WATER CHANGE =',F8.2,'%')
2681 ! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2682 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
2683 !
2684 !-----------------------------------------------------------------------
2685 !--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
2686 !-----------------------------------------------------------------------
2687 !
2688 CUTOP(I,J)=REAL(LTOP)
2689 CUBOT(I,J)=REAL(LCL)
2690 !
2691 !-----------------------------------------------------------------------
2692 END SUBROUTINE KF_eta_PARA
2693 !********************************************************************
2694 ! ***********************************************************************
2695 SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
2696 !
2697 ! Lookup table variables:
2698 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2699 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2700 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2701 ! REAL, SAVE, DIMENSION(1:200) :: ALU
2702 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2703 ! End of Lookup table variables:
2704 !-----------------------------------------------------------------------
2705 IMPLICIT NONE
2706 !-----------------------------------------------------------------------
2707 REAL, INTENT(IN ) :: P,THES,XLV1,XLV0
2708 REAL, INTENT(OUT ) :: QNEWLQ,QNEWIC
2709 REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
2710 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11, &
2711 TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
2712 INTEGER :: IPTB,ITHTB
2713 !-----------------------------------------------------------------------
2715 !c******** LOOKUP TABLE VARIABLES... ****************************
2716 ! parameter(kfnt=250,kfnp=220)
2717 !c
2718 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
2719 ! * alu(200),rdpr,rdthk,plutop
2720 !C***************************************************************
2721 !c
2722 !c***********************************************************************
2723 !c scaling pressure and tt table index
2724 !c***********************************************************************
2725 !c
2726 tp=(p-plutop)*rdpr
2727 qq=tp-aint(tp)
2728 iptb=int(tp)+1
2730 !
2731 !***********************************************************************
2732 ! base and scaling factor for the
2733 !***********************************************************************
2734 !
2735 ! scaling the and tt table index
2736 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2737 tth=(thes-bth)*rdthk
2738 pp =tth-aint(tth)
2739 ithtb=int(tth)+1
2740 IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
2741 write(98,*)'**** OUT OF BOUNDS *********'
2742 ! flush(98)
2743 ENDIF
2744 !
2745 t00=ttab(ithtb ,iptb )
2746 t10=ttab(ithtb+1,iptb )
2747 t01=ttab(ithtb ,iptb+1)
2748 t11=ttab(ithtb+1,iptb+1)
2749 !
2750 q00=qstab(ithtb ,iptb )
2751 q10=qstab(ithtb+1,iptb )
2752 q01=qstab(ithtb ,iptb+1)
2753 q11=qstab(ithtb+1,iptb+1)
2754 !
2755 !***********************************************************************
2756 ! parcel temperature
2757 !***********************************************************************
2758 !
2759 temp=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
2760 !
2761 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
2762 !
2763 DQ=QS-QU
2764 IF(DQ.LE.0.)THEN
2765 QNEW=QU-QS
2766 QU=QS
2767 ELSE
2768 !
2769 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2770 ! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
2771 !
2772 QNEW=0.
2773 QTOT=QLIQ+QICE
2774 !
2775 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
2776 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
2777 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
2778 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
2779 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
2780 !
2781 !...subsaturated values only occur in calculations involving various mixtures of
2782 !...updraft and environmental air for estimation of entrainment and detrainment.
2783 !...For these purposes, assume that reasonable estimates can be given using
2784 !...liquid water saturation calculations only - i.e., ignore the effect of the
2785 !...ice phase in this process only...will not affect conservative properties...
2786 !
2787 IF(QTOT.GE.DQ)THEN
2788 qliq=qliq-dq*qliq/(qtot+1.e-10)
2789 qice=qice-dq*qice/(qtot+1.e-10)
2790 QU=QS
2791 ELSE
2792 RLL=XLV0-XLV1*TEMP
2793 CPP=1004.5*(1.+0.89*QU)
2794 IF(QTOT.LT.1.E-10)THEN
2795 !
2796 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2797 TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
2798 ELSE
2799 !
2800 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
2801 ! THE TEMPERATURE IS GIVEN BY:
2802 !
2803 TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
2804 QU=QU+QTOT
2805 QTOT=0.
2806 QLIQ=0.
2807 QICE=0.
2808 ENDIF
2809 ENDIF
2810 ENDIF
2811 TU=TEMP
2812 qnewlq=qnew
2813 qnewic=0.
2814 !
2815 END SUBROUTINE TPMIX2
2816 !******************************************************************************
2817 SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
2818 !-----------------------------------------------------------------------
2819 IMPLICIT NONE
2820 !-----------------------------------------------------------------------
2821 REAL, INTENT(IN ) :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
2822 REAL, INTENT(INOUT) :: TU,THTEU,QU,QICE
2823 REAL :: RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
2824 !-----------------------------------------------------------------------
2825 !
2826 !...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
2827 !...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
2828 !...TTFRZ TO TBFRZ...
2829 !...FOR COLDER TEMPERATURES, FREEZE ALL LIQUID WATER...
2830 !...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
2831 !...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
2832 !
2833 RLC=2.5E6-2369.276*(TU-273.16)
2834 RLS=2833922.-259.532*(TU-273.16)
2835 RLF=RLS-RLC
2836 CPP=1004.5*(1.+0.89*QU)
2837 !
2838 ! A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
2839 ! FOR SATURATION VAPOR PRESSURE...
2840 !
2841 A=(CLIQ-BLIQ*DLIQ)/((TU-DLIQ)*(TU-DLIQ))
2842 DTFRZ = RLF*QFRZ/(CPP+RLS*QU*A)
2843 TU = TU+DTFRZ
2845 ES = ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2846 QS = ES*0.622/(P-ES)
2847 !
2848 !...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
2849 !...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
2850 !...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
2851 !...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
2852 !...TEMPERATURE TO THE SATURATION VALUE...
2853 !
2854 DQEVAP = QS-QU
2855 QICE = QICE-DQEVAP
2856 QU = QU+DQEVAP
2857 PII=(1.E5/P)**(0.2854*(1.-0.28*QU))
2858 THTEU=TU*PII*EXP((3374.6525/TU-2.5403)*QU*(1.+0.81*QU))
2859 !
2860 END SUBROUTINE DTFRZNEW
2861 ! --------------------------------------------------------------------------------
2863 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
2864 QNEWIC,QLQOUT,QICOUT,G)
2866 !-----------------------------------------------------------------------
2867 IMPLICIT NONE
2868 !-----------------------------------------------------------------------
2869 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2870 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2871 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
2872 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2873 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
2875 REAL, INTENT(IN ) :: G
2876 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
2877 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2878 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2879 !
2880 QTOT=QLIQ+QICE
2881 QNEW=QNEWLQ+QNEWIC
2882 !
2883 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
2884 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2885 ! LEVELS...
2886 !
2887 QEST=0.5*(QTOT+QNEW)
2888 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
2889 IF(G1.LT.0.0)G1=0.
2890 WAVG=0.5*(SQRT(WTW)+SQRT(G1))
2891 CONV=RATE*DZ/WAVG ! KF90 Eq. 9
2892 !
2893 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2894 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2895 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2896 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2897 !
2898 RATIO3=QNEWLQ/(QNEW+1.E-8)
2899 ! OLDQ=QTOT
2900 QTOT=QTOT+0.6*QNEW
2901 OLDQ=QTOT
2902 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
2903 QTOT=QTOT*EXP(-CONV) ! KF90 Eq. 9
2904 !
2905 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2906 ! PARCEL AT THIS LEVEL...
2907 !
2908 DQ=OLDQ-QTOT
2909 QLQOUT=RATIO4*DQ
2910 QICOUT=(1.-RATIO4)*DQ
2911 !
2912 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2913 ! LATE VERTICAL VELOCITY
2914 !
2915 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
2916 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
2917 IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
2918 !
2919 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2920 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2921 !
2922 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
2923 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
2924 QNEWLQ=0.
2925 QNEWIC=0.
2927 END SUBROUTINE CONDLOAD
2929 ! ----------------------------------------------------------------------
2930 SUBROUTINE PROF5(EQ,EE,UD)
2931 !
2932 !***********************************************************************
2933 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
2934 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2935 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2936 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
2937 ! "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
2938 ! ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
2939 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2940 ! JACK KAIN
2941 ! 7/6/89
2942 ! Solves for KF90 Eq. 2
2943 !
2944 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2945 !-----------------------------------------------------------------------
2946 IMPLICIT NONE
2947 !-----------------------------------------------------------------------
2948 REAL, INTENT(IN ) :: EQ
2949 REAL, INTENT(INOUT) :: EE,UD
2950 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2952 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2953 0.9372980,0.33267,0.166666667,0.202765151/
2954 X=(EQ-0.5)/SIGMA
2955 Y=6.*EQ-3.
2956 EY=EXP(Y*Y/(-2))
2957 E45=EXP(-4.5)
2958 T2=1./(1.+P*ABS(Y))
2959 T1=0.500498
2960 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
2961 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
2962 IF(Y.GE.0.)THEN
2963 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2964 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
2965 EQ)
2966 ELSE
2967 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2968 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
2969 EQ/2.-EQ)
2970 ENDIF
2971 EE=EE/FE
2972 UD=UD/FE
2974 END SUBROUTINE PROF5
2976 ! ------------------------------------------------------------------------
2977 SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
2978 !
2979 ! Lookup table variables:
2980 ! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2981 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2982 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2983 ! REAL, SAVE, DIMENSION(1:200) :: ALU
2984 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2985 ! End of Lookup table variables:
2986 !-----------------------------------------------------------------------
2987 IMPLICIT NONE
2988 !-----------------------------------------------------------------------
2989 REAL, INTENT(IN ) :: P,THES
2990 REAL, INTENT(INOUT) :: TS,QS
2991 INTEGER, INTENT(IN ) :: i,j ! avail for debugging
2992 REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
2993 INTEGER :: IPTB,ITHTB
2994 CHARACTER*256 :: MESS
2995 !-----------------------------------------------------------------------
2997 !
2998 !******** LOOKUP TABLE VARIABLES (F77 format)... ****************************
2999 ! parameter(kfnt=250,kfnp=220)
3000 !
3001 ! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
3002 ! alu(200),rdpr,rdthk,plutop
3003 !***************************************************************
3004 !
3005 !***********************************************************************
3006 ! scaling pressure and tt table index
3007 !***********************************************************************
3008 !
3009 tp=(p-plutop)*rdpr
3010 qq=tp-aint(tp)
3011 iptb=int(tp)+1
3012 !
3013 !***********************************************************************
3014 ! base and scaling factor for the
3015 !***********************************************************************
3016 !
3017 ! scaling the and tt table index
3018 bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
3019 tth=(thes-bth)*rdthk
3020 pp =tth-aint(tth)
3021 ithtb=int(tth)+1
3022 !
3023 t00=ttab(ithtb ,iptb )
3024 t10=ttab(ithtb+1,iptb )
3025 t01=ttab(ithtb ,iptb+1)
3026 t11=ttab(ithtb+1,iptb+1)
3027 !
3028 q00=qstab(ithtb ,iptb )
3029 q10=qstab(ithtb+1,iptb )
3030 q01=qstab(ithtb ,iptb+1)
3031 q11=qstab(ithtb+1,iptb+1)
3032 !
3033 !***********************************************************************
3034 ! parcel temperature and saturation mixing ratio
3035 !***********************************************************************
3036 !
3037 ts=(t00+(t10-t00)*pp+(t01-t00)*qq+(t00-t10-t01+t11)*pp*qq)
3038 !
3039 qs=(q00+(q10-q00)*pp+(q01-q00)*qq+(q00-q10-q01+q11)*pp*qq)
3040 !
3041 END SUBROUTINE TPMIX2DD
3043 ! -----------------------------------------------------------------------
3044 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
3045 !
3046 !-----------------------------------------------------------------------
3047 IMPLICIT NONE
3048 !-----------------------------------------------------------------------
3049 REAL, INTENT(IN ) :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
3050 REAL, INTENT(INOUT) :: THT1
3051 REAL :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT, &
3052 T00,P00,C1,C2,C3,C4,C5
3053 INTEGER :: INDLU
3054 !-----------------------------------------------------------------------
3055 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
3056 0.278296,1.0723E-3/
3057 !
3058 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
3059 !
3060 ! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
3061 ! For example, KF90 Eq. 10 no longer used
3062 !
3063 EE=Q1*P1/(0.622+Q1)
3064 ! TLOG=ALOG(EE/ALIQ)
3065 ! ...calculate LOG term using lookup table...
3066 !
3067 astrt=1.e-3
3068 ainc=0.075
3069 a1=ee/aliq
3070 tp=(a1-astrt)/ainc
3071 indlu=int(tp)+1
3072 value=(indlu-1)*ainc+astrt
3073 aintrp=(a1-value)/ainc
3074 tlog=aintrp*alu(indlu+1)+(1-aintrp)*alu(indlu)
3075 !
3076 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
3077 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
3078 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
3079 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
3080 !
3081 END SUBROUTINE ENVIRTHT
3082 ! ***********************************************************************
3083 !====================================================================
3084 SUBROUTINE kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
3085 RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
3086 SVP1,SVP2,SVP3,SVPT0, &
3087 P_FIRST_SCALAR,restart,allowed_to_read, &
3088 ids, ide, jds, jde, kds, kde, &
3089 ims, ime, jms, jme, kms, kme, &
3090 its, ite, jts, jte, kts, kte )
3091 !--------------------------------------------------------------------
3092 IMPLICIT NONE
3093 !--------------------------------------------------------------------
3094 LOGICAL , INTENT(IN) :: restart,allowed_to_read
3095 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
3096 ims, ime, jms, jme, kms, kme, &
3097 its, ite, jts, jte, kts, kte
3098 INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
3100 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
3101 RTHCUTEN, &
3102 RQVCUTEN, &
3103 RQCCUTEN, &
3104 RQRCUTEN, &
3105 RQICUTEN, &
3106 RQSCUTEN
3108 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
3110 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
3112 INTEGER :: i, j, k, itf, jtf, ktf
3113 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
3115 jtf=min0(jte,jde-1)
3116 ktf=min0(kte,kde-1)
3117 itf=min0(ite,ide-1)
3119 IF(.not.restart)THEN
3121 DO j=jts,jtf
3122 DO k=kts,ktf
3123 DO i=its,itf
3124 RTHCUTEN(i,k,j)=0.
3125 RQVCUTEN(i,k,j)=0.
3126 RQCCUTEN(i,k,j)=0.
3127 RQRCUTEN(i,k,j)=0.
3128 ENDDO
3129 ENDDO
3130 ENDDO
3132 IF (P_QI .ge. P_FIRST_SCALAR) THEN
3133 DO j=jts,jtf
3134 DO k=kts,ktf
3135 DO i=its,itf
3136 RQICUTEN(i,k,j)=0.
3137 ENDDO
3138 ENDDO
3139 ENDDO
3140 ENDIF
3142 IF (P_QS .ge. P_FIRST_SCALAR) THEN
3143 DO j=jts,jtf
3144 DO k=kts,ktf
3145 DO i=its,itf
3146 RQSCUTEN(i,k,j)=0.
3147 ENDDO
3148 ENDDO
3149 ENDDO
3150 ENDIF
3152 DO j=jts,jtf
3153 DO i=its,itf
3154 NCA(i,j)=-100.
3155 ENDDO
3156 ENDDO
3158 DO j=jts,jtf
3159 DO k=kts,ktf
3160 DO i=its,itf
3161 W0AVG(i,k,j)=0.
3162 ENDDO
3163 ENDDO
3164 ENDDO
3166 endif
3168 CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
3170 END SUBROUTINE kf_eta_init
3172 !-------------------------------------------------------
3174 subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
3175 !
3176 ! This subroutine is a lookup table.
3177 ! Given a series of series of saturation equivalent potential
3178 ! temperatures, the temperature is calculated.
3179 !
3180 !--------------------------------------------------------------------
3181 IMPLICIT NONE
3182 !--------------------------------------------------------------------
3183 ! Lookup table variables
3184 ! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
3185 ! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3186 ! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3187 ! REAL, SAVE, DIMENSION(1:200) :: ALU
3188 ! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3189 ! End of Lookup table variables
3191 INTEGER :: KP,IT,ITCNT,I
3192 REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
3193 TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
3194 ASTRT,AINC,A1,THTGS
3195 ! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
3196 REAL :: ALIQ,BLIQ,CLIQ,DLIQ
3197 REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
3198 !
3199 ! equivalent potential temperature increment
3200 data dth/1./
3201 ! minimum starting temp
3202 data tmin/150./
3203 ! tolerance for accuracy of temperature
3204 data toler/0.001/
3205 ! top pressure (pascals)
3206 plutop=5000.0
3207 ! bottom pressure (pascals)
3208 pbot=110000.0
3210 ALIQ = SVP1*1000.
3211 BLIQ = SVP2
3212 CLIQ = SVP2*SVPT0
3213 DLIQ = SVP3
3215 !
3216 ! compute parameters
3217 !
3218 ! 1._over_(sat. equiv. theta increment)
3219 rdthk=1./dth
3220 ! pressure increment
3221 !
3222 DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
3223 ! dpr=(pbot-plutop)/REAL(kfnp-1)
3224 ! 1._over_(pressure increment)
3225 rdpr=1./dpr
3226 ! compute the spread of thes
3227 ! thespd=dth*(kfnt-1)
3228 !
3229 ! calculate the starting sat. equiv. theta
3230 !
3231 temp=tmin
3232 p=plutop-dpr
3233 do kp=1,kfnp
3234 p=p+dpr
3235 es=aliq*exp((bliq*temp-cliq)/(temp-dliq))
3236 qs=0.622*es/(p-es)
3237 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
3238 the0k(kp)=temp*pi*exp((3374.6525/temp-2.5403)*qs* &
3239 (1.+0.81*qs))
3240 enddo
3241 !
3242 ! compute temperatures for each sat. equiv. potential temp.
3243 !
3244 p=plutop-dpr
3245 do kp=1,kfnp
3246 thes=the0k(kp)-dth
3247 p=p+dpr
3248 do it=1,kfnt
3249 ! define sat. equiv. pot. temp.
3250 thes=thes+dth
3251 ! iterate to find temperature
3252 ! find initial guess
3253 if(it.eq.1) then
3254 tgues=tmin
3255 else
3256 tgues=ttab(it-1,kp)
3257 endif
3258 es=aliq*exp((bliq*tgues-cliq)/(tgues-dliq))
3259 qs=0.622*es/(p-es)
3260 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
3261 thgues=tgues*pi*exp((3374.6525/tgues-2.5403)*qs* &
3262 (1.+0.81*qs))
3263 f0=thgues-thes
3264 t1=tgues-0.5*f0
3265 t0=tgues
3266 itcnt=0
3267 ! iteration loop
3268 do itcnt=1,11
3269 es=aliq*exp((bliq*t1-cliq)/(t1-dliq))
3270 qs=0.622*es/(p-es)
3271 pi=(1.e5/p)**(0.2854*(1.-0.28*qs))
3272 thtgs=t1*pi*exp((3374.6525/t1-2.5403)*qs*(1.+0.81*qs))
3273 f1=thtgs-thes
3274 if(abs(f1).lt.toler)then
3275 exit
3276 endif
3277 ! itcnt=itcnt+1
3278 dt=f1*(t1-t0)/(f1-f0)
3279 t0=t1
3280 f0=f1
3281 t1=t1-dt
3282 enddo
3283 ttab(it,kp)=t1
3284 qstab(it,kp)=qs
3285 enddo
3286 enddo
3287 !
3288 ! lookup table for tlog(emix/aliq)
3289 !
3290 ! set up intial values for lookup tables
3291 !
3292 astrt=1.e-3
3293 ainc=0.075
3294 !
3295 a1=astrt-ainc
3296 do i=1,200
3297 a1=a1+ainc
3298 alu(i)=alog(a1)
3299 enddo
3300 !
3301 END SUBROUTINE KF_LUTAB
3303 END MODULE module_cu_kfeta