module_cu_sas.F
References to this file elsewhere.
1 !
2 MODULE module_cu_sas
3
4 #ifndef IFORT_KLUDGE
5 USE module_configure, ONLY: nl_get_start_year, nl_get_start_month, &
6 nl_get_start_day, nl_get_start_hour
7 #endif
8 CONTAINS
9
10 !-----------------------------------------------------------------
11 SUBROUTINE CU_SAS( &
12 DT,ITIMESTEP,STEPCU &
13 ,RAINCV,HTOP,HBOT &
14 ,U3D,V3D,W,T3D,QV3D,QC3D,QI3D,PI3D,RHO3D &
15 ,DZ8W,PCPS,P8W,XLAND,CU_ACT_FLAG &
16 ,ids,ide, jds,jde, kds,kde &
17 ,ims,ime, jms,jme, kms,kme &
18 ,its,ite, jts,jte, kts,kte &
19 ,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
20 ,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN &
21 )
22
23 !-------------------------------------------------------------------
24 USE MODULE_GFS_MACHINE , ONLY : kind_phys, kind_evod
25 USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
26 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
27 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
28 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
29 &, EPS => con_eps, EPSM1 => con_epsm1 &
30 &, ROVCP => con_rocp, RD => con_rd
31 !-------------------------------------------------------------------
32 IMPLICIT NONE
33 !-------------------------------------------------------------------
34 !-- U3D 3D u-velocity interpolated to theta points (m/s)
35 !-- V3D 3D v-velocity interpolated to theta points (m/s)
36 !-- TH3D 3D potential temperature (K)
37 !-- T3D temperature (K)
38 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
39 !-- QC3D 3D cloud mixing ratio (Kg/Kg)
40 !-- QI3D 3D ice mixing ratio (Kg/Kg)
41 !-- P8w 3D pressure at full levels (Pa)
42 !-- Pcps 3D pressure (Pa)
43 !-- PI3D 3D exner function (dimensionless)
44 !-- rr3D 3D dry air density (kg/m^3)
45 !-- RUBLTEN U tendency due to
46 ! PBL parameterization (m/s^2)
47 !-- RVBLTEN V tendency due to
48 ! PBL parameterization (m/s^2)
49 !-- RTHBLTEN Theta tendency due to
50 ! PBL parameterization (K/s)
51 !-- RQVBLTEN Qv tendency due to
52 ! PBL parameterization (kg/kg/s)
53 !-- RQCBLTEN Qc tendency due to
54 ! PBL parameterization (kg/kg/s)
55 !-- RQIBLTEN Qi tendency due to
56 ! PBL parameterization (kg/kg/s)
57 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
58 !-- GRAV acceleration due to gravity (m/s^2)
59 !-- ROVCP R/CP
60 !-- RD gas constant for dry air (J/kg/K)
61 !-- ROVG R/G
62 !-- P_QI species index for cloud ice
63 !-- dz8w dz between full levels (m)
64 !-- z height above sea level (m)
65 !-- PSFC pressure at the surface (Pa)
66 !-- UST u* in similarity theory (m/s)
67 !-- PBL PBL height (m)
68 !-- PSIM similarity stability function for momentum
69 !-- PSIH similarity stability function for heat
70 !-- HFX upward heat flux at the surface (W/m^2)
71 !-- QFX upward moisture flux at the surface (kg/m^2/s)
72 !-- TSK surface temperature (K)
73 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
74 !-- WSPD wind speed at lowest model level (m/s)
75 !-- BR bulk Richardson number in surface layer
76 !-- DT time step (s)
77 !-- rvovrd R_v divided by R_d (dimensionless)
78 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
79 !-- KARMAN Von Karman constant
80 !-- ids start index for i in domain
81 !-- ide end index for i in domain
82 !-- jds start index for j in domain
83 !-- jde end index for j in domain
84 !-- kds start index for k in domain
85 !-- kde end index for k in domain
86 !-- ims start index for i in memory
87 !-- ime end index for i in memory
88 !-- jms start index for j in memory
89 !-- jme end index for j in memory
90 !-- kms start index for k in memory
91 !-- kme end index for k in memory
92 !-- its start index for i in tile
93 !-- ite end index for i in tile
94 !-- jts start index for j in tile
95 !-- jte end index for j in tile
96 !-- kts start index for k in tile
97 !-- kte end index for k in tile
98 !-------------------------------------------------------------------
99
100 INTEGER :: ICLDCK
101
102 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
103 ims,ime, jms,jme, kms,kme, &
104 its,ite, jts,jte, kts,kte, &
105 ITIMESTEP, &
106 STEPCU
107
108 REAL, INTENT(IN) :: &
109 DT
110
111
112 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
113 XLAND
114
115 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
116 RAINCV
117
118 REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
119 HBOT, &
120 HTOP
121
122 LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) :: &
123 CU_ACT_FLAG
124
125
126 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
127 DZ8W, &
128 P8w, &
129 Pcps, &
130 PI3D, &
131 QC3D, &
132 QI3D, &
133 QV3D, &
134 RHO3D, &
135 T3D, &
136 U3D, &
137 V3D, &
138 W
139
140 !--------------------------- OPTIONAL VARS ----------------------------
141
142 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), &
143 OPTIONAL, INTENT(INOUT) :: &
144 RQCCUTEN, &
145 RQICUTEN, &
146 RQVCUTEN, &
147 RTHCUTEN
148
149 !
150 ! Flags relating to the optional tendency arrays declared above
151 ! Models that carry the optional tendencies will provdide the
152 ! optional arguments at compile time; these flags all the model
153 ! to determine at run-time whether a particular tracer is in
154 ! use or not.
155 !
156 LOGICAL, OPTIONAL :: &
157 F_QV &
158 ,F_QC &
159 ,F_QR &
160 ,F_QI &
161 ,F_QS
162
163
164 !--------------------------- LOCAL VARS ------------------------------
165
166 REAL, DIMENSION(ims:ime, jms:jme) :: &
167 PSFC
168
169 REAL (kind=kind_evod),save :: seed0
170 ! REAL (kind=kind_evod) :: seed0
171 REAL (kind=kind_evod) :: wrk
172
173 REAL (kind=kind_phys) :: &
174 DELT, &
175 DPSHC, &
176 RDELT, &
177 RSEED
178
179 REAL (kind=kind_phys), DIMENSION(ids:ide,jds:jde) :: &
180 RANNUM
181
182 REAL (kind=kind_phys), DIMENSION(its:ite) :: &
183 CLDWRK, &
184 PS, &
185 RCS, &
186 RN, &
187 SLIMSK, &
188 XKT2
189
190 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
191 PRSI
192
193 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
194 DEL, &
195 DOT, &
196 PHIL, &
197 PRSL, &
198 PRSLK, &
199 Q1, &
200 T1, &
201 U1, &
202 V1, &
203 ZI, &
204 ZL
205
206 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) :: &
207 QL
208
209 INTEGER, DIMENSION(its:ite) :: &
210 KBOT, &
211 KTOP, &
212 KUO
213
214 INTEGER :: &
215 I, &
216 ! IGPVS, &
217 IM, &
218 J, &
219 JCAP, &
220 K, &
221 KM, &
222 KP, &
223 KX, &
224 NCLOUD
225
226 INTEGER :: start_year,start_month,start_day,start_hour
227
228 integer :: iseed
229 ! integer, save :: krsize
230 integer :: krsize
231 integer, allocatable :: nrnd(:)
232 real :: fsec
233
234 ! DATA IGPVS/0/
235
236 !-----------------------------------------------------------------------
237 !
238 !*** CHECK TO SEE IF THIS IS A CONVECTION TIMESTEP
239 !
240 ICLDCK=MOD(ITIMESTEP,STEPCU)
241 !-----------------------------------------------------------------------
242
243
244 IF(ICLDCK.EQ.0.OR.ITIMESTEP.EQ.0)THEN
245
246 DO J=JTS,JTE
247 DO I=ITS,ITE
248 CU_ACT_FLAG(I,J)=.TRUE.
249 ENDDO
250 ENDDO
251
252 IM=ITE-ITS+1
253 KX=KTE-KTS+1
254 JCAP=126
255 DPSHC=30_kind_phys
256 DELT=DT*STEPCU
257 RDELT=1./DELT
258 NCLOUD=1
259
260
261 DO J=jms,jme
262 DO I=ims,ime
263 PSFC(i,j)=P8w(i,kms,j)
264 ENDDO
265 ENDDO
266
267 if(itimestep.eq.0) then
268 CALL GFUNCPHYS
269
270 CALL nl_get_start_year(1,start_year)
271 CALL nl_get_start_month(1,start_month)
272 CALL nl_get_start_day(1,start_day)
273 CALL nl_get_start_hour(1,start_hour)
274
275 call random_seed(size=krsize)
276 if (.not. allocated (nrnd)) allocate (nrnd(krsize))
277
278 seed0 = start_year + start_month + start_day + start_hour
279 nrnd = start_hour + start_day*24
280 call random_seed
281 call random_seed(put=nrnd)
282 call random_number(wrk)
283 seed0 = seed0 + nint(wrk*1000.0)
284
285 endif
286
287 fsec = ITIMESTEP*DT
288 iseed = mod(100.0*sqrt(fsec),1.0e9) + 1 + seed0
289 call random_seed(size=krsize)
290 if (.not. allocated (nrnd)) allocate (nrnd(krsize))
291 nrnd = iseed
292 call random_seed
293 call random_seed(put=nrnd)
294 call random_number(rannum)
295
296 ! igpvs=1
297
298 !------------- J LOOP (OUTER) --------------------------------------------------
299
300 DO J=jts,jte
301
302 ! --------------- compute zi and zl -----------------------------------------
303 DO i=its,ite
304 ZI(I,KTS)=0.0
305 ENDDO
306
307 DO k=kts+1,kte
308 KM=K-1
309 DO i=its,ite
310 ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
311 ENDDO
312 ENDDO
313
314 DO k=kts+1,kte
315 KM=K-1
316 DO i=its,ite
317 ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
318 ENDDO
319 ENDDO
320
321 DO i=its,ite
322 ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
323 ENDDO
324
325 ! --------------- end compute zi and zl -------------------------------------
326
327
328 ! call random_number(XKT2)
329 DO i=its,ite
330 xkt2(i)=rannum(i,j)
331 PS(i)=PSFC(i,j)*.001
332 RCS(i)=1.
333 SLIMSK(i)=ABS(XLAND(i,j)-2.)
334 ENDDO
335
336 DO i=its,ite
337 PRSI(i,kts)=PS(i)
338 ENDDO
339
340 DO k=kts,kte
341 kp=k+1
342 DO i=its,ite
343 PRSL(I,K)=Pcps(i,k,j)*.001
344 PHIL(I,K)=ZL(I,K)*GRAV
345 DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
346 ENDDO
347 ENDDO
348
349 DO k=kts,kte
350 DO i=its,ite
351 DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
352 U1(i,k)=U3D(i,k,j)
353 V1(i,k)=V3D(i,k,j)
354 Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
355 T1(i,k)=T3D(i,k,j)
356 QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
357 QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
358 PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
359 ENDDO
360 ENDDO
361
362 DO k=kts+1,kte+1
363 km=k-1
364 DO i=its,ite
365 PRSI(i,k)=PRSI(i,km)-del(i,km)
366 ENDDO
367 ENDDO
368
369
370 CALL SASCNV(IM,IM,KX,JCAP,DELT,DEL,PRSL,PS,PHIL, &
371 QL,Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT, &
372 KTOP,KUO,SLIMSK,DOT,XKT2,NCLOUD)
373
374 CALL SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KUO,Q1,T1,DPSHC)
375
376 DO I=ITS,ITE
377 RAINCV(I,J)=RN(I)*1000./STEPCU
378 HBOT(I,J)=KBOT(I)
379 HTOP(I,J)=KTOP(I)
380 ENDDO
381
382 DO K=KTS,KTE
383 DO I=ITS,ITE
384 RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
385 RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
386 ENDDO
387 ENDDO
388
389 IF(PRESENT(RQCCUTEN))THEN
390 IF ( F_QC ) THEN
391 DO K=KTS,KTE
392 DO I=ITS,ITE
393 RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
394 ENDDO
395 ENDDO
396 ENDIF
397 ENDIF
398
399 IF(PRESENT(RQICUTEN))THEN
400 IF ( F_QI ) THEN
401 DO K=KTS,KTE
402 DO I=ITS,ITE
403 RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
404 ENDDO
405 ENDDO
406 ENDIF
407 ENDIF
408
409
410 ENDDO
411
412 ENDIF
413
414 END SUBROUTINE CU_SAS
415
416 !====================================================================
417 SUBROUTINE sasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
418 RESTART,P_QC,P_QI,P_FIRST_SCALAR, &
419 allowed_to_read, &
420 ids, ide, jds, jde, kds, kde, &
421 ims, ime, jms, jme, kms, kme, &
422 its, ite, jts, jte, kts, kte)
423 !--------------------------------------------------------------------
424 IMPLICIT NONE
425 !--------------------------------------------------------------------
426 LOGICAL , INTENT(IN) :: allowed_to_read,restart
427 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
428 ims, ime, jms, jme, kms, kme, &
429 its, ite, jts, jte, kts, kte
430 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
431
432 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
433 RTHCUTEN, &
434 RQVCUTEN, &
435 RQCCUTEN, &
436 RQICUTEN
437
438 INTEGER :: i, j, k, itf, jtf, ktf
439
440 jtf=min0(jte,jde-1)
441 ktf=min0(kte,kde-1)
442 itf=min0(ite,ide-1)
443
444 IF(.not.restart)THEN
445 DO j=jts,jtf
446 DO k=kts,ktf
447 DO i=its,itf
448 RTHCUTEN(i,k,j)=0.
449 RQVCUTEN(i,k,j)=0.
450 ENDDO
451 ENDDO
452 ENDDO
453
454 IF (P_QC .ge. P_FIRST_SCALAR) THEN
455 DO j=jts,jtf
456 DO k=kts,ktf
457 DO i=its,itf
458 RQCCUTEN(i,k,j)=0.
459 ENDDO
460 ENDDO
461 ENDDO
462 ENDIF
463
464 IF (P_QI .ge. P_FIRST_SCALAR) THEN
465 DO j=jts,jtf
466 DO k=kts,ktf
467 DO i=its,itf
468 RQICUTEN(i,k,j)=0.
469 ENDDO
470 ENDDO
471 ENDDO
472 ENDIF
473 ENDIF
474
475 END SUBROUTINE sasinit
476
477 ! ------------------------------------------------------------------------
478
479 SUBROUTINE SASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PS,PHIL,QL, &
480 ! SUBROUTINE SASCNV(IM,IX,KM,JCAP,DLT,DEL,PRSL,PHIL,QL, &
481 & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK, &
482 & DOT,XKT2,ncloud)
483 ! for cloud water version
484 ! parameter(ncloud=0)
485 ! SUBROUTINE SASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
486 ! & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
487 ! & DOT,xkt2,ncloud)
488 !
489 USE MODULE_GFS_MACHINE , ONLY : kind_phys,kind_evod
490 USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
491 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
492 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
493 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
494 &, EPS => con_eps, EPSM1 => con_epsm1
495
496 implicit none
497 !
498 ! include 'constant.h'
499 !
500 integer IM, IX, KM, JCAP, ncloud, &
501 & KBOT(IM), KTOP(IM), KUO(IM)
502 real(kind=kind_phys) DELT
503 real(kind=kind_phys) PS(IM), DEL(IX,KM), PRSL(IX,KM), &
504 ! real(kind=kind_phys) DEL(IX,KM), PRSL(IX,KM),
505 & QL(IX,KM,2), Q1(IX,KM), T1(IX,KM), &
506 & U1(IX,KM), V1(IX,KM), RCS(IM), &
507 & CLDWRK(IM), RN(IM), SLIMSK(IM), &
508 & DOT(IX,KM), XKT2(IM), PHIL(IX,KM)
509 !
510 integer I, INDX, jmn, k, knumb, latd, lond, km1
511 !
512 real(kind=kind_phys) adw, alpha, alphal, alphas, &
513 & aup, beta, betal, betas, &
514 & c0, cpoel, dellat, delta, &
515 & desdt, deta, detad, dg, &
516 & dh, dhh, dlnsig, dp, &
517 & dq, dqsdp, dqsdt, dt, &
518 & dt2, dtmax, dtmin, dv1, &
519 & dv1q, dv2, dv2q, dv1u, &
520 & dv1v, dv2u, dv2v, dv3u, &
521 & dv3v, dv3, dv3q, dvq1, &
522 & dz, dz1, e1, edtmax, &
523 & edtmaxl, edtmaxs, el2orc, elocp, &
524 & es, etah, &
525 & evef, evfact, evfactl, fact1, &
526 & fact2, factor, fjcap, fkm, &
527 & fuv, g, gamma, onemf, &
528 & onemfu, pdetrn, pdpdwn, pprime, &
529 & qc, qlk, qrch, qs, &
530 & rain, rfact, shear, tem1, &
531 & tem2, terr, val, val1, &
532 & val2, w1, w1l, w1s, &
533 & w2, w2l, w2s, w3, &
534 & w3l, w3s, w4, w4l, &
535 & w4s, xdby, xpw, xpwd, &
536 & xqc, xqrch, xlambu, mbdt, &
537 & tem
538 !
539 !
540 integer JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM), &
541 & KT2(IM), KTCON(IM), LMIN(IM), &
542 & kbm(IM), kbmax(IM), kmax(IM)
543 !
544 real(kind=kind_phys) AA1(IM), ACRT(IM), ACRTFCT(IM), &
545 & DELHBAR(IM), DELQ(IM), DELQ2(IM), &
546 & DELQBAR(IM), DELQEV(IM), DELTBAR(IM), &
547 & DELTV(IM), DTCONV(IM), EDT(IM), &
548 & EDTO(IM), EDTX(IM), FLD(IM), &
549 & HCDO(IM), HKBO(IM), HMAX(IM), &
550 & HMIN(IM), HSBAR(IM), UCDO(IM), &
551 & UKBO(IM), VCDO(IM), VKBO(IM), &
552 & PBCDIF(IM), PDOT(IM), PO(IM,KM), &
553 & PWAVO(IM), PWEVO(IM), &
554 ! & PSFC(IM), PWAVO(IM), PWEVO(IM), &
555 & QCDO(IM), QCOND(IM), QEVAP(IM), &
556 & QKBO(IM), RNTOT(IM), VSHEAR(IM), &
557 & XAA0(IM), XHCD(IM), XHKB(IM), &
558 & XK(IM), XLAMB(IM), XLAMD(IM), &
559 & XMB(IM), XMBMAX(IM), XPWAV(IM), &
560 & XPWEV(IM), XQCD(IM), XQKB(IM)
561 !
562 ! PHYSICAL PARAMETERS
563 PARAMETER(G=grav)
564 PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP, &
565 & EL2ORC=HVAP*HVAP/(RV*CP))
566 PARAMETER(TERR=0.,C0=.002,DELTA=fv)
567 PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
568 ! LOCAL VARIABLES AND ARRAYS
569 real(kind=kind_phys) PFLD(IM,KM), TO(IM,KM), QO(IM,KM), &
570 & UO(IM,KM), VO(IM,KM), QESO(IM,KM)
571 ! cloud water
572 real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM), TVO(IM,KM), &
573 & DBYO(IM,KM), ZO(IM,KM), SUMZ(IM,KM), &
574 & SUMH(IM,KM), HEO(IM,KM), HESO(IM,KM), &
575 & QRCD(IM,KM), DELLAH(IM,KM), DELLAQ(IM,KM),&
576 & DELLAU(IM,KM), DELLAV(IM,KM), HCKO(IM,KM), &
577 & UCKO(IM,KM), VCKO(IM,KM), QCKO(IM,KM), &
578 & ETA(IM,KM), ETAU(IM,KM), ETAD(IM,KM), &
579 & QRCDO(IM,KM), PWO(IM,KM), PWDO(IM,KM), &
580 & RHBAR(IM), TX1(IM)
581 !
582 LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
583 !
584 real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
585 ! SAVE PCRIT, ACRITT
586 DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
587 & 350.,300.,250.,200.,150./
588 DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
589 & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
590 ! GDAS DERIVED ACRIT
591 ! DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
592 ! & .743,.813,.886,.947,1.138,1.377,1.896/
593 !
594 real(kind=kind_phys) TF, TCR, TCRF, RZERO, RONE
595 parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF))
596 parameter (RZERO=0.0,RONE=1.0)
597 !-----------------------------------------------------------------------
598 !
599 km1 = km - 1
600 ! INITIALIZE ARRAYS
601 !
602 DO I=1,IM
603 RN(I)=0.
604 KBOT(I)=KM+1
605 KTOP(I)=0
606 KUO(I)=0
607 CNVFLG(I) = .TRUE.
608 DTCONV(I) = 3600.
609 CLDWRK(I) = 0.
610 PDOT(I) = 0.
611 KT2(I) = 0
612 QLKO_KTCON(I) = 0.
613 DELLAL(I) = 0.
614 ENDDO
615 !!
616 DO K = 1, 15
617 ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
618 ENDDO
619 DT2 = DELT
620 !cmr dtmin = max(dt2,1200.)
621 val = 1200.
622 dtmin = max(dt2, val )
623 !cmr dtmax = max(dt2,3600.)
624 val = 3600.
625 dtmax = max(dt2, val )
626 ! MODEL TUNABLE PARAMETERS ARE ALL HERE
627 MBDT = 10.
628 EDTMAXl = .3
629 EDTMAXs = .3
630 ALPHAl = .5
631 ALPHAs = .5
632 BETAl = .15
633 betas = .15
634 BETAl = .05
635 betas = .05
636 ! EVEF = 0.07
637 evfact = 0.3
638 evfactl = 0.3
639 PDPDWN = 0.
640 PDETRN = 200.
641 xlambu = 1.e-4
642 fjcap = (float(jcap) / 126.) ** 2
643 !cmr fjcap = max(fjcap,1.)
644 val = 1.
645 fjcap = max(fjcap,val)
646 fkm = (float(km) / 28.) ** 2
647 !cmr fkm = max(fkm,1.)
648 fkm = max(fkm,val)
649 W1l = -8.E-3
650 W2l = -4.E-2
651 W3l = -5.E-3
652 W4l = -5.E-4
653 W1s = -2.E-4
654 W2s = -2.E-3
655 W3s = -1.E-3
656 W4s = -2.E-5
657 !CCCC IF(IM.EQ.384) THEN
658 LATD = 92
659 lond = 189
660 !CCCC ELSEIF(IM.EQ.768) THEN
661 !CCCC LATD = 80
662 !CCCC ELSE
663 !CCCC LATD = 0
664 !CCCC ENDIF
665 !
666 ! DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
667 ! AND THE MAXIMUM THETAE FOR UPDRAFT
668 !
669 DO I=1,IM
670 KBMAX(I) = KM
671 KBM(I) = KM
672 KMAX(I) = KM
673 TX1(I) = 1.0 / PS(I)
674 ENDDO
675 !
676 DO K = 1, KM
677 DO I=1,IM
678 IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
679 IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I) = K + 1
680 IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
681 ENDDO
682 ENDDO
683 DO I=1,IM
684 KBMAX(I) = MIN(KBMAX(I),KMAX(I))
685 KBM(I) = MIN(KBM(I),KMAX(I))
686 ENDDO
687 !
688 ! CONVERT SURFACE PRESSURE TO MB FROM CB
689 !
690 !!
691 DO K = 1, KM
692 DO I=1,IM
693 if (K .le. kmax(i)) then
694 PFLD(I,k) = PRSL(I,K) * 10.0
695 PWO(I,k) = 0.
696 PWDO(I,k) = 0.
697 TO(I,k) = T1(I,k)
698 QO(I,k) = Q1(I,k)
699 UO(I,k) = U1(I,k)
700 VO(I,k) = V1(I,k)
701 DBYO(I,k) = 0.
702 SUMZ(I,k) = 0.
703 SUMH(I,k) = 0.
704 endif
705 ENDDO
706 ENDDO
707
708 !
709 ! COLUMN VARIABLES
710 ! P IS PRESSURE OF THE LAYER (MB)
711 ! T IS TEMPERATURE AT T-DT (K)..TN
712 ! Q IS MIXING RATIO AT T-DT (KG/KG)..QN
713 ! TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
714 ! QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
715 !
716 DO K = 1, KM
717 DO I=1,IM
718 if (k .le. kmax(i)) then
719 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
720 !
721 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
722 !
723 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
724 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
725 val1 = 1.E-8
726 QESO(I,k) = MAX(QESO(I,k), val1)
727 !cmr QO(I,k) = max(QO(I,k),1.e-10)
728 val2 = 1.e-10
729 QO(I,k) = max(QO(I,k), val2 )
730 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
731 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
732 endif
733 ENDDO
734 ENDDO
735
736 !
737 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
738 !
739 DO K = 1, KM
740 DO I=1,IM
741 ZO(I,k) = PHIL(I,k) / G
742 ENDDO
743 ENDDO
744 ! COMPUTE MOIST STATIC ENERGY
745 DO K = 1, KM
746 DO I=1,IM
747 if (K .le. kmax(i)) then
748 ! tem = G * ZO(I,k) + CP * TO(I,k)
749 tem = PHIL(I,k) + CP * TO(I,k)
750 HEO(I,k) = tem + HVAP * QO(I,k)
751 HESO(I,k) = tem + HVAP * QESO(I,k)
752 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
753 endif
754 ENDDO
755 ENDDO
756 !
757 ! DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
758 ! THIS IS THE LEVEL WHERE UPDRAFT STARTS
759 !
760 DO I=1,IM
761 HMAX(I) = HEO(I,1)
762 KB(I) = 1
763 ENDDO
764 !!
765 DO K = 2, KM
766 DO I=1,IM
767 if (k .le. kbm(i)) then
768 IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
769 KB(I) = K
770 HMAX(I) = HEO(I,k)
771 ENDIF
772 endif
773 ENDDO
774 ENDDO
775 ! DO K = 1, KMAX - 1
776 ! TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
777 ! QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
778 ! QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
779 ! HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
780 ! HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
781 ! ENDDO
782 DO K = 1, KM1
783 DO I=1,IM
784 if (k .le. kmax(i)-1) then
785 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
786 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
787 !jfe ES = 10. * FPVS(TO(I,k+1))
788 !
789 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
790 !
791 PPRIME = PFLD(I,k+1) + EPSM1 * ES
792 QS = EPS * ES / PPRIME
793 DQSDP = - QS / PPRIME
794 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
795 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
796 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
797 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
798 DQ = DQSDT * DT + DQSDP * DP
799 TO(I,k) = TO(I,k+1) + DT
800 QO(I,k) = QO(I,k+1) + DQ
801 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
802 endif
803 ENDDO
804 ENDDO
805 !
806 DO K = 1, KM1
807 DO I=1,IM
808 if (k .le. kmax(I)-1) then
809 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
810 !
811 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
812 !
813 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
814 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
815 val1 = 1.E-8
816 QESO(I,k) = MAX(QESO(I,k), val1)
817 !cmr QO(I,k) = max(QO(I,k),1.e-10)
818 val2 = 1.e-10
819 QO(I,k) = max(QO(I,k), val2 )
820 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
821 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
822 & CP * TO(I,k) + HVAP * QO(I,k)
823 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
824 & CP * TO(I,k) + HVAP * QESO(I,k)
825 UO(I,k) = .5 * (UO(I,k) + UO(I,k+1))
826 VO(I,k) = .5 * (VO(I,k) + VO(I,k+1))
827 endif
828 ENDDO
829 ENDDO
830 ! k = kmax
831 ! HEO(I,k) = HEO(I,k)
832 ! hesol(k) = HESO(I,k)
833 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
834 ! PRINT *, ' HEO ='
835 ! PRINT 6001, (HEO(I,K),K=1,KMAX)
836 ! PRINT *, ' HESO ='
837 ! PRINT 6001, (HESO(I,K),K=1,KMAX)
838 ! PRINT *, ' TO ='
839 ! PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
840 ! PRINT *, ' QO ='
841 ! PRINT 6003, (QO(I,K),K=1,KMAX)
842 ! PRINT *, ' QSO ='
843 ! PRINT 6003, (QESO(I,K),K=1,KMAX)
844 ! ENDIF
845 !
846 ! LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
847 !
848 DO I=1,IM
849 IF(CNVFLG(I)) THEN
850 INDX = KB(I)
851 HKBO(I) = HEO(I,INDX)
852 QKBO(I) = QO(I,INDX)
853 UKBO(I) = UO(I,INDX)
854 VKBO(I) = VO(I,INDX)
855 ENDIF
856 FLG(I) = CNVFLG(I)
857 KBCON(I) = KMAX(I)
858 ENDDO
859 !!
860 DO K = 1, KM
861 DO I=1,IM
862 if (k .le. kbmax(i)) then
863 IF(FLG(I).AND.K.GT.KB(I)) THEN
864 HSBAR(I) = HESO(I,k)
865 IF(HKBO(I).GT.HSBAR(I)) THEN
866 FLG(I) = .FALSE.
867 KBCON(I) = K
868 ENDIF
869 ENDIF
870 endif
871 ENDDO
872 ENDDO
873 DO I=1,IM
874 IF(CNVFLG(I)) THEN
875 PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
876 PDOT(I) = 10.* DOT(I,KBCON(I))
877 IF(PBCDIF(I).GT.150.) CNVFLG(I) = .FALSE.
878 IF(KBCON(I).EQ.KMAX(I)) CNVFLG(I) = .FALSE.
879 ENDIF
880 ENDDO
881 !!
882 TOTFLG = .TRUE.
883 DO I=1,IM
884 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
885 ENDDO
886 IF(TOTFLG) RETURN
887 ! FOUND LFC, CAN DEFINE REST OF VARIABLES
888 6001 FORMAT(2X,-2P10F12.2)
889 6002 FORMAT(2X,10F12.2)
890 6003 FORMAT(2X,3P10F12.2)
891
892 !
893 ! DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
894 !
895 DO I = 1, IM
896 alpha = alphas
897 if(SLIMSK(I).eq.1.) alpha = alphal
898 IF(CNVFLG(I)) THEN
899 IF(KB(I).EQ.1) THEN
900 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
901 ELSE
902 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) &
903 & - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
904 ENDIF
905 IF(KBCON(I).NE.KB(I)) THEN
906 !cmr XLAMB(I) = -ALOG(ALPHA) / DZ
907 XLAMB(I) = - LOG(ALPHA) / DZ
908 ELSE
909 XLAMB(I) = 0.
910 ENDIF
911 ENDIF
912 ENDDO
913 ! DETERMINE UPDRAFT MASS FLUX
914 DO K = 1, KM
915 DO I = 1, IM
916 if (k .le. kmax(i) .and. CNVFLG(I)) then
917 ETA(I,k) = 1.
918 ETAU(I,k) = 1.
919 ENDIF
920 ENDDO
921 ENDDO
922 DO K = KM1, 2, -1
923 DO I = 1, IM
924 if (k .le. kbmax(i)) then
925 IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
926 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
927 ETA(I,k) = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
928 ETAU(I,k) = ETA(I,k)
929 ENDIF
930 endif
931 ENDDO
932 ENDDO
933 DO I = 1, IM
934 IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
935 DZ = .5 * (ZO(I,2) - ZO(I,1))
936 ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
937 ETAU(I,1) = ETA(I,1)
938 ENDIF
939 ENDDO
940 !
941 ! WORK UP UPDRAFT CLOUD PROPERTIES
942 !
943 DO I = 1, IM
944 IF(CNVFLG(I)) THEN
945 INDX = KB(I)
946 HCKO(I,INDX) = HKBO(I)
947 QCKO(I,INDX) = QKBO(I)
948 UCKO(I,INDX) = UKBO(I)
949 VCKO(I,INDX) = VKBO(I)
950 PWAVO(I) = 0.
951 ENDIF
952 ENDDO
953 !
954 ! CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
955 !
956 DO K = 2, KM1
957 DO I = 1, IM
958 if (k .le. kmax(i)-1) then
959 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
960 FACTOR = ETA(I,k-1) / ETA(I,k)
961 ONEMF = 1. - FACTOR
962 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
963 & .5 * (HEO(I,k) + HEO(I,k+1))
964 UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF * &
965 & .5 * (UO(I,k) + UO(I,k+1))
966 VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF * &
967 & .5 * (VO(I,k) + VO(I,k+1))
968 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
969 ENDIF
970 IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
971 HCKO(I,k) = HCKO(I,k-1)
972 UCKO(I,k) = UCKO(I,k-1)
973 VCKO(I,k) = VCKO(I,k-1)
974 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
975 ENDIF
976 endif
977 ENDDO
978 ENDDO
979 ! DETERMINE CLOUD TOP
980 DO I = 1, IM
981 FLG(I) = CNVFLG(I)
982 KTCON(I) = 1
983 ENDDO
984 ! DO K = 2, KMAX
985 ! KK = KMAX - K + 1
986 ! IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
987 ! KTCON(I) = KK + 1
988 ! FLG(I) = .FALSE.
989 ! ENDIF
990 ! ENDDO
991 DO K = 2, KM
992 DO I = 1, IM
993 if (k .le. kmax(i)) then
994 IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
995 KTCON(I) = K
996 FLG(I) = .FALSE.
997 ENDIF
998 endif
999 ENDDO
1000 ENDDO
1001 DO I = 1, IM
1002 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.) &
1003 & CNVFLG(I) = .FALSE.
1004 ENDDO
1005 TOTFLG = .TRUE.
1006 DO I = 1, IM
1007 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1008 ENDDO
1009 IF(TOTFLG) RETURN
1010 !
1011 ! SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
1012 !
1013 DO I = 1, IM
1014 HMIN(I) = HEO(I,KBCON(I))
1015 LMIN(I) = KBMAX(I)
1016 JMIN(I) = KBMAX(I)
1017 ENDDO
1018 DO I = 1, IM
1019 DO K = KBCON(I), KBMAX(I)
1020 IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
1021 LMIN(I) = K + 1
1022 HMIN(I) = HEO(I,k)
1023 ENDIF
1024 ENDDO
1025 ENDDO
1026 !
1027 ! Make sure that JMIN(I) is within the cloud
1028 !
1029 DO I = 1, IM
1030 IF(CNVFLG(I)) THEN
1031 JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
1032 XMBMAX(I) = .1
1033 JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
1034 ENDIF
1035 ENDDO
1036 !
1037 ! ENTRAINING CLOUD
1038 !
1039 do k = 2, km1
1040 DO I = 1, IM
1041 if (k .le. kmax(i)-1) then
1042 if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
1043 SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
1044 SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1)) &
1045 & * HEO(I,k)
1046 ENDIF
1047 endif
1048 enddo
1049 enddo
1050 !!
1051 DO I = 1, IM
1052 IF(CNVFLG(I)) THEN
1053 ! call random_number(XKT2)
1054 ! call srand(fhour)
1055 ! XKT2(I) = rand()
1056 KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
1057 ! KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1058 ! KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1059 tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1060 tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1061 if (abs(tem2) .gt. 0.000001) THEN
1062 XLAMB(I) = tem1 / tem2
1063 else
1064 CNVFLG(I) = .false.
1065 ENDIF
1066 ! XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1067 ! & / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1068 XLAMB(I) = max(XLAMB(I),RZERO)
1069 XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
1070 ENDIF
1071 ENDDO
1072 !!
1073 DO I = 1, IM
1074 DWNFLG(I) = CNVFLG(I)
1075 DWNFLG2(I) = CNVFLG(I)
1076 IF(CNVFLG(I)) THEN
1077 if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
1078 if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)&
1079 & DWNFLG(I) = .false.
1080 do k = JMIN(I), KT2(I)
1081 if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
1082 enddo
1083 ! IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
1084 ! & DWNFLG(I)=.FALSE.
1085 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN) &
1086 & DWNFLG2(I)=.FALSE.
1087 ENDIF
1088 ENDDO
1089 !!
1090 DO K = 2, KM1
1091 DO I = 1, IM
1092 if (k .le. kmax(i)-1) then
1093 IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1094 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1095 ! ETA(I,k) = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
1096 ! to simplify matter, we will take the linear approach here
1097 !
1098 ETA(I,k) = ETA(I,k-1) * (1. + XLAMB(I) * dz)
1099 ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
1100 ENDIF
1101 endif
1102 ENDDO
1103 ENDDO
1104 !!
1105 DO K = 2, KM1
1106 DO I = 1, IM
1107 if (k .le. kmax(i)-1) then
1108 ! IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1109 IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
1110 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1111 ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
1112 ENDIF
1113 endif
1114 ENDDO
1115 ENDDO
1116 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1117 ! PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
1118 ! PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
1119 ! ENDIF
1120 ! IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
1121 ! print *, ' xlamb =', xlamb
1122 ! print *, ' eta =', (eta(k),k=1,KT2(I))
1123 ! print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
1124 ! print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
1125 ! print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
1126 ! print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
1127 ! ENDIF
1128 DO I = 1, IM
1129 if(DWNFLG(I)) THEN
1130 KTCON(I) = KT2(I)
1131 ENDIF
1132 ENDDO
1133 !
1134 ! CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
1135 !
1136 DO K = 2, KM1
1137 DO I = 1, IM
1138 if (k .le. kmax(i)-1) then
1139 !jfe
1140 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1141 !jfe IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1142 FACTOR = ETA(I,k-1) / ETA(I,k)
1143 ONEMF = 1. - FACTOR
1144 fuv = ETAU(I,k-1) / ETAU(I,k)
1145 onemfu = 1. - fuv
1146 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1147 & .5 * (HEO(I,k) + HEO(I,k+1))
1148 UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu * &
1149 & .5 * (UO(I,k) + UO(I,k+1))
1150 VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu * &
1151 & .5 * (VO(I,k) + VO(I,k+1))
1152 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1153 ENDIF
1154 endif
1155 ENDDO
1156 ENDDO
1157 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1158 ! PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
1159 ! PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
1160 ! ENDIF
1161 DO I = 1, IM
1162 if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I)) &
1163 & THEN
1164 CNVFLG(I) = .false.
1165 DWNFLG(I) = .false.
1166 DWNFLG2(I) = .false.
1167 ENDIF
1168 ENDDO
1169 !!
1170 TOTFLG = .TRUE.
1171 DO I = 1, IM
1172 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1173 ENDDO
1174 IF(TOTFLG) RETURN
1175 !!
1176 !
1177 ! COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
1178 !
1179 DO I = 1, IM
1180 AA1(I) = 0.
1181 RHBAR(I) = 0.
1182 ENDDO
1183 DO K = 1, KM
1184 DO I = 1, IM
1185 if (k .le. kmax(i)) then
1186 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1187 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1188 DZ1 = (ZO(I,k) - ZO(I,k-1))
1189 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1190 QRCH = QESO(I,k) &
1191 & + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
1192 FACTOR = ETA(I,k-1) / ETA(I,k)
1193 ONEMF = 1. - FACTOR
1194 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1195 & .5 * (QO(I,k) + QO(I,k+1))
1196 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
1197 RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
1198 !
1199 ! BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1200 !
1201 IF(DQ.GT.0.) THEN
1202 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1203 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1204 AA1(I) = AA1(I) - DZ1 * G * QLK
1205 QC = QLK + QRCH
1206 PWO(I,k) = ETAH * C0 * DZ * QLK
1207 QCKO(I,k) = QC
1208 PWAVO(I) = PWAVO(I) + PWO(I,k)
1209 ENDIF
1210 ENDIF
1211 endif
1212 ENDDO
1213 ENDDO
1214 DO I = 1, IM
1215 RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
1216 ENDDO
1217 !
1218 ! this section is ready for cloud water
1219 !
1220 if(ncloud.gt.0) THEN
1221 !
1222 ! compute liquid and vapor separation at cloud top
1223 !
1224 DO I = 1, IM
1225 k = KTCON(I)
1226 IF(CNVFLG(I)) THEN
1227 GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
1228 QRCH = QESO(I,K) &
1229 & + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
1230 DQ = QCKO(I,K-1) - QRCH
1231 !
1232 ! CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1233 !
1234 IF(DQ.GT.0.) THEN
1235 QLKO_KTCON(I) = dq
1236 QCKO(I,K-1) = QRCH
1237 ENDIF
1238 ENDIF
1239 ENDDO
1240 ENDIF
1241 !
1242 ! CALCULATE CLOUD WORK FUNCTION AT T+DT
1243 !
1244 DO K = 1, KM
1245 DO I = 1, IM
1246 if (k .le. kmax(i)) then
1247 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1248 DZ1 = ZO(I,k) - ZO(I,k-1)
1249 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1250 RFACT = 1. + DELTA * CP * GAMMA &
1251 & * TO(I,k-1) / HVAP
1252 AA1(I) = AA1(I) + &
1253 & DZ1 * (G / (CP * TO(I,k-1))) &
1254 & * DBYO(I,k-1) / (1. + GAMMA) &
1255 & * RFACT
1256 val = 0.
1257 AA1(I)=AA1(I)+ &
1258 & DZ1 * G * DELTA * &
1259 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1260 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1261 ENDIF
1262 endif
1263 ENDDO
1264 ENDDO
1265 DO I = 1, IM
1266 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1267 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1268 IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1269 ENDDO
1270 !!
1271 TOTFLG = .TRUE.
1272 DO I = 1, IM
1273 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1274 ENDDO
1275 IF(TOTFLG) RETURN
1276 !!
1277 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1278 !cccc PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
1279 !cccc ENDIF
1280 !
1281 !------- DOWNDRAFT CALCULATIONS
1282 !
1283 !
1284 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1285 !
1286 DO I = 1, IM
1287 IF(CNVFLG(I)) THEN
1288 VSHEAR(I) = 0.
1289 ENDIF
1290 ENDDO
1291 DO K = 1, KM
1292 DO I = 1, IM
1293 if (k .le. kmax(i)) then
1294 IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
1295 shear=rcs(I) * sqrt((UO(I,k+1)-UO(I,k)) ** 2 &
1296 & + (VO(I,k+1)-VO(I,k)) ** 2)
1297 VSHEAR(I) = VSHEAR(I) + SHEAR
1298 ENDIF
1299 endif
1300 ENDDO
1301 ENDDO
1302 DO I = 1, IM
1303 EDT(I) = 0.
1304 IF(CNVFLG(I)) THEN
1305 KNUMB = KTCON(I) - KB(I) + 1
1306 KNUMB = MAX(KNUMB,1)
1307 VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
1308 E1=1.591-.639*VSHEAR(I) &
1309 & +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
1310 EDT(I)=1.-E1
1311 !cmr EDT(I) = MIN(EDT(I),.9)
1312 val = .9
1313 EDT(I) = MIN(EDT(I),val)
1314 !cmr EDT(I) = MAX(EDT(I),.0)
1315 val = .0
1316 EDT(I) = MAX(EDT(I),val)
1317 EDTO(I)=EDT(I)
1318 EDTX(I)=EDT(I)
1319 ENDIF
1320 ENDDO
1321 ! DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
1322 DO I = 1, IM
1323 KBDTR(I) = KBCON(I)
1324 beta = betas
1325 if(SLIMSK(I).eq.1.) beta = betal
1326 IF(CNVFLG(I)) THEN
1327 KBDTR(I) = KBCON(I)
1328 KBDTR(I) = MAX(KBDTR(I),1)
1329 XLAMD(I) = 0.
1330 IF(KBDTR(I).GT.1) THEN
1331 DZ = .5 * ZO(I,KBDTR(I)) + .5 * ZO(I,KBDTR(I)-1) &
1332 & - ZO(I,1)
1333 XLAMD(I) = LOG(BETA) / DZ
1334 ENDIF
1335 ENDIF
1336 ENDDO
1337 ! DETERMINE DOWNDRAFT MASS FLUX
1338 DO K = 1, KM
1339 DO I = 1, IM
1340 IF(k .le. kmax(i)) then
1341 IF(CNVFLG(I)) THEN
1342 ETAD(I,k) = 1.
1343 ENDIF
1344 QRCDO(I,k) = 0.
1345 endif
1346 ENDDO
1347 ENDDO
1348 DO K = KM1, 2, -1
1349 DO I = 1, IM
1350 if (k .le. kbmax(i)) then
1351 IF(CNVFLG(I).AND.K.LT.KBDTR(I)) THEN
1352 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1353 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1354 ENDIF
1355 endif
1356 ENDDO
1357 ENDDO
1358 K = 1
1359 DO I = 1, IM
1360 IF(CNVFLG(I).AND.KBDTR(I).GT.1) THEN
1361 DZ = .5 * (ZO(I,2) - ZO(I,1))
1362 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1363 ENDIF
1364 ENDDO
1365 !
1366 !--- DOWNDRAFT MOISTURE PROPERTIES
1367 !
1368 DO I = 1, IM
1369 PWEVO(I) = 0.
1370 FLG(I) = CNVFLG(I)
1371 ENDDO
1372 DO I = 1, IM
1373 IF(CNVFLG(I)) THEN
1374 JMN = JMIN(I)
1375 HCDO(I) = HEO(I,JMN)
1376 QCDO(I) = QO(I,JMN)
1377 QRCDO(I,JMN) = QESO(I,JMN)
1378 UCDO(I) = UO(I,JMN)
1379 VCDO(I) = VO(I,JMN)
1380 ENDIF
1381 ENDDO
1382 DO K = KM1, 1, -1
1383 DO I = 1, IM
1384 if (k .le. kmax(i)-1) then
1385 IF(CNVFLG(I).AND.K.LT.JMIN(I)) THEN
1386 DQ = QESO(I,k)
1387 DT = TO(I,k)
1388 GAMMA = EL2ORC * DQ / DT**2
1389 DH = HCDO(I) - HESO(I,k)
1390 QRCDO(I,k) = DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1391 DETAD = ETAD(I,k+1) - ETAD(I,k)
1392 PWDO(I,k) = ETAD(I,k+1) * QCDO(I) - &
1393 & ETAD(I,k) * QRCDO(I,k)
1394 PWDO(I,k) = PWDO(I,k) - DETAD * &
1395 & .5 * (QRCDO(I,k) + QRCDO(I,k+1))
1396 QCDO(I) = QRCDO(I,k)
1397 PWEVO(I) = PWEVO(I) + PWDO(I,k)
1398 ENDIF
1399 endif
1400 ENDDO
1401 ENDDO
1402 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG(I)) THEN
1403 ! PRINT *, ' PWAVO(I), PWEVO(I) =', PWAVO(I), PWEVO(I)
1404 ! ENDIF
1405 !
1406 !--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
1407 !--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
1408 !--- EVAPORATE (PWEV)
1409 !
1410 DO I = 1, IM
1411 edtmax = edtmaxl
1412 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1413 IF(DWNFLG2(I)) THEN
1414 IF(PWEVO(I).LT.0.) THEN
1415 EDTO(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
1416 EDTO(I) = MIN(EDTO(I),EDTMAX)
1417 ELSE
1418 EDTO(I) = 0.
1419 ENDIF
1420 ELSE
1421 EDTO(I) = 0.
1422 ENDIF
1423 ENDDO
1424 !
1425 !
1426 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1427 !
1428 !
1429 DO K = KM1, 1, -1
1430 DO I = 1, IM
1431 if (k .le. kmax(i)-1) then
1432 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1433 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1434 DHH=HCDO(I)
1435 DT=TO(I,k+1)
1436 DG=GAMMA
1437 DH=HESO(I,k+1)
1438 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1439 AA1(I)=AA1(I)+EDTO(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1440 & *(1.+DELTA*CP*DG*DT/HVAP)
1441 val=0.
1442 AA1(I)=AA1(I)+EDTO(I)* &
1443 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
1444 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1445 ENDIF
1446 endif
1447 ENDDO
1448 ENDDO
1449 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1450 !cccc PRINT *, ' AA1(I) AFTER DWNDRFT =', AA1(I)
1451 !cccc ENDIF
1452 DO I = 1, IM
1453 IF(AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1454 IF(AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1455 IF(AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1456 ENDDO
1457 !!
1458 TOTFLG = .TRUE.
1459 DO I = 1, IM
1460 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1461 ENDDO
1462 IF(TOTFLG) RETURN
1463 !!
1464 !
1465 !
1466 !--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
1467 !--- WILL DO TO THE ENVIRONMENT?
1468 !
1469 DO K = 1, KM
1470 DO I = 1, IM
1471 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1472 DELLAH(I,k) = 0.
1473 DELLAQ(I,k) = 0.
1474 DELLAU(I,k) = 0.
1475 DELLAV(I,k) = 0.
1476 ENDIF
1477 ENDDO
1478 ENDDO
1479 DO I = 1, IM
1480 IF(CNVFLG(I)) THEN
1481 DP = 1000. * DEL(I,1)
1482 DELLAH(I,1) = EDTO(I) * ETAD(I,1) * (HCDO(I) &
1483 & - HEO(I,1)) * G / DP
1484 DELLAQ(I,1) = EDTO(I) * ETAD(I,1) * (QCDO(I) &
1485 & - QO(I,1)) * G / DP
1486 DELLAU(I,1) = EDTO(I) * ETAD(I,1) * (UCDO(I) &
1487 & - UO(I,1)) * G / DP
1488 DELLAV(I,1) = EDTO(I) * ETAD(I,1) * (VCDO(I) &
1489 & - VO(I,1)) * G / DP
1490 ENDIF
1491 ENDDO
1492 !
1493 !--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1494 !
1495 DO K = 2, KM1
1496 DO I = 1, IM
1497 if (k .le. kmax(i)-1) then
1498 IF(CNVFLG(I).AND.K.LT.KTCON(I)) THEN
1499 AUP = 1.
1500 IF(K.LE.KB(I)) AUP = 0.
1501 ADW = 1.
1502 IF(K.GT.JMIN(I)) ADW = 0.
1503 DV1= HEO(I,k)
1504 DV2 = .5 * (HEO(I,k) + HEO(I,k+1))
1505 DV3= HEO(I,k-1)
1506 DV1Q= QO(I,k)
1507 DV2Q = .5 * (QO(I,k) + QO(I,k+1))
1508 DV3Q= QO(I,k-1)
1509 DV1U= UO(I,k)
1510 DV2U = .5 * (UO(I,k) + UO(I,k+1))
1511 DV3U= UO(I,k-1)
1512 DV1V= VO(I,k)
1513 DV2V = .5 * (VO(I,k) + VO(I,k+1))
1514 DV3V= VO(I,k-1)
1515 DP = 1000. * DEL(I,K)
1516 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1517 DETA = ETA(I,k) - ETA(I,k-1)
1518 DETAD = ETAD(I,k) - ETAD(I,k-1)
1519 DELLAH(I,k) = DELLAH(I,k) + &
1520 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1 &
1521 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3 &
1522 & - AUP * DETA * DV2 &
1523 & + ADW * EDTO(I) * DETAD * HCDO(I)) * G / DP
1524 DELLAQ(I,k) = DELLAQ(I,k) + &
1525 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1Q &
1526 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3Q &
1527 & - AUP * DETA * DV2Q &
1528 & +ADW*EDTO(I)*DETAD*.5*(QRCDO(I,k)+QRCDO(I,k-1))) * G / DP
1529 DELLAU(I,k) = DELLAU(I,k) + &
1530 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1U &
1531 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3U &
1532 & - AUP * DETA * DV2U &
1533 & + ADW * EDTO(I) * DETAD * UCDO(I) &
1534 & ) * G / DP
1535 DELLAV(I,k) = DELLAV(I,k) + &
1536 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1V &
1537 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3V &
1538 & - AUP * DETA * DV2V &
1539 & + ADW * EDTO(I) * DETAD * VCDO(I) &
1540 & ) * G / DP
1541 ENDIF
1542 endif
1543 ENDDO
1544 ENDDO
1545 !
1546 !------- CLOUD TOP
1547 !
1548 DO I = 1, IM
1549 IF(CNVFLG(I)) THEN
1550 INDX = KTCON(I)
1551 DP = 1000. * DEL(I,INDX)
1552 DV1 = HEO(I,INDX-1)
1553 DELLAH(I,INDX) = ETA(I,INDX-1) * &
1554 & (HCKO(I,INDX-1) - DV1) * G / DP
1555 DVQ1 = QO(I,INDX-1)
1556 DELLAQ(I,INDX) = ETA(I,INDX-1) * &
1557 & (QCKO(I,INDX-1) - DVQ1) * G / DP
1558 DV1U = UO(I,INDX-1)
1559 DELLAU(I,INDX) = ETA(I,INDX-1) * &
1560 & (UCKO(I,INDX-1) - DV1U) * G / DP
1561 DV1V = VO(I,INDX-1)
1562 DELLAV(I,INDX) = ETA(I,INDX-1) * &
1563 & (VCKO(I,INDX-1) - DV1V) * G / DP
1564 !
1565 ! cloud water
1566 !
1567 DELLAL(I) = ETA(I,INDX-1) * QLKO_KTCON(I) * g / dp
1568 ENDIF
1569 ENDDO
1570 !
1571 !------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
1572 !
1573 DO K = 1, KM
1574 DO I = 1, IM
1575 if (k .le. kmax(i)) then
1576 IF(CNVFLG(I).and.k.gt.KTCON(I)) THEN
1577 QO(I,k) = Q1(I,k)
1578 TO(I,k) = T1(I,k)
1579 UO(I,k) = U1(I,k)
1580 VO(I,k) = V1(I,k)
1581 ENDIF
1582 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1583 QO(I,k) = DELLAQ(I,k) * MBDT + Q1(I,k)
1584 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
1585 TO(I,k) = DELLAT * MBDT + T1(I,k)
1586 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1587 val = 1.e-10
1588 QO(I,k) = max(QO(I,k), val )
1589 ENDIF
1590 endif
1591 ENDDO
1592 ENDDO
1593 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1594 !
1595 !--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
1596 !--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
1597 !--- WOULD HAVE ON THE STABILITY,
1598 !--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
1599 !--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
1600 !--- DESTABILIZATION.
1601 !
1602 !--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
1603 !
1604 DO K = 1, KM
1605 DO I = 1, IM
1606 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1607 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1608 !
1609 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1610 !
1611 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k)+EPSM1*QESO(I,k))
1612 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1613 val = 1.E-8
1614 QESO(I,k) = MAX(QESO(I,k), val )
1615 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
1616 ENDIF
1617 ENDDO
1618 ENDDO
1619 DO I = 1, IM
1620 IF(CNVFLG(I)) THEN
1621 XAA0(I) = 0.
1622 XPWAV(I) = 0.
1623 ENDIF
1624 ENDDO
1625 !
1626 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
1627 !
1628 ! DO I = 1, IM
1629 ! IF(CNVFLG(I)) THEN
1630 ! DLNSIG = LOG(PRSL(I,1)/PS(I))
1631 ! ZO(I,1) = TERR - DLNSIG * RD / G * TVO(I,1)
1632 ! ENDIF
1633 ! ENDDO
1634 ! DO K = 2, KM
1635 ! DO I = 1, IM
1636 ! IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1637 ! DLNSIG = LOG(PRSL(I,K) / PRSL(I,K-1))
1638 ! ZO(I,k) = ZO(I,k-1) - DLNSIG * RD / G
1639 ! & * .5 * (TVO(I,k) + TVO(I,k-1))
1640 ! ENDIF
1641 ! ENDDO
1642 ! ENDDO
1643 !
1644 !--- MOIST STATIC ENERGY
1645 !
1646 DO K = 1, KM1
1647 DO I = 1, IM
1648 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1649 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
1650 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
1651 !jfe ES = 10. * FPVS(TO(I,k+1))
1652 !
1653 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
1654 !
1655 PPRIME = PFLD(I,k+1) + EPSM1 * ES
1656 QS = EPS * ES / PPRIME
1657 DQSDP = - QS / PPRIME
1658 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
1659 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
1660 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
1661 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
1662 DQ = DQSDT * DT + DQSDP * DP
1663 TO(I,k) = TO(I,k+1) + DT
1664 QO(I,k) = QO(I,k+1) + DQ
1665 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
1666 ENDIF
1667 ENDDO
1668 ENDDO
1669 DO K = 1, KM1
1670 DO I = 1, IM
1671 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1672 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1673 !
1674 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1675 !
1676 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1 * QESO(I,k))
1677 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1678 val1 = 1.E-8
1679 QESO(I,k) = MAX(QESO(I,k), val1)
1680 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1681 val2 = 1.e-10
1682 QO(I,k) = max(QO(I,k), val2 )
1683 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
1684 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1685 & CP * TO(I,k) + HVAP * QO(I,k)
1686 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1687 & CP * TO(I,k) + HVAP * QESO(I,k)
1688 ENDIF
1689 ENDDO
1690 ENDDO
1691 DO I = 1, IM
1692 k = kmax(i)
1693 IF(CNVFLG(I)) THEN
1694 HEO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QO(I,k)
1695 HESO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QESO(I,k)
1696 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
1697 ENDIF
1698 ENDDO
1699 DO I = 1, IM
1700 IF(CNVFLG(I)) THEN
1701 INDX = KB(I)
1702 XHKB(I) = HEO(I,INDX)
1703 XQKB(I) = QO(I,INDX)
1704 HCKO(I,INDX) = XHKB(I)
1705 QCKO(I,INDX) = XQKB(I)
1706 ENDIF
1707 ENDDO
1708 !
1709 !
1710 !**************************** STATIC CONTROL
1711 !
1712 !
1713 !------- MOISTURE AND CLOUD WORK FUNCTIONS
1714 !
1715 DO K = 2, KM1
1716 DO I = 1, IM
1717 if (k .le. kmax(i)-1) then
1718 ! IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1719 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KTCON(I)) THEN
1720 FACTOR = ETA(I,k-1) / ETA(I,k)
1721 ONEMF = 1. - FACTOR
1722 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1723 & .5 * (HEO(I,k) + HEO(I,k+1))
1724 ENDIF
1725 ! IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1726 ! HEO(I,k) = HEO(I,k-1)
1727 ! ENDIF
1728 endif
1729 ENDDO
1730 ENDDO
1731 DO K = 2, KM1
1732 DO I = 1, IM
1733 if (k .le. kmax(i)-1) then
1734 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1735 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1736 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1737 XDBY = HCKO(I,k) - HESO(I,k)
1738 !cmr XDBY = MAX(XDBY,0.)
1739 val = 0.
1740 XDBY = MAX(XDBY,val)
1741 XQRCH = QESO(I,k) &
1742 & + GAMMA * XDBY / (HVAP * (1. + GAMMA))
1743 FACTOR = ETA(I,k-1) / ETA(I,k)
1744 ONEMF = 1. - FACTOR
1745 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1746 & .5 * (QO(I,k) + QO(I,k+1))
1747 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * XQRCH
1748 IF(DQ.GT.0.) THEN
1749 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1750 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1751 XAA0(I) = XAA0(I) - (ZO(I,k) - ZO(I,k-1)) * G * QLK
1752 XQC = QLK + XQRCH
1753 XPW = ETAH * C0 * DZ * QLK
1754 QCKO(I,k) = XQC
1755 XPWAV(I) = XPWAV(I) + XPW
1756 ENDIF
1757 ENDIF
1758 ! IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LT.KTCON(I)) THEN
1759 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1760 DZ1 = ZO(I,k) - ZO(I,k-1)
1761 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1762 RFACT = 1. + DELTA * CP * GAMMA &
1763 & * TO(I,k-1) / HVAP
1764 XDBY = HCKO(I,k-1) - HESO(I,k-1)
1765 XAA0(I) = XAA0(I) &
1766 & + DZ1 * (G / (CP * TO(I,k-1))) &
1767 & * XDBY / (1. + GAMMA) &
1768 & * RFACT
1769 val=0.
1770 XAA0(I)=XAA0(I)+ &
1771 & DZ1 * G * DELTA * &
1772 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1773 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1774 ENDIF
1775 endif
1776 ENDDO
1777 ENDDO
1778 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1779 !cccc PRINT *, ' XAA BEFORE DWNDRFT =', XAA0(I)
1780 !cccc ENDIF
1781 !
1782 !------- DOWNDRAFT CALCULATIONS
1783 !
1784 !
1785 !--- DOWNDRAFT MOISTURE PROPERTIES
1786 !
1787 DO I = 1, IM
1788 XPWEV(I) = 0.
1789 ENDDO
1790 DO I = 1, IM
1791 IF(DWNFLG2(I)) THEN
1792 JMN = JMIN(I)
1793 XHCD(I) = HEO(I,JMN)
1794 XQCD(I) = QO(I,JMN)
1795 QRCD(I,JMN) = QESO(I,JMN)
1796 ENDIF
1797 ENDDO
1798 DO K = KM1, 1, -1
1799 DO I = 1, IM
1800 if (k .le. kmax(i)-1) then
1801 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1802 DQ = QESO(I,k)
1803 DT = TO(I,k)
1804 GAMMA = EL2ORC * DQ / DT**2
1805 DH = XHCD(I) - HESO(I,k)
1806 QRCD(I,k)=DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1807 DETAD = ETAD(I,k+1) - ETAD(I,k)
1808 XPWD = ETAD(I,k+1) * QRCD(I,k+1) - &
1809 & ETAD(I,k) * QRCD(I,k)
1810 XPWD = XPWD - DETAD * &
1811 & .5 * (QRCD(I,k) + QRCD(I,k+1))
1812 XPWEV(I) = XPWEV(I) + XPWD
1813 ENDIF
1814 endif
1815 ENDDO
1816 ENDDO
1817 !
1818 DO I = 1, IM
1819 edtmax = edtmaxl
1820 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1821 IF(DWNFLG2(I)) THEN
1822 IF(XPWEV(I).GE.0.) THEN
1823 EDTX(I) = 0.
1824 ELSE
1825 EDTX(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
1826 EDTX(I) = MIN(EDTX(I),EDTMAX)
1827 ENDIF
1828 ELSE
1829 EDTX(I) = 0.
1830 ENDIF
1831 ENDDO
1832 !
1833 !
1834 !
1835 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1836 !
1837 !
1838 DO K = KM1, 1, -1
1839 DO I = 1, IM
1840 if (k .le. kmax(i)-1) then
1841 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1842 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1843 DHH=XHCD(I)
1844 DT= TO(I,k+1)
1845 DG= GAMMA
1846 DH= HESO(I,k+1)
1847 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1848 XAA0(I)=XAA0(I)+EDTX(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1849 & *(1.+DELTA*CP*DG*DT/HVAP)
1850 val=0.
1851 XAA0(I)=XAA0(I)+EDTX(I)* &
1852 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
1853 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1854 ENDIF
1855 endif
1856 ENDDO
1857 ENDDO
1858 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1859 !cccc PRINT *, ' XAA AFTER DWNDRFT =', XAA0(I)
1860 !cccc ENDIF
1861 !
1862 ! CALCULATE CRITICAL CLOUD WORK FUNCTION
1863 !
1864 DO I = 1, IM
1865 ACRT(I) = 0.
1866 IF(CNVFLG(I)) THEN
1867 ! IF(CNVFLG(I).AND.SLIMSK(I).NE.1.) THEN
1868 IF(PFLD(I,KTCON(I)).LT.PCRIT(15))THEN
1869 ACRT(I)=ACRIT(15)*(975.-PFLD(I,KTCON(I))) &
1870 & /(975.-PCRIT(15))
1871 ELSE IF(PFLD(I,KTCON(I)).GT.PCRIT(1))THEN
1872 ACRT(I)=ACRIT(1)
1873 ELSE
1874 !cmr K = IFIX((850. - PFLD(I,KTCON(I)))/50.) + 2
1875 K = int((850. - PFLD(I,KTCON(I)))/50.) + 2
1876 K = MIN(K,15)
1877 K = MAX(K,2)
1878 ACRT(I)=ACRIT(K)+(ACRIT(K-1)-ACRIT(K))* &
1879 & (PFLD(I,KTCON(I))-PCRIT(K))/(PCRIT(K-1)-PCRIT(K))
1880 ENDIF
1881 ! ELSE
1882 ! ACRT(I) = .5 * (PFLD(I,KBCON(I)) - PFLD(I,KTCON(I)))
1883 ENDIF
1884 ENDDO
1885 DO I = 1, IM
1886 ACRTFCT(I) = 1.
1887 IF(CNVFLG(I)) THEN
1888 if(SLIMSK(I).eq.1.) THEN
1889 w1 = w1l
1890 w2 = w2l
1891 w3 = w3l
1892 w4 = w4l
1893 else
1894 w1 = w1s
1895 w2 = w2s
1896 w3 = w3s
1897 w4 = w4s
1898 ENDIF
1899 !C IF(CNVFLG(I).AND.SLIMSK(I).EQ.1.) THEN
1900 ! ACRTFCT(I) = PDOT(I) / W3
1901 !
1902 ! modify critical cloud workfunction by cloud base vertical velocity
1903 !
1904 IF(PDOT(I).LE.W4) THEN
1905 ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
1906 ELSEIF(PDOT(I).GE.-W4) THEN
1907 ACRTFCT(I) = - (PDOT(I) + W4) / (W4 - W3)
1908 ELSE
1909 ACRTFCT(I) = 0.
1910 ENDIF
1911 !cmr ACRTFCT(I) = MAX(ACRTFCT(I),-1.)
1912 val1 = -1.
1913 ACRTFCT(I) = MAX(ACRTFCT(I),val1)
1914 !cmr ACRTFCT(I) = MIN(ACRTFCT(I),1.)
1915 val2 = 1.
1916 ACRTFCT(I) = MIN(ACRTFCT(I),val2)
1917 ACRTFCT(I) = 1. - ACRTFCT(I)
1918 !
1919 ! modify ACRTFCT(I) by colume mean rh if RHBAR(I) is greater than 80 percent
1920 !
1921 ! if(RHBAR(I).ge..8) THEN
1922 ! ACRTFCT(I) = ACRTFCT(I) * (.9 - min(RHBAR(I),.9)) * 10.
1923 ! ENDIF
1924 !
1925 ! modify adjustment time scale by cloud base vertical velocity
1926 !
1927 DTCONV(I) = DT2 + max((1800. - DT2),RZERO) * &
1928 & (PDOT(I) - W2) / (W1 - W2)
1929 ! DTCONV(I) = MAX(DTCONV(I), DT2)
1930 ! DTCONV(I) = 1800. * (PDOT(I) - w2) / (w1 - w2)
1931 DTCONV(I) = max(DTCONV(I),dtmin)
1932 DTCONV(I) = min(DTCONV(I),dtmax)
1933
1934 ENDIF
1935 ENDDO
1936 !
1937 !--- LARGE SCALE FORCING
1938 !
1939 DO I= 1, IM
1940 FLG(I) = CNVFLG(I)
1941 IF(CNVFLG(I)) THEN
1942 ! F = AA1(I) / DTCONV(I)
1943 FLD(I) = (AA1(I) - ACRT(I) * ACRTFCT(I)) / DTCONV(I)
1944 IF(FLD(I).LE.0.) FLG(I) = .FALSE.
1945 ENDIF
1946 CNVFLG(I) = FLG(I)
1947 IF(CNVFLG(I)) THEN
1948 ! XAA0(I) = MAX(XAA0(I),0.)
1949 XK(I) = (XAA0(I) - AA1(I)) / MBDT
1950 IF(XK(I).GE.0.) FLG(I) = .FALSE.
1951 ENDIF
1952 !
1953 !--- KERNEL, CLOUD BASE MASS FLUX
1954 !
1955 CNVFLG(I) = FLG(I)
1956 IF(CNVFLG(I)) THEN
1957 XMB(I) = -FLD(I) / XK(I)
1958 XMB(I) = MIN(XMB(I),XMBMAX(I))
1959 ENDIF
1960 ENDDO
1961 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1962 ! print *, ' RHBAR(I), ACRTFCT(I) =', RHBAR(I), ACRTFCT(I)
1963 ! PRINT *, ' A1, XA =', AA1(I), XAA0(I)
1964 ! PRINT *, ' XMB(I), ACRT =', XMB(I), ACRT
1965 ! ENDIF
1966 TOTFLG = .TRUE.
1967 DO I = 1, IM
1968 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1969 ENDDO
1970 IF(TOTFLG) RETURN
1971 !
1972 ! restore t0 and QO to t1 and q1 in case convection stops
1973 !
1974 do k = 1, km
1975 DO I = 1, IM
1976 if (k .le. kmax(i)) then
1977 TO(I,k) = T1(I,k)
1978 QO(I,k) = Q1(I,k)
1979 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
1980 !
1981 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
1982 !
1983 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
1984 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1985 val = 1.E-8
1986 QESO(I,k) = MAX(QESO(I,k), val )
1987 endif
1988 enddo
1989 enddo
1990 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1991 !
1992 !--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
1993 !--- MULTIPLIED BY THE MASS FLUX NECESSARY TO KEEP THE
1994 !--- EQUILIBRIUM WITH THE LARGER-SCALE.
1995 !
1996 DO I = 1, IM
1997 DELHBAR(I) = 0.
1998 DELQBAR(I) = 0.
1999 DELTBAR(I) = 0.
2000 QCOND(I) = 0.
2001 ENDDO
2002 DO K = 1, KM
2003 DO I = 1, IM
2004 if (k .le. kmax(i)) then
2005 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2006 AUP = 1.
2007 IF(K.Le.KB(I)) AUP = 0.
2008 ADW = 1.
2009 IF(K.GT.JMIN(I)) ADW = 0.
2010 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
2011 T1(I,k) = T1(I,k) + DELLAT * XMB(I) * DT2
2012 Q1(I,k) = Q1(I,k) + DELLAQ(I,k) * XMB(I) * DT2
2013 U1(I,k) = U1(I,k) + DELLAU(I,k) * XMB(I) * DT2
2014 V1(I,k) = V1(I,k) + DELLAV(I,k) * XMB(I) * DT2
2015 DP = 1000. * DEL(I,K)
2016 DELHBAR(I) = DELHBAR(I) + DELLAH(I,k)*XMB(I)*DP/G
2017 DELQBAR(I) = DELQBAR(I) + DELLAQ(I,k)*XMB(I)*DP/G
2018 DELTBAR(I) = DELTBAR(I) + DELLAT*XMB(I)*DP/G
2019 ENDIF
2020 endif
2021 ENDDO
2022 ENDDO
2023 DO K = 1, KM
2024 DO I = 1, IM
2025 if (k .le. kmax(i)) then
2026 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2027 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
2028 !
2029 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
2030 !
2031 QESO(I,k) = EPS * QESO(I,k)/(PFLD(I,k) + EPSM1*QESO(I,k))
2032 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
2033 val = 1.E-8
2034 QESO(I,k) = MAX(QESO(I,k), val )
2035 !
2036 ! cloud water
2037 !
2038 if(ncloud.gt.0.and.CNVFLG(I).and.k.eq.KTCON(I)) THEN
2039 tem = DELLAL(I) * XMB(I) * dt2
2040 tem1 = MAX(RZERO, MIN(RONE, (TCR-t1(I,K))*TCRF))
2041 if (QL(I,k,2) .gt. -999.0) then
2042 QL(I,k,1) = QL(I,k,1) + tem * tem1 ! Ice
2043 QL(I,k,2) = QL(I,k,2) + tem *(1.0-tem1) ! Water
2044 else
2045 tem2 = QL(I,k,1) + tem
2046 QL(I,k,1) = tem2 * tem1 ! Ice
2047 QL(I,k,2) = tem2 - QL(I,k,1) ! Water
2048 endif
2049 ! QL(I,k) = QL(I,k) + DELLAL(I) * XMB(I) * dt2
2050 dp = 1000. * del(i,k)
2051 DELLAL(I) = DELLAL(I) * XMB(I) * dp / g
2052 ENDIF
2053 ENDIF
2054 endif
2055 ENDDO
2056 ENDDO
2057 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2058 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2059 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2060 ! PRINT *, ' DELLBAR ='
2061 ! PRINT 6003, HVAP*DELLbar
2062 ! PRINT *, ' DELLAQ ='
2063 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2064 ! PRINT *, ' DELLAT ='
2065 ! PRINT 6003, (DELLAH(i,k)*XMB(I)-HVAP*DELLAQ(I,k)*XMB(I), &
2066 ! & K=1,KMAX)
2067 ! ENDIF
2068 DO I = 1, IM
2069 RNTOT(I) = 0.
2070 DELQEV(I) = 0.
2071 DELQ2(I) = 0.
2072 FLG(I) = CNVFLG(I)
2073 ENDDO
2074 DO K = KM, 1, -1
2075 DO I = 1, IM
2076 if (k .le. kmax(i)) then
2077 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2078 AUP = 1.
2079 IF(K.Le.KB(I)) AUP = 0.
2080 ADW = 1.
2081 IF(K.GT.JMIN(I)) ADW = 0.
2082 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2083 RNTOT(I) = RNTOT(I) + rain * XMB(I) * .001 * dt2
2084 ENDIF
2085 endif
2086 ENDDO
2087 ENDDO
2088 DO K = KM, 1, -1
2089 DO I = 1, IM
2090 if (k .le. kmax(i)) then
2091 DELTV(I) = 0.
2092 DELQ(I) = 0.
2093 QEVAP(I) = 0.
2094 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2095 AUP = 1.
2096 IF(K.Le.KB(I)) AUP = 0.
2097 ADW = 1.
2098 IF(K.GT.JMIN(I)) ADW = 0.
2099 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2100 RN(I) = RN(I) + rain * XMB(I) * .001 * dt2
2101 ENDIF
2102 IF(FLG(I).AND.K.LE.KTCON(I)) THEN
2103 evef = EDT(I) * evfact
2104 if(SLIMSK(I).eq.1.) evef=EDT(I) * evfactl
2105 ! if(SLIMSK(I).eq.1.) evef=.07
2106 ! if(SLIMSK(I).ne.1.) evef = 0.
2107 QCOND(I) = EVEF * (Q1(I,k) - QESO(I,k)) &
2108 & / (1. + EL2ORC * QESO(I,k) / T1(I,k)**2)
2109 DP = 1000. * DEL(I,K)
2110 IF(RN(I).GT.0..AND.QCOND(I).LT.0.) THEN
2111 QEVAP(I) = -QCOND(I) * (1.-EXP(-.32*SQRT(DT2*RN(I))))
2112 QEVAP(I) = MIN(QEVAP(I), RN(I)*1000.*G/DP)
2113 DELQ2(I) = DELQEV(I) + .001 * QEVAP(I) * dp / g
2114 ENDIF
2115 if(RN(I).gt.0..and.QCOND(I).LT.0..and. &
2116 & DELQ2(I).gt.RNTOT(I)) THEN
2117 QEVAP(I) = 1000.* g * (RNTOT(I) - DELQEV(I)) / dp
2118 FLG(I) = .false.
2119 ENDIF
2120 IF(RN(I).GT.0..AND.QEVAP(I).gt.0.) THEN
2121 Q1(I,k) = Q1(I,k) + QEVAP(I)
2122 T1(I,k) = T1(I,k) - ELOCP * QEVAP(I)
2123 RN(I) = RN(I) - .001 * QEVAP(I) * DP / G
2124 DELTV(I) = - ELOCP*QEVAP(I)/DT2
2125 DELQ(I) = + QEVAP(I)/DT2
2126 DELQEV(I) = DELQEV(I) + .001*dp*QEVAP(I)/g
2127 ENDIF
2128 DELLAQ(I,k) = DELLAQ(I,k) + DELQ(I) / XMB(I)
2129 DELQBAR(I) = DELQBAR(I) + DELQ(I)*DP/G
2130 DELTBAR(I) = DELTBAR(I) + DELTV(I)*DP/G
2131 ENDIF
2132 endif
2133 ENDDO
2134 ENDDO
2135 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2136 ! PRINT *, ' DELLAH ='
2137 ! PRINT 6003, (DELLAH(k)*XMB(I),K=1,KMAX)
2138 ! PRINT *, ' DELLAQ ='
2139 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2140 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2141 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2142 ! PRINT *, ' PRECIP =', HVAP*RN(I)*1000./DT2
2143 !CCCC PRINT *, ' DELLBAR ='
2144 !CCCC PRINT *, HVAP*DELLbar
2145 ! ENDIF
2146 !
2147 ! PRECIPITATION RATE CONVERTED TO ACTUAL PRECIP
2148 ! IN UNIT OF M INSTEAD OF KG
2149 !
2150 DO I = 1, IM
2151 IF(CNVFLG(I)) THEN
2152 !
2153 ! IN THE EVENT OF UPPER LEVEL RAIN EVAPORATION AND LOWER LEVEL DOWNDRAF
2154 ! MOISTENING, RN CAN BECOME NEGATIVE, IN THIS CASE, WE BACK OUT OF TH
2155 ! HEATING AND THE MOISTENING
2156 !
2157 if(RN(I).lt.0..and..not.FLG(I)) RN(I) = 0.
2158 IF(RN(I).LE.0.) THEN
2159 RN(I) = 0.
2160 ELSE
2161 KTOP(I) = KTCON(I)
2162 KBOT(I) = KBCON(I)
2163 KUO(I) = 1
2164 CLDWRK(I) = AA1(I)
2165 ENDIF
2166 ENDIF
2167 ENDDO
2168 DO K = 1, KM
2169 DO I = 1, IM
2170 if (k .le. kmax(i)) then
2171 IF(CNVFLG(I).AND.RN(I).LE.0.) THEN
2172 T1(I,k) = TO(I,k)
2173 Q1(I,k) = QO(I,k)
2174 ENDIF
2175 endif
2176 ENDDO
2177 ENDDO
2178 !!
2179 RETURN
2180 END SUBROUTINE SASCNV
2181
2182 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2183
2184 SUBROUTINE SHALCV(IM,IX,KM,DT,DEL,PRSI,PRSL,PRSLK,KUO,Q,T,DPSHC)
2185 !
2186 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2187 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
2188 &, RD => con_RD
2189
2190 implicit none
2191 !
2192 ! include 'constant.h'
2193 !
2194 integer IM, IX, KM, KUO(IM)
2195 real(kind=kind_phys) DEL(IX,KM), PRSI(IX,KM+1), PRSL(IX,KM), &
2196 & PRSLK(IX,KM), &
2197 & Q(IX,KM), T(IX,KM), DT, DPSHC
2198 !
2199 ! Locals
2200 !
2201 real(kind=kind_phys) ck, cpdt, dmse, dsdz1, dsdz2, &
2202 & dsig, dtodsl, dtodsu, eldq, g, &
2203 & gocp, rtdls
2204 !
2205 integer k,k1,k2,kliftl,kliftu,kt,N2,I,iku,ik1,ik,ii
2206 integer INDEX2(IM), KLCL(IM), KBOT(IM), KTOP(IM),kk &
2207 &, KTOPM(IM)
2208 !!
2209 ! PHYSICAL PARAMETERS
2210 PARAMETER(G=GRAV, GOCP=G/CP)
2211 ! BOUNDS OF PARCEL ORIGIN
2212 PARAMETER(KLIFTL=2,KLIFTU=2)
2213 LOGICAL LSHC(IM)
2214 real(kind=kind_phys) Q2(IM*KM), T2(IM*KM), &
2215 & PRSL2(IM*KM), PRSLK2(IM*KM), &
2216 & AL(IM*(KM-1)), AD(IM*KM), AU(IM*(KM-1))
2217 !-----------------------------------------------------------------------
2218 ! COMPRESS FIELDS TO POINTS WITH NO DEEP CONVECTION
2219 ! AND MOIST STATIC INSTABILITY.
2220 DO I=1,IM
2221 LSHC(I)=.FALSE.
2222 ENDDO
2223 DO K=1,KM-1
2224 DO I=1,IM
2225 IF(KUO(I).EQ.0) THEN
2226 ELDQ = HVAP*(Q(I,K)-Q(I,K+1))
2227 CPDT = CP*(T(I,K)-T(I,K+1))
2228 RTDLS = (PRSL(I,K)-PRSL(I,K+1)) / &
2229 & PRSI(I,K+1)*RD*0.5*(T(I,K)+T(I,K+1))
2230 DMSE = ELDQ+CPDT-RTDLS
2231 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2232 ENDIF
2233 ENDDO
2234 ENDDO
2235 N2 = 0
2236 DO I=1,IM
2237 IF(LSHC(I)) THEN
2238 N2 = N2 + 1
2239 INDEX2(N2) = I
2240 ENDIF
2241 ENDDO
2242 IF(N2.EQ.0) RETURN
2243 DO K=1,KM
2244 KK = (K-1)*N2
2245 DO I=1,N2
2246 IK = KK + I
2247 ii = index2(i)
2248 Q2(IK) = Q(II,K)
2249 T2(IK) = T(II,K)
2250 PRSL2(IK) = PRSL(II,K)
2251 PRSLK2(IK) = PRSLK(II,K)
2252 ENDDO
2253 ENDDO
2254 do i=1,N2
2255 ktopm(i) = KM
2256 enddo
2257 do k=2,KM
2258 do i=1,N2
2259 ii = index2(i)
2260 if (prsi(ii,1)-prsi(ii,k) .le. dpshc) ktopm(i) = k
2261 enddo
2262 enddo
2263
2264 !-----------------------------------------------------------------------
2265 ! COMPUTE MOIST ADIABAT AND DETERMINE LIMITS OF SHALLOW CONVECTION.
2266 ! CHECK FOR MOIST STATIC INSTABILITY AGAIN WITHIN CLOUD.
2267 CALL MSTADBT3(N2,KM-1,KLIFTL,KLIFTU,PRSL2,PRSLK2,T2,Q2, &
2268 & KLCL,KBOT,KTOP,AL,AU)
2269 DO I=1,N2
2270 KBOT(I) = min(KLCL(I)-1, ktopm(i)-1)
2271 KTOP(I) = min(KTOP(I)+1, ktopm(i))
2272 LSHC(I) = .FALSE.
2273 ENDDO
2274 DO K=1,KM-1
2275 KK = (K-1)*N2
2276 DO I=1,N2
2277 IF(K.GE.KBOT(I).AND.K.LT.KTOP(I)) THEN
2278 IK = KK + I
2279 IKU = IK + N2
2280 ELDQ = HVAP * (Q2(IK)-Q2(IKU))
2281 CPDT = CP * (T2(IK)-T2(IKU))
2282 RTDLS = (PRSL2(IK)-PRSL2(IKU)) / &
2283 & PRSI(index2(i),K+1)*RD*0.5*(T2(IK)+T2(IKU))
2284 DMSE = ELDQ + CPDT - RTDLS
2285 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2286 AU(IK) = G/RTDLS
2287 ENDIF
2288 ENDDO
2289 ENDDO
2290 K1=KM+1
2291 K2=0
2292 DO I=1,N2
2293 IF(.NOT.LSHC(I)) THEN
2294 KBOT(I) = KM+1
2295 KTOP(I) = 0
2296 ENDIF
2297 K1 = MIN(K1,KBOT(I))
2298 K2 = MAX(K2,KTOP(I))
2299 ENDDO
2300 KT = K2-K1+1
2301 IF(KT.LT.2) RETURN
2302 !-----------------------------------------------------------------------
2303 ! SET EDDY VISCOSITY COEFFICIENT CKU AT SIGMA INTERFACES.
2304 ! COMPUTE DIAGONALS AND RHS FOR TRIDIAGONAL MATRIX SOLVER.
2305 ! EXPAND FINAL FIELDS.
2306 KK = (K1-1) * N2
2307 DO I=1,N2
2308 IK = KK + I
2309 AD(IK) = 1.
2310 ENDDO
2311 !
2312 ! DTODSU=DT/DEL(K1)
2313 DO K=K1,K2-1
2314 ! DTODSL=DTODSU
2315 ! DTODSU= DT/DEL(K+1)
2316 ! DSIG=SL(K)-SL(K+1)
2317 KK = (K-1) * N2
2318 DO I=1,N2
2319 ii = index2(i)
2320 DTODSL = DT/DEL(II,K)
2321 DTODSU = DT/DEL(II,K+1)
2322 DSIG = PRSL(II,K) - PRSL(II,K+1)
2323 IK = KK + I
2324 IKU = IK + N2
2325 IF(K.EQ.KBOT(I)) THEN
2326 CK=1.5
2327 ELSEIF(K.EQ.KTOP(I)-1) THEN
2328 CK=1.
2329 ELSEIF(K.EQ.KTOP(I)-2) THEN
2330 CK=3.
2331 ELSEIF(K.GT.KBOT(I).AND.K.LT.KTOP(I)-2) THEN
2332 CK=5.
2333 ELSE
2334 CK=0.
2335 ENDIF
2336 DSDZ1 = CK*DSIG*AU(IK)*GOCP
2337 DSDZ2 = CK*DSIG*AU(IK)*AU(IK)
2338 AU(IK) = -DTODSL*DSDZ2
2339 AL(IK) = -DTODSU*DSDZ2
2340 AD(IK) = AD(IK)-AU(IK)
2341 AD(IKU) = 1.-AL(IK)
2342 T2(IK) = T2(IK)+DTODSL*DSDZ1
2343 T2(IKU) = T2(IKU)-DTODSU*DSDZ1
2344 ENDDO
2345 ENDDO
2346 IK1=(K1-1)*N2+1
2347 CALL TRIDI2T3(N2,KT,AL(IK1),AD(IK1),AU(IK1),Q2(IK1),T2(IK1), &
2348 & AU(IK1),Q2(IK1),T2(IK1))
2349 DO K=K1,K2
2350 KK = (K-1)*N2
2351 DO I=1,N2
2352 IK = KK + I
2353 Q(INDEX2(I),K) = Q2(IK)
2354 T(INDEX2(I),K) = T2(IK)
2355 ENDDO
2356 ENDDO
2357 !-----------------------------------------------------------------------
2358 RETURN
2359 END SUBROUTINE SHALCV
2360 !-----------------------------------------------------------------------
2361 SUBROUTINE TRIDI2T3(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
2362 !yt INCLUDE DBTRIDI2;
2363 !!
2364 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2365 implicit none
2366 integer k,n,l,i
2367 real(kind=kind_phys) fk
2368 !!
2369 real(kind=kind_phys) &
2370 & CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
2371 & AU(L,N-1),A1(L,N),A2(L,N)
2372 !-----------------------------------------------------------------------
2373 DO I=1,L
2374 FK=1./CM(I,1)
2375 AU(I,1)=FK*CU(I,1)
2376 A1(I,1)=FK*R1(I,1)
2377 A2(I,1)=FK*R2(I,1)
2378 ENDDO
2379 DO K=2,N-1
2380 DO I=1,L
2381 FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
2382 AU(I,K)=FK*CU(I,K)
2383 A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
2384 A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
2385 ENDDO
2386 ENDDO
2387 DO I=1,L
2388 FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
2389 A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
2390 A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
2391 ENDDO
2392 DO K=N-1,1,-1
2393 DO I=1,L
2394 A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
2395 A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
2396 ENDDO
2397 ENDDO
2398 !-----------------------------------------------------------------------
2399 RETURN
2400 END SUBROUTINE TRIDI2T3
2401 !-----------------------------------------------------------------------
2402
2403 SUBROUTINE MSTADBT3(IM,KM,K1,K2,PRSL,PRSLK,TENV,QENV, &
2404 & KLCL,KBOT,KTOP,TCLD,QCLD)
2405 !yt INCLUDE DBMSTADB;
2406 !!
2407 USE MODULE_GFS_MACHINE, ONLY : kind_phys
2408 USE MODULE_GFS_FUNCPHYS, ONLY : FTDP, FTHE, FTLCL, STMA
2409 USE MODULE_GFS_PHYSCONS, EPS => con_eps, EPSM1 => con_epsm1, FV => con_FVirt
2410
2411 implicit none
2412 !!
2413 ! include 'constant.h'
2414 !!
2415 integer k,k1,k2,km,i,im
2416 real(kind=kind_phys) pv,qma,slklcl,tdpd,thelcl,tlcl
2417 real(kind=kind_phys) tma,tvcld,tvenv
2418 !!
2419 real(kind=kind_phys) PRSL(IM,KM), PRSLK(IM,KM), TENV(IM,KM), &
2420 & QENV(IM,KM), TCLD(IM,KM), QCLD(IM,KM)
2421 INTEGER KLCL(IM), KBOT(IM), KTOP(IM)
2422 ! LOCAL ARRAYS
2423 real(kind=kind_phys) SLKMA(IM), THEMA(IM)
2424 !-----------------------------------------------------------------------
2425 ! DETERMINE WARMEST POTENTIAL WET-BULB TEMPERATURE BETWEEN K1 AND K2.
2426 ! COMPUTE ITS LIFTING CONDENSATION LEVEL.
2427 !
2428 DO I=1,IM
2429 SLKMA(I) = 0.
2430 THEMA(I) = 0.
2431 ENDDO
2432 DO K=K1,K2
2433 DO I=1,IM
2434 PV = 1000.0 * PRSL(I,K)*QENV(I,K)/(EPS-EPSM1*QENV(I,K))
2435 TDPD = TENV(I,K)-FTDP(PV)
2436 IF(TDPD.GT.0.) THEN
2437 TLCL = FTLCL(TENV(I,K),TDPD)
2438 SLKLCL = PRSLK(I,K)*TLCL/TENV(I,K)
2439 ELSE
2440 TLCL = TENV(I,K)
2441 SLKLCL = PRSLK(I,K)
2442 ENDIF
2443 THELCL=FTHE(TLCL,SLKLCL)
2444 IF(THELCL.GT.THEMA(I)) THEN
2445 SLKMA(I) = SLKLCL
2446 THEMA(I) = THELCL
2447 ENDIF
2448 ENDDO
2449 ENDDO
2450 !-----------------------------------------------------------------------
2451 ! SET CLOUD TEMPERATURES AND HUMIDITIES WHEREVER THE PARCEL LIFTED UP
2452 ! THE MOIST ADIABAT IS BUOYANT WITH RESPECT TO THE ENVIRONMENT.
2453 DO I=1,IM
2454 KLCL(I)=KM+1
2455 KBOT(I)=KM+1
2456 KTOP(I)=0
2457 ENDDO
2458 DO K=1,KM
2459 DO I=1,IM
2460 TCLD(I,K)=0.
2461 QCLD(I,K)=0.
2462 ENDDO
2463 ENDDO
2464 DO K=K1,KM
2465 DO I=1,IM
2466 IF(PRSLK(I,K).LE.SLKMA(I)) THEN
2467 KLCL(I)=MIN(KLCL(I),K)
2468 CALL STMA(THEMA(I),PRSLK(I,K),TMA,QMA)
2469 ! TMA=FTMA(THEMA(I),PRSLK(I,K),QMA)
2470 TVCLD=TMA*(1.+FV*QMA)
2471 TVENV=TENV(I,K)*(1.+FV*QENV(I,K))
2472 IF(TVCLD.GT.TVENV) THEN
2473 KBOT(I)=MIN(KBOT(I),K)
2474 KTOP(I)=MAX(KTOP(I),K)
2475 TCLD(I,K)=TMA-TENV(I,K)
2476 QCLD(I,K)=QMA-QENV(I,K)
2477 ENDIF
2478 ENDIF
2479 ENDDO
2480 ENDDO
2481 !-----------------------------------------------------------------------
2482 RETURN
2483 END SUBROUTINE MSTADBT3
2484
2485 END MODULE module_cu_sas