1 !WRF:MODEL_LAYER:PHYSICS
2 !
3
4 MODULE module_cu_kf
5
6 USE module_wrf_error
7
8 REAL , PARAMETER :: RAD = 1500.
9
10 CONTAINS
11
12 !-------------------------------------------------------------
13 SUBROUTINE KFCPS( &
14 ids,ide, jds,jde, kds,kde &
15 ,ims,ime, jms,jme, kms,kme &
16 ,its,ite, jts,jte, kts,kte &
17 ,DT,KTAU,DX &
18 ,rho &
19 ,RAINCV,NCA &
20 ,U,V,TH,T,W,QV,dz8w,Pcps,pi &
21 ,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
22 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
23 ,STEPCU,CU_ACT_FLAG,warm_rain &
24 ! optional arguments
25 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
26 ,RQVCUTEN,RQCCUTEN,RQRCUTEN,RQICUTEN,RQSCUTEN &
27 ,RTHCUTEN &
28 )
29 !-------------------------------------------------------------
30 IMPLICIT NONE
31 !-------------------------------------------------------------
32 INTEGER, INTENT(IN ) :: &
33 ids,ide, jds,jde, kds,kde, &
34 ims,ime, jms,jme, kms,kme, &
35 its,ite, jts,jte, kts,kte
36
37 INTEGER, INTENT(IN ) :: STEPCU
38 LOGICAL, INTENT(IN ) :: warm_rain
39
40 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
41 REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
42 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
43
44 INTEGER, INTENT(IN ) :: KTAU
45
46 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
47 INTENT(IN ) :: &
48 U, &
49 V, &
50 W, &
51 TH, &
52 QV, &
53 T, &
54 dz8w, &
55 Pcps, &
56 rho, &
57 pi
58 !
59 REAL, INTENT(IN ) :: DT, DX
60
61 REAL, DIMENSION( ims:ime , jms:jme ), &
62 INTENT(INOUT) :: &
63 RAINCV &
64 , NCA
65
66 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
67 INTENT(INOUT) :: &
68 W0AVG
69
70 LOGICAL, DIMENSION( ims:ime , jms:jme ), &
71 INTENT(INOUT) :: CU_ACT_FLAG
72
73 !
74 ! Optional arguments
75 !
76
77 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
78 OPTIONAL, &
79 INTENT(INOUT) :: &
80 RTHCUTEN &
81 ,RQVCUTEN &
82 ,RQCCUTEN &
83 ,RQRCUTEN &
84 ,RQICUTEN &
85 ,RQSCUTEN
86
87 !
88 ! Flags relating to the optional tendency arrays declared above
89 ! Models that carry the optional tendencies will provdide the
90 ! optional arguments at compile time; these flags all the model
91 ! to determine at run-time whether a particular tracer is in
92 ! use or not.
93 !
94
95 LOGICAL, OPTIONAL :: &
96 F_QV &
97 ,F_QC &
98 ,F_QR &
99 ,F_QI &
100 ,F_QS
101
102
103
104 ! LOCAL VARS
105
106 REAL, DIMENSION( kts:kte ) :: &
107 U1D, &
108 V1D, &
109 T1D, &
110 DZ1D, &
111 QV1D, &
112 P1D, &
113 RHO1D, &
114 W0AVG1D
115
116 REAL, DIMENSION( kts:kte ):: &
117 DQDT, &
118 DQIDT, &
119 DQCDT, &
120 DQRDT, &
121 DQSDT, &
122 DTDT
123
124 REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp,RTHCUMAX
125
126 INTEGER :: i,j,k,NTST,ICLDCK
127
128 LOGICAL :: qi_flag , qs_flag
129
130 !----------------------------------------------------------------------
131
132 !--- CALL CUMULUS PARAMETERIZATION
133 !
134 !...TST IS THE NUMBER OF TIME STEPS IN 10 MINUTES...W0AVG IS CLOSE TO A
135 !...RUNNING MEAN VERTICAL VELOCITY...NOTE THAT IF YOU CHANGE TST, IT WIL
136 !...CHANGE THE FREQUENCY OF THE CONVECTIVE INTITIATION CHECK (SEE BELOW)
137 !...NOTE THAT THE ORDERING OF VERTICAL LAYERS MUST BE REVERSED FOR W0AVG
138 !...BECAUSE THE ORDERING IS REVERSED IN KFPARA...
139 !
140 DXSQ=DX*DX
141 qi_flag = .FALSE.
142 qs_flag = .FALSE.
143 IF ( PRESENT( F_QI ) ) qi_flag = f_qi
144 IF ( PRESENT( F_QS ) ) qs_flag = f_qs
145
146 !----------------------
147 NTST=STEPCU
148 TST=float(NTST*2)
149 !----------------------
150 ! NTST=NINT(1200./(DT*2.))
151 ! TST=float(NTST)
152 ! NTST=NINT(0.5*TST)
153 ! NTST=MAX0(NTST,1)
154 !----------------------
155 ICLDCK=MOD(KTAU,NTST)
156 !----------------------
157
158 DO J = jts,jte
159 DO K=kts,kte
160 DO I= its,ite
161 SCR1=-5.0E-4*G*rho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
162 TV=T(I,K,J)*(1.+EP1*QV(I,K,J))
163 RHOE=Pcps(I,K,J)/(R*TV)
164 W0=-101.9368*SCR1/RHOE
165 W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
166 ENDDO
167 ENDDO
168 ENDDO
169 !
170 !...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
171 !
172 RTHCUMAX=0.
173
174 IF(ICLDCK.EQ.0 .or. KTAU .eq. 1) then
175
176 DO J = jts,jte
177 DO I= its,ite
178 CU_ACT_FLAG(i,j) = .true.
179 ENDDO
180 ENDDO
181
182 DO J = jts,jte
183 DO I=its,ite
184 IF ( NINT(NCA(I,J)) .gt. 0 ) then
185 CU_ACT_FLAG(i,j) = .false.
186 ELSE
187
188 DO k=kts,kte
189 DQDT(k)=0.
190 DQIDT(k)=0.
191 DQCDT(k)=0.
192 DQRDT(k)=0.
193 DQSDT(k)=0.
194 DTDT(k)=0.
195 ENDDO
196 RAINCV(I,J)=0.
197 !
198 ! assign vars from 3D to 1D
199
200 DO K=kts,kte
201 U1D(K) =U(I,K,J)
202 V1D(K) =V(I,K,J)
203 T1D(K) =T(I,K,J)
204 RHO1D(K) =rho(I,K,J)
205 QV1D(K)=QV(I,K,J)
206 P1D(K) =Pcps(I,K,J)
207 W0AVG1D(K) =W0AVG(I,K,J)
208 DZ1D(k)=dz8w(I,K,J)
209 ENDDO
210
211 !
212 CALL KFPARA(I, J, &
213 U1D,V1D,T1D,QV1D,P1D,DZ1D, &
214 W0AVG1D,DT,DX,DXSQ,RHO1D, &
215 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
216 EP2,SVP1,SVP2,SVP3,SVPT0, &
217 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
218 RAINCV,NCA, &
219 warm_rain,qi_flag,qs_flag, &
220 ids,ide, jds,jde, kds,kde, &
221 ims,ime, jms,jme, kms,kme, &
222 its,ite, jts,jte, kts,kte )
223
224 IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN ) ) THEN
225 DO K=kts,kte
226 RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
227 RQVCUTEN(I,K,J)=DQDT(K)
228 ENDDO
229 ENDIF
230
231 IF( PRESENT(RQRCUTEN) .AND. PRESENT(RQCCUTEN) .AND. &
232 PRESENT(F_QR) ) THEN
233 IF ( F_QR ) THEN
234 DO K=kts,kte
235 RQRCUTEN(I,K,J)=DQRDT(K)
236 RQCCUTEN(I,K,J)=DQCDT(K)
237 ENDDO
238 ELSE
239 ! This is the case for Eta microphysics without 3d rain field
240 DO K=kts,kte
241 RQRCUTEN(I,K,J)=0.
242 RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
243 ENDDO
244 ENDIF
245 ENDIF
246
247 !...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
248
249 IF( PRESENT( RQICUTEN ) .AND. qi_flag )THEN
250 DO K=kts,kte
251 RQICUTEN(I,K,J)=DQIDT(K)
252 ENDDO
253 ENDIF
254
255 IF( PRESENT ( RQSCUTEN ) .AND. qs_flag )THEN
256 DO K=kts,kte
257 RQSCUTEN(I,K,J)=DQSDT(K)
258 ENDDO
259 ENDIF
260 !
261 ENDIF
262 ENDDO
263 ENDDO
264
265 ENDIF
266
267 END SUBROUTINE KFCPS
268
269 !-----------------------------------------------------------
270 SUBROUTINE KFPARA (I, J, &
271 U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
272 DT,DX,DXSQ,rho, &
273 XLV0,XLV1,XLS0,XLS1,CP,R,G, &
274 EP2,SVP1,SVP2,SVP3,SVPT0, &
275 DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
276 RAINCV,NCA, &
277 warm_rain,qi_flag,qs_flag, &
278 ids,ide, jds,jde, kds,kde, &
279 ims,ime, jms,jme, kms,kme, &
280 its,ite, jts,jte, kts,kte )
281 !-----------------------------------------------------------
282 IMPLICIT NONE
283 !-----------------------------------------------------------
284 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
285 ims,ime, jms,jme, kms,kme, &
286 its,ite, jts,jte, kts,kte, &
287 I,J
288 LOGICAL, INTENT(IN ) :: warm_rain
289 LOGICAL :: qi_flag, qs_flag
290
291 !
292 REAL, DIMENSION( kts:kte ), &
293 INTENT(IN ) :: U0, &
294 V0, &
295 T0, &
296 QV0, &
297 P0, &
298 rho, &
299 DZQ, &
300 W0AVG1D
301 !
302 REAL, INTENT(IN ) :: DT,DX,DXSQ
303 !
304
305 REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
306 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
307 !
308 REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
309 DQDT, &
310 DQIDT, &
311 DQCDT, &
312 DQRDT, &
313 DQSDT, &
314 DTDT
315
316 REAL, DIMENSION( ims:ime , jms:jme ), &
317 INTENT(INOUT) :: RAINCV, &
318 NCA
319 !
320 !...DEFINE LOCAL VARIABLES...
321 !
322 REAL, DIMENSION( kts:kte ) :: &
323 Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
324 QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
325 UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
326 UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
327 THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
328 QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
329 DETLQ2,DETIC2,RATIO,RATIO2
330
331 REAL, DIMENSION( kts:kte ) :: &
332 DOMGDP,EXN,RHOE,TVQU,DP,RH,EQFRC,WSPD, &
333 QDT,FXM,THTAG,THTESG,THPA,THFXTOP, &
334 THFXBOT,QPA,QFXTOP,QFXBOT,QLPA,QLFXIN, &
335 QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
336 QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
337 QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
338
339 REAL, DIMENSION( kts:kte+1 ) :: OMG
340 REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
341
342 ! LOCAL VARS
343
344 REAL :: P00,T00,CV,B61,RLF,RHIC,RHBC,PIE, &
345 TTFRZ,TBFRZ,C5,RATE
346 REAL :: GDRY,ROCP,ALIQ,BLIQ, &
347 CLIQ,DLIQ,AICE,BICE,CICE,DICE
348 REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
349 ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
350 CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
351 ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
352 TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
353 UPNEW,ABE,WKLCL,THTUDL,TUDL,TTEMP,FRC1, &
354 QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
355 DZZ,WSQ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
356 THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
357 UD1,CLDHGT,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
358 THTA,P150,USR,VCONV,TIMEC,SHSIGN,VWS,PEF, &
359 CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
360 DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
361 DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
362 UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
363 DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
364 AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
365 DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
366 TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
367 UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
368 RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
369 DDFRC,TDC,DEFRC
370
371 INTEGER :: KX,K,KL
372 !
373 INTEGER :: ISTOP,ML,L5,L4,KMIX,LOW, &
374 LC,MXLAYR,LLFC,NLAYRS,NK, &
375 KPBL,KLCL,LCL,LET,IFLAG, &
376 KFRZ,NK1,LTOP,NJ,LTOP1, &
377 LTOPM1,LVF,KSTART,KMIN,LFS, &
378 ND,NIC,LDB,LDT,ND1,NDK, &
379 NM,LMAX,NCOUNT,NOITR, &
380 NSTEP,NTC
381 !
382 DATA P00,T00/1.E5,273.16/
383 DATA CV,B61,RLF/717.,0.608,3.339E5/
384 DATA RHIC,RHBC/1.,0.90/
385 DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
386 DATA RATE/0.01/
387 !-----------------------------------------------------------
388 GDRY=-G/CP
389 ROCP=R/CP
390 KL=kte
391 KX=kte
392 !
393 ! ALIQ = 613.3
394 ! BLIQ = 17.502
395 ! CLIQ = 4780.8
396 ! DLIQ = 32.19
397 ALIQ = SVP1*1000.
398 BLIQ = SVP2
399 CLIQ = SVP2*SVPT0
400 DLIQ = SVP3
401 AICE = 613.2
402 BICE = 22.452
403 CICE = 6133.0
404 DICE = 0.61
405 !
406
407 !...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER
408 !...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)
409 !...FIELD. 'FBFRC' IS THE FRACTION OF AVAILABLE
410 !...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...
411 !
412 FBFRC=0.0
413 !
414 !...SCHEME IS CALLED ONCE ON EACH NORTH-SOUTH SLICE, THE LOOP BELOW
415 !...CHECKS FOR THE POSSIBILITY OF INITIATING PARAMETERIZED
416 !...CONVECTION AT EACH POINT WITHIN THE SLICE
417 !
418 !...SEE IF IT IS NECESSARY TO CHECK FOR CONVECTIVE TRIGGERING AT THIS
419 !...GRID POINT. IF NCA>0, CONVECTION IS ALREADY ACTIVE AT THIS POINT,
420 !...JUST FEED BACK THE TENDENCIES SAVED FROM THE TIME WHEN CONVECTION
421 !...WAS INITIATED. IF NCA<0, CONVECTION IS NOT ACTIVE
422 !...AND YOU MAY WANT TO CHECK TO SEE IF IT CAN BE ACTIVATED FOR THE
423 !...CURRENT CONDITIONS. IN PREVIOUS APLICATIONS OF THIS SCHEME,
424 !...THE VARIABLE ICLDCK WAS USED BELOW TO SAVE TIME BY ONLY CHECKING
425 !...FOR THE POSSIBILITY OF CONVECTIVE INITIATION EVERY 5 OR 10
426 !...MINUTES...
427 !
428
429 ! 10 CONTINUE
430 !SUE P300=1000.*(PSB(I,J)*A(KL)+PTOP-30.)+PP3D(I,J,KL)
431
432 P300=P0(1)-30000.
433 !
434 !...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
435 !...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
436 !...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
437 !
438 !...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
439 !...FROM BOTTOM-UP IN THE KF SCHEME...
440 !
441 ML=0
442 !SUE tmprpsb=1./PSB(I,J)
443 !SUE CELL=PTOP*tmprpsb
444
445 DO 15 K=1,KX
446 !SUE P0(K)=1.E3*(A(NK)*PSB(I,J)+PTOP)+PP3D(I,J,NK)
447 !
448 !...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
449 !
450 ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
451 QES(K)=EP2*ES/(P0(K)-ES)
452 Q0(K)=AMIN1(QES(K),QV0(K))
453 Q0(K)=AMAX1(0.000001,Q0(K))
454 QL0(K)=0.
455 QI0(K)=0.
456 QR0(K)=0.
457 QS0(K)=0.
458
459 TV0(K)=T0(K)*(1.+B61*Q0(K))
460 RHOE(K)=P0(K)/(R*TV0(K))
461
462 DP(K)=rho(k)*g*DZQ(k)
463 !
464 !...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
465 ! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
466 !
467 IF(P0(K).GE.500E2)L5=K
468 IF(P0(K).GE.400E2)L4=K
469 IF(P0(K).GE.P300)LLFC=K
470 IF(T0(K).GT.T00)ML=K
471 15 CONTINUE
472
473 Z0(1)=.5*DZQ(1)
474 DO 20 K=2,KL
475 Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
476 DZA(K-1)=Z0(K)-Z0(K-1)
477 20 CONTINUE
478 DZA(KL)=0.
479 KMIX=1
480 25 LOW=KMIX
481
482 IF(LOW.GT.LLFC)GOTO 325
483
484 LC=LOW
485 MXLAYR=0
486 !
487 !...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
488 !...UNSTABLE AIR 50 TO 100 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
489 !...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
490 !...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 60 mb..
491 !
492 NLAYRS=0
493 DPTHMX=0.
494 DO 63 NK=LC,KX
495 DPTHMX=DPTHMX+DP(NK)
496 NLAYRS=NLAYRS+1
497 63 IF(DPTHMX.GT.6.E3)GOTO 64
498 GOTO 325
499 64 KPBL=LC+NLAYRS-1
500 KMIX=LC+1
501 18 THMIX=0.
502 QMIX=0.
503 ZMIX=0.
504 PMIX=0.
505 DPTHMX=0.
506 !
507 !...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
508 !...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
509 !...LAYERS...
510 !
511 DO 17 NK=LC,KPBL
512 DPTHMX=DPTHMX+DP(NK)
513 ROCPQ=0.2854*(1.-0.28*Q0(NK))
514 THMIX=THMIX+DP(NK)*T0(NK)*(P00/P0(NK))**ROCPQ
515 QMIX=QMIX+DP(NK)*Q0(NK)
516 ZMIX=ZMIX+DP(NK)*Z0(NK)
517 17 PMIX=PMIX+DP(NK)*P0(NK)
518 THMIX=THMIX/DPTHMX
519 QMIX=QMIX/DPTHMX
520 ZMIX=ZMIX/DPTHMX
521 PMIX=PMIX/DPTHMX
522 ROCPQ=0.2854*(1.-0.28*QMIX)
523 TMIX=THMIX*(PMIX/P00)**ROCPQ
524 EMIX=QMIX*PMIX/(EP2+QMIX)
525 !
526 !...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL, PRESSURE
527 !...LEVEL OF LCL...
528 !
529 TLOG=ALOG(EMIX/ALIQ)
530 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
531 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
532 TDPT)
533 TLCL=AMIN1(TLCL,TMIX)
534 TVLCL=TLCL*(1.+0.608*QMIX)
535 CPORQ=1./ROCPQ
536 PLCL=P00*(TLCL/THMIX)**CPORQ
537 DO 29 NK=LC,KL
538 KLCL=NK
539 IF(PLCL.GE.P0(NK))GOTO 35
540 29 CONTINUE
541 GOTO 325
542 35 K=KLCL-1
543 DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
544 !
545 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
546 !
547 TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
548 QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
549 TVEN=TENV*(1.+0.608*QENV)
550 TVBAR=0.5*(TV0(K)+TVEN)
551 ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
552 ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
553 !
554 !...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
555 !...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0AVG IS AN
556 !...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
557 !...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
558 !...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
559 !...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
560 !...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
561 !...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
562 !...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
563 !
564 WKLCL=0.02*ZLCL/2.5E3
565 WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))*DLP)*DX/25.E3- &
566 WKLCL
567 WABS=ABS(WKL)+1.E-10
568 WSIGNE=WKL/WABS
569 DTLCL=4.64*WSIGNE*WABS**0.33
570 GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
571 WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
572 IF(TLCL+DTLCL.GT.TENV)GOTO 45
573 IF(KPBL.GE.LLFC)GOTO 325
574 GOTO 25
575 !
576 !...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
577 !...EQUIVALENT POTENTIAL TEMPERATURE
578 !...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
579 !
580 45 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
581 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
582 ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
583 TVAVG=0.5*(TV0(KLCL)+TENV*(1.+0.608*QENV))
584 PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
585 QESE=EP2*ES/(PLCL-ES)
586 GDT=G*DTLCL*(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
587 WLCL=1.+.5*WSIGNE*SQRT(ABS(GDT)+1.E-10)
588 THTES(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
589 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
590 WTW=WLCL*WLCL
591 IF(WLCL.LT.0.)GOTO 25
592 TVLCL=TLCL*(1.+0.608*QMIX)
593 RHOLCL=PLCL/(R*TVLCL)
594 !
595 LCL=KLCL
596 LET=LCL
597 !
598 !*******************************************************************
599 ! *
600 ! COMPUTE UPDRAFT PROPERTIES *
601 ! *
602 !*******************************************************************
603 !
604 !
605 !...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
606 !
607 WU(K)=WLCL
608 AU0=PIE*RAD*RAD
609 UMF(K)=RHOLCL*AU0
610 VMFLCL=UMF(K)
611 UPOLD=VMFLCL
612 UPNEW=UPOLD
613 !
614 !...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
615 !...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE BUOYANT ENERGY,
616 ! TRPPT IS THE TOTAL RATE OF PRECIPITATION PRODUCTION...
617 !
618 RATIO2(K)=0.
619 UER(K)=0.
620 ABE=0.
621 TRPPT=0.
622 TU(K)=TLCL
623 TVU(K)=TVLCL
624 QU(K)=QMIX
625 EQFRC(K)=1.
626 QLIQ(K)=0.
627 QICE(K)=0.
628 QLQOUT(K)=0.
629 QICOUT(K)=0.
630 DETLQ(K)=0.
631 DETIC(K)=0.
632 PPTLIQ(K)=0.
633 PPTICE(K)=0.
634 IFLAG=0
635 KFRZ=LC
636 !
637 !...THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH
638 ! RESPECT TO UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILUTE
639 ! PARCEL IS THTUDL, UNDILUTE TEMPERATURE IS GIVEN BY TUDL...
640 !
641 THTUDL=THETEU(K)
642 TUDL=TLCL
643 !
644 !...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
645 ! PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
646 ! FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
647 ! INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
648 ! PREVIOUS MODEL LEVEL...
649 !
650 TTEMP=TTFRZ
651 !
652 !...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
653 ! MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
654 ! MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
655 !
656 DO 60 NK=K,KL-1
657 NK1=NK+1
658 RATIO2(NK1)=RATIO2(NK)
659 !
660 !...UPDATE UPDRAFT PROPERTIES AT THE NEXT MODEL LVL TO REFLECT
661 ! ENTRAINMENT OF ENVIRONMENTAL AIR...
662 !
663 FRC1=0.
664 TU(NK1)=T0(NK1)
665 THETEU(NK1)=THETEU(NK)
666 QU(NK1)=QU(NK)
667 QLIQ(NK1)=QLIQ(NK)
668 QICE(NK1)=QICE(NK)
669
670 CALL TPMIX(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),QLIQ(NK1), &
671 QICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0, &
672 XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
673 TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
674 !
675 !...CHECK TO SEE IF UPDRAFT TEMP IS WITHIN THE FREEZING INTERVAL,
676 ! IF IT IS, CALCULATE THE FRACTIONAL CONVERSION TO GLACIATION
677 ! AND ADJUST QNEWLQ TO REFLECT THE GRADUAL CHANGE IN THETAU
678 ! SINCE THE LAST MODEL LEVEL...THE GLACIATION EFFECTS WILL BE
679 ! DETERMINED AFTER THE AMOUNT OF CONDENSATE AVAILABLE AFTER
680 ! PRECIP FALLOUT IS DETERMINED...TTFRZ IS THE TEMP AT WHICH
681 ! GLACIATION BEGINS, TBFRZ THE TEMP AT WHICH IT ENDS...
682 !
683 IF(TU(NK1).LE.TTFRZ.AND.IFLAG.LT.1)THEN
684 IF(TU(NK1).GT.TBFRZ)THEN
685 IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
686 FRC1=(TTEMP-TU(NK1))/(TTFRZ-TBFRZ)
687 R1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
688 ELSE
689 FRC1=(TTEMP-TBFRZ)/(TTFRZ-TBFRZ)
690 R1=1.
691 IFLAG=1
692 ENDIF
693 QNWFRZ=QNEWLQ
694 QNEWIC=QNEWIC+QNEWLQ*R1*0.5
695 QNEWLQ=QNEWLQ-QNEWLQ*R1*0.5
696 EFFQ=(TTFRZ-TBFRZ)/(TTEMP-TBFRZ)
697 TTEMP=TU(NK1)
698 ENDIF
699 !
700 ! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
701 !
702 IF(NK.EQ.K)THEN
703 BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
704 BOTERM=2.*(Z0(NK1)-ZLCL)*G*BE/1.5
705 ENTERM=0.
706 DZZ=Z0(NK1)-ZLCL
707 ELSE
708 BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
709 BOTERM=2.*DZA(NK)*G*BE/1.5
710 ENTERM=2.*UER(NK)*WTW/UPOLD
711 DZZ=DZA(NK)
712 ENDIF
713 WSQ=WTW
714 CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,RATE, &
715 QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1), G)
716
717 !...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
718 ! IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
719 !
720 IF(WTW.LE.0.)GOTO 65
721 WABS=SQRT(ABS(WTW))
722 WU(NK1)=WTW/WABS
723 !
724 ! UPDATE THE ABE FOR UNDILUTE ASCENT...
725 !
726 THTES(NK1)=T0(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QES(NK1))) &
727 * &
728 EXP((3374.6525/T0(NK1)-2.5403)*QES(NK1)*(1.+0.81* &
729 QES(NK1)))
730 UDLBE=((2.*THTUDL)/(THTES(NK)+THTES(NK1))-1.)*DZZ
731 IF(UDLBE.GT.0.)ABE=ABE+UDLBE*G
732 !
733 ! DETERMINE THE EFFECTS OF CLOUD GLACIATION IF WITHIN THE SPECIFIED
734 ! TEMP INTERVAL...
735 !
736 IF(FRC1.GT.1.E-6)THEN
737 CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QLIQ(NK1), &
738 QICE(NK1),RATIO2(NK1),TTFRZ,TBFRZ,QNWFRZ,RL,FRC1,EFFQ, &
739 IFLAG,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE &
740 ,CICE,DICE)
741 ENDIF
742 !
743 ! CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMP.
744 ! WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED WITH RESPECT TO
745 ! SAME DEGREE OF GLACIATION FOR ALL ENTRAINING AIR...
746 !
747 CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),RATIO2(NK1), &
748 RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
749
750 !...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
751 !
752 REI=VMFLCL*DP(NK1)*0.03/RAD
753 TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
754 !
755 !...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, NO
756 ! ENTRAINMENT IS ALLOWED AT THIS LEVEL...
757 !
758 IF(TVQU(NK1).LE.TV0(NK1))THEN
759 UER(NK1)=0.0
760 UDR(NK1)=REI
761 EE2=0.
762 UD2=1.
763 EQFRC(NK1)=0.
764 GOTO 55
765 ENDIF
766 LET=NK1
767 TTMP=TVQU(NK1)
768 !
769 !...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL
770 ! AIR FOR ESTIMATION OF ENTRAINMENT AND DETRAINMENT RATES...
771 !
772 F1=0.95
773 F2=1.-F1
774 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
775 QTMP=F1*Q0(NK1)+F2*QU(NK1)
776 TMPLIQ=F2*QLIQ(NK1)
777 TMPICE=F2*QICE(NK1)
778 CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
779 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
780 DLIQ,AICE,BICE,CICE,DICE)
781 TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
782 IF(TU95.GT.TV0(NK1))THEN
783 EE2=1.
784 UD2=0.
785 EQFRC(NK1)=1.0
786 GOTO 50
787 ENDIF
788 F1=0.10
789 F2=1.-F1
790 THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
791 QTMP=F1*Q0(NK1)+F2*QU(NK1)
792 TMPLIQ=F2*QLIQ(NK1)
793 TMPICE=F2*QICE(NK1)
794 CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
795 QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
796 DLIQ,AICE,BICE,CICE,DICE)
797 TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
798 IF(TU10.EQ.TVQU(NK1))THEN
799 EE2=1.
800 UD2=0.
801 EQFRC(NK1)=1.0
802 GOTO 50
803 ENDIF
804 EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
805 EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
806 EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
807 IF(EQFRC(NK1).EQ.1)THEN
808 EE2=1.
809 UD2=0.
810 GOTO 50
811 ELSEIF(EQFRC(NK1).EQ.0.)THEN
812 EE2=0.
813 UD2=1.
814 GOTO 50
815 ELSE
816 !
817 !...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
818 ! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
819 !
820 CALL PROF5(EQFRC(NK1),EE2,UD2)
821 ENDIF
822 !
823 50 IF(NK.EQ.K)THEN
824 EE1=1.
825 UD1=0.
826 ENDIF
827 !
828 !...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE
829 ! FRACTIONAL VALUES IN THE LAYER...
830 !
831 UER(NK1)=0.5*REI*(EE1+EE2)
832 UDR(NK1)=0.5*REI*(UD1+UD2)
833 !
834 !...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
835 ! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATION
836 !
837 55 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
838 !
839 !...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL
840 ! UPDRAFT FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE
841 ! PREVIOUS MODEL
842 !
843 IF(UDLBE.GT.0.)ABE=ABE-UDLBE*G
844 LET=NK
845 ! WRITE(98,1015)P0(NK1)/100.
846 GOTO 65
847 ENDIF
848 EE1=EE2
849 UD1=UD2
850 UPOLD=UMF(NK)-UDR(NK1)
851 UPNEW=UPOLD+UER(NK1)
852 UMF(NK1)=UPNEW
853 !
854 !...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND ICE IN
855 ! THE DETRAINING UPDRAFT MASS...
856 !
857 DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
858 DETIC(NK1)=QICE(NK1)*UDR(NK1)
859 QDT(NK1)=QU(NK1)
860 QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
861 THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
862 QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
863 QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
864 !
865 !...KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE IS
866 ! GENERATING PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF LIQUID
867 ! PRECIP AT A GIVING MODEL LVL, PPTICE THE SAME FOR ICE, TRPPT IS
868 ! THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE CURRENT MODEL LEVEL
869 !
870 IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
871 PPTLIQ(NK1)=QLQOUT(NK1)*(UMF(NK)-UDR(NK1))
872 PPTICE(NK1)=QICOUT(NK1)*(UMF(NK)-UDR(NK1))
873 TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
874 IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
875 60 CONTINUE
876 !
877 !...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
878 ! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
879 ! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE
880 ! BETWEEN THE LET AND CLOUD TOP...
881 !
882 !...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL
883 ! VELOCITY FIRST BECOMES NEGATIVE...
884 !
885 65 LTOP=NK
886 CLDHGT=Z0(LTOP)-ZLCL
887 !
888 !...IF CLOUD TOP HGT IS LESS THAN SPECIFIED MINIMUM HEIGHT, GO BACK AND
889 ! THE NEXT HIGHEST 60MB LAYER TO SEE IF A BIGGER CLOUD CAN BE OBTAINED
890 ! THAT SOURCE AIR...
891 !
892 ! IF(CLDHGT.LT.4.E3.OR.ABE.LT.1.)THEN
893 IF(CLDHGT.LT.3.E3.OR.ABE.LT.1.)THEN
894 DO 70 NK=K,LTOP
895 UMF(NK)=0.
896 UDR(NK)=0.
897 UER(NK)=0.
898 DETLQ(NK)=0.
899 DETIC(NK)=0.
900 PPTLIQ(NK)=0.
901 70 PPTICE(NK)=0.
902 GOTO 25
903 ENDIF
904 !
905 !...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS
906 ! FLUX THIS LEVEL...
907 !
908 IF(LET.EQ.LTOP)THEN
909 UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
910 DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
911 DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
912 TRPPT=TRPPT-(PPTLIQ(LTOP)+PPTICE(LTOP))
913 UER(LTOP)=0.
914 UMF(LTOP)=0.
915 GOTO 85
916 ENDIF
917 !
918 ! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
919 !
920 DPTT=0.
921 DO 71 NJ=LET+1,LTOP
922 71 DPTT=DPTT+DP(NJ)
923 DUMFDP=UMF(LET)/DPTT
924 !
925 !...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
926 ! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
927 ! PTOP
928 !
929 DO 75 NK=LET+1,LTOP
930 UDR(NK)=DP(NK)*DUMFDP
931 UMF(NK)=UMF(NK-1)-UDR(NK)
932 DETLQ(NK)=QLIQ(NK)*UDR(NK)
933 DETIC(NK)=QICE(NK)*UDR(NK)
934 TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
935 PPTLIQ(NK)=(UMF(NK-1)-UDR(NK))*QLQOUT(NK)
936 PPTICE(NK)=(UMF(NK-1)-UDR(NK))*QICOUT(NK)
937 TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
938 75 CONTINUE
939 !
940 !...SEND UPDRAFT CHARACTERISTICS TO OUTPUT FILES...
941 !
942 85 CONTINUE
943 !
944 !...EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR
945 ! THE UPDRAFT AIR...ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL...
946 !
947 DO 90 NK=1,K
948 IF(NK.GE.LC)THEN
949 IF(NK.EQ.LC)THEN
950 UMF(NK)=VMFLCL*DP(NK)/DPTHMX
951 UER(NK)=VMFLCL*DP(NK)/DPTHMX
952 ELSEIF(NK.LE.KPBL)THEN
953 UER(NK)=VMFLCL*DP(NK)/DPTHMX
954 UMF(NK)=UMF(NK-1)+UER(NK)
955 ELSE
956 UMF(NK)=VMFLCL
957 UER(NK)=0.
958 ENDIF
959 TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
960 QU(NK)=QMIX
961 WU(NK)=WLCL
962 ELSE
963 TU(NK)=0.
964 QU(NK)=0.
965 UMF(NK)=0.
966 WU(NK)=0.
967 UER(NK)=0.
968 ENDIF
969 UDR(NK)=0.
970 QDT(NK)=0.
971 QLIQ(NK)=0.
972 QICE(NK)=0.
973 QLQOUT(NK)=0.
974 QICOUT(NK)=0.
975 PPTLIQ(NK)=0.
976 PPTICE(NK)=0.
977 DETLQ(NK)=0.
978 DETIC(NK)=0.
979 RATIO2(NK)=0.
980 EE=Q0(NK)*P0(NK)/(EP2+Q0(NK))
981 TLOG=ALOG(EE/ALIQ)
982 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
983 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T0(NK)-T00))*( &
984 T0(NK)-TDPT)
985 THTA=T0(NK)*(1.E5/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
986 THETEE(NK)=THTA* &
987 EXP((3374.6525/TSAT-2.5403)*Q0(NK)*(1.+0.81*Q0(NK)) &
988 )
989 THTES(NK)=THTA* &
990 EXP((3374.6525/T0(NK)-2.5403)*QES(NK)*(1.+0.81* &
991 QES(NK)))
992 EQFRC(NK)=1.0
993 90 CONTINUE
994 !
995 LTOP1=LTOP+1
996 LTOPM1=LTOP-1
997 !
998 !...DEFINE VARIABLES ABOVE CLOUD TOP...
999 !
1000 DO 95 NK=LTOP1,KX
1001 UMF(NK)=0.
1002 UDR(NK)=0.
1003 UER(NK)=0.
1004 QDT(NK)=0.
1005 QLIQ(NK)=0.
1006 QICE(NK)=0.
1007 QLQOUT(NK)=0.
1008 QICOUT(NK)=0.
1009 DETLQ(NK)=0.
1010 DETIC(NK)=0.
1011 PPTLIQ(NK)=0.
1012 PPTICE(NK)=0.
1013 IF(NK.GT.LTOP1)THEN
1014 TU(NK)=0.
1015 QU(NK)=0.
1016 WU(NK)=0.
1017 ENDIF
1018 THTA0(NK)=0.
1019 THTAU(NK)=0.
1020 EMS(NK)=DP(NK)*DXSQ/G
1021 EMSD(NK)=1./EMS(NK)
1022 TG(NK)=T0(NK)
1023 QG(NK)=Q0(NK)
1024 QLG(NK)=0.
1025 QIG(NK)=0.
1026 QRG(NK)=0.
1027 QSG(NK)=0.
1028 95 OMG(NK)=0.
1029 OMG(KL+1)=0.
1030 P150=P0(KLCL)-1.50E4
1031 DO 100 NK=1,LTOP
1032 THTAD(NK)=0.
1033 EMS(NK)=DP(NK)*DXSQ/G
1034 EMSD(NK)=1./EMS(NK)
1035 !
1036 !...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION
1037 ! SCHEME
1038 !
1039 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
1040 THTAU(NK)=TU(NK)*EXN(NK)
1041 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
1042 THTA0(NK)=T0(NK)*EXN(NK)
1043 !
1044 !...LVF IS THE LEVEL AT WHICH MOISTURE FLUX IS ESTIMATED AS THE BASIS
1045 !...FOR PRECIPITATION EFFICIENCY CALCULATIONS...
1046 !
1047 IF(P0(NK).GT.P150)LVF=NK
1048 100 OMG(NK)=0.
1049 LVF=MIN0(LVF,LET)
1050 USR=UMF(LVF+1)*(QU(LVF+1)+QLIQ(LVF+1)+QICE(LVF+1))
1051 USR=AMIN1(USR,TRPPT)
1052 IF(USR.LT.1.E-8)USR=TRPPT
1053 !
1054 ! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1055 ! * TMIX-T00,PMIX,QMIX,ABE
1056 ! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1057 ! * WLCL,CLDHGT
1058 !
1059 !...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1060 !...AND MIDTROPOSPHERE IS USED.
1061 !
1062 WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
1063 WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
1064 WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
1065 VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1066 if (VCONV .gt. 0.) then
1067 TIMEC=DX/VCONV
1068 else
1069 TIMEC=3600.
1070 endif
1071 ! TIMEC=DX/VCONV
1072 TADVEC=TIMEC
1073 TIMEC=AMAX1(1800.,TIMEC)
1074 TIMEC=AMIN1(3600.,TIMEC)
1075 NIC=NINT(TIMEC/DT)
1076 TIMEC=FLOAT(NIC)*DT
1077 !
1078 !...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1079 !
1080 ! SHSIGN = CVMGT(1.,-1.,WSPD(LTOP).GT.WSPD(KLCL))
1081 IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1082 SHSIGN=1.
1083 ELSE
1084 SHSIGN=-1.
1085 ENDIF
1086 VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))* &
1087 (V0(LTOP)-V0(KLCL))
1088 VWS=1.E3*SHSIGN*SQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1089 PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
1090 PEF=AMAX1(PEF,.2)
1091 PEF=AMIN1(PEF,.9)
1092 !
1093 !...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1094 !
1095 CBH=(ZLCL-Z0(1))*3.281E-3
1096 IF(CBH.LT.3.)THEN
1097 RCBH=.02
1098 ELSE
1099 RCBH=.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*(- &
1100 1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
1101 ENDIF
1102 IF(CBH.GT.25)RCBH=2.4
1103 PEFCBH=1./(1.+RCBH)
1104 PEFCBH=AMIN1(PEFCBH,.9)
1105 !
1106 !... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1107 !
1108 PEFF=.5*(PEF+PEFCBH)
1109 PEFF2=PEFF
1110 ! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1111 !
1112 !*****************************************************************
1113 ! *
1114 ! COMPUTE DOWNDRAFT PROPERTIES *
1115 ! *
1116 !*****************************************************************
1117 !
1118 !...LET DOWNDRAFT ORIGINATE AT THE LEVEL OF MINIMUM SATURATION EQUIVALEN
1119 !...POTENTIAL TEMPERATURE (SEQT) IN THE CLOUD LAYER, EXTEND DOWNWARD TO
1120 !...SURFACE, OR TO THE LAYER BELOW CLOUD BASE AT WHICH ENVIR SEQT IS LES
1121 !...THAN MIN SEQT IN THE CLOUD LAYER...LET DOWNDRAFT DETRAIN OVER A LAYE
1122 !...OF SPECIFIED PRESSURE-DEPTH (DPDD)...
1123 !
1124 TDER=0.
1125 KSTART=MAX0(KPBL,KLCL)
1126 THTMIN=THTES(KSTART+1)
1127 KMIN=KSTART+1
1128 DO 104 NK=KSTART+2,LTOP-1
1129 THTMIN=AMIN1(THTMIN,THTES(NK))
1130 IF(THTMIN.EQ.THTES(NK))KMIN=NK
1131 104 CONTINUE
1132 LFS=KMIN
1133 IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT(P0(LFS),T0(LFS),Q0(LFS), &
1134 THETEE(LFS),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1135 EQFRC(LFS)=(THTES(LFS)-THETEU(LFS))/(THETEE(LFS)-THETEU(LFS))
1136 EQFRC(LFS)=AMAX1(EQFRC(LFS),0.)
1137 EQFRC(LFS)=AMIN1(EQFRC(LFS),1.)
1138 THETED(LFS)=THTES(LFS)
1139 !
1140 !...ESTIMATE THE EFFECT OF MELTING PRECIPITATION IN THE DOWNDRAFT...
1141 !
1142 IF(ML.GT.0)THEN
1143 DTMLTD=0.5*(QU(KLCL)-QU(LTOP))*RLF/CP
1144 ELSE
1145 DTMLTD=0.
1146 ENDIF
1147 TZ(LFS)=T0(LFS)-DTMLTD
1148 ES=ALIQ*EXP((TZ(LFS)*BLIQ-CLIQ)/(TZ(LFS)-DLIQ))
1149 QS=EP2*ES/(P0(LFS)-ES)
1150 QD(LFS)=EQFRC(LFS)*Q0(LFS)+(1.-EQFRC(LFS))*QU(LFS)
1151 THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QD(LFS)))
1152 IF(QD(LFS).GE.QS)THEN
1153 THETED(LFS)=THTAD(LFS)* &
1154 EXP((3374.6525/TZ(LFS)-2.5403)*QS*(1.+0.81*QS))
1155 ELSE
1156 CALL ENVIRTHT(P0(LFS),TZ(LFS),QD(LFS),THETED(LFS),0.,RL,EP2,ALIQ, &
1157 BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1158 ENDIF
1159 DO 107 NK=1,LFS
1160 ND=LFS-NK
1161 IF(THETED(LFS).GT.THTES(ND).OR.ND.EQ.1)THEN
1162 LDB=ND
1163 !
1164 !...IF DOWNDRAFT NEVER BECOMES NEGATIVELY BUOYANT OR IF IT
1165 !...IS SHALLOWER 50 mb, DON'T ALLOW IT TO OCCUR AT ALL...
1166 !
1167 IF(NK.EQ.1.OR.(P0(LDB)-P0(LFS)).LT.50.E2)GOTO 141
1168 ! testing ---- no downdraft
1169 ! GOTO 141
1170 GOTO 110
1171 ENDIF
1172 107 CONTINUE
1173 !
1174 !...ALLOW DOWNDRAFT TO DETRAIN IN A SINGLE LAYER, BUT WITH DOWNDRAFT AIR
1175 !...TYPICALLY FLUSHED UP INTO HIGHER LAYERS AS ALLOWED IN THE TOTAL
1176 !...VERTICAL ADVECTION CALCULATIONS FARTHER DOWN IN THE CODE...
1177 !
1178 110 DPDD=DP(LDB)
1179 LDT=LDB
1180 FRC=1.
1181 DPT=0.
1182 ! DO 115 NK=LDB,LFS
1183 ! DPT=DPT+DP(NK)
1184 ! IF(DPT.GT.DPDD)THEN
1185 ! LDT=NK
1186 ! FRC=(DPDD+DP(NK)-DPT)/DP(NK)
1187 ! GOTO 120
1188 ! ENDIF
1189 ! IF(NK.EQ.LFS-1)THEN
1190 ! LDT=NK
1191 ! FRC=1.
1192 ! DPDD=DPT
1193 ! GOTO 120
1194 ! ENDIF
1195 !115 CONTINUE
1196 120 CONTINUE
1197 !
1198 !...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX..
1199 !
1200 TVD(LFS)=T0(LFS)*(1.+0.608*QES(LFS))
1201 RDD=P0(LFS)/(R*TVD(LFS))
1202 A1=(1.-PEFF)*AU0
1203 DMF(LFS)=-A1*RDD
1204 DER(LFS)=EQFRC(LFS)*DMF(LFS)
1205 DDR(LFS)=0.
1206 DO 140 ND=LFS-1,LDB,-1
1207 ND1=ND+1
1208 IF(ND.LE.LDT)THEN
1209 DER(ND)=0.
1210 DDR(ND)=-DMF(LDT+1)*DP(ND)*FRC/DPDD
1211 DMF(ND)=DMF(ND1)+DDR(ND)
1212 FRC=1.
1213 THETED(ND)=THETED(ND1)
1214 QD(ND)=QD(ND1)
1215 ELSE
1216 DER(ND)=DMF(LFS)*0.03*DP(ND)/RAD
1217 DDR(ND)=0.
1218 DMF(ND)=DMF(ND1)+DER(ND)
1219 IF(RATIO2(ND).GT.0.)CALL ENVIRTHT(P0(ND),T0(ND),Q0(ND), &
1220 THETEE(ND),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1221 THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
1222 QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)
1223 ENDIF
1224 140 CONTINUE
1225 TDER=0.
1226 !
1227 !...CALCULATION AN EVAPORATION RATE FOR GIVEN MASS FLUX...
1228 !
1229 DO 135 ND=LDB,LDT
1230 TZ(ND)= &
1231 TPDD(P0(ND),THETED(LDT),T0(ND),QS,QD(ND),1.0,XLV0,XLV1, &
1232 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1233 ES=ALIQ*EXP((TZ(ND)*BLIQ-CLIQ)/(TZ(ND)-DLIQ))
1234 QS=EP2*ES/(P0(ND)-ES)
1235 DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
1236 RL=XLV0-XLV1*TZ(ND)
1237 DTMP=RL*QS*(1.-RHBC)/(CP+RL*RHBC*QS*DSSDT)
1238 T1RH=TZ(ND)+DTMP
1239 ES=RHBC*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1240 QSRH=EP2*ES/(P0(ND)-ES)
1241 !
1242 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1243 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1244 !
1245 IF(QSRH.LT.QD(ND))THEN
1246 QSRH=QD(ND)
1247 ! T1RH=T1+(QS-QSRH)*RL/CP
1248 T1RH=TZ(ND)
1249 ENDIF
1250 TZ(ND)=T1RH
1251 QS=QSRH
1252 TDER=TDER+(QS-QD(ND))*DDR(ND)
1253 QD(ND)=QS
1254 135 THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
1255 !
1256 !...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1257 !...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1258 !
1259 141 IF(TDER.LT.1.)THEN
1260 ! WRITE(98,3004)I,J
1261 3004 FORMAT(' ','I=',I3,2X,'J=',I3)
1262 PPTFLX=TRPPT
1263 CPR=TRPPT
1264 TDER=0.
1265 CNDTNF=0.
1266 UPDINC=1.
1267 LDB=LFS
1268 DO 117 NDK=1,LTOP
1269 DMF(NDK)=0.
1270 DER(NDK)=0.
1271 DDR(NDK)=0.
1272 THTAD(NDK)=0.
1273 WD(NDK)=0.
1274 TZ(NDK)=0.
1275 117 QD(NDK)=0.
1276 AINCM2=100.
1277 GOTO 165
1278 ENDIF
1279 !
1280 !...ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT IS
1281 !...CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP...
1282 !
1283 DEVDMF=TDER/DMF(LFS)
1284 PPR=0.
1285 PPTFLX=PEFF*USR
1286 RCED=TRPPT-PPTFLX
1287 !
1288 !...PPR IS THE TOTAL AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE
1289 !...UPDRAFT FROM CLOUD BASE TO THE LFS...UPDRAFT MASS FLUX WILL BE
1290 !...INCREASED UP TO THE LFS TO ACCOUNT FOR UPDRAFT AIR MIXING WITH
1291 !...ENVIRONMENTAL AIR TO THE UPDRAFT, SO PPR WILL INCREASE
1292 !...PROPORTIONATELY...
1293 !
1294 DO 132 NM=KLCL,LFS
1295 132 PPR=PPR+PPTLIQ(NM)+PPTICE(NM)
1296 IF(LFS.GE.KLCL)THEN
1297 DPPTDF=(1.-PEFF)*PPR*(1.-EQFRC(LFS))/UMF(LFS)
1298 ELSE
1299 DPPTDF=0.
1300 ENDIF
1301 !
1302 !...CNDTNF IS THE AMOUNT OF CONDENSATE TRANSFERRED ALONG WITH UPDRAFT
1303 !...MASS THE DOWNDRAFT AT THE LFS...
1304 !
1305 CNDTNF=(QLIQ(LFS)+QICE(LFS))*(1.-EQFRC(LFS))
1306 DMFLFS=RCED/(DEVDMF+DPPTDF+CNDTNF)
1307 IF(DMFLFS.GT.0.)THEN
1308 TDER=0.
1309 GOTO 141
1310 ENDIF
1311 !
1312 !...DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT
1313 !...MASS FLUX TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP, UPDINC IS T
1314 !...WHICH TO INCREASE THE UPDRAFT MASS FLUX BELOW THE LFS TO ACCOUNT FOR
1315 !...TRANSFER OF MASS FROM UPDRAFT TO DOWNDRAFT...
1316 !
1317 ! DDINC=DMFLFS/DMF(LFS)
1318 IF(LFS.GE.KLCL)THEN
1319 UPDINC=(UMF(LFS)-(1.-EQFRC(LFS))*DMFLFS)/UMF(LFS)
1320 !
1321 !...LIMIT UPDINC TO LESS THAN OR EQUAL TO 1.5...
1322 !
1323 IF(UPDINC.GT.1.5)THEN
1324 UPDINC=1.5
1325 DMFLFS2=UMF(LFS)*(UPDINC-1.)/(EQFRC(LFS)-1.)
1326 RCED2=DMFLFS2*(DEVDMF+DPPTDF+CNDTNF)
1327 PPTFLX=PPTFLX+(RCED-RCED2)
1328 PEFF2=PPTFLX/USR
1329 RCED=RCED2
1330 DMFLFS=DMFLFS2
1331 ENDIF
1332 ELSE
1333 UPDINC=1.
1334 ENDIF
1335 DDINC=DMFLFS/DMF(LFS)
1336 DO 149 NK=LDB,LFS
1337 DMF(NK)=DMF(NK)*DDINC
1338 DER(NK)=DER(NK)*DDINC
1339 DDR(NK)=DDR(NK)*DDINC
1340 149 CONTINUE
1341 CPR=TRPPT+PPR*(UPDINC-1.)
1342 PPTFLX=PPTFLX+PEFF*PPR*(UPDINC-1.)
1343 PEFF=PEFF2
1344 TDER=TDER*DDINC
1345 !
1346 !...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1347 ! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1348 ! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1349 !
1350 DO 155 NK=LC,LFS
1351 UMF(NK)=UMF(NK)*UPDINC
1352 UDR(NK)=UDR(NK)*UPDINC
1353 UER(NK)=UER(NK)*UPDINC
1354 PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1355 PPTICE(NK)=PPTICE(NK)*UPDINC
1356 DETLQ(NK)=DETLQ(NK)*UPDINC
1357 155 DETIC(NK)=DETIC(NK)*UPDINC
1358 !
1359 !...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1360 !...DOWNDRAFT...
1361 !
1362 IF(LDB.GT.1)THEN
1363 DO 156 NK=1,LDB-1
1364 DMF(NK)=0.
1365 DER(NK)=0.
1366 DDR(NK)=0.
1367 WD(NK)=0.
1368 TZ(NK)=0.
1369 QD(NK)=0.
1370 THTAD(NK)=0.
1371 156 CONTINUE
1372 ENDIF
1373 DO 157 NK=LFS+1,KX
1374 DMF(NK)=0.
1375 DER(NK)=0.
1376 DDR(NK)=0.
1377 WD(NK)=0.
1378 TZ(NK)=0.
1379 QD(NK)=0.
1380 THTAD(NK)=0.
1381 157 CONTINUE
1382 DO 158 NK=LDT+1,LFS-1
1383 TZ(NK)=0.
1384 QD(NK)=0.
1385 158 CONTINUE
1386 !
1387 !
1388 !...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE
1389 ! INFLOW INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN
1390 ! IS AVAILABLE IN THAT LAYER INITIALLY...
1391 !
1392 165 AINCMX=1000.
1393 LMAX=MAX0(KLCL,LFS)
1394 DO 166 NK=LC,LMAX
1395 IF((UER(NK)-DER(NK)).GT.0.)AINCM1=EMS(NK)/((UER(NK)-DER(NK))* &
1396 TIMEC)
1397 AINCMX=AMIN1(AINCMX,AINCM1)
1398 166 CONTINUE
1399 AINC=1.
1400 IF(AINCMX.LT.AINC)AINC=AINCMX
1401 !
1402 !...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT...THEY
1403 !...WILL ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE
1404 !...STABILIZATION CLOSURE...
1405 !
1406 NCOUNT=0
1407 TDER2=TDER
1408 PPTFL2=PPTFLX
1409 DO 170 NK=1,LTOP
1410 DETLQ2(NK)=DETLQ(NK)
1411 DETIC2(NK)=DETIC(NK)
1412 UDR2(NK)=UDR(NK)
1413 UER2(NK)=UER(NK)
1414 DDR2(NK)=DDR(NK)
1415 DER2(NK)=DER(NK)
1416 UMF2(NK)=UMF(NK)
1417 DMF2(NK)=DMF(NK)
1418 170 CONTINUE
1419 FABE=1.
1420 STAB=0.95
1421 NOITR=0
1422 IF(AINC/AINCMX.GT.0.999)THEN
1423 NCOUNT=0
1424 GOTO 255
1425 ENDIF
1426 ISTOP=0
1427 175 NCOUNT=NCOUNT+1
1428 !
1429 !*****************************************************************
1430 ! *
1431 ! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1432 ! *
1433 !*****************************************************************
1434 !
1435 !...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1436 !...SATISFY MASS CONTINUITY...
1437 !
1438 185 CONTINUE
1439 DTT=TIMEC
1440 DO 200 NK=1,LTOP
1441 DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1442 IF(NK.GT.1)THEN
1443 OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1444 DTT1=0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
1445 DTT=AMIN1(DTT,DTT1)
1446 ENDIF
1447 200 CONTINUE
1448 DO 488 NK=1,LTOP
1449 THPA(NK)=THTA0(NK)
1450 QPA(NK)=Q0(NK)
1451 NSTEP=NINT(TIMEC/DTT+1)
1452 DTIME=TIMEC/FLOAT(NSTEP)
1453 FXM(NK)=OMG(NK)*DXSQ/G
1454 488 CONTINUE
1455 !
1456 !...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1457 !
1458 DO 495 NTC=1,NSTEP
1459 !
1460 !...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
1461 !...SIGN OF OMEGA...
1462 !
1463 DO 493 NK=1,LTOP
1464 THFXTOP(NK)=0.
1465 THFXBOT(NK)=0.
1466 QFXTOP(NK)=0.
1467 493 QFXBOT(NK)=0.
1468 DO 494 NK=2,LTOP
1469 IF(OMG(NK).LE.0.)THEN
1470 THFXBOT(NK)=-FXM(NK)*THPA(NK-1)
1471 QFXBOT(NK)=-FXM(NK)*QPA(NK-1)
1472 THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1473 QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1474 ELSE
1475 THFXBOT(NK)=-FXM(NK)*THPA(NK)
1476 QFXBOT(NK)=-FXM(NK)*QPA(NK)
1477 THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1478 QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1479 ENDIF
1480 494 CONTINUE
1481 !
1482 !...UPDATE THE THETA AND QV VALUES AT EACH LEVEL..
1483 !
1484 DO 492 NK=1,LTOP
1485 THPA(NK)=THPA(NK)+(THFXBOT(NK)+UDR(NK)*THTAU(NK)+DDR(NK)* &
1486 THTAD(NK)+THFXTOP(NK)-(UER(NK)-DER(NK))*THTA0(NK))* &
1487 DTIME*EMSD(NK)
1488 QPA(NK)=QPA(NK)+(QFXBOT(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)+ &
1489 QFXTOP(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
1490
1491 492 CONTINUE
1492 495 CONTINUE
1493 DO 498 NK=1,LTOP
1494 THTAG(NK)=THPA(NK)
1495 QG(NK)=QPA(NK)
1496 498 CONTINUE
1497 !
1498 !...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO,
1499 !...BORROW MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO.
1500 !
1501 DO 499 NK=1,LTOP
1502 IF(QG(NK).LT.0.)THEN
1503 IF(NK.EQ.1)THEN
1504 CALL wrf_error_fatal ( 'module_cu_kf.F: problem with kf scheme: qg = 0 at the surface' )
1505 ENDIF
1506 NK1=NK+1
1507 IF(NK.EQ.LTOP)NK1=KLCL
1508 TMA=QG(NK1)*EMS(NK1)
1509 TMB=QG(NK-1)*EMS(NK-1)
1510 TMM=(QG(NK)-1.E-9)*EMS(NK)
1511 BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1512 ACOEFF=BCOEFF*TMA/TMB
1513 TMB=TMB*(1.-BCOEFF)
1514 TMA=TMA*(1.-ACOEFF)
1515 IF(NK.EQ.LTOP)THEN
1516 QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
1517 IF(ABS(QVDIFF).GT.1.)THEN
1518 PRINT *,'--WARNING-- CLOUD BASE WATER VAPOR CHANGES BY ', &
1519 QVDIFF, &
1520 ' PERCENT WHEN MOISTURE IS BORROWED TO PREVENT NEG VALUES', &
1521 ' IN KAIN-FRITSCH'
1522 ENDIF
1523 ENDIF
1524 QG(NK)=1.E-9
1525 QG(NK1)=TMA*EMSD(NK1)
1526 QG(NK-1)=TMB*EMSD(NK-1)
1527 ENDIF
1528 499 CONTINUE
1529 TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
1530 IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1531 ! WRITE(98,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;'
1532 ! * ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1533 WRITE(6,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;' &
1534 ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1535 ISTOP=1
1536 GOTO 265
1537 ENDIF
1538 !
1539 !...CONVERT THETA TO T...
1540 !
1541 ! PAY ATTENTION ...
1542 !
1543 DO 230 NK=1,LTOP
1544 EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
1545 TG(NK)=THTAG(NK)/EXN(NK)
1546 TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
1547 230 CONTINUE
1548 !
1549 !*******************************************************************
1550 ! *
1551 ! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
1552 ! *
1553 !*******************************************************************
1554 !
1555 !...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1556 !
1557 THMIX=0.
1558 QMIX=0.
1559 PMIX=0.
1560 DO 217 NK=LC,KPBL
1561 ROCPQ=0.2854*(1.-0.28*QG(NK))
1562 THMIX=THMIX+DP(NK)*TG(NK)*(P00/P0(NK))**ROCPQ
1563 QMIX=QMIX+DP(NK)*QG(NK)
1564 217 PMIX=PMIX+DP(NK)*P0(NK)
1565 THMIX=THMIX/DPTHMX
1566 QMIX=QMIX/DPTHMX
1567 PMIX=PMIX/DPTHMX
1568 ROCPQ=0.2854*(1.-0.28*QMIX)
1569 TMIX=THMIX*(PMIX/P00)**ROCPQ
1570 ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
1571 QS=EP2*ES/(PMIX-ES)
1572 !
1573 !...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1574 !
1575 IF(QMIX.GT.QS)THEN
1576 RL=XLV0-XLV1*TMIX
1577 CPM=CP*(1.+0.887*QMIX)
1578 DSSDT=QS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1579 DQ=(QMIX-QS)/(1.+RL*DSSDT/CPM)
1580 TMIX=TMIX+RL/CP*DQ
1581 QMIX=QMIX-DQ
1582 ROCPQ=0.2854*(1.-0.28*QMIX)
1583 THMIX=TMIX*(P00/PMIX)**ROCPQ
1584 TLCL=TMIX
1585 PLCL=PMIX
1586 ELSE
1587 QMIX=AMAX1(QMIX,0.)
1588 EMIX=QMIX*PMIX/(EP2+QMIX)
1589 TLOG=ALOG(EMIX/ALIQ)
1590 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1591 TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX- &
1592 TDPT)
1593 TLCL=AMIN1(TLCL,TMIX)
1594 CPORQ=1./ROCPQ
1595 PLCL=P00*(TLCL/THMIX)**CPORQ
1596 ENDIF
1597 TVLCL=TLCL*(1.+0.608*QMIX)
1598 DO 235 NK=LC,KL
1599 KLCL=NK
1600 235 IF(PLCL.GE.P0(NK))GOTO 240
1601 240 K=KLCL-1
1602 DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
1603 !
1604 !...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1605 !
1606 TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
1607 QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
1608 TVEN=TENV*(1.+0.608*QENV)
1609 TVBAR=0.5*(TVG(K)+TVEN)
1610 ! ZLCL=Z0(K)+R*TVBAR*ALOG(P0(K)/PLCL)/G
1611 ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
1612 TVAVG=0.5*(TVEN+TG(KLCL)*(1.+0.608*QG(KLCL)))
1613 PLCL=P0(KLCL)*EXP(G/(R*TVAVG)*(Z0(KLCL)-ZLCL))
1614 THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))* &
1615 EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
1616 ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
1617 QESE=EP2*ES/(PLCL-ES)
1618 THTESG(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))* &
1619 EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
1620 !
1621 !...COMPUTE ADJUSTED ABE(ABEG).
1622 !
1623 ABEG=0.
1624 THTUDL=THETEU(K)
1625 DO 245 NK=K,LTOPM1
1626 NK1=NK+1
1627 ES=ALIQ*EXP((TG(NK1)*BLIQ-CLIQ)/(TG(NK1)-DLIQ))
1628 QESE=EP2*ES/(P0(NK1)-ES)
1629 THTESG(NK1)=TG(NK1)*(1.E5/P0(NK1))**(0.2854*(1.-0.28*QESE))* &
1630 EXP((3374.6525/TG(NK1)-2.5403)*QESE*(1.+0.81*QESE) &
1631 )
1632 ! DZZ=CVMGT(Z0(KLCL)-ZLCL,DZA(NK),NK.EQ.K)
1633 IF(NK.EQ.K)THEN
1634 DZZ=Z0(KLCL)-ZLCL
1635 ELSE
1636 DZZ=DZA(NK)
1637 ENDIF
1638 BE=((2.*THTUDL)/(THTESG(NK1)+THTESG(NK))-1.)*DZZ
1639 245 IF(BE.GT.0.)ABEG=ABEG+BE*G
1640 !
1641 !...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
1642 !...THE PERIOD TIMEC...
1643 !
1644 IF(NOITR.EQ.1)THEN
1645 ! WRITE(98,1060)FABE
1646 GOTO 265
1647 ENDIF
1648 DABE=AMAX1(ABE-ABEG,0.1*ABE)
1649 FABE=ABEG/(ABE+1.E-8)
1650 IF(FABE.GT.1.)THEN
1651 ! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS '
1652 ! *,'GRID POINT; NO CONVECTION ALLOWED!'
1653 GOTO 325
1654 ENDIF
1655 IF(NCOUNT.NE.1)THEN
1656 DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
1657 IF(DFDA.GT.0.)THEN
1658 NOITR=1
1659 AINC=AINCOLD
1660 GOTO 255
1661 ENDIF
1662 ENDIF
1663 AINCOLD=AINC
1664 FABEOLD=FABE
1665 IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
1666 ! WRITE(98,1055)FABE
1667 GOTO 265
1668 ENDIF
1669 IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB)GOTO 265
1670 IF(NCOUNT.GT.10)THEN
1671 ! WRITE(98,1060)FABE
1672 GOTO 265
1673 ENDIF
1674 !
1675 !...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE
1676 !...CONVECTIVE MASS FLUX BY THE FACTOR AINC:
1677 !
1678 IF(FABE.EQ.0.)THEN
1679 AINC=AINC*0.5
1680 ELSE
1681 AINC=AINC*STAB*ABE/(DABE+1.E-8)
1682 ENDIF
1683 255 AINC=AMIN1(AINCMX,AINC)
1684 !...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION
1685 !...WILL BE MINIMAL SO JUST IGNORE IT...
1686 IF(AINC.LT.0.05)GOTO 325
1687 ! AINC=AMAX1(AINC,0.05)
1688 TDER=TDER2*AINC
1689 PPTFLX=PPTFL2*AINC
1690 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABEOLD,AINCOLD
1691 DO 260 NK=1,LTOP
1692 UMF(NK)=UMF2(NK)*AINC
1693 DMF(NK)=DMF2(NK)*AINC
1694 DETLQ(NK)=DETLQ2(NK)*AINC
1695 DETIC(NK)=DETIC2(NK)*AINC
1696 UDR(NK)=UDR2(NK)*AINC
1697 UER(NK)=UER2(NK)*AINC
1698 DER(NK)=DER2(NK)*AINC
1699 DDR(NK)=DDR2(NK)*AINC
1700 260 CONTINUE
1701 !
1702 !...GO BACK UP FOR ANOTHER ITERATION...
1703 !
1704 GOTO 175
1705 265 CONTINUE
1706 !
1707 !...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
1708 !...GRID POINT...
1709 !
1710 !...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
1711 !
1712 !...FRC2 IS THE FRACTION OF TOTAL CONDENSATE
1713 !...GENERATED THAT GOES INTO PRECIPITIATION
1714 FRC2=PPTFLX/(CPR*AINC)
1715 DO 270 NK=1,LTOP
1716 QLPA(NK)=QL0(NK)
1717 QIPA(NK)=QI0(NK)
1718 QRPA(NK)=QR0(NK)
1719 QSPA(NK)=QS0(NK)
1720 RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2
1721 SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2
1722 270 CONTINUE
1723 DO 290 NTC=1,NSTEP
1724 !
1725 !...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH
1726 !...LAYER BASED ON THE SIGN OF OMEGA...
1727 !
1728 DO 275 NK=1,LTOP
1729 QLFXIN(NK)=0.
1730 QLFXOUT(NK)=0.
1731 QIFXIN(NK)=0.
1732 QIFXOUT(NK)=0.
1733 QRFXIN(NK)=0.
1734 QRFXOUT(NK)=0.
1735 QSFXIN(NK)=0.
1736 QSFXOUT(NK)=0.
1737 275 CONTINUE
1738 DO 280 NK=2,LTOP
1739 IF(OMG(NK).LE.0.)THEN
1740 QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
1741 QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
1742 QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
1743 QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
1744 QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
1745 QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
1746 QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
1747 QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
1748 ELSE
1749 QLFXOUT(NK)=FXM(NK)*QLPA(NK)
1750 QIFXOUT(NK)=FXM(NK)*QIPA(NK)
1751 QRFXOUT(NK)=FXM(NK)*QRPA(NK)
1752 QSFXOUT(NK)=FXM(NK)*QSPA(NK)
1753 QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
1754 QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
1755 QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
1756 QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
1757 ENDIF
1758 280 CONTINUE
1759 !
1760 !...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
1761 !
1762 DO 285 NK=1,LTOP
1763 QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME* &
1764 EMSD(NK)
1765 QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME* &
1766 EMSD(NK)
1767 QRPA(NK)=QRPA(NK)+(QRFXIN(NK)+QLQOUT(NK)*UDR(NK)-QRFXOUT(NK) &
1768 +RAINFB(NK))*DTIME*EMSD(NK)
1769 QSPA(NK)=QSPA(NK)+(QSFXIN(NK)+QICOUT(NK)*UDR(NK)-QSFXOUT(NK) &
1770 +SNOWFB(NK))*DTIME*EMSD(NK)
1771 285 CONTINUE
1772 290 CONTINUE
1773 DO 295 NK=1,LTOP
1774 QLG(NK)=QLPA(NK)
1775 QIG(NK)=QIPA(NK)
1776 QRG(NK)=QRPA(NK)
1777 QSG(NK)=QSPA(NK)
1778 295 CONTINUE
1779 ! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
1780 !
1781 !...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
1782 !
1783 IF(ISTOP.EQ.1)THEN
1784 WRITE(6,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
1785 ' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
1786 ,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
1787 DO 300 K=LTOP,1,-1
1788 DTT=(TG(K)-T0(K))*86400./TIMEC
1789 RL=XLV0-XLV1*TG(K)
1790 DR=-(QG(K)-Q0(K))*RL*86400./(TIMEC*CP)
1791 UDFRC=UDR(K)*TIMEC*EMSD(K)
1792 UEFRC=UER(K)*TIMEC*EMSD(K)
1793 DDFRC=DDR(K)*TIMEC*EMSD(K)
1794 DEFRC=-DER(K)*TIMEC*EMSD(K)
1795 WRITE (6,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)* &
1796 1.E4,UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC &
1797 ,DDFRC,EMS(K)/1.E11,W0AVG1D(K)*1.E2,DETLQ(K) &
1798 *TIMEC*EMSD(K)*1.E3,DETIC(K)*TIMEC*EMSD(K)* &
1799 1.E3
1800 300 CONTINUE
1801 WRITE(6,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
1802 'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
1803 DO 305 K=KX,1,-1
1804 DTT=TG(K)-T0(K)
1805 TUC=TU(K)-T00
1806 IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
1807 TDC=TZ(K)-T00
1808 IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
1809 ES=ALIQ*EXP((BLIQ*TG(K)-CLIQ)/(TG(K)-DLIQ))
1810 QGS=ES*EP2/(P0(K)-ES)
1811 RH0=Q0(K)/QES(K)
1812 RHG=QG(K)/QGS
1813 WRITE (6,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC &
1814 ,TDC,Q0(K)*1000.,QG(K)*1000.,(QG(K)-Q0(K))* &
1815 1000.,QU(K)*1000.,QD(K)*1000.,QLG(K)*1000., &
1816 QIG(K)*1000.,QRG(K)*1000.,QSG(K)*1000.,RH0,RHG
1817 305 CONTINUE
1818 !
1819 !...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
1820 !...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
1821 !
1822 IF(ISTOP.EQ.1)THEN
1823 DO 310 K=1,KX
1824 WRITE ( wrf_err_message , 1115 ) &
1825 Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
1826 U0(K),V0(K),DP(K)/100.,W0AVG1D(K)
1827 CALL wrf_message ( TRIM( wrf_err_message ) )
1828 310 CONTINUE
1829 CALL wrf_error_fatal ( 'module_cu_kf.F: KAIN-FRITSCH' )
1830 ENDIF
1831 ENDIF
1832 CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
1833 ! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
1834 !
1835 ! EVALUATE MOISTURE BUDGET...
1836 !
1837 QINIT=0.
1838 QFNL=0.
1839 DPT=0.
1840 DO 315 NK=1,LTOP
1841 DPT=DPT+DP(NK)
1842 QINIT=QINIT+Q0(NK)*EMS(NK)
1843 QFNL=QFNL+QG(NK)*EMS(NK)
1844 QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
1845 315 CONTINUE
1846 QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC)
1847 ERR2=(QFNL-QINIT)*100./QINIT
1848 ! WRITE(98,1110)QINIT,QFNL,ERR2
1849 ! IF(ABS(ERR2).GT.0.05)STOP 'QVERR'
1850 IF(ABS(ERR2).GT.0.05)CALL wrf_error_fatal( 'module_cu_kf.F: QVERR' )
1851 RELERR=ERR2*QINIT/(PPTFLX*TIMEC+1.E-10)
1852 ! WRITE(98,1120)RELERR
1853 ! WRITE(98,*)'TDER, CPR, USR, TRPPT =',
1854 ! *TDER,CPR*AINC,USR*AINC,TRPPT*AINC
1855 !
1856 !...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
1857 !
1858 !...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
1859 !...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
1860 !
1861 IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
1862 NCA(I,J)=FLOAT(NIC)
1863 DO 320 K=1,KX
1864 ! IF(IMOIST.NE.2)THEN
1865 !
1866 !...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
1867 !...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
1868 !...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
1869 !...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
1870 !...OF QG...
1871 !
1872 ! RLC=XLV0-XLV1*TG(K)
1873 ! RLS=XLS0-XLS1*TG(K)
1874 ! CPM=CP*(1.+0.887*QG(K))
1875 ! TG(K)=TG(K)-(RLC*(QLG(K)+QRG(K))+RLS*(QIG(K)+QSG(K)))/CPM
1876 ! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
1877 ! DQCDT(K)=0.
1878 ! DQIDT(K)=0.
1879 ! DQRDT(K)=0.
1880 ! DQSDT(K)=0.
1881 ! ELSE
1882 IF(.NOT. qi_flag .and. warm_rain)THEN
1883 !
1884 !...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
1885 !
1886 CPM=CP*(1.+0.887*QG(K))
1887 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
1888 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
1889 DQIDT(K)=0.
1890 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
1891 DQSDT(K)=0.
1892 ELSEIF(.NOT. qi_flag .and. .not. warm_rain)THEN
1893 !
1894 !...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
1895 !...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
1896 !
1897 CPM=CP*(1.+0.887*QG(K))
1898 IF(K.LE.ML)THEN
1899 TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
1900 ELSEIF(K.GT.ML)THEN
1901 TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
1902 ENDIF
1903 DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
1904 DQIDT(K)=0.
1905 DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
1906 DQSDT(K)=0.
1907 ELSEIF(qi_flag) THEN
1908 !
1909 !...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE
1910 !...TENDENCY OF HYDROMETEORS DIRECTLY...
1911 !
1912 DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
1913 DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
1914 DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
1915 IF (qs_flag ) THEN
1916 DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
1917 ELSE
1918 DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
1919 ENDIF
1920 ELSE
1921 CALL wrf_error_fatal ( 'module_cu_kf: THIS COMBINATION OF IMOIST, IICE NOT ALLOWED' )
1922 ENDIF
1923 ! ENDIF
1924 DTDT(K)=(TG(K)-T0(K))/TIMEC
1925 DQDT(K)=(QG(K)-Q0(K))/TIMEC
1926 320 CONTINUE
1927
1928 ! RAINCV is in the unit of mm
1929
1930 RAINCV(I,J)=DT*PPTFLX*(1.-FBFRC)/DXSQ
1931 ! RNC=0.1*TIMEC*PPTFLX/DXSQ
1932 RNC=RAINCV(I,J)*NIC
1933 ! WRITE(98,909)RNC
1934 909 FORMAT(' CONVECTIVE RAINFALL =',F8.4,' CM')
1935
1936 325 CONTINUE
1937
1938 1000 FORMAT(' ',10A8)
1939 1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
1940 1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
1941 1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
1942 1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
1943 ' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
1944 I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
1945 ' CAPE=',0PF7.1)
1946 1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
1947 E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
1948 F8.1)
1949 1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
1950 ,F6.3,'VWS=',F5.2)
1951 1040 FORMAT(' ','PRECIP EFF = 100%, ENVIR CANNOT SUPPORT DOWND' &
1952 ,'RAFTS')
1953 !1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10...PPTFLX' &
1954 ! ' IS DIFFERENT FROM THAT GIVEN BY PRECIP EFF RELATION')
1955 ! FLIC HAS TROUBLE WITH THIS ONE.
1956 1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10')
1957 1050 FORMAT(' ','LCOUNT= ',I3,' PPTFLX/CPR, PEFF= ',F5.3,1X,F5.3, &
1958 'DMF(LFS)/UMF(LCL)= ',F5.3)
1959 1055 FORMAT(/'*** DEGREE OF STABILIZATION =',F5.3,', NO MORE MASS F' &
1960 ,'LUX IS ALLOWED')
1961 !1060 FORMAT(/' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED ' &
1962 ! 'DEGREE OF STABILIZATION! FABE= ',F6.4)
1963 1060 FORMAT(/' ITERATION DOES NOT CONVERGE. FABE= ',F6.4)
1964 1070 FORMAT (16A8)
1965 1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
1966 1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, NSTEP=',F5.0,I3, &
1967 'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
1968 1085 FORMAT (A3,16A7,2A8)
1969 1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
1970 1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ', &
1971 F10.0)
1972 1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =', &
1973 E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'PERCENT')
1974 1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
1975 ' TOTAL WATER CHANGE =',F8.2,'PERCENT')
1976 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4 &
1977 )
1978 1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3, &
1979 'PERCENT')
1980
1981 END SUBROUTINE KFPARA
1982
1983 !-----------------------------------------------------------------------
1984 SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
1985 QNEWIC,QLQOUT,QICOUT,G)
1986 !-----------------------------------------------------------------------
1987 IMPLICIT NONE
1988 !-----------------------------------------------------------------------
1989 ! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
1990 ! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
1991 ! CULATED USING THE EQUATION DQ=-RATE*Q*DT, BUT TO SIMULATE A QUASI-
1992 ! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
1993 ! RESOLUTION THIS IS EXPRESSED AS Q=Q*EXP(-RATE*DZ).
1994
1995 REAL, INTENT(IN ) :: G
1996 REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
1997 REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
1998
1999 REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2000
2001 QTOT=QLIQ+QICE
2002 QNEW=QNEWLQ+QNEWIC
2003 !
2004 ! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY C
2005 ! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2006 ! LEVELS...
2007 !
2008 QEST=0.5*(QTOT+QNEW)
2009 G1=WTW+BOTERM-ENTERM-2.*G*DZ*QEST/1.5
2010 IF(G1.LT.0.0)G1=0.
2011 WAVG=(SQRT(WTW)+SQRT(G1))/2.
2012 CONV=RATE*DZ/WAVG
2013 !
2014 ! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2015 ! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2016 ! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2017 ! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2018 !
2019 RATIO3=QNEWLQ/(QNEW+1.E-10)
2020 ! OLDQ=QTOT
2021 QTOT=QTOT+0.6*QNEW
2022 OLDQ=QTOT
2023 RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-10)
2024 QTOT=QTOT*EXP(-CONV)
2025 !
2026 ! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2027 ! PARCEL AT THIS LEVEL...
2028 !
2029 DQ=OLDQ-QTOT
2030 QLQOUT=RATIO4*DQ
2031 QICOUT=(1.-RATIO4)*DQ
2032 !
2033 ! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2034 ! LATE VERTICAL VELOCITY
2035 !
2036 PPTDRG=0.5*(OLDQ+QTOT-0.2*QNEW)
2037 WTW=WTW+BOTERM-ENTERM-2.*G*DZ*PPTDRG/1.5
2038 !
2039 ! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2040 ! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2041 !
2042 QLIQ=RATIO4*QTOT+RATIO3*0.4*QNEW
2043 QICE=(1.-RATIO4)*QTOT+(1.-RATIO3)*0.4*QNEW
2044 QNEWLQ=0.
2045 QNEWIC=0.
2046
2047 END SUBROUTINE CONDLOAD
2048
2049 !-----------------------------------------------------------------------
2050 SUBROUTINE DTFRZNEW(TU,P,THTEU,QVAP,QLIQ,QICE,RATIO2,TTFRZ,TBFRZ, &
2051 QNWFRZ,RL,FRC1,EFFQ,IFLAG,XLV0,XLV1,XLS0,XLS1, &
2052 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2053 !-----------------------------------------------------------------------
2054 IMPLICIT NONE
2055 !-----------------------------------------------------------------------
2056 REAL, INTENT(IN ) :: XLV0,XLV1
2057 REAL, INTENT(IN ) :: P,TTFRZ,TBFRZ,EFFQ,XLS0,XLS1,EP2,ALIQ, &
2058 BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE
2059 REAL, INTENT(INOUT) :: TU,THTEU,QVAP,QLIQ,QICE,RATIO2, &
2060 FRC1,RL,QNWFRZ
2061 INTEGER, INTENT(INOUT) :: IFLAG
2062
2063 REAL :: CCP,RV,C5,QLQFRZ,QNEW,ESLIQ,ESICE,RLC,RLS,PI,ES,RLF,A, &
2064 B,C,DQVAP,DTFRZ,TU1,QVAP1
2065 !-----------------------------------------------------------------------
2066 !
2067 !...ALLOW GLACIATION OF THE UPDRAFT TO OCCUR AS AN APPROXIMATELY LINEAR
2068 ! FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE TTFRZ TO TBFRZ...
2069 !
2070
2071 RV=461.5
2072 C5=1.0723E-3
2073 !
2074 !...ADJUST THE LIQUID WATER CONCENTRATIONS FROM FRESH CONDENSATE AND THA
2075 ! BROUGHT UP FROM LOWER LEVELS TO AN AMOUNT THAT WOULD BE PRESENT IF N
2076 ! LIQUID WATER HAD FROZEN THUS FAR...THIS IS NECESSARY BECAUSE THE
2077 ! EXPRESSION FOR TEMP CHANGE IS MULTIPLIED BY THE FRACTION EQUAL TO TH
2078 ! PARCEL TEMP DECREASE SINCE THE LAST MODEL LEVEL DIVIDED BY THE TOTAL
2079 ! GLACIATION INTERVAL, SO THAT EFFECTIVELY THIS APPROXIMATELY ALLOWS A
2080 ! AMOUNT OF LIQUID WATER TO FREEZE WHICH IS EQUAL TO THIS SAME FRACTIO
2081 ! OF THE LIQUID WATER THAT WAS PRESENT BEFORE THE GLACIATION PROCESS W
2082 ! INITIATED...ALSO, TO ALLOW THETAU TO CONVERT APPROXIMATELY LINEARLY
2083 ! ITS VALUE WITH RESPECT TO ICE, WE NEED TO ALLOW A PORTION OF THE FRE
2084 ! CONDENSATE TO CONTRIBUTE TO THE GLACIATION PROCESS; THE FRACTIONAL
2085 ! AMOUNT THAT APPLIES TO THIS PORTION IS 1/2 OF THE FRACTIONAL AMOUNT
2086 ! FROZEN OF THE "OLD" CONDENSATE BECAUSE THIS FRESH CONDENSATE IS ONLY
2087 ! PRODUCED GRADUALLY OVER THE LAYER...NOTE THAT IN TERMS OF THE DYNAMI
2088 ! OF THE PRECIPITATION PROCESS, IE. PRECIPITATION FALLOUT, THIS FRACTI
2089 ! AMNT OF FRESH CONDENSATE HAS ALREADY BEEN INCLUDED IN THE ICE CATEGO
2090 !
2091 QLQFRZ=QLIQ*EFFQ
2092 QNEW=QNWFRZ*EFFQ*0.5
2093 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2094 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2095 RLC=2.5E6-2369.276*(TU-273.16)
2096 RLS=2833922.-259.532*(TU-273.16)
2097 RLF=RLS-RLC
2098 CCP=1005.7*(1.+0.89*QVAP)
2099 !
2100 ! A = D(ES)/DT IS THAT CALCULATED FROM BUCK`S (1981) EMPIRICAL FORMULAS
2101 ! FOR SATURATION VAPOR PRESSURE...
2102 !
2103 A=(CICE-BICE*DICE)/((TU-DICE)*(TU-DICE))
2104 B=RLS*EP2/P
2105 C=A*B*ESICE/CCP
2106 DQVAP=B*(ESLIQ-ESICE)/(RLS+RLS*C)-RLF*(QLQFRZ+QNEW)/(RLS+RLS/C)
2107 DTFRZ=(RLF*(QLQFRZ+QNEW)+B*(ESLIQ-ESICE))/(CCP+A*B*ESICE)
2108 TU1=TU
2109 QVAP1=QVAP
2110 TU=TU+FRC1*DTFRZ
2111 QVAP=QVAP-FRC1*DQVAP
2112 ES=QVAP*P/(EP2+QVAP)
2113 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2114 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2115 RATIO2=(ESLIQ-ES)/(ESLIQ-ESICE)
2116 !
2117 ! TYPICALLY, RATIO2 IS VERY CLOSE TO (TTFRZ-TU)/(TTFRZ-TBFRZ), USUALLY
2118 ! WITHIN 1% (USING TU BEFORE GALCIATION EFFECTS ARE APPLIED); IF THE
2119 ! INITIAL UPDRAFT TEMP IS BELOW TBFRZ AND RATIO2 IS STILL LESS THAN 1,
2120 ! AN ADJUSTMENT TO FRC1 AND RATIO2 IS INTRODUCED SO THAT GLACIATION
2121 ! EFFECTS ARE NOT UNDERESTIMATED; CONVERSELY, IF RATIO2 IS GREATER THAN
2122 ! FRC1 IS ADJUSTED SO THAT GLACIATION EFFECTS ARE NOT OVERESTIMATED...
2123 !
2124 IF(IFLAG.GT.0.AND.RATIO2.LT.1)THEN
2125 FRC1=FRC1+(1.-RATIO2)
2126 TU=TU1+FRC1*DTFRZ
2127 QVAP=QVAP1-FRC1*DQVAP
2128 RATIO2=1.
2129 IFLAG=1
2130 GOTO 20
2131 ENDIF
2132 IF(RATIO2.GT.1.)THEN
2133 FRC1=FRC1-(RATIO2-1.)
2134 FRC1=AMAX1(0.0,FRC1)
2135 TU=TU1+FRC1*DTFRZ
2136 QVAP=QVAP1-FRC1*DQVAP
2137 RATIO2=1.
2138 IFLAG=1
2139 ENDIF
2140 !
2141 ! CALCULATE A HYBRID VALUE OF THETAU, ASSUMING THAT THE LATENT HEAT OF
2142 ! VAPORIZATION/SUBLIMATION CAN BE ESTIMATED USING THE SAME WEIGHTING
2143 ! FUNCTION AS THAT USED TO CALCULATE SATURATION VAPOR PRESSURE, CALCU-
2144 ! LATE NEW LIQUID WATER AND ICE CONCENTRATIONS...
2145 !
2146 20 RLC=XLV0-XLV1*TU
2147 RLS=XLS0-XLS1*TU
2148 RL=RATIO2*RLS+(1.-RATIO2)*RLC
2149 PI=(1.E5/P)**(0.2854*(1.-0.28*QVAP))
2150 THTEU=TU*PI*EXP(RL*QVAP*C5/TU*(1.+0.81*QVAP))
2151 IF(IFLAG.EQ.1)THEN
2152 QICE=QICE+FRC1*DQVAP+QLIQ
2153 QLIQ=0.
2154 ELSE
2155 QICE=QICE+FRC1*(DQVAP+QLQFRZ)
2156 QLIQ=QLIQ-FRC1*QLQFRZ
2157 ENDIF
2158 QNWFRZ=0.
2159
2160 END SUBROUTINE DTFRZNEW
2161
2162 !-----------------------------------------------------------------------
2163 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2164 ! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2165 ! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN F
2166 ! HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMA
2167 ! TABLES ED. BY ABRAMOWITZ AND STEGUN, NAT L BUREAU OF STANDARDS APPLI
2168 ! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2169 ! JACK KAIN
2170 ! 7/6/89
2171 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2172 !***********************************************************************
2173 !***** GAUSSIAN TYPE MIXING PROFILE....******************************
2174 SUBROUTINE PROF5(EQ,EE,UD)
2175 !-----------------------------------------------------------------------
2176 IMPLICIT NONE
2177 !-----------------------------------------------------------------------
2178 REAL, INTENT(IN ) :: EQ
2179 REAL, INTENT(INOUT) :: EE,UD
2180 REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2181
2182 DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2183 0.9372980,0.33267,0.166666667,0.202765151/
2184 X=(EQ-0.5)/SIGMA
2185 Y=6.*EQ-3.
2186 EY=EXP(Y*Y/(-2))
2187 E45=EXP(-4.5)
2188 T2=1./(1.+P*ABS(Y))
2189 T1=0.500498
2190 C1=A1*T1+A2*T1*T1+A3*T1*T1*T1
2191 C2=A1*T2+A2*T2*T2+A3*T2*T2*T2
2192 IF(Y.GE.0.)THEN
2193 EE=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2194 UD=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*(0.5+EQ*EQ/2.- &
2195 EQ)
2196 ELSE
2197 EE=SIGMA*(0.5*(EY*C2-E45*C1)+SIGMA*(E45-EY))-E45*EQ*EQ/2.
2198 UD=SIGMA*(0.5*(SQRT2P-E45*C1-EY*C2)+SIGMA*(E45-EY))-E45*(0.5+EQ* &
2199 EQ/2.-EQ)
2200 ENDIF
2201 EE=EE/FE
2202 UD=UD/FE
2203
2204 END SUBROUTINE PROF5
2205
2206 !-----------------------------------------------------------------------
2207 SUBROUTINE TPMIX(P,THTU,TU,QU,QLIQ,QICE,QNEWLQ,QNEWIC,RATIO2,RL, &
2208 XLV0,XLV1,XLS0,XLS1, &
2209 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2210 !-----------------------------------------------------------------------
2211 IMPLICIT NONE
2212 !-----------------------------------------------------------------------
2213 REAL, INTENT(IN ) :: XLV0,XLV1
2214 REAL, INTENT(IN ) :: P,THTU,RATIO2,RL,XLS0, &
2215 XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,&
2216 CICE,DICE
2217 REAL, INTENT(INOUT) :: QU,QLIQ,QICE,TU,QNEWLQ,QNEWIC
2218 REAL :: ES,QS,PI,THTGS,F0,T1,T0,C5,RV,ESLIQ,ESICE,F1,DT,QNEW, &
2219 DQ, QTOT,DQICE,DQLIQ,RLL,CCP
2220 INTEGER :: ITCNT
2221 !-----------------------------------------------------------------------
2222 !
2223 !...THIS SUBROUTINE ITERATIVELY EXTRACTS WET-BULB TEMPERATURE FROM EQUIV
2224 ! POTENTIAL TEMPERATURE, THEN CHECKS TO SEE IF SUFFICIENT MOISTURE IS
2225 ! AVAILABLE TO ACHIEVE SATURATION...IF NOT, TEMPERATURE IS ADJUSTED
2226 ! ACCORDINGLY, IF SO, THE RESIDUAL LIQUID WATER/ICE CONCENTRATION IS
2227 ! DETERMINED...
2228 !
2229 C5=1.0723E-3
2230 RV=461.5
2231 !
2232 ! ITERATE TO FIND WET BULB TEMPERATURE AS A FUNCTION OF EQUIVALENT POT
2233 ! TEMP AND PRS, ASSUMING SATURATION VAPOR PRESSURE...RATIO2 IS THE DEG
2234 ! OF GLACIATION...
2235 !
2236 IF(RATIO2.LT.1.E-6)THEN
2237 ES=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2238 QS=EP2*ES/(P-ES)
2239 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2240 THTGS=TU*PI*EXP((3374.6525/TU-2.5403)*QS*(1.+0.81*QS))
2241 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2242 ES=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2243 QS=EP2*ES/(P-ES)
2244 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2245 THTGS=TU*PI*EXP((3114.834/TU-0.278296)*QS*(1.+0.81*QS))
2246 ELSE
2247 ESLIQ=ALIQ*EXP((BLIQ*TU-CLIQ)/(TU-DLIQ))
2248 ESICE=AICE*EXP((BICE*TU-CICE)/(TU-DICE))
2249 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
2250 QS=EP2*ES/(P-ES)
2251 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2252 THTGS=TU*PI*EXP(RL*QS*C5/TU*(1.+0.81*QS))
2253 ENDIF
2254 F0=THTGS-THTU
2255 T1=TU-0.5*F0
2256 T0=TU
2257 ITCNT=0
2258 90 IF(RATIO2.LT.1.E-6)THEN
2259 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2260 QS=EP2*ES/(P-ES)
2261 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2262 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*QS*(1.+0.81*QS))
2263 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2264 ES=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
2265 QS=EP2*ES/(P-ES)
2266 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2267 THTGS=T1*PI*EXP((3114.834/T1-0.278296)*QS*(1.+0.81*QS))
2268 ELSE
2269 ESLIQ=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2270 ESICE=AICE*EXP((BICE*T1-CICE)/(T1-DICE))
2271 ES=(1.-RATIO2)*ESLIQ+RATIO2*ESICE
2272 QS=EP2*ES/(P-ES)
2273 PI=(1.E5/P)**(0.2854*(1.-0.28*QS))
2274 THTGS=T1*PI*EXP(RL*QS*C5/T1*(1.+0.81*QS))
2275 ENDIF
2276 F1=THTGS-THTU
2277 IF(ABS(F1).LT.0.01)GOTO 50
2278 ITCNT=ITCNT+1
2279 IF(ITCNT.GT.10)GOTO 50
2280 DT=F1*(T1-T0)/(F1-F0)
2281 T0=T1
2282 F0=F1
2283 T1=T1-DT
2284 GOTO 90
2285 !
2286 ! IF THE PARCEL IS SUPERSATURATED, CALCULATE CONCENTRATION OF FRESH
2287 ! CONDENSATE...
2288 !
2289 50 IF(QS.LE.QU)THEN
2290 QNEW=QU-QS
2291 QU=QS
2292 GOTO 96
2293 ENDIF
2294 !
2295 ! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2296 ! ADJUSTED...IF LIQUID WATER OR ICE IS PRESENT, IT IS ALLOWED TO EVAPO
2297 ! SUBLIMATE.
2298 !
2299 QNEW=0.
2300 DQ=QS-QU
2301 QTOT=QLIQ+QICE
2302 !
2303 ! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS
2304 ! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MI
2305 ! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURA
2306 ! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPR
2307 ! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE
2308 !
2309 !...NOTE THAT THE LIQ AND ICE MAY BE PRESENT IN PROPORTIONS SLIGHTLY DIF
2310 ! THAN SUGGESTED BY THE VALUE OF RATIO2...CHECK TO MAKE SURE THAT LIQ
2311 ! ICE CONCENTRATIONS ARE NOT REDUCED TO BELOW ZERO WHEN EVAPORATION/
2312 ! SUBLIMATION OCCURS...
2313 !
2314 IF(QTOT.GE.DQ)THEN
2315 DQICE=0.0
2316 DQLIQ=0.0
2317 QLIQ=QLIQ-(1.-RATIO2)*DQ
2318 IF(QLIQ.LT.0.)THEN
2319 DQICE=0.0-QLIQ
2320 QLIQ=0.0
2321 ENDIF
2322 QICE=QICE-RATIO2*DQ+DQICE
2323 IF(QICE.LT.0.)THEN
2324 DQLIQ=0.0-QICE
2325 QICE=0.0
2326 ENDIF
2327 QLIQ=QLIQ+DQLIQ
2328 QU=QS
2329 GOTO 96
2330 ELSE
2331 IF(RATIO2.LT.1.E-6)THEN
2332 RLL=XLV0-XLV1*T1
2333 ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2334 RLL=XLS0-XLS1*T1
2335 ELSE
2336 RLL=RL
2337 ENDIF
2338 CCP=1005.7*(1.+0.89*QU)
2339 IF(QTOT.LT.1.E-10)THEN
2340 !
2341 !...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2342 T1=T1+RLL*(DQ/(1.+DQ))/CCP
2343 GOTO 96
2344 ELSE
2345 !
2346 !...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURA
2347 ! THE TEMPERATURE IS GIVEN BY:
2348 T1=T1+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CCP
2349 QU=QU+QTOT
2350 QTOT=0.
2351 ENDIF
2352 QLIQ=0
2353 QICE=0.
2354 ENDIF
2355 96 TU=T1
2356 QNEWLQ=(1.-RATIO2)*QNEW
2357 QNEWIC=RATIO2*QNEW
2358 IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPMIX =', &
2359 ITCNT
2360
2361 END SUBROUTINE TPMIX
2362 !-----------------------------------------------------------------------
2363 SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,R1,RL, &
2364 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2365 !-----------------------------------------------------------------------
2366 IMPLICIT NONE
2367 !-----------------------------------------------------------------------
2368 REAL, INTENT(IN ) :: P1,T1,Q1,R1,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,&
2369 BICE,CICE,DICE
2370 REAL, INTENT(INOUT) :: THT1
2371 REAL:: T00,P00,C1,C2,C3,C4,C5,EE,TLOG,TDPT,TSAT,THT,TFPT,TLOGIC, &
2372 TSATLQ,TSATIC
2373
2374 DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834,&
2375 0.278296,1.0723E-3/
2376 !
2377 ! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
2378 !
2379
2380 IF(R1.LT.1.E-6)THEN
2381 EE=Q1*P1/(EP2+Q1)
2382 TLOG=ALOG(EE/ALIQ)
2383 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2384 TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2385 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2386 THT1=THT*EXP((C1/TSAT-C2)*Q1*(1.+0.81*Q1))
2387 ELSEIF(ABS(R1-1.).LT.1.E-6)THEN
2388 EE=Q1*P1/(EP2+Q1)
2389 TLOG=ALOG(EE/AICE)
2390 TFPT=(CICE-DICE*TLOG)/(BICE-TLOG)
2391 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2392 TSAT=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
2393 THT1=THT*EXP((C3/TSAT-C4)*Q1*(1.+0.81*Q1))
2394 ELSE
2395 EE=Q1*P1/(EP2+Q1)
2396 TLOG=ALOG(EE/ALIQ)
2397 TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2398 TLOGIC=ALOG(EE/AICE)
2399 TFPT=(CICE-DICE*TLOGIC)/(BICE-TLOGIC)
2400 THT=T1*(P00/P1)**(0.2854*(1.-0.28*Q1))
2401 TSATLQ=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T1-T00))*(T1-TDPT)
2402 TSATIC=TFPT-(.182+1.13E-3*(TFPT-T00)-3.58E-4*(T1-T00))*(T1-TFPT)
2403 TSAT=R1*TSATIC+(1.-R1)*TSATLQ
2404 THT1=THT*EXP(RL*Q1*C5/TSAT*(1.+0.81*Q1))
2405 ENDIF
2406
2407 END SUBROUTINE ENVIRTHT
2408
2409 !-----------------------------------------------------------------------
2410 !************************* TPDD.FOR ************************************
2411 ! THIS SUBROUTINE ITERATIVELY EXTRACTS TEMPERATURE FROM EQUIVALENT *
2412 ! POTENTIAL TEMP. IT IS DESIGNED FOR USE WITH DOWNDRAFT CALCULATIONS.
2413 ! IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, PARCEL *
2414 ! TEMP, SPECIFIC HUMIDITY, AND LIQUID WATER CONTENT ARE ITERATIVELY *
2415 ! CALCULATED. *
2416 !***********************************************************************
2417 FUNCTION TPDD(P,THTED,TGS,RS,RD,RH,XLV0,XLV1, &
2418 EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2419 !-----------------------------------------------------------------------
2420 IMPLICIT NONE
2421 !-----------------------------------------------------------------------
2422 REAL, INTENT(IN ) :: XLV0,XLV1
2423 REAL, INTENT(IN ) :: P,THTED,TGS,RD,RH,EP2,ALIQ,BLIQ, &
2424 CLIQ,DLIQ,AICE,BICE,CICE,DICE
2425 REAL, INTENT(INOUT) :: RS
2426 REAL :: TPDD,ES,PI,THTGS,F0,T1,T0,CCP,F1,DT,RL,DSSDT,T1RH,RSRH
2427 INTEGER :: ITCNT
2428 !-----------------------------------------------------------------------
2429 ES=ALIQ*EXP((BLIQ*TGS-CLIQ)/(TGS-DLIQ))
2430 RS=EP2*ES/(P-ES)
2431 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
2432 THTGS=TGS*PI*EXP((3374.6525/TGS-2.5403)*RS*(1.+0.81*RS))
2433 F0=THTGS-THTED
2434 T1=TGS-0.5*F0
2435 T0=TGS
2436 CCP=1005.7
2437 !
2438 !...ITERATE TO FIND WET-BULB TEMPERATURE...
2439 !
2440 ITCNT=0
2441 90 ES=ALIQ*EXP((BLIQ*T1-CLIQ)/(T1-DLIQ))
2442 RS=EP2*ES/(P-ES)
2443 PI=(1.E5/P)**(0.2854*(1.-0.28*RS))
2444 THTGS=T1*PI*EXP((3374.6525/T1-2.5403)*RS*(1.+0.81*RS))
2445 F1=THTGS-THTED
2446 IF(ABS(F1).LT.0.05)GOTO 50
2447 ITCNT=ITCNT+1
2448 IF(ITCNT.GT.10)GOTO 50
2449 DT=F1*(T1-T0)/(F1-F0)
2450 T0=T1
2451 F0=F1
2452 T1=T1-DT
2453 GOTO 90
2454 50 RL=XLV0-XLV1*T1
2455 !
2456 !...IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, ESTIMATE THE
2457 ! TEMPERATURE AND MIXING RATIO WHICH WILL YIELD THE APPROPRIATE VALUE.
2458 !
2459 IF(RH.EQ.1.)GOTO 110
2460 DSSDT=(CLIQ-BLIQ*DLIQ)/((T1-DLIQ)*(T1-DLIQ))
2461 DT=RL*RS*(1.-RH)/(CCP+RL*RH*RS*DSSDT)
2462 T1RH=T1+DT
2463 ES=RH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
2464 RSRH=EP2*ES/(P-ES)
2465 !
2466 !...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
2467 !...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
2468 !
2469 IF(RSRH.LT.RD)THEN
2470 RSRH=RD
2471 T1RH=T1+(RS-RSRH)*RL/CCP
2472 ENDIF
2473 T1=T1RH
2474 RS=RSRH
2475 110 TPDD=T1
2476 IF(ITCNT.GT.10)PRINT*,'***** NUMBER OF ITERATIONS IN TPDD = ', &
2477 ITCNT
2478
2479 END FUNCTION TPDD
2480
2481 !====================================================================
2482 SUBROUTINE kfinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2483 RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
2484 P_FIRST_SCALAR,restart,allowed_to_read, &
2485 ids, ide, jds, jde, kds, kde, &
2486 ims, ime, jms, jme, kms, kme, &
2487 its, ite, jts, jte, kts, kte )
2488 !--------------------------------------------------------------------
2489 IMPLICIT NONE
2490 !--------------------------------------------------------------------
2491 LOGICAL , INTENT(IN) :: restart, allowed_to_read
2492 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2493 ims, ime, jms, jme, kms, kme, &
2494 its, ite, jts, jte, kts, kte
2495 INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
2496
2497 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
2498 RTHCUTEN, &
2499 RQVCUTEN, &
2500 RQCCUTEN, &
2501 RQRCUTEN, &
2502 RQICUTEN, &
2503 RQSCUTEN
2504
2505 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2506
2507 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2508
2509 INTEGER :: i, j, k, itf, jtf, ktf
2510
2511 jtf=min0(jte,jde-1)
2512 ktf=min0(kte,kde-1)
2513 itf=min0(ite,ide-1)
2514
2515 IF(.not.restart)THEN
2516 DO j=jts,jtf
2517 DO k=kts,ktf
2518 DO i=its,itf
2519 RTHCUTEN(i,k,j)=0.
2520 RQVCUTEN(i,k,j)=0.
2521 RQCCUTEN(i,k,j)=0.
2522 RQRCUTEN(i,k,j)=0.
2523 ENDDO
2524 ENDDO
2525 ENDDO
2526
2527 IF (P_QI .ge. P_FIRST_SCALAR) THEN
2528 DO j=jts,jtf
2529 DO k=kts,ktf
2530 DO i=its,itf
2531 RQICUTEN(i,k,j)=0.
2532 ENDDO
2533 ENDDO
2534 ENDDO
2535 ENDIF
2536
2537 IF (P_QS .ge. P_FIRST_SCALAR) THEN
2538 DO j=jts,jtf
2539 DO k=kts,ktf
2540 DO i=its,itf
2541 RQSCUTEN(i,k,j)=0.
2542 ENDDO
2543 ENDDO
2544 ENDDO
2545 ENDIF
2546
2547 DO j=jts,jtf
2548 DO i=its,itf
2549 NCA(i,j)=-100.
2550 ENDDO
2551 ENDDO
2552
2553 DO j=jts,jtf
2554 DO k=kts,ktf
2555 DO i=its,itf
2556 W0AVG(i,k,j)=0.
2557 ENDDO
2558 ENDDO
2559 ENDDO
2560
2561 ENDIF
2562
2563 END SUBROUTINE kfinit
2564
2565 !-------------------------------------------------------
2566
2567 END MODULE module_cu_kf
2568