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