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