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