| blob | d92bf2c4674af7b1c7895d030441a19756759596 |
1 !!
2 MODULE module_cu_sas
4 CONTAINS
6 !-----------------------------------------------------------------
7 SUBROUTINE CU_SAS(DT,ITIMESTEP,STEPCU, &
8 RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
9 RUCUTEN,RVCUTEN, &
10 RAINCV,PRATEC,HTOP,HBOT, &
11 U3D,V3D,W,T3D,QV3D,QC3D,QI3D,PI3D,RHO3D, &
12 DZ8W,PCPS,P8W,XLAND,CU_ACT_FLAG, &
13 P_QC, &
14 MOMMIX, & ! gopal's doing
15 PGCON,sas_mass_flux, &
16 pert_sas, ens_random_seed, ens_sasamp, &
17 shalconv,shal_pgcon, &
18 HPBL2D,EVAP2D,HEAT2D, & !Kwon for shallow convection
19 P_QI,P_FIRST_SCALAR, &
20 ids,ide, jds,jde, kds,kde, &
21 ims,ime, jms,jme, kms,kme, &
22 its,ite, jts,jte, kts,kte )
24 !-------------------------------------------------------------------
25 USE MODULE_GFS_MACHINE , ONLY : kind_phys
26 USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
27 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
28 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
29 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
30 &, EPS => con_eps, EPSM1 => con_epsm1 &
31 &, ROVCP => con_rocp, RD => con_rd
32 !-------------------------------------------------------------------
33 IMPLICIT NONE
34 !-------------------------------------------------------------------
35 !-- U3D 3D u-velocity interpolated to theta points (m/s)
36 !-- V3D 3D v-velocity interpolated to theta points (m/s)
37 !-- TH3D 3D potential temperature (K)
38 !-- T3D temperature (K)
39 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
40 !-- QC3D 3D cloud mixing ratio (Kg/Kg)
41 !-- QI3D 3D ice mixing ratio (Kg/Kg)
42 !-- P8w 3D pressure at full levels (Pa)
43 !-- Pcps 3D pressure (Pa)
44 !-- PI3D 3D exner function (dimensionless)
45 !-- rr3D 3D dry air density (kg/m^3)
46 !-- RUBLTEN U tendency due to
47 ! PBL parameterization (m/s^2)
48 !-- RVBLTEN V tendency due to
49 ! PBL parameterization (m/s^2)
50 !-- RTHBLTEN Theta tendency due to
51 ! PBL parameterization (K/s)
52 !-- RQVBLTEN Qv tendency due to
53 ! PBL parameterization (kg/kg/s)
54 !-- RQCBLTEN Qc tendency due to
55 ! PBL parameterization (kg/kg/s)
56 !-- RQIBLTEN Qi tendency due to
57 ! PBL parameterization (kg/kg/s)
58 !
59 !-- MOMMIX MOMENTUM MIXING COEFFICIENT (can be set in the namelist)
60 !-- RUCUTEN U tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
61 !-- RVCUTEN V tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
62 !
63 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
64 !-- GRAV acceleration due to gravity (m/s^2)
65 !-- ROVCP R/CP
66 !-- RD gas constant for dry air (J/kg/K)
67 !-- ROVG R/G
68 !-- P_QI species index for cloud ice
69 !-- dz8w dz between full levels (m)
70 !-- z height above sea level (m)
71 !-- PSFC pressure at the surface (Pa)
72 !-- UST u* in similarity theory (m/s)
73 !-- PBL PBL height (m)
74 !-- PSIM similarity stability function for momentum
75 !-- PSIH similarity stability function for heat
76 !-- HFX upward heat flux at the surface (W/m^2)
77 !-- QFX upward moisture flux at the surface (kg/m^2/s)
78 !-- TSK surface temperature (K)
79 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
80 !-- WSPD wind speed at lowest model level (m/s)
81 !-- BR bulk Richardson number in surface layer
82 !-- DT time step (s)
83 !-- rvovrd R_v divided by R_d (dimensionless)
84 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
85 !-- KARMAN Von Karman constant
86 !-- ids start index for i in domain
87 !-- ide end index for i in domain
88 !-- jds start index for j in domain
89 !-- jde end index for j in domain
90 !-- kds start index for k in domain
91 !-- kde end index for k in domain
92 !-- ims start index for i in memory
93 !-- ime end index for i in memory
94 !-- jms start index for j in memory
95 !-- jme end index for j in memory
96 !-- kms start index for k in memory
97 !-- kme end index for k in memory
98 !-- its start index for i in tile
99 !-- ite end index for i in tile
100 !-- jts start index for j in tile
101 !-- jte end index for j in tile
102 !-- kts start index for k in tile
103 !-- kte end index for k in tile
104 !-------------------------------------------------------------------
106 INTEGER :: ICLDCK
108 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
109 ims,ime, jms,jme, kms,kme, &
110 its,ite, jts,jte, kts,kte, &
111 ITIMESTEP, & !NSTD
112 P_FIRST_SCALAR, &
113 P_QC, &
114 P_QI, &
115 STEPCU
117 REAL, INTENT(IN) :: &
118 DT
120 REAL, OPTIONAL, INTENT(IN) :: PGCON,sas_mass_flux,shal_pgcon
121 INTEGER, OPTIONAL, INTENT(IN) :: shalconv
122 REAL(kind=kind_phys) :: PGCON_USE,SHAL_PGCON_USE,massf
123 logical,intent(in) :: pert_sas
124 integer,intent(in) :: ens_random_seed
125 real,intent(in) :: ens_sasamp
126 INTEGER :: shalconv_use
127 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
128 RQCCUTEN, &
129 RQICUTEN, &
130 RQVCUTEN, &
131 RTHCUTEN
132 REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(INOUT) :: &
133 RUCUTEN, &
134 RVCUTEN
135 REAL, OPTIONAL, INTENT(IN) :: MOMMIX
137 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
138 INTENT(IN) :: HPBL2D,EVAP2D,HEAT2D !Kwon for sha
140 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
141 XLAND
143 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
144 RAINCV, PRATEC
146 REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
147 HBOT, &
148 HTOP
150 LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) :: &
151 CU_ACT_FLAG
154 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
155 DZ8W, &
156 P8w, &
157 Pcps, &
158 PI3D, &
159 QC3D, &
160 QI3D, &
161 QV3D, &
162 RHO3D, &
163 T3D, &
164 U3D, &
165 V3D, &
166 W
168 !--------------------------- LOCAL VARS ------------------------------
170 REAL, DIMENSION(ims:ime, jms:jme) :: &
171 PSFC
173 REAL, DIMENSION(its:ite, jts:jte) :: &
174 RAINCV1, PRATEC1
175 REAL, DIMENSION(its:ite, jts:jte) :: &
176 RAINCV2, PRATEC2
178 REAL (kind=kind_phys) :: &
179 DELT, &
180 DPSHC, &
181 RDELT, &
182 RSEED
184 REAL (kind=kind_phys), DIMENSION(its:ite) :: &
185 CLDWRK, &
186 PS, &
187 RCS, &
188 RN, &
189 SLIMSK, &
190 HPBL,EVAP,HEAT !Kwon for shallow convection
192 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
193 PRSI
195 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
196 DEL, &
197 DOT, &
198 PHIL, &
199 PRSL, &
200 PRSLK, &
201 Q1, &
202 T1, &
203 U1, &
204 V1, &
205 ZI, &
206 ZL
208 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) :: &
209 QL
211 INTEGER, DIMENSION(its:ite) :: &
212 KBOT, &
213 KTOP, &
214 KCNV
216 INTEGER :: &
217 I, &
218 IGPVS, &
219 IM, &
220 J, &
221 JCAP, &
222 K, &
223 KM, &
224 KP, &
225 KX, &
226 NCLOUD
228 DATA IGPVS/0/
230 !-----------------------------------------------------------------------
231 !
233 if(present(shalconv)) then
234 shalconv_use=shalconv
235 else
236 #if (NMM_CORE==1)
237 shalconv_use=0
238 #else
239 #if (EM_CORE==1)
240 shalconv_use=1
241 #else
242 shalconv_use=0
243 #endif
244 #endif
245 endif
247 if(present(pgcon)) then
248 pgcon_use = pgcon
249 else
250 ! pgcon_use = 0.7 ! Gregory et al. (1997, QJRMS)
251 pgcon_use = 0.55 ! Zhang & Wu (2003,JAS), used in GFS (25km res spectral)
252 ! pgcon_use = 0.2 ! HWRF, for model tuning purposes
253 ! pgcon_use = 0.3 ! GFDL, or so I am told
255 ! For those attempting to tune pgcon:
257 ! The value of 0.55 comes from an observational study of
258 ! synoptic-scale deep convection and 0.7 came from an
259 ! incorrect fit to the same data. That value is likely
260 ! correct for deep convection at gridscales near that of GFS,
261 ! but is questionable in shallow convection, or for scales
262 ! much finer than synoptic scales.
264 ! Then again, the assumptions of SAS break down when the
265 ! gridscale is near the convection scale anyway. In a large
266 ! storm such as a hurricane, there is often no environment to
267 ! detrain into since adjancent gridsquares are also undergoing
268 ! active convection. Each gridsquare will no longer have many
269 ! updrafts and downdrafts. At sub-convective timescales, you
270 ! will find unstable columns for many (say, 5 second length)
271 ! timesteps in a real atmosphere during a convection cell's
272 ! lifetime, so forcing it to be neutrally stable is unphysical.
274 ! Hence, in scales near the convection scale (cells have
275 ! ~0.5-4km diameter in hurricanes), this parameter is more of a
276 ! tuning parameter to get a scheme that is inappropriate for
277 ! that resolution to do a reasonable job.
279 ! Your mileage might vary.
281 ! - Sam Trahan
282 endif
284 if(present(sas_mass_flux)) then
285 massf=sas_mass_flux
286 ! Use this to reduce the fluxes added by SAS to prevent
287 ! computational instability as a result of large fluxes.
288 else
289 massf=9e9 ! large number to disable check
290 endif
292 if(present(shal_pgcon)) then
293 if(shal_pgcon>=0) then
294 shal_pgcon_use = shal_pgcon
295 else
296 ! shal_pgcon<0 means use deep pgcon
297 shal_pgcon_use = pgcon_use
298 endif
299 else
300 ! Default: Same as deep convection pgcon
301 shal_pgcon_use = pgcon_use
302 ! Read the warning above though. It may be advisable for
303 ! these to be different.
304 endif
306 DO J=JTS,JTE
307 DO I=ITS,ITE
308 CU_ACT_FLAG(I,J)=.TRUE.
309 ENDDO
310 ENDDO
312 IM=ITE-ITS+1
313 KX=KTE-KTS+1
314 JCAP=126
315 DPSHC=30_kind_phys
316 DELT=DT*STEPCU
317 RDELT=1./DELT
318 NCLOUD=1
321 DO J=jms,jme
322 DO I=ims,ime
323 PSFC(i,j)=P8w(i,kms,j)
324 ENDDO
325 ENDDO
327 if(igpvs.eq.0) CALL GFUNCPHYS
328 igpvs=1
330 !------------- J LOOP (OUTER) --------------------------------------------------
332 big_outer_j_loop: DO J=jts,jte
334 ! --------------- compute zi and zl -----------------------------------------
335 DO i=its,ite
336 ZI(I,KTS)=0.0
337 ENDDO
339 DO k=kts+1,kte
340 KM=K-1
341 DO i=its,ite
342 ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
343 ENDDO
344 ENDDO
346 DO k=kts+1,kte
347 KM=K-1
348 DO i=its,ite
349 ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
350 ENDDO
351 ENDDO
353 DO i=its,ite
354 ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
355 ENDDO
357 ! --------------- end compute zi and zl -------------------------------------
359 DO i=its,ite
360 PS(i)=PSFC(i,j)*.001
361 RCS(i)=1.
362 SLIMSK(i)=ABS(XLAND(i,j)-2.)
363 ENDDO
365 #if (NMM_CORE == 1)
366 if(shalconv_use==1) then
367 DO i=its,ite
368 HPBL(I) = HPBL2D(I,J) !kwon for shallow convection
369 EVAP(I) = EVAP2D(I,J) !kwon for shallow convection
370 HEAT(I) = HEAT2D(I,J) !kwon for shallow convection
371 ENDDO
372 endif
373 #endif
375 DO i=its,ite
376 PRSI(i,kts)=PS(i)
377 ENDDO
379 DO k=kts,kte
380 kp=k+1
381 DO i=its,ite
382 PRSL(I,K)=Pcps(i,k,j)*.001
383 PHIL(I,K)=ZL(I,K)*GRAV
384 DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
385 ENDDO
386 ENDDO
388 DO k=kts,kte
389 DO i=its,ite
390 DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
391 U1(i,k)=U3D(i,k,j)
392 V1(i,k)=V3D(i,k,j)
393 Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
394 T1(i,k)=T3D(i,k,j)
395 QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
396 QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
397 PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
398 ENDDO
399 ENDDO
401 DO k=kts+1,kte+1
402 km=k-1
403 DO i=its,ite
404 PRSI(i,k)=PRSI(i,km)-del(i,km)
405 ENDDO
406 ENDDO
408 CALL SASCNVN(IM,IM,KX,JCAP,DELT,DEL,PRSL,PS,PHIL, &
409 QL,Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT, &
410 KTOP,KCNV,SLIMSK,DOT,NCLOUD,PGCON_USE,massf, &
411 pert_sas, ens_random_seed, ens_sasamp)
412 !
413 do i=its,ite
414 RAINCV1(I,J)=RN(I)*1000./STEPCU
415 PRATEC1(I,J)=RN(I)*1000./(STEPCU * DT)
416 enddo
417 !
418 do i=its,ite
419 RAINCV2(I,J)=0.
420 PRATEC2(I,J)=0.
421 enddo
422 !
424 if_shallow_conv: if(shalconv_use==1) then
425 #if (NMM_CORE == 1)
426 ! NMM calls the new shallow convection developed by J Han
427 ! (Added to WRF by Y.Kwon)
428 call shalcnv(im,im,kx,jcap,delt,del,prsl,ps,phil,ql, &
429 & q1,t1,u1,v1,rcs,rn,kbot,ktop,kcnv,slimsk, &
430 & dot,ncloud,hpbl,heat,evap,shal_pgcon_use)
431 !
432 DO I=ITS,ITE
433 RAINCV2(I,J)=RN(I)*1000./STEPCU
434 PRATEC2(I,J)=RN(I)*1000./(STEPCU * DT)
435 ENDDO
436 !
437 #else
438 #if (EM_CORE == 1)
439 ! NOTE: ARW should be able to call the new shalcnv here, but
440 ! they need to add the three new variables, so I'm leaving the
441 ! old shallow convection call here - Sam Trahan
442 CALL OLD_ARW_SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KCNV,Q1,T1,DPSHC)
443 #else
444 ! Shallow convection is untested for other cores.
445 #endif
446 #endif
447 endif if_shallow_conv
449 DO I=ITS,ITE
450 RAINCV(I,J)= RAINCV1(I,J) + RAINCV2(I,J)
451 PRATEC(I,J)= PRATEC1(I,J) + PRATEC2(I,J)
452 HBOT(I,J)=KBOT(I)
453 HTOP(I,J)=KTOP(I)
454 ENDDO
456 DO K=KTS,KTE
457 DO I=ITS,ITE
458 RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
459 RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
460 ENDDO
461 ENDDO
463 !===============================================================================
464 ! ADD MOMENTUM MIXING TERM AS TENDENCIES. This is gopal's doing for SAS
465 ! MOMMIX is the reduction factor set to 0.7 by default. Because NMM has
466 ! divergence damping term, a reducion factor for cumulum mixing may be
467 ! required otherwise storms were too weak.
468 !===============================================================================
469 !
470 #if (NMM_CORE == 1)
471 DO K=KTS,KTE
472 DO I=ITS,ITE
473 ! RUCUTEN(I,J,K)=MOMMIX*(U1(I,K)-U3D(I,K,J))*RDELT
474 ! RVCUTEN(I,J,K)=MOMMIX*(V1(I,K)-V3D(I,K,J))*RDELT
475 RUCUTEN(I,J,K)=(U1(I,K)-U3D(I,K,J))*RDELT
476 RVCUTEN(I,J,K)=(V1(I,K)-V3D(I,K,J))*RDELT
477 ENDDO
478 ENDDO
479 #endif
482 IF(P_QC .ge. P_FIRST_SCALAR)THEN
483 DO K=KTS,KTE
484 DO I=ITS,ITE
485 RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
486 ENDDO
487 ENDDO
488 ENDIF
490 IF(P_QI .ge. P_FIRST_SCALAR)THEN
491 DO K=KTS,KTE
492 DO I=ITS,ITE
493 RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
494 ENDDO
495 ENDDO
496 ENDIF
498 ENDDO big_outer_j_loop ! Outer most J loop
500 END SUBROUTINE CU_SAS
502 !====================================================================
503 SUBROUTINE sasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
504 RUCUTEN,RVCUTEN, &
505 RESTART,P_QC,P_QI,P_FIRST_SCALAR, &
506 allowed_to_read, &
507 ids, ide, jds, jde, kds, kde, &
508 ims, ime, jms, jme, kms, kme, &
509 its, ite, jts, jte, kts, kte )
510 !--------------------------------------------------------------------
511 IMPLICIT NONE
512 !--------------------------------------------------------------------
513 LOGICAL , INTENT(IN) :: allowed_to_read,restart
514 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
515 ims, ime, jms, jme, kms, kme, &
516 its, ite, jts, jte, kts, kte
517 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
519 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
520 RTHCUTEN, &
521 RQVCUTEN, &
522 RQCCUTEN, &
523 RQICUTEN
524 REAL, DIMENSION( ims:ime , jms:jme , kms:kme ) , INTENT(OUT) :: &
525 RUCUTEN, & ! gopal's doing for SAS
526 RVCUTEN
528 INTEGER :: i, j, k, itf, jtf, ktf
530 jtf=min0(jte,jde-1)
531 ktf=min0(kte,kde-1)
532 itf=min0(ite,ide-1)
534 #if ( HWRF == 1 )
535 !zhang's doing
536 IF(.not.restart .or. .not.allowed_to_read)THEN
537 !end of zhang's doing
538 #else
539 IF(.not.restart)THEN
540 #endif
541 DO j=jts,jtf
542 DO k=kts,ktf
543 DO i=its,itf
544 RTHCUTEN(i,k,j)=0.
545 RQVCUTEN(i,k,j)=0.
546 RUCUTEN(i,j,k)=0.
547 RVCUTEN(i,j,k)=0.
548 ENDDO
549 ENDDO
550 ENDDO
552 IF (P_QC .ge. P_FIRST_SCALAR) THEN
553 DO j=jts,jtf
554 DO k=kts,ktf
555 DO i=its,itf
556 RQCCUTEN(i,k,j)=0.
557 ENDDO
558 ENDDO
559 ENDDO
560 ENDIF
562 IF (P_QI .ge. P_FIRST_SCALAR) THEN
563 DO j=jts,jtf
564 DO k=kts,ktf
565 DO i=its,itf
566 RQICUTEN(i,k,j)=0.
567 ENDDO
568 ENDDO
569 ENDDO
570 ENDIF
571 ENDIF
573 END SUBROUTINE sasinit
575 ! ------------------------------------------------------------------------
577 SUBROUTINE SASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PS,PHIL,QL, &
578 ! SUBROUTINE SASCNV(IM,IX,KM,JCAP,DLT,DEL,PRSL,PHIL,QL, &
579 & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK, &
580 & DOT,XKT2,ncloud)
581 ! for cloud water version
582 ! parameter(ncloud=0)
583 ! SUBROUTINE SASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
584 ! & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
585 ! & DOT,xkt2,ncloud)
586 !
587 USE MODULE_GFS_MACHINE , ONLY : kind_phys
588 USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
589 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
590 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
591 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
592 &, EPS => con_eps, EPSM1 => con_epsm1
594 implicit none
595 !
596 ! include 'constant.h'
597 !
598 integer IM, IX, KM, JCAP, ncloud, &
599 & KBOT(IM), KTOP(IM), KUO(IM), J
600 real(kind=kind_phys) DELT
601 real(kind=kind_phys) PS(IM), DEL(IX,KM), PRSL(IX,KM), &
602 ! real(kind=kind_phys) DEL(IX,KM), PRSL(IX,KM),
603 & QL(IX,KM,2), Q1(IX,KM), T1(IX,KM), &
604 & U1(IX,KM), V1(IX,KM), RCS(IM), &
605 & CLDWRK(IM), RN(IM), SLIMSK(IM), &
606 & DOT(IX,KM), XKT2(IM), PHIL(IX,KM)
607 !
608 integer I, INDX, jmn, k, knumb, latd, lond, km1
609 !
610 real(kind=kind_phys) adw, alpha, alphal, alphas, &
611 & aup, beta, betal, betas, &
612 & c0, cpoel, dellat, delta, &
613 & desdt, deta, detad, dg, &
614 & dh, dhh, dlnsig, dp, &
615 & dq, dqsdp, dqsdt, dt, &
616 & dt2, dtmax, dtmin, dv1, &
617 & dv1q, dv2, dv2q, dv1u, &
618 & dv1v, dv2u, dv2v, dv3u, &
619 & dv3v, dv3, dv3q, dvq1, &
620 & dz, dz1, e1, edtmax, &
621 & edtmaxl, edtmaxs, el2orc, elocp, &
622 & es, etah, &
623 & evef, evfact, evfactl, fact1, &
624 & fact2, factor, fjcap, fkm, &
625 & fuv, g, gamma, onemf, &
626 & onemfu, pdetrn, pdpdwn, pprime, &
627 & qc, qlk, qrch, qs, &
628 & rain, rfact, shear, tem1, &
629 & tem2, terr, val, val1, &
630 & val2, w1, w1l, w1s, &
631 & w2, w2l, w2s, w3, &
632 & w3l, w3s, w4, w4l, &
633 & w4s, xdby, xpw, xpwd, &
634 & xqc, xqrch, xlambu, mbdt, &
635 & tem
636 !
637 !
638 integer JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM), &
639 & KT2(IM), KTCON(IM), LMIN(IM), &
640 & kbm(IM), kbmax(IM), kmax(IM)
641 !
642 real(kind=kind_phys) AA1(IM), ACRT(IM), ACRTFCT(IM), &
643 & DELHBAR(IM), DELQ(IM), DELQ2(IM), &
644 & DELQBAR(IM), DELQEV(IM), DELTBAR(IM), &
645 & DELTV(IM), DTCONV(IM), EDT(IM), &
646 & EDTO(IM), EDTX(IM), FLD(IM), &
647 & HCDO(IM), HKBO(IM), HMAX(IM), &
648 & HMIN(IM), HSBAR(IM), UCDO(IM), &
649 & UKBO(IM), VCDO(IM), VKBO(IM), &
650 & PBCDIF(IM), PDOT(IM), PO(IM,KM), &
651 & PWAVO(IM), PWEVO(IM), &
652 ! & PSFC(IM), PWAVO(IM), PWEVO(IM), &
653 & QCDO(IM), QCOND(IM), QEVAP(IM), &
654 & QKBO(IM), RNTOT(IM), VSHEAR(IM), &
655 & XAA0(IM), XHCD(IM), XHKB(IM), &
656 & XK(IM), XLAMB(IM), XLAMD(IM), &
657 & XMB(IM), XMBMAX(IM), XPWAV(IM), &
658 & XPWEV(IM), XQCD(IM), XQKB(IM)
659 !
660 ! PHYSICAL PARAMETERS
661 PARAMETER(G=grav)
662 PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP, &
663 & EL2ORC=HVAP*HVAP/(RV*CP))
664 PARAMETER(TERR=0.,C0=.002,DELTA=fv)
665 PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
666 ! LOCAL VARIABLES AND ARRAYS
667 real(kind=kind_phys) PFLD(IM,KM), TO(IM,KM), QO(IM,KM), &
668 & UO(IM,KM), VO(IM,KM), QESO(IM,KM)
669 ! cloud water
670 real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM), TVO(IM,KM), &
671 & DBYO(IM,KM), ZO(IM,KM), SUMZ(IM,KM), &
672 & SUMH(IM,KM), HEO(IM,KM), HESO(IM,KM), &
673 & QRCD(IM,KM), DELLAH(IM,KM), DELLAQ(IM,KM),&
674 & DELLAU(IM,KM), DELLAV(IM,KM), HCKO(IM,KM), &
675 & UCKO(IM,KM), VCKO(IM,KM), QCKO(IM,KM), &
676 & ETA(IM,KM), ETAU(IM,KM), ETAD(IM,KM), &
677 & QRCDO(IM,KM), PWO(IM,KM), PWDO(IM,KM), &
678 & RHBAR(IM), TX1(IM)
679 !
680 LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
681 !
682 real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
683 ! SAVE PCRIT, ACRITT
684 DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
685 & 350.,300.,250.,200.,150./
686 DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
687 & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
688 ! GDAS DERIVED ACRIT
689 ! DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
690 ! & .743,.813,.886,.947,1.138,1.377,1.896/
691 !
692 real(kind=kind_phys) TF, TCR, TCRF, RZERO, RONE
693 parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF))
694 parameter (RZERO=0.0,RONE=1.0)
695 !-----------------------------------------------------------------------
696 !
697 km1 = km - 1
698 ! INITIALIZE ARRAYS
699 !
700 DO I=1,IM
701 RN(I)=0.
702 KBOT(I)=KM+1
703 KTOP(I)=0
704 KUO(I)=0
705 CNVFLG(I) = .TRUE.
706 DTCONV(I) = 3600.
707 CLDWRK(I) = 0.
708 PDOT(I) = 0.
709 KT2(I) = 0
710 QLKO_KTCON(I) = 0.
711 DELLAL(I) = 0.
712 ENDDO
713 !!
714 DO K = 1, 15
715 ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
716 ENDDO
717 DT2 = DELT
718 !cmr dtmin = max(dt2,1200.)
719 val = 1200.
720 dtmin = max(dt2, val )
721 !cmr dtmax = max(dt2,3600.)
722 val = 3600.
723 dtmax = max(dt2, val )
724 ! MODEL TUNABLE PARAMETERS ARE ALL HERE
725 MBDT = 10.
726 EDTMAXl = .3
727 EDTMAXs = .3
728 ALPHAl = .5
729 ALPHAs = .5
730 BETAl = .15
731 betas = .15
732 BETAl = .05
733 betas = .05
734 ! change for hurricane model
735 BETAl = .5
736 betas = .5
737 ! EVEF = 0.07
738 evfact = 0.3
739 evfactl = 0.3
740 ! change for hurricane model
741 evfact = 0.6
742 evfactl = .6
743 #if ( EM_CORE == 1 )
744 ! HAWAII TEST - ZCX
745 ALPHAl = .5
746 ALPHAs = .75
747 BETAl = .05
748 betas = .05
749 evfact = 0.5
750 evfactl = 0.5
751 #endif
752 PDPDWN = 0.
753 PDETRN = 200.
754 xlambu = 1.e-4
755 fjcap = (float(jcap) / 126.) ** 2
756 !cmr fjcap = max(fjcap,1.)
757 val = 1.
758 fjcap = max(fjcap,val)
759 fkm = (float(km) / 28.) ** 2
760 !cmr fkm = max(fkm,1.)
761 fkm = max(fkm,val)
762 W1l = -8.E-3
763 W2l = -4.E-2
764 W3l = -5.E-3
765 W4l = -5.E-4
766 W1s = -2.E-4
767 W2s = -2.E-3
768 W3s = -1.E-3
769 W4s = -2.E-5
770 !CCCC IF(IM.EQ.384) THEN
771 LATD = 92
772 lond = 189
773 !CCCC ELSEIF(IM.EQ.768) THEN
774 !CCCC LATD = 80
775 !CCCC ELSE
776 !CCCC LATD = 0
777 !CCCC ENDIF
778 !
779 ! DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
780 ! AND THE MAXIMUM THETAE FOR UPDRAFT
781 !
782 DO I=1,IM
783 KBMAX(I) = KM
784 KBM(I) = KM
785 KMAX(I) = KM
786 TX1(I) = 1.0 / PS(I)
787 ENDDO
788 !
789 DO K = 1, KM
790 DO I=1,IM
791 IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
792 IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I) = K + 1
793 IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
794 ENDDO
795 ENDDO
796 DO I=1,IM
797 KBMAX(I) = MIN(KBMAX(I),KMAX(I))
798 KBM(I) = MIN(KBM(I),KMAX(I))
799 ENDDO
800 !
801 ! CONVERT SURFACE PRESSURE TO MB FROM CB
802 !
803 !!
804 DO K = 1, KM
805 DO I=1,IM
806 if (K .le. kmax(i)) then
807 PFLD(I,k) = PRSL(I,K) * 10.0
808 PWO(I,k) = 0.
809 PWDO(I,k) = 0.
810 TO(I,k) = T1(I,k)
811 QO(I,k) = Q1(I,k)
812 UO(I,k) = U1(I,k)
813 VO(I,k) = V1(I,k)
814 DBYO(I,k) = 0.
815 SUMZ(I,k) = 0.
816 SUMH(I,k) = 0.
817 endif
818 ENDDO
819 ENDDO
821 !
822 ! COLUMN VARIABLES
823 ! P IS PRESSURE OF THE LAYER (MB)
824 ! T IS TEMPERATURE AT T-DT (K)..TN
825 ! Q IS MIXING RATIO AT T-DT (KG/KG)..QN
826 ! TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
827 ! QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
828 !
829 DO K = 1, KM
830 DO I=1,IM
831 if (k .le. kmax(i)) then
832 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
833 !
834 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
835 !
836 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
837 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
838 val1 = 1.E-8
839 QESO(I,k) = MAX(QESO(I,k), val1)
840 !cmr QO(I,k) = max(QO(I,k),1.e-10)
841 val2 = 1.e-10
842 QO(I,k) = max(QO(I,k), val2 )
843 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
844 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
845 endif
846 ENDDO
847 ENDDO
849 !
850 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
851 !
852 DO K = 1, KM
853 DO I=1,IM
854 ZO(I,k) = PHIL(I,k) / G
855 ENDDO
856 ENDDO
857 ! COMPUTE MOIST STATIC ENERGY
858 DO K = 1, KM
859 DO I=1,IM
860 if (K .le. kmax(i)) then
861 ! tem = G * ZO(I,k) + CP * TO(I,k)
862 tem = PHIL(I,k) + CP * TO(I,k)
863 HEO(I,k) = tem + HVAP * QO(I,k)
864 HESO(I,k) = tem + HVAP * QESO(I,k)
865 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
866 endif
867 ENDDO
868 ENDDO
869 !
870 ! DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
871 ! THIS IS THE LEVEL WHERE UPDRAFT STARTS
872 !
873 DO I=1,IM
874 HMAX(I) = HEO(I,1)
875 KB(I) = 1
876 ENDDO
877 !!
878 DO K = 2, KM
879 DO I=1,IM
880 if (k .le. kbm(i)) then
881 IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
882 KB(I) = K
883 HMAX(I) = HEO(I,k)
884 ENDIF
885 endif
886 ENDDO
887 ENDDO
888 ! DO K = 1, KMAX - 1
889 ! TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
890 ! QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
891 ! QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
892 ! HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
893 ! HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
894 ! ENDDO
895 DO K = 1, KM1
896 DO I=1,IM
897 if (k .le. kmax(i)-1) then
898 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
899 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
900 !jfe ES = 10. * FPVS(TO(I,k+1))
901 !
902 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
903 !
904 PPRIME = PFLD(I,k+1) + EPSM1 * ES
905 QS = EPS * ES / PPRIME
906 DQSDP = - QS / PPRIME
907 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
908 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
909 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
910 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
911 DQ = DQSDT * DT + DQSDP * DP
912 TO(I,k) = TO(I,k+1) + DT
913 QO(I,k) = QO(I,k+1) + DQ
914 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
915 endif
916 ENDDO
917 ENDDO
918 !
919 DO K = 1, KM1
920 DO I=1,IM
921 if (k .le. kmax(I)-1) then
922 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
923 !
924 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
925 !
926 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
927 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
928 val1 = 1.E-8
929 QESO(I,k) = MAX(QESO(I,k), val1)
930 !cmr QO(I,k) = max(QO(I,k),1.e-10)
931 val2 = 1.e-10
932 QO(I,k) = max(QO(I,k), val2 )
933 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
934 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
935 & CP * TO(I,k) + HVAP * QO(I,k)
936 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
937 & CP * TO(I,k) + HVAP * QESO(I,k)
938 UO(I,k) = .5 * (UO(I,k) + UO(I,k+1))
939 VO(I,k) = .5 * (VO(I,k) + VO(I,k+1))
940 endif
941 ENDDO
942 ENDDO
943 ! k = kmax
944 ! HEO(I,k) = HEO(I,k)
945 ! hesol(k) = HESO(I,k)
946 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
947 ! PRINT *, ' HEO ='
948 ! PRINT 6001, (HEO(I,K),K=1,KMAX)
949 ! PRINT *, ' HESO ='
950 ! PRINT 6001, (HESO(I,K),K=1,KMAX)
951 ! PRINT *, ' TO ='
952 ! PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
953 ! PRINT *, ' QO ='
954 ! PRINT 6003, (QO(I,K),K=1,KMAX)
955 ! PRINT *, ' QSO ='
956 ! PRINT 6003, (QESO(I,K),K=1,KMAX)
957 ! ENDIF
958 !
959 ! LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
960 !
961 DO I=1,IM
962 IF(CNVFLG(I)) THEN
963 INDX = KB(I)
964 HKBO(I) = HEO(I,INDX)
965 QKBO(I) = QO(I,INDX)
966 UKBO(I) = UO(I,INDX)
967 VKBO(I) = VO(I,INDX)
968 ENDIF
969 FLG(I) = CNVFLG(I)
970 KBCON(I) = KMAX(I)
971 ENDDO
972 !!
973 DO K = 1, KM
974 DO I=1,IM
975 if (k .le. kbmax(i)) then
976 IF(FLG(I).AND.K.GT.KB(I)) THEN
977 HSBAR(I) = HESO(I,k)
978 IF(HKBO(I).GT.HSBAR(I)) THEN
979 FLG(I) = .FALSE.
980 KBCON(I) = K
981 ENDIF
982 ENDIF
983 endif
984 ENDDO
985 ENDDO
986 DO I=1,IM
987 IF(CNVFLG(I)) THEN
988 PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
989 PDOT(I) = 10.* DOT(I,KBCON(I))
990 IF(PBCDIF(I).GT.150.) CNVFLG(I) = .FALSE.
991 IF(KBCON(I).EQ.KMAX(I)) CNVFLG(I) = .FALSE.
992 ENDIF
993 ENDDO
994 !!
995 TOTFLG = .TRUE.
996 DO I=1,IM
997 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
998 ENDDO
999 IF(TOTFLG) RETURN
1000 ! FOUND LFC, CAN DEFINE REST OF VARIABLES
1001 6001 FORMAT(2X,-2P10F12.2)
1002 6002 FORMAT(2X,10F12.2)
1003 6003 FORMAT(2X,3P10F12.2)
1005 !
1006 ! DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
1007 !
1008 DO I = 1, IM
1009 alpha = alphas
1010 if(SLIMSK(I).eq.1.) alpha = alphal
1011 IF(CNVFLG(I)) THEN
1012 IF(KB(I).EQ.1) THEN
1013 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
1014 ELSE
1015 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) &
1016 & - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
1017 ENDIF
1018 IF(KBCON(I).NE.KB(I)) THEN
1019 !cmr XLAMB(I) = -ALOG(ALPHA) / DZ
1020 XLAMB(I) = - LOG(ALPHA) / DZ
1021 ELSE
1022 XLAMB(I) = 0.
1023 ENDIF
1024 ENDIF
1025 ENDDO
1026 ! DETERMINE UPDRAFT MASS FLUX
1027 DO K = 1, KM
1028 DO I = 1, IM
1029 if (k .le. kmax(i) .and. CNVFLG(I)) then
1030 ETA(I,k) = 1.
1031 ETAU(I,k) = 1.
1032 ENDIF
1033 ENDDO
1034 ENDDO
1035 DO K = KM1, 2, -1
1036 DO I = 1, IM
1037 if (k .le. kbmax(i)) then
1038 IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
1039 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1040 ETA(I,k) = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
1041 ETAU(I,k) = ETA(I,k)
1042 ENDIF
1043 endif
1044 ENDDO
1045 ENDDO
1046 DO I = 1, IM
1047 IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
1048 DZ = .5 * (ZO(I,2) - ZO(I,1))
1049 ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
1050 ETAU(I,1) = ETA(I,1)
1051 ENDIF
1052 ENDDO
1053 !
1054 ! WORK UP UPDRAFT CLOUD PROPERTIES
1055 !
1056 DO I = 1, IM
1057 IF(CNVFLG(I)) THEN
1058 INDX = KB(I)
1059 HCKO(I,INDX) = HKBO(I)
1060 QCKO(I,INDX) = QKBO(I)
1061 UCKO(I,INDX) = UKBO(I)
1062 VCKO(I,INDX) = VKBO(I)
1063 PWAVO(I) = 0.
1064 ENDIF
1065 ENDDO
1066 !
1067 ! CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
1068 !
1069 DO K = 2, KM1
1070 DO I = 1, IM
1071 if (k .le. kmax(i)-1) then
1072 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1073 FACTOR = ETA(I,k-1) / ETA(I,k)
1074 ONEMF = 1. - FACTOR
1075 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1076 & .5 * (HEO(I,k) + HEO(I,k+1))
1077 UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF * &
1078 & .5 * (UO(I,k) + UO(I,k+1))
1079 VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF * &
1080 & .5 * (VO(I,k) + VO(I,k+1))
1081 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1082 ENDIF
1083 IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1084 HCKO(I,k) = HCKO(I,k-1)
1085 UCKO(I,k) = UCKO(I,k-1)
1086 VCKO(I,k) = VCKO(I,k-1)
1087 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1088 ENDIF
1089 endif
1090 ENDDO
1091 ENDDO
1092 ! DETERMINE CLOUD TOP
1093 DO I = 1, IM
1094 FLG(I) = CNVFLG(I)
1095 KTCON(I) = 1
1096 ENDDO
1097 ! DO K = 2, KMAX
1098 ! KK = KMAX - K + 1
1099 ! IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
1100 ! KTCON(I) = KK + 1
1101 ! FLG(I) = .FALSE.
1102 ! ENDIF
1103 ! ENDDO
1104 DO K = 2, KM
1105 DO I = 1, IM
1106 if (k .le. kmax(i)) then
1107 IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
1108 KTCON(I) = K
1109 FLG(I) = .FALSE.
1110 ENDIF
1111 endif
1112 ENDDO
1113 ENDDO
1114 DO I = 1, IM
1115 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.) &
1116 & CNVFLG(I) = .FALSE.
1117 ENDDO
1118 TOTFLG = .TRUE.
1119 DO I = 1, IM
1120 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1121 ENDDO
1122 IF(TOTFLG) RETURN
1123 !
1124 ! SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
1125 !
1126 DO I = 1, IM
1127 HMIN(I) = HEO(I,KBCON(I))
1128 LMIN(I) = KBMAX(I)
1129 JMIN(I) = KBMAX(I)
1130 ENDDO
1131 DO I = 1, IM
1132 DO K = KBCON(I), KBMAX(I)
1133 IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
1134 LMIN(I) = K + 1
1135 HMIN(I) = HEO(I,k)
1136 ENDIF
1137 ENDDO
1138 ENDDO
1139 !
1140 ! Make sure that JMIN(I) is within the cloud
1141 !
1142 DO I = 1, IM
1143 IF(CNVFLG(I)) THEN
1144 JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
1145 XMBMAX(I) = .1
1146 JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
1147 ENDIF
1148 ENDDO
1149 !
1150 ! ENTRAINING CLOUD
1151 !
1152 do k = 2, km1
1153 DO I = 1, IM
1154 if (k .le. kmax(i)-1) then
1155 if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
1156 SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
1157 SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1)) &
1158 & * HEO(I,k)
1159 ENDIF
1160 endif
1161 enddo
1162 enddo
1163 !!
1164 DO I = 1, IM
1165 IF(CNVFLG(I)) THEN
1166 ! call random_number(XKT2)
1167 ! call srand(fhour)
1168 ! XKT2(I) = rand()
1169 KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
1170 ! KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1171 ! KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1172 tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1173 tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1174 if (abs(tem2) .gt. 0.000001) THEN
1175 XLAMB(I) = tem1 / tem2
1176 else
1177 CNVFLG(I) = .false.
1178 ENDIF
1179 ! XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1180 ! & / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1181 XLAMB(I) = max(XLAMB(I),RZERO)
1182 XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
1183 ENDIF
1184 ENDDO
1185 !!
1186 DO I = 1, IM
1187 DWNFLG(I) = CNVFLG(I)
1188 DWNFLG2(I) = CNVFLG(I)
1189 IF(CNVFLG(I)) THEN
1190 if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
1191 if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)&
1192 & DWNFLG(I) = .false.
1193 do k = JMIN(I), KT2(I)
1194 if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
1195 enddo
1196 ! IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
1197 ! & DWNFLG(I)=.FALSE.
1198 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN) &
1199 & DWNFLG2(I)=.FALSE.
1200 ENDIF
1201 ENDDO
1202 !!
1203 DO K = 2, KM1
1204 DO I = 1, IM
1205 if (k .le. kmax(i)-1) then
1206 IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1207 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1208 ! ETA(I,k) = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
1209 ! to simplify matter, we will take the linear approach here
1210 !
1211 ETA(I,k) = ETA(I,k-1) * (1. + XLAMB(I) * dz)
1212 ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
1213 ENDIF
1214 endif
1215 ENDDO
1216 ENDDO
1217 !!
1218 DO K = 2, KM1
1219 DO I = 1, IM
1220 if (k .le. kmax(i)-1) then
1221 ! IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1222 IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
1223 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1224 ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
1225 ENDIF
1226 endif
1227 ENDDO
1228 ENDDO
1229 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1230 ! PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
1231 ! PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
1232 ! ENDIF
1233 ! IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
1234 ! print *, ' xlamb =', xlamb
1235 ! print *, ' eta =', (eta(k),k=1,KT2(I))
1236 ! print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
1237 ! print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
1238 ! print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
1239 ! print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
1240 ! ENDIF
1241 DO I = 1, IM
1242 if(DWNFLG(I)) THEN
1243 KTCON(I) = KT2(I)
1244 ENDIF
1245 ENDDO
1246 !
1247 ! CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
1248 !
1249 DO K = 2, KM1
1250 DO I = 1, IM
1251 if (k .le. kmax(i)-1) then
1252 !jfe
1253 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1254 !jfe IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1255 FACTOR = ETA(I,k-1) / ETA(I,k)
1256 ONEMF = 1. - FACTOR
1257 fuv = ETAU(I,k-1) / ETAU(I,k)
1258 onemfu = 1. - fuv
1259 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1260 & .5 * (HEO(I,k) + HEO(I,k+1))
1261 UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu * &
1262 & .5 * (UO(I,k) + UO(I,k+1))
1263 VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu * &
1264 & .5 * (VO(I,k) + VO(I,k+1))
1265 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1266 ENDIF
1267 endif
1268 ENDDO
1269 ENDDO
1270 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1271 ! PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
1272 ! PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
1273 ! ENDIF
1274 DO I = 1, IM
1275 if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I)) &
1276 & THEN
1277 CNVFLG(I) = .false.
1278 DWNFLG(I) = .false.
1279 DWNFLG2(I) = .false.
1280 ENDIF
1281 ENDDO
1282 !!
1283 TOTFLG = .TRUE.
1284 DO I = 1, IM
1285 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1286 ENDDO
1287 IF(TOTFLG) RETURN
1288 !!
1289 !
1290 ! COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
1291 !
1292 DO I = 1, IM
1293 AA1(I) = 0.
1294 RHBAR(I) = 0.
1295 ENDDO
1296 DO K = 1, KM
1297 DO I = 1, IM
1298 if (k .le. kmax(i)) then
1299 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1300 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1301 DZ1 = (ZO(I,k) - ZO(I,k-1))
1302 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1303 QRCH = QESO(I,k) &
1304 & + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
1305 FACTOR = ETA(I,k-1) / ETA(I,k)
1306 ONEMF = 1. - FACTOR
1307 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1308 & .5 * (QO(I,k) + QO(I,k+1))
1309 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
1310 RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
1311 !
1312 ! BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1313 !
1314 IF(DQ.GT.0.) THEN
1315 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1316 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1317 AA1(I) = AA1(I) - DZ1 * G * QLK
1318 QC = QLK + QRCH
1319 PWO(I,k) = ETAH * C0 * DZ * QLK
1320 QCKO(I,k) = QC
1321 PWAVO(I) = PWAVO(I) + PWO(I,k)
1322 ENDIF
1323 ENDIF
1324 endif
1325 ENDDO
1326 ENDDO
1327 DO I = 1, IM
1328 RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
1329 ENDDO
1330 !
1331 ! this section is ready for cloud water
1332 !
1333 if(ncloud.gt.0) THEN
1334 !
1335 ! compute liquid and vapor separation at cloud top
1336 !
1337 DO I = 1, IM
1338 k = KTCON(I)
1339 IF(CNVFLG(I)) THEN
1340 GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
1341 QRCH = QESO(I,K) &
1342 & + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
1343 DQ = QCKO(I,K-1) - QRCH
1344 !
1345 ! CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1346 !
1347 IF(DQ.GT.0.) THEN
1348 QLKO_KTCON(I) = dq
1349 QCKO(I,K-1) = QRCH
1350 ENDIF
1351 ENDIF
1352 ENDDO
1353 ENDIF
1354 !
1355 ! CALCULATE CLOUD WORK FUNCTION AT T+DT
1356 !
1357 DO K = 1, KM
1358 DO I = 1, IM
1359 if (k .le. kmax(i)) then
1360 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1361 DZ1 = ZO(I,k) - ZO(I,k-1)
1362 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1363 RFACT = 1. + DELTA * CP * GAMMA &
1364 & * TO(I,k-1) / HVAP
1365 AA1(I) = AA1(I) + &
1366 & DZ1 * (G / (CP * TO(I,k-1))) &
1367 & * DBYO(I,k-1) / (1. + GAMMA) &
1368 & * RFACT
1369 val = 0.
1370 AA1(I)=AA1(I)+ &
1371 & DZ1 * G * DELTA * &
1372 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1373 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1374 ENDIF
1375 endif
1376 ENDDO
1377 ENDDO
1378 DO I = 1, IM
1379 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1380 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1381 IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1382 ENDDO
1383 !!
1384 TOTFLG = .TRUE.
1385 DO I = 1, IM
1386 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1387 ENDDO
1388 IF(TOTFLG) RETURN
1389 !!
1390 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1391 !cccc PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
1392 !cccc ENDIF
1393 !
1394 !------- DOWNDRAFT CALCULATIONS
1395 !
1396 !
1397 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1398 !
1399 DO I = 1, IM
1400 IF(CNVFLG(I)) THEN
1401 VSHEAR(I) = 0.
1402 ENDIF
1403 ENDDO
1404 DO K = 1, KM
1405 DO I = 1, IM
1406 if (k .le. kmax(i)) then
1407 IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
1408 shear=rcs(I) * sqrt((UO(I,k+1)-UO(I,k)) ** 2 &
1409 & + (VO(I,k+1)-VO(I,k)) ** 2)
1410 VSHEAR(I) = VSHEAR(I) + SHEAR
1411 ENDIF
1412 endif
1413 ENDDO
1414 ENDDO
1415 DO I = 1, IM
1416 EDT(I) = 0.
1417 IF(CNVFLG(I)) THEN
1418 KNUMB = KTCON(I) - KB(I) + 1
1419 KNUMB = MAX(KNUMB,1)
1420 VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
1421 E1=1.591-.639*VSHEAR(I) &
1422 & +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
1423 EDT(I)=1.-E1
1424 !cmr EDT(I) = MIN(EDT(I),.9)
1425 val = .9
1426 EDT(I) = MIN(EDT(I),val)
1427 !cmr EDT(I) = MAX(EDT(I),.0)
1428 val = .0
1429 EDT(I) = MAX(EDT(I),val)
1430 EDTO(I)=EDT(I)
1431 EDTX(I)=EDT(I)
1432 ENDIF
1433 ENDDO
1434 ! DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
1435 DO I = 1, IM
1436 KBDTR(I) = KBCON(I)
1437 beta = betas
1438 if(SLIMSK(I).eq.1.) beta = betal
1439 IF(CNVFLG(I)) THEN
1440 KBDTR(I) = KBCON(I)
1441 KBDTR(I) = MAX(KBDTR(I),1)
1442 XLAMD(I) = 0.
1443 IF(KBDTR(I).GT.1) THEN
1444 DZ = .5 * ZO(I,KBDTR(I)) + .5 * ZO(I,KBDTR(I)-1) &
1445 & - ZO(I,1)
1446 XLAMD(I) = LOG(BETA) / DZ
1447 ENDIF
1448 ENDIF
1449 ENDDO
1450 ! DETERMINE DOWNDRAFT MASS FLUX
1451 DO K = 1, KM
1452 DO I = 1, IM
1453 IF(k .le. kmax(i)) then
1454 IF(CNVFLG(I)) THEN
1455 ETAD(I,k) = 1.
1456 ENDIF
1457 QRCDO(I,k) = 0.
1458 endif
1459 ENDDO
1460 ENDDO
1461 DO K = KM1, 2, -1
1462 DO I = 1, IM
1463 if (k .le. kbmax(i)) then
1464 IF(CNVFLG(I).AND.K.LT.KBDTR(I)) THEN
1465 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1466 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1467 ENDIF
1468 endif
1469 ENDDO
1470 ENDDO
1471 K = 1
1472 DO I = 1, IM
1473 IF(CNVFLG(I).AND.KBDTR(I).GT.1) THEN
1474 DZ = .5 * (ZO(I,2) - ZO(I,1))
1475 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1476 ENDIF
1477 ENDDO
1478 !
1479 !--- DOWNDRAFT MOISTURE PROPERTIES
1480 !
1481 DO I = 1, IM
1482 PWEVO(I) = 0.
1483 FLG(I) = CNVFLG(I)
1484 ENDDO
1485 DO I = 1, IM
1486 IF(CNVFLG(I)) THEN
1487 JMN = JMIN(I)
1488 HCDO(I) = HEO(I,JMN)
1489 QCDO(I) = QO(I,JMN)
1490 QRCDO(I,JMN) = QESO(I,JMN)
1491 UCDO(I) = UO(I,JMN)
1492 VCDO(I) = VO(I,JMN)
1493 ENDIF
1494 ENDDO
1495 DO K = KM1, 1, -1
1496 DO I = 1, IM
1497 if (k .le. kmax(i)-1) then
1498 IF(CNVFLG(I).AND.K.LT.JMIN(I)) THEN
1499 DQ = QESO(I,k)
1500 DT = TO(I,k)
1501 GAMMA = EL2ORC * DQ / DT**2
1502 DH = HCDO(I) - HESO(I,k)
1503 QRCDO(I,k) = DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1504 DETAD = ETAD(I,k+1) - ETAD(I,k)
1505 PWDO(I,k) = ETAD(I,k+1) * QCDO(I) - &
1506 & ETAD(I,k) * QRCDO(I,k)
1507 PWDO(I,k) = PWDO(I,k) - DETAD * &
1508 & .5 * (QRCDO(I,k) + QRCDO(I,k+1))
1509 QCDO(I) = QRCDO(I,k)
1510 PWEVO(I) = PWEVO(I) + PWDO(I,k)
1511 ENDIF
1512 endif
1513 ENDDO
1514 ENDDO
1515 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG(I)) THEN
1516 ! PRINT *, ' PWAVO(I), PWEVO(I) =', PWAVO(I), PWEVO(I)
1517 ! ENDIF
1518 !
1519 !--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
1520 !--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
1521 !--- EVAPORATE (PWEV)
1522 !
1523 DO I = 1, IM
1524 edtmax = edtmaxl
1525 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1526 IF(DWNFLG2(I)) THEN
1527 IF(PWEVO(I).LT.0.) THEN
1528 EDTO(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
1529 EDTO(I) = MIN(EDTO(I),EDTMAX)
1530 ELSE
1531 EDTO(I) = 0.
1532 ENDIF
1533 ELSE
1534 EDTO(I) = 0.
1535 ENDIF
1536 ENDDO
1537 !
1538 !
1539 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1540 !
1541 !
1542 DO K = KM1, 1, -1
1543 DO I = 1, IM
1544 if (k .le. kmax(i)-1) then
1545 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1546 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1547 DHH=HCDO(I)
1548 DT=TO(I,k+1)
1549 DG=GAMMA
1550 DH=HESO(I,k+1)
1551 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1552 AA1(I)=AA1(I)+EDTO(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1553 & *(1.+DELTA*CP*DG*DT/HVAP)
1554 val=0.
1555 AA1(I)=AA1(I)+EDTO(I)* &
1556 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
1557 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1558 ENDIF
1559 endif
1560 ENDDO
1561 ENDDO
1562 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1563 !cccc PRINT *, ' AA1(I) AFTER DWNDRFT =', AA1(I)
1564 !cccc ENDIF
1565 DO I = 1, IM
1566 IF(AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1567 IF(AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1568 IF(AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1569 ENDDO
1570 !!
1571 TOTFLG = .TRUE.
1572 DO I = 1, IM
1573 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1574 ENDDO
1575 IF(TOTFLG) RETURN
1576 !!
1577 !
1578 !
1579 !--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
1580 !--- WILL DO TO THE ENVIRONMENT?
1581 !
1582 DO K = 1, KM
1583 DO I = 1, IM
1584 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1585 DELLAH(I,k) = 0.
1586 DELLAQ(I,k) = 0.
1587 DELLAU(I,k) = 0.
1588 DELLAV(I,k) = 0.
1589 ENDIF
1590 ENDDO
1591 ENDDO
1592 DO I = 1, IM
1593 IF(CNVFLG(I)) THEN
1594 DP = 1000. * DEL(I,1)
1595 DELLAH(I,1) = EDTO(I) * ETAD(I,1) * (HCDO(I) &
1596 & - HEO(I,1)) * G / DP
1597 DELLAQ(I,1) = EDTO(I) * ETAD(I,1) * (QCDO(I) &
1598 & - QO(I,1)) * G / DP
1599 DELLAU(I,1) = EDTO(I) * ETAD(I,1) * (UCDO(I) &
1600 & - UO(I,1)) * G / DP
1601 DELLAV(I,1) = EDTO(I) * ETAD(I,1) * (VCDO(I) &
1602 & - VO(I,1)) * G / DP
1603 ENDIF
1604 ENDDO
1605 !
1606 !--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1607 !
1608 DO K = 2, KM1
1609 DO I = 1, IM
1610 if (k .le. kmax(i)-1) then
1611 IF(CNVFLG(I).AND.K.LT.KTCON(I)) THEN
1612 AUP = 1.
1613 IF(K.LE.KB(I)) AUP = 0.
1614 ADW = 1.
1615 IF(K.GT.JMIN(I)) ADW = 0.
1616 DV1= HEO(I,k)
1617 DV2 = .5 * (HEO(I,k) + HEO(I,k+1))
1618 DV3= HEO(I,k-1)
1619 DV1Q= QO(I,k)
1620 DV2Q = .5 * (QO(I,k) + QO(I,k+1))
1621 DV3Q= QO(I,k-1)
1622 DV1U= UO(I,k)
1623 DV2U = .5 * (UO(I,k) + UO(I,k+1))
1624 DV3U= UO(I,k-1)
1625 DV1V= VO(I,k)
1626 DV2V = .5 * (VO(I,k) + VO(I,k+1))
1627 DV3V= VO(I,k-1)
1628 DP = 1000. * DEL(I,K)
1629 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1630 DETA = ETA(I,k) - ETA(I,k-1)
1631 DETAD = ETAD(I,k) - ETAD(I,k-1)
1632 DELLAH(I,k) = DELLAH(I,k) + &
1633 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1 &
1634 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3 &
1635 & - AUP * DETA * DV2 &
1636 & + ADW * EDTO(I) * DETAD * HCDO(I)) * G / DP
1637 DELLAQ(I,k) = DELLAQ(I,k) + &
1638 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1Q &
1639 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3Q &
1640 & - AUP * DETA * DV2Q &
1641 & +ADW*EDTO(I)*DETAD*.5*(QRCDO(I,k)+QRCDO(I,k-1))) * G / DP
1642 DELLAU(I,k) = DELLAU(I,k) + &
1643 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1U &
1644 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3U &
1645 & - AUP * DETA * DV2U &
1646 & + ADW * EDTO(I) * DETAD * UCDO(I) &
1647 & ) * G / DP
1648 DELLAV(I,k) = DELLAV(I,k) + &
1649 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1V &
1650 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3V &
1651 & - AUP * DETA * DV2V &
1652 & + ADW * EDTO(I) * DETAD * VCDO(I) &
1653 & ) * G / DP
1654 ENDIF
1655 endif
1656 ENDDO
1657 ENDDO
1658 !
1659 !------- CLOUD TOP
1660 !
1661 DO I = 1, IM
1662 IF(CNVFLG(I)) THEN
1663 INDX = KTCON(I)
1664 DP = 1000. * DEL(I,INDX)
1665 DV1 = HEO(I,INDX-1)
1666 DELLAH(I,INDX) = ETA(I,INDX-1) * &
1667 & (HCKO(I,INDX-1) - DV1) * G / DP
1668 DVQ1 = QO(I,INDX-1)
1669 DELLAQ(I,INDX) = ETA(I,INDX-1) * &
1670 & (QCKO(I,INDX-1) - DVQ1) * G / DP
1671 DV1U = UO(I,INDX-1)
1672 DELLAU(I,INDX) = ETA(I,INDX-1) * &
1673 & (UCKO(I,INDX-1) - DV1U) * G / DP
1674 DV1V = VO(I,INDX-1)
1675 DELLAV(I,INDX) = ETA(I,INDX-1) * &
1676 & (VCKO(I,INDX-1) - DV1V) * G / DP
1677 !
1678 ! cloud water
1679 !
1680 DELLAL(I) = ETA(I,INDX-1) * QLKO_KTCON(I) * g / dp
1681 ENDIF
1682 ENDDO
1683 !
1684 !------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
1685 !
1686 DO K = 1, KM
1687 DO I = 1, IM
1688 if (k .le. kmax(i)) then
1689 IF(CNVFLG(I).and.k.gt.KTCON(I)) THEN
1690 QO(I,k) = Q1(I,k)
1691 TO(I,k) = T1(I,k)
1692 UO(I,k) = U1(I,k)
1693 VO(I,k) = V1(I,k)
1694 ENDIF
1695 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1696 QO(I,k) = DELLAQ(I,k) * MBDT + Q1(I,k)
1697 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
1698 TO(I,k) = DELLAT * MBDT + T1(I,k)
1699 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1700 val = 1.e-10
1701 QO(I,k) = max(QO(I,k), val )
1702 ENDIF
1703 endif
1704 ENDDO
1705 ENDDO
1706 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1707 !
1708 !--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
1709 !--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
1710 !--- WOULD HAVE ON THE STABILITY,
1711 !--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
1712 !--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
1713 !--- DESTABILIZATION.
1714 !
1715 !--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
1716 !
1717 DO K = 1, KM
1718 DO I = 1, IM
1719 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1720 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1721 !
1722 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1723 !
1724 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k)+EPSM1*QESO(I,k))
1725 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1726 val = 1.E-8
1727 QESO(I,k) = MAX(QESO(I,k), val )
1728 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
1729 ENDIF
1730 ENDDO
1731 ENDDO
1732 DO I = 1, IM
1733 IF(CNVFLG(I)) THEN
1734 XAA0(I) = 0.
1735 XPWAV(I) = 0.
1736 ENDIF
1737 ENDDO
1738 !
1739 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
1740 !
1741 ! DO I = 1, IM
1742 ! IF(CNVFLG(I)) THEN
1743 ! DLNSIG = LOG(PRSL(I,1)/PS(I))
1744 ! ZO(I,1) = TERR - DLNSIG * RD / G * TVO(I,1)
1745 ! ENDIF
1746 ! ENDDO
1747 ! DO K = 2, KM
1748 ! DO I = 1, IM
1749 ! IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1750 ! DLNSIG = LOG(PRSL(I,K) / PRSL(I,K-1))
1751 ! ZO(I,k) = ZO(I,k-1) - DLNSIG * RD / G
1752 ! & * .5 * (TVO(I,k) + TVO(I,k-1))
1753 ! ENDIF
1754 ! ENDDO
1755 ! ENDDO
1756 !
1757 !--- MOIST STATIC ENERGY
1758 !
1759 DO K = 1, KM1
1760 DO I = 1, IM
1761 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1762 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
1763 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
1764 !jfe ES = 10. * FPVS(TO(I,k+1))
1765 !
1766 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
1767 !
1768 PPRIME = PFLD(I,k+1) + EPSM1 * ES
1769 QS = EPS * ES / PPRIME
1770 DQSDP = - QS / PPRIME
1771 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
1772 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
1773 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
1774 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
1775 DQ = DQSDT * DT + DQSDP * DP
1776 TO(I,k) = TO(I,k+1) + DT
1777 QO(I,k) = QO(I,k+1) + DQ
1778 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
1779 ENDIF
1780 ENDDO
1781 ENDDO
1782 DO K = 1, KM1
1783 DO I = 1, IM
1784 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1785 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1786 !
1787 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1788 !
1789 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1 * QESO(I,k))
1790 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1791 val1 = 1.E-8
1792 QESO(I,k) = MAX(QESO(I,k), val1)
1793 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1794 val2 = 1.e-10
1795 QO(I,k) = max(QO(I,k), val2 )
1796 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
1797 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1798 & CP * TO(I,k) + HVAP * QO(I,k)
1799 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1800 & CP * TO(I,k) + HVAP * QESO(I,k)
1801 ENDIF
1802 ENDDO
1803 ENDDO
1804 DO I = 1, IM
1805 k = kmax(i)
1806 IF(CNVFLG(I)) THEN
1807 HEO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QO(I,k)
1808 HESO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QESO(I,k)
1809 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
1810 ENDIF
1811 ENDDO
1812 DO I = 1, IM
1813 IF(CNVFLG(I)) THEN
1814 INDX = KB(I)
1815 XHKB(I) = HEO(I,INDX)
1816 XQKB(I) = QO(I,INDX)
1817 HCKO(I,INDX) = XHKB(I)
1818 QCKO(I,INDX) = XQKB(I)
1819 ENDIF
1820 ENDDO
1821 !
1822 !
1823 !**************************** STATIC CONTROL
1824 !
1825 !
1826 !------- MOISTURE AND CLOUD WORK FUNCTIONS
1827 !
1828 DO K = 2, KM1
1829 DO I = 1, IM
1830 if (k .le. kmax(i)-1) then
1831 ! IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1832 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KTCON(I)) THEN
1833 FACTOR = ETA(I,k-1) / ETA(I,k)
1834 ONEMF = 1. - FACTOR
1835 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1836 & .5 * (HEO(I,k) + HEO(I,k+1))
1837 ENDIF
1838 ! IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1839 ! HEO(I,k) = HEO(I,k-1)
1840 ! ENDIF
1841 endif
1842 ENDDO
1843 ENDDO
1844 DO K = 2, KM1
1845 DO I = 1, IM
1846 if (k .le. kmax(i)-1) then
1847 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1848 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1849 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1850 XDBY = HCKO(I,k) - HESO(I,k)
1851 !cmr XDBY = MAX(XDBY,0.)
1852 val = 0.
1853 XDBY = MAX(XDBY,val)
1854 XQRCH = QESO(I,k) &
1855 & + GAMMA * XDBY / (HVAP * (1. + GAMMA))
1856 FACTOR = ETA(I,k-1) / ETA(I,k)
1857 ONEMF = 1. - FACTOR
1858 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1859 & .5 * (QO(I,k) + QO(I,k+1))
1860 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * XQRCH
1861 IF(DQ.GT.0.) THEN
1862 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1863 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1864 XAA0(I) = XAA0(I) - (ZO(I,k) - ZO(I,k-1)) * G * QLK
1865 XQC = QLK + XQRCH
1866 XPW = ETAH * C0 * DZ * QLK
1867 QCKO(I,k) = XQC
1868 XPWAV(I) = XPWAV(I) + XPW
1869 ENDIF
1870 ENDIF
1871 ! IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LT.KTCON(I)) THEN
1872 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1873 DZ1 = ZO(I,k) - ZO(I,k-1)
1874 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1875 RFACT = 1. + DELTA * CP * GAMMA &
1876 & * TO(I,k-1) / HVAP
1877 XDBY = HCKO(I,k-1) - HESO(I,k-1)
1878 XAA0(I) = XAA0(I) &
1879 & + DZ1 * (G / (CP * TO(I,k-1))) &
1880 & * XDBY / (1. + GAMMA) &
1881 & * RFACT
1882 val=0.
1883 XAA0(I)=XAA0(I)+ &
1884 & DZ1 * G * DELTA * &
1885 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1886 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1887 ENDIF
1888 endif
1889 ENDDO
1890 ENDDO
1891 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1892 !cccc PRINT *, ' XAA BEFORE DWNDRFT =', XAA0(I)
1893 !cccc ENDIF
1894 !
1895 !------- DOWNDRAFT CALCULATIONS
1896 !
1897 !
1898 !--- DOWNDRAFT MOISTURE PROPERTIES
1899 !
1900 DO I = 1, IM
1901 XPWEV(I) = 0.
1902 ENDDO
1903 DO I = 1, IM
1904 IF(DWNFLG2(I)) THEN
1905 JMN = JMIN(I)
1906 XHCD(I) = HEO(I,JMN)
1907 XQCD(I) = QO(I,JMN)
1908 QRCD(I,JMN) = QESO(I,JMN)
1909 ENDIF
1910 ENDDO
1911 DO K = KM1, 1, -1
1912 DO I = 1, IM
1913 if (k .le. kmax(i)-1) then
1914 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1915 DQ = QESO(I,k)
1916 DT = TO(I,k)
1917 GAMMA = EL2ORC * DQ / DT**2
1918 DH = XHCD(I) - HESO(I,k)
1919 QRCD(I,k)=DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1920 DETAD = ETAD(I,k+1) - ETAD(I,k)
1921 XPWD = ETAD(I,k+1) * QRCD(I,k+1) - &
1922 & ETAD(I,k) * QRCD(I,k)
1923 XPWD = XPWD - DETAD * &
1924 & .5 * (QRCD(I,k) + QRCD(I,k+1))
1925 XPWEV(I) = XPWEV(I) + XPWD
1926 ENDIF
1927 endif
1928 ENDDO
1929 ENDDO
1930 !
1931 DO I = 1, IM
1932 edtmax = edtmaxl
1933 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1934 IF(DWNFLG2(I)) THEN
1935 IF(XPWEV(I).GE.0.) THEN
1936 EDTX(I) = 0.
1937 ELSE
1938 EDTX(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
1939 EDTX(I) = MIN(EDTX(I),EDTMAX)
1940 ENDIF
1941 ELSE
1942 EDTX(I) = 0.
1943 ENDIF
1944 ENDDO
1945 !
1946 !
1947 !
1948 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1949 !
1950 !
1951 DO K = KM1, 1, -1
1952 DO I = 1, IM
1953 if (k .le. kmax(i)-1) then
1954 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1955 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1956 DHH=XHCD(I)
1957 DT= TO(I,k+1)
1958 DG= GAMMA
1959 DH= HESO(I,k+1)
1960 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1961 XAA0(I)=XAA0(I)+EDTX(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1962 & *(1.+DELTA*CP*DG*DT/HVAP)
1963 val=0.
1964 XAA0(I)=XAA0(I)+EDTX(I)* &
1965 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
1966 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1967 ENDIF
1968 endif
1969 ENDDO
1970 ENDDO
1971 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1972 !cccc PRINT *, ' XAA AFTER DWNDRFT =', XAA0(I)
1973 !cccc ENDIF
1974 !
1975 ! CALCULATE CRITICAL CLOUD WORK FUNCTION
1976 !
1977 DO I = 1, IM
1978 ACRT(I) = 0.
1979 IF(CNVFLG(I)) THEN
1980 ! IF(CNVFLG(I).AND.SLIMSK(I).NE.1.) THEN
1981 IF(PFLD(I,KTCON(I)).LT.PCRIT(15))THEN
1982 ACRT(I)=ACRIT(15)*(975.-PFLD(I,KTCON(I))) &
1983 & /(975.-PCRIT(15))
1984 ELSE IF(PFLD(I,KTCON(I)).GT.PCRIT(1))THEN
1985 ACRT(I)=ACRIT(1)
1986 ELSE
1987 !cmr K = IFIX((850. - PFLD(I,KTCON(I)))/50.) + 2
1988 K = int((850. - PFLD(I,KTCON(I)))/50.) + 2
1989 K = MIN(K,15)
1990 K = MAX(K,2)
1991 ACRT(I)=ACRIT(K)+(ACRIT(K-1)-ACRIT(K))* &
1992 & (PFLD(I,KTCON(I))-PCRIT(K))/(PCRIT(K-1)-PCRIT(K))
1993 ENDIF
1994 ! ELSE
1995 ! ACRT(I) = .5 * (PFLD(I,KBCON(I)) - PFLD(I,KTCON(I)))
1996 ENDIF
1997 ENDDO
1998 DO I = 1, IM
1999 ACRTFCT(I) = 1.
2000 IF(CNVFLG(I)) THEN
2001 if(SLIMSK(I).eq.1.) THEN
2002 w1 = w1l
2003 w2 = w2l
2004 w3 = w3l
2005 w4 = w4l
2006 else
2007 w1 = w1s
2008 w2 = w2s
2009 w3 = w3s
2010 w4 = w4s
2011 ENDIF
2012 !C IF(CNVFLG(I).AND.SLIMSK(I).EQ.1.) THEN
2013 ! ACRTFCT(I) = PDOT(I) / W3
2014 !
2015 ! modify critical cloud workfunction by cloud base vertical velocity
2016 !
2017 IF(PDOT(I).LE.W4) THEN
2018 ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
2019 ELSEIF(PDOT(I).GE.-W4) THEN
2020 ACRTFCT(I) = - (PDOT(I) + W4) / (W4 - W3)
2021 ELSE
2022 ACRTFCT(I) = 0.
2023 ENDIF
2024 !cmr ACRTFCT(I) = MAX(ACRTFCT(I),-1.)
2025 val1 = -1.
2026 ACRTFCT(I) = MAX(ACRTFCT(I),val1)
2027 !cmr ACRTFCT(I) = MIN(ACRTFCT(I),1.)
2028 val2 = 1.
2029 ACRTFCT(I) = MIN(ACRTFCT(I),val2)
2030 ACRTFCT(I) = 1. - ACRTFCT(I)
2031 !
2032 ! modify ACRTFCT(I) by colume mean rh if RHBAR(I) is greater than 80 percent
2033 !
2034 ! if(RHBAR(I).ge..8) THEN
2035 ! ACRTFCT(I) = ACRTFCT(I) * (.9 - min(RHBAR(I),.9)) * 10.
2036 ! ENDIF
2037 !
2038 ! modify adjustment time scale by cloud base vertical velocity
2039 !
2040 DTCONV(I) = DT2 + max((1800. - DT2),RZERO) * &
2041 & (PDOT(I) - W2) / (W1 - W2)
2042 ! DTCONV(I) = MAX(DTCONV(I), DT2)
2043 ! DTCONV(I) = 1800. * (PDOT(I) - w2) / (w1 - w2)
2044 DTCONV(I) = max(DTCONV(I),dtmin)
2045 DTCONV(I) = min(DTCONV(I),dtmax)
2047 ENDIF
2048 ENDDO
2049 !
2050 !--- LARGE SCALE FORCING
2051 !
2052 DO I= 1, IM
2053 FLG(I) = CNVFLG(I)
2054 IF(CNVFLG(I)) THEN
2055 ! F = AA1(I) / DTCONV(I)
2056 FLD(I) = (AA1(I) - ACRT(I) * ACRTFCT(I)) / DTCONV(I)
2057 IF(FLD(I).LE.0.) FLG(I) = .FALSE.
2058 ENDIF
2059 CNVFLG(I) = FLG(I)
2060 IF(CNVFLG(I)) THEN
2061 ! XAA0(I) = MAX(XAA0(I),0.)
2062 XK(I) = (XAA0(I) - AA1(I)) / MBDT
2063 IF(XK(I).GE.0.) FLG(I) = .FALSE.
2064 ENDIF
2065 !
2066 !--- KERNEL, CLOUD BASE MASS FLUX
2067 !
2068 CNVFLG(I) = FLG(I)
2069 IF(CNVFLG(I)) THEN
2070 XMB(I) = -FLD(I) / XK(I)
2071 XMB(I) = MIN(XMB(I),XMBMAX(I))
2072 ENDIF
2073 ENDDO
2074 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
2075 ! print *, ' RHBAR(I), ACRTFCT(I) =', RHBAR(I), ACRTFCT(I)
2076 ! PRINT *, ' A1, XA =', AA1(I), XAA0(I)
2077 ! PRINT *, ' XMB(I), ACRT =', XMB(I), ACRT
2078 ! ENDIF
2079 TOTFLG = .TRUE.
2080 DO I = 1, IM
2081 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
2082 ENDDO
2083 IF(TOTFLG) RETURN
2084 !
2085 ! restore t0 and QO to t1 and q1 in case convection stops
2086 !
2087 do k = 1, km
2088 DO I = 1, IM
2089 if (k .le. kmax(i)) then
2090 TO(I,k) = T1(I,k)
2091 QO(I,k) = Q1(I,k)
2092 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
2093 !
2094 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
2095 !
2096 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
2097 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
2098 val = 1.E-8
2099 QESO(I,k) = MAX(QESO(I,k), val )
2100 endif
2101 enddo
2102 enddo
2103 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2104 !
2105 !--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
2106 !--- MULTIPLIED BY THE MASS FLUX NECESSARY TO KEEP THE
2107 !--- EQUILIBRIUM WITH THE LARGER-SCALE.
2108 !
2109 DO I = 1, IM
2110 DELHBAR(I) = 0.
2111 DELQBAR(I) = 0.
2112 DELTBAR(I) = 0.
2113 QCOND(I) = 0.
2114 ENDDO
2115 DO K = 1, KM
2116 DO I = 1, IM
2117 if (k .le. kmax(i)) then
2118 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2119 AUP = 1.
2120 IF(K.Le.KB(I)) AUP = 0.
2121 ADW = 1.
2122 IF(K.GT.JMIN(I)) ADW = 0.
2123 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
2124 T1(I,k) = T1(I,k) + DELLAT * XMB(I) * DT2
2125 Q1(I,k) = Q1(I,k) + DELLAQ(I,k) * XMB(I) * DT2
2126 U1(I,k) = U1(I,k) + DELLAU(I,k) * XMB(I) * DT2
2127 V1(I,k) = V1(I,k) + DELLAV(I,k) * XMB(I) * DT2
2128 DP = 1000. * DEL(I,K)
2129 DELHBAR(I) = DELHBAR(I) + DELLAH(I,k)*XMB(I)*DP/G
2130 DELQBAR(I) = DELQBAR(I) + DELLAQ(I,k)*XMB(I)*DP/G
2131 DELTBAR(I) = DELTBAR(I) + DELLAT*XMB(I)*DP/G
2132 ENDIF
2133 endif
2134 ENDDO
2135 ENDDO
2136 DO K = 1, KM
2137 DO I = 1, IM
2138 if (k .le. kmax(i)) then
2139 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2140 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
2141 !
2142 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
2143 !
2144 QESO(I,k) = EPS * QESO(I,k)/(PFLD(I,k) + EPSM1*QESO(I,k))
2145 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
2146 val = 1.E-8
2147 QESO(I,k) = MAX(QESO(I,k), val )
2148 !
2149 ! cloud water
2150 !
2151 if(ncloud.gt.0.and.CNVFLG(I).and.k.eq.KTCON(I)) THEN
2152 tem = DELLAL(I) * XMB(I) * dt2
2153 tem1 = MAX(RZERO, MIN(RONE, (TCR-t1(I,K))*TCRF))
2154 if (QL(I,k,2) .gt. -999.0) then
2155 QL(I,k,1) = QL(I,k,1) + tem * tem1 ! Ice
2156 QL(I,k,2) = QL(I,k,2) + tem *(1.0-tem1) ! Water
2157 else
2158 tem2 = QL(I,k,1) + tem
2159 QL(I,k,1) = tem2 * tem1 ! Ice
2160 QL(I,k,2) = tem2 - QL(I,k,1) ! Water
2161 endif
2162 ! QL(I,k) = QL(I,k) + DELLAL(I) * XMB(I) * dt2
2163 dp = 1000. * del(i,k)
2164 DELLAL(I) = DELLAL(I) * XMB(I) * dp / g
2165 ENDIF
2166 ENDIF
2167 endif
2168 ENDDO
2169 ENDDO
2170 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2171 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2172 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2173 ! PRINT *, ' DELLBAR ='
2174 ! PRINT 6003, HVAP*DELLbar
2175 ! PRINT *, ' DELLAQ ='
2176 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2177 ! PRINT *, ' DELLAT ='
2178 ! PRINT 6003, (DELLAH(i,k)*XMB(I)-HVAP*DELLAQ(I,k)*XMB(I), &
2179 ! & K=1,KMAX)
2180 ! ENDIF
2181 DO I = 1, IM
2182 RNTOT(I) = 0.
2183 DELQEV(I) = 0.
2184 DELQ2(I) = 0.
2185 FLG(I) = CNVFLG(I)
2186 ENDDO
2187 DO K = KM, 1, -1
2188 DO I = 1, IM
2189 if (k .le. kmax(i)) then
2190 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2191 AUP = 1.
2192 IF(K.Le.KB(I)) AUP = 0.
2193 ADW = 1.
2194 IF(K.GT.JMIN(I)) ADW = 0.
2195 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2196 RNTOT(I) = RNTOT(I) + rain * XMB(I) * .001 * dt2
2197 ENDIF
2198 endif
2199 ENDDO
2200 ENDDO
2201 DO K = KM, 1, -1
2202 DO I = 1, IM
2203 if (k .le. kmax(i)) then
2204 DELTV(I) = 0.
2205 DELQ(I) = 0.
2206 QEVAP(I) = 0.
2207 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2208 AUP = 1.
2209 IF(K.Le.KB(I)) AUP = 0.
2210 ADW = 1.
2211 IF(K.GT.JMIN(I)) ADW = 0.
2212 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2213 RN(I) = RN(I) + rain * XMB(I) * .001 * dt2
2214 ENDIF
2215 IF(FLG(I).AND.K.LE.KTCON(I)) THEN
2216 evef = EDT(I) * evfact
2217 if(SLIMSK(I).eq.1.) evef=EDT(I) * evfactl
2218 ! if(SLIMSK(I).eq.1.) evef=.07
2219 ! if(SLIMSK(I).ne.1.) evef = 0.
2220 QCOND(I) = EVEF * (Q1(I,k) - QESO(I,k)) &
2221 & / (1. + EL2ORC * QESO(I,k) / T1(I,k)**2)
2222 DP = 1000. * DEL(I,K)
2223 IF(RN(I).GT.0..AND.QCOND(I).LT.0.) THEN
2224 QEVAP(I) = -QCOND(I) * (1.-EXP(-.32*SQRT(DT2*RN(I))))
2225 QEVAP(I) = MIN(QEVAP(I), RN(I)*1000.*G/DP)
2226 DELQ2(I) = DELQEV(I) + .001 * QEVAP(I) * dp / g
2227 ENDIF
2228 if(RN(I).gt.0..and.QCOND(I).LT.0..and. &
2229 & DELQ2(I).gt.RNTOT(I)) THEN
2230 QEVAP(I) = 1000.* g * (RNTOT(I) - DELQEV(I)) / dp
2231 FLG(I) = .false.
2232 ENDIF
2233 IF(RN(I).GT.0..AND.QEVAP(I).gt.0.) THEN
2234 Q1(I,k) = Q1(I,k) + QEVAP(I)
2235 T1(I,k) = T1(I,k) - ELOCP * QEVAP(I)
2236 RN(I) = RN(I) - .001 * QEVAP(I) * DP / G
2237 DELTV(I) = - ELOCP*QEVAP(I)/DT2
2238 DELQ(I) = + QEVAP(I)/DT2
2239 DELQEV(I) = DELQEV(I) + .001*dp*QEVAP(I)/g
2240 ENDIF
2241 DELLAQ(I,k) = DELLAQ(I,k) + DELQ(I) / XMB(I)
2242 DELQBAR(I) = DELQBAR(I) + DELQ(I)*DP/G
2243 DELTBAR(I) = DELTBAR(I) + DELTV(I)*DP/G
2244 ENDIF
2245 endif
2246 ENDDO
2247 ENDDO
2248 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2249 ! PRINT *, ' DELLAH ='
2250 ! PRINT 6003, (DELLAH(k)*XMB(I),K=1,KMAX)
2251 ! PRINT *, ' DELLAQ ='
2252 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2253 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2254 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2255 ! PRINT *, ' PRECIP =', HVAP*RN(I)*1000./DT2
2256 !CCCC PRINT *, ' DELLBAR ='
2257 !CCCC PRINT *, HVAP*DELLbar
2258 ! ENDIF
2259 !
2260 ! PRECIPITATION RATE CONVERTED TO ACTUAL PRECIP
2261 ! IN UNIT OF M INSTEAD OF KG
2262 !
2263 DO I = 1, IM
2264 IF(CNVFLG(I)) THEN
2265 !
2266 ! IN THE EVENT OF UPPER LEVEL RAIN EVAPORATION AND LOWER LEVEL DOWNDRAF
2267 ! MOISTENING, RN CAN BECOME NEGATIVE, IN THIS CASE, WE BACK OUT OF TH
2268 ! HEATING AND THE MOISTENING
2269 !
2270 if(RN(I).lt.0..and..not.FLG(I)) RN(I) = 0.
2271 IF(RN(I).LE.0.) THEN
2272 RN(I) = 0.
2273 ELSE
2274 KTOP(I) = KTCON(I)
2275 KBOT(I) = KBCON(I)
2276 KUO(I) = 1
2277 CLDWRK(I) = AA1(I)
2278 ENDIF
2279 ENDIF
2280 ENDDO
2281 DO K = 1, KM
2282 DO I = 1, IM
2283 if (k .le. kmax(i)) then
2284 IF(CNVFLG(I).AND.RN(I).LE.0.) THEN
2285 T1(I,k) = TO(I,k)
2286 Q1(I,k) = QO(I,k)
2287 ENDIF
2288 endif
2289 ENDDO
2290 ENDDO
2291 !!
2292 RETURN
2293 END SUBROUTINE SASCNV
2295 ! ------------------------------------------------------------------------
2297 SUBROUTINE OLD_ARW_SHALCV(IM,IX,KM,DT,DEL,PRSI,PRSL,PRSLK,KUO,Q,T,DPSHC)
2298 !
2299 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2300 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
2301 &, RD => con_RD
2303 implicit none
2304 !
2305 ! include 'constant.h'
2306 !
2307 integer IM, IX, KM, KUO(IM)
2308 real(kind=kind_phys) DEL(IX,KM), PRSI(IX,KM+1), PRSL(IX,KM), &
2309 & PRSLK(IX,KM), &
2310 & Q(IX,KM), T(IX,KM), DT, DPSHC
2311 !
2312 ! Locals
2313 !
2314 real(kind=kind_phys) ck, cpdt, dmse, dsdz1, dsdz2, &
2315 & dsig, dtodsl, dtodsu, eldq, g, &
2316 & gocp, rtdls
2317 !
2318 integer k,k1,k2,kliftl,kliftu,kt,N2,I,iku,ik1,ik,ii
2319 integer INDEX2(IM), KLCL(IM), KBOT(IM), KTOP(IM),kk &
2320 &, KTOPM(IM)
2321 !!
2322 ! PHYSICAL PARAMETERS
2323 PARAMETER(G=GRAV, GOCP=G/CP)
2324 ! BOUNDS OF PARCEL ORIGIN
2325 PARAMETER(KLIFTL=2,KLIFTU=2)
2326 LOGICAL LSHC(IM)
2327 real(kind=kind_phys) Q2(IM*KM), T2(IM*KM), &
2328 & PRSL2(IM*KM), PRSLK2(IM*KM), &
2329 & AL(IM*(KM-1)), AD(IM*KM), AU(IM*(KM-1))
2330 !-----------------------------------------------------------------------
2331 ! COMPRESS FIELDS TO POINTS WITH NO DEEP CONVECTION
2332 ! AND MOIST STATIC INSTABILITY.
2333 DO I=1,IM
2334 LSHC(I)=.FALSE.
2335 ENDDO
2336 DO K=1,KM-1
2337 DO I=1,IM
2338 IF(KUO(I).EQ.0) THEN
2339 ELDQ = HVAP*(Q(I,K)-Q(I,K+1))
2340 CPDT = CP*(T(I,K)-T(I,K+1))
2341 RTDLS = (PRSL(I,K)-PRSL(I,K+1)) / &
2342 & PRSI(I,K+1)*RD*0.5*(T(I,K)+T(I,K+1))
2343 DMSE = ELDQ+CPDT-RTDLS
2344 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2345 ENDIF
2346 ENDDO
2347 ENDDO
2348 N2 = 0
2349 DO I=1,IM
2350 IF(LSHC(I)) THEN
2351 N2 = N2 + 1
2352 INDEX2(N2) = I
2353 ENDIF
2354 ENDDO
2355 IF(N2.EQ.0) RETURN
2356 DO K=1,KM
2357 KK = (K-1)*N2
2358 DO I=1,N2
2359 IK = KK + I
2360 ii = index2(i)
2361 Q2(IK) = Q(II,K)
2362 T2(IK) = T(II,K)
2363 PRSL2(IK) = PRSL(II,K)
2364 PRSLK2(IK) = PRSLK(II,K)
2365 ENDDO
2366 ENDDO
2367 do i=1,N2
2368 ktopm(i) = KM
2369 enddo
2370 do k=2,KM
2371 do i=1,N2
2372 ii = index2(i)
2373 if (prsi(ii,1)-prsi(ii,k) .le. dpshc) ktopm(i) = k
2374 enddo
2375 enddo
2377 !-----------------------------------------------------------------------
2378 ! COMPUTE MOIST ADIABAT AND DETERMINE LIMITS OF SHALLOW CONVECTION.
2379 ! CHECK FOR MOIST STATIC INSTABILITY AGAIN WITHIN CLOUD.
2380 CALL MSTADBT3(N2,KM-1,KLIFTL,KLIFTU,PRSL2,PRSLK2,T2,Q2, &
2381 & KLCL,KBOT,KTOP,AL,AU)
2382 DO I=1,N2
2383 KBOT(I) = min(KLCL(I)-1, ktopm(i)-1)
2384 KTOP(I) = min(KTOP(I)+1, ktopm(i))
2385 LSHC(I) = .FALSE.
2386 ENDDO
2387 DO K=1,KM-1
2388 KK = (K-1)*N2
2389 DO I=1,N2
2390 IF(K.GE.KBOT(I).AND.K.LT.KTOP(I)) THEN
2391 IK = KK + I
2392 IKU = IK + N2
2393 ELDQ = HVAP * (Q2(IK)-Q2(IKU))
2394 CPDT = CP * (T2(IK)-T2(IKU))
2395 RTDLS = (PRSL2(IK)-PRSL2(IKU)) / &
2396 & PRSI(index2(i),K+1)*RD*0.5*(T2(IK)+T2(IKU))
2397 DMSE = ELDQ + CPDT - RTDLS
2398 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2399 AU(IK) = G/RTDLS
2400 ENDIF
2401 ENDDO
2402 ENDDO
2403 K1=KM+1
2404 K2=0
2405 DO I=1,N2
2406 IF(.NOT.LSHC(I)) THEN
2407 KBOT(I) = KM+1
2408 KTOP(I) = 0
2409 ENDIF
2410 K1 = MIN(K1,KBOT(I))
2411 K2 = MAX(K2,KTOP(I))
2412 ENDDO
2413 KT = K2-K1+1
2414 IF(KT.LT.2) RETURN
2415 !-----------------------------------------------------------------------
2416 ! SET EDDY VISCOSITY COEFFICIENT CKU AT SIGMA INTERFACES.
2417 ! COMPUTE DIAGONALS AND RHS FOR TRIDIAGONAL MATRIX SOLVER.
2418 ! EXPAND FINAL FIELDS.
2419 KK = (K1-1) * N2
2420 DO I=1,N2
2421 IK = KK + I
2422 AD(IK) = 1.
2423 ENDDO
2424 !
2425 ! DTODSU=DT/DEL(K1)
2426 DO K=K1,K2-1
2427 ! DTODSL=DTODSU
2428 ! DTODSU= DT/DEL(K+1)
2429 ! DSIG=SL(K)-SL(K+1)
2430 KK = (K-1) * N2
2431 DO I=1,N2
2432 ii = index2(i)
2433 DTODSL = DT/DEL(II,K)
2434 DTODSU = DT/DEL(II,K+1)
2435 DSIG = PRSL(II,K) - PRSL(II,K+1)
2436 IK = KK + I
2437 IKU = IK + N2
2438 IF(K.EQ.KBOT(I)) THEN
2439 CK=1.5
2440 ELSEIF(K.EQ.KTOP(I)-1) THEN
2441 CK=1.
2442 ELSEIF(K.EQ.KTOP(I)-2) THEN
2443 CK=3.
2444 ELSEIF(K.GT.KBOT(I).AND.K.LT.KTOP(I)-2) THEN
2445 CK=5.
2446 ELSE
2447 CK=0.
2448 ENDIF
2449 DSDZ1 = CK*DSIG*AU(IK)*GOCP
2450 DSDZ2 = CK*DSIG*AU(IK)*AU(IK)
2451 AU(IK) = -DTODSL*DSDZ2
2452 AL(IK) = -DTODSU*DSDZ2
2453 AD(IK) = AD(IK)-AU(IK)
2454 AD(IKU) = 1.-AL(IK)
2455 T2(IK) = T2(IK)+DTODSL*DSDZ1
2456 T2(IKU) = T2(IKU)-DTODSU*DSDZ1
2457 ENDDO
2458 ENDDO
2459 IK1=(K1-1)*N2+1
2460 CALL TRIDI2T3(N2,KT,AL(IK1),AD(IK1),AU(IK1),Q2(IK1),T2(IK1), &
2461 & AU(IK1),Q2(IK1),T2(IK1))
2462 DO K=K1,K2
2463 KK = (K-1)*N2
2464 DO I=1,N2
2465 IK = KK + I
2466 Q(INDEX2(I),K) = Q2(IK)
2467 T(INDEX2(I),K) = T2(IK)
2468 ENDDO
2469 ENDDO
2470 !-----------------------------------------------------------------------
2471 RETURN
2472 END SUBROUTINE OLD_ARW_SHALCV
2474 !-----------------------------------------------------------------------
2475 SUBROUTINE TRIDI2T3(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
2476 !yt INCLUDE DBTRIDI2;
2477 !!
2478 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2479 implicit none
2480 integer k,n,l,i
2481 real(kind=kind_phys) fk
2482 !!
2483 real(kind=kind_phys) &
2484 & CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
2485 & AU(L,N-1),A1(L,N),A2(L,N)
2486 !-----------------------------------------------------------------------
2487 DO I=1,L
2488 FK=1./CM(I,1)
2489 AU(I,1)=FK*CU(I,1)
2490 A1(I,1)=FK*R1(I,1)
2491 A2(I,1)=FK*R2(I,1)
2492 ENDDO
2493 DO K=2,N-1
2494 DO I=1,L
2495 FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
2496 AU(I,K)=FK*CU(I,K)
2497 A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
2498 A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
2499 ENDDO
2500 ENDDO
2501 DO I=1,L
2502 FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
2503 A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
2504 A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
2505 ENDDO
2506 DO K=N-1,1,-1
2507 DO I=1,L
2508 A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
2509 A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
2510 ENDDO
2511 ENDDO
2512 !-----------------------------------------------------------------------
2513 RETURN
2514 END SUBROUTINE TRIDI2T3
2515 !-----------------------------------------------------------------------
2517 SUBROUTINE MSTADBT3(IM,KM,K1,K2,PRSL,PRSLK,TENV,QENV, &
2518 & KLCL,KBOT,KTOP,TCLD,QCLD)
2519 !yt INCLUDE DBMSTADB;
2520 !!
2521 USE MODULE_GFS_MACHINE, ONLY : kind_phys
2522 USE MODULE_GFS_FUNCPHYS, ONLY : FTDP, FTHE, FTLCL, STMA
2523 USE MODULE_GFS_PHYSCONS, EPS => con_eps, EPSM1 => con_epsm1, FV => con_FVirt
2525 implicit none
2526 !!
2527 ! include 'constant.h'
2528 !!
2529 integer k,k1,k2,km,i,im
2530 real(kind=kind_phys) pv,qma,slklcl,tdpd,thelcl,tlcl
2531 real(kind=kind_phys) tma,tvcld,tvenv
2532 !!
2533 real(kind=kind_phys) PRSL(IM,KM), PRSLK(IM,KM), TENV(IM,KM), &
2534 & QENV(IM,KM), TCLD(IM,KM), QCLD(IM,KM)
2535 INTEGER KLCL(IM), KBOT(IM), KTOP(IM)
2536 ! LOCAL ARRAYS
2537 real(kind=kind_phys) SLKMA(IM), THEMA(IM)
2538 !-----------------------------------------------------------------------
2539 ! DETERMINE WARMEST POTENTIAL WET-BULB TEMPERATURE BETWEEN K1 AND K2.
2540 ! COMPUTE ITS LIFTING CONDENSATION LEVEL.
2541 !
2542 DO I=1,IM
2543 SLKMA(I) = 0.
2544 THEMA(I) = 0.
2545 ENDDO
2546 DO K=K1,K2
2547 DO I=1,IM
2548 PV = 1000.0 * PRSL(I,K)*QENV(I,K)/(EPS-EPSM1*QENV(I,K))
2549 TDPD = TENV(I,K)-FTDP(PV)
2550 IF(TDPD.GT.0.) THEN
2551 TLCL = FTLCL(TENV(I,K),TDPD)
2552 SLKLCL = PRSLK(I,K)*TLCL/TENV(I,K)
2553 ELSE
2554 TLCL = TENV(I,K)
2555 SLKLCL = PRSLK(I,K)
2556 ENDIF
2557 THELCL=FTHE(TLCL,SLKLCL)
2558 IF(THELCL.GT.THEMA(I)) THEN
2559 SLKMA(I) = SLKLCL
2560 THEMA(I) = THELCL
2561 ENDIF
2562 ENDDO
2563 ENDDO
2564 !-----------------------------------------------------------------------
2565 ! SET CLOUD TEMPERATURES AND HUMIDITIES WHEREVER THE PARCEL LIFTED UP
2566 ! THE MOIST ADIABAT IS BUOYANT WITH RESPECT TO THE ENVIRONMENT.
2567 DO I=1,IM
2568 KLCL(I)=KM+1
2569 KBOT(I)=KM+1
2570 KTOP(I)=0
2571 ENDDO
2572 DO K=1,KM
2573 DO I=1,IM
2574 TCLD(I,K)=0.
2575 QCLD(I,K)=0.
2576 ENDDO
2577 ENDDO
2578 DO K=K1,KM
2579 DO I=1,IM
2580 IF(PRSLK(I,K).LE.SLKMA(I)) THEN
2581 KLCL(I)=MIN(KLCL(I),K)
2582 CALL STMA(THEMA(I),PRSLK(I,K),TMA,QMA)
2583 ! TMA=FTMA(THEMA(I),PRSLK(I,K),QMA)
2584 TVCLD=TMA*(1.+FV*QMA)
2585 TVENV=TENV(I,K)*(1.+FV*QENV(I,K))
2586 IF(TVCLD.GT.TVENV) THEN
2587 KBOT(I)=MIN(KBOT(I),K)
2588 KTOP(I)=MAX(KTOP(I),K)
2589 TCLD(I,K)=TMA-TENV(I,K)
2590 QCLD(I,K)=QMA-QENV(I,K)
2591 ENDIF
2592 ENDIF
2593 ENDDO
2594 ENDDO
2595 !-----------------------------------------------------------------------
2596 RETURN
2597 END SUBROUTINE MSTADBT3
2599 subroutine sascnvn(im,ix,km,jcap,delt,del,prsl,ps,phil,ql, &
2600 & q1,t1,u1,v1,rcs,cldwrk,rn,kbot,ktop,kcnv,slimsk, &
2601 & dot,ncloud,pgcon,sas_mass_flux, &
2602 & pert_sas, ens_random_seed, ens_sasamp)
2603 ! & dot,ncloud,ud_mf,dd_mf,dt_mf)
2604 ! & dot,ncloud,ud_mf,dd_mf,dt_mf,me)
2605 !
2606 ! use machine , only : kind_phys
2607 ! use funcphys , only : fpvs
2608 ! use physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
2609 USE MODULE_GFS_MACHINE, ONLY : kind_phys
2610 USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
2611 USE MODULE_GFS_PHYSCONS, grav => con_g, cp => con_cp &
2612 &, hvap => con_hvap &
2613 &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
2614 &, cvap => con_cvap, cliq => con_cliq &
2615 &, eps => con_eps, epsm1 => con_epsm1
2616 implicit none
2617 !
2618 integer im, ix, km, jcap, ncloud, &
2619 & kbot(im), ktop(im), kcnv(im)
2620 ! &, me
2621 real(kind=kind_phys) delt,sas_mass_flux
2622 real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
2623 & ql(ix,km,2),q1(ix,km), t1(ix,km), &
2624 & u1(ix,km), v1(ix,km), rcs(im), &
2625 & cldwrk(im), rn(im), slimsk(im), &
2626 & dot(ix,km), phil(ix,km)
2627 ! hchuang code change mass flux output
2628 ! &, ud_mf(im,km),dd_mf(im,km),dt_mf(im,km)
2629 !
2630 integer i, j, indx, jmn, k, kk, latd, lond, km1
2631 !
2632 real(kind=kind_phys) clam, cxlamu, xlamde, xlamdd
2633 !
2634 real(kind=kind_phys) adw, aup, aafac, &
2635 & beta, betal, betas, &
2636 & c0, cpoel, dellat, delta, &
2637 & desdt, deta, detad, dg, &
2638 & dh, dhh, dlnsig, dp, &
2639 & dq, dqsdp, dqsdt, dt, &
2640 & dt2, dtmax, dtmin, dv1h, &
2641 & dv1q, dv2h, dv2q, dv1u, &
2642 & dv1v, dv2u, dv2v, dv3q, &
2643 & dv3h, dv3u, dv3v, &
2644 & dz, dz1, e1, edtmax, &
2645 & edtmaxl, edtmaxs, el2orc, elocp, &
2646 & es, etah, cthk, dthk, &
2647 & evef, evfact, evfactl, fact1, &
2648 & fact2, factor, fjcap, fkm, &
2649 & g, gamma, pprime, &
2650 & qlk, qrch, qs, c1, &
2651 & rain, rfact, shear, tem1, &
2652 & tem2, terr, val, val1, &
2653 & val2, w1, w1l, w1s, &
2654 & w2, w2l, w2s, w3, &
2655 & w3l, w3s, w4, w4l, &
2656 & w4s, xdby, xpw, xpwd, &
2657 & xqrch, mbdt, tem, &
2658 & ptem, ptem1
2659 !
2660 real(kind=kind_phys), intent(in) :: pgcon
2661 logical,intent(in) :: pert_sas
2662 integer,intent(in) :: ens_random_seed
2663 real,intent(in) :: ens_sasamp
2665 integer kb(im), kbcon(im), kbcon1(im), &
2666 & ktcon(im), ktcon1(im), &
2667 & jmin(im), lmin(im), kbmax(im), &
2668 & kbm(im), kmax(im)
2669 !
2670 real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im), &
2671 & delhbar(im), delq(im), delq2(im), &
2672 & delqbar(im), delqev(im), deltbar(im), &
2673 & deltv(im), dtconv(im), edt(im), &
2674 & edto(im), edtx(im), fld(im), &
2675 & hcdo(im,km), hmax(im), hmin(im), &
2676 & ucdo(im,km), vcdo(im,km),aa2(im), &
2677 & pbcdif(im), pdot(im), po(im,km), &
2678 & pwavo(im), pwevo(im), xlamud(im), &
2679 & qcdo(im,km), qcond(im), qevap(im), &
2680 & rntot(im), vshear(im), xaa0(im), &
2681 & xk(im), xlamd(im), &
2682 & xmb(im), xmbmax(im), xpwav(im), &
2683 & xpwev(im), delubar(im),delvbar(im)
2684 !cj
2685 real(kind=kind_phys) cincr, cincrmax, cincrmin
2686 real(kind=kind_phys) xmbmx1
2687 !cj
2688 !c physical parameters
2689 parameter(g=grav)
2690 parameter(cpoel=cp/hvap,elocp=hvap/cp, &
2691 & el2orc=hvap*hvap/(rv*cp))
2692 parameter(terr=0.,c0=.002,c1=.002,delta=fv)
2693 parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
2694 parameter(cthk=150.,cincrmax=180.,cincrmin=120.,dthk=25.)
2695 !c local variables and arrays
2696 real(kind=kind_phys) pfld(im,km),to(im,km), qo(im,km), &
2697 & uo(im,km), vo(im,km), qeso(im,km)
2698 !c cloud water
2699 real(kind=kind_phys)qlko_ktcon(im),dellal(im,km),tvo(im,km), &
2700 & dbyo(im,km), zo(im,km), xlamue(im,km), &
2701 & fent1(im,km),fent2(im,km), frh(im,km), &
2702 & heo(im,km), heso(im,km), &
2703 & qrcd(im,km), dellah(im,km), dellaq(im,km), &
2704 & dellau(im,km),dellav(im,km), hcko(im,km), &
2705 & ucko(im,km), vcko(im,km), qcko(im,km), &
2706 & eta(im,km), etad(im,km), zi(im,km), &
2707 & qrcdo(im,km),pwo(im,km), pwdo(im,km), &
2708 & tx1(im), sumx(im)
2709 ! &, rhbar(im)
2710 !
2711 logical totflg, cnvflg(im), flg(im)
2712 !
2713 real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
2714 ! save pcrit, acritt
2715 data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,&
2716 & 350.,300.,250.,200.,150./
2717 data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
2718 & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
2719 !c gdas derived acrit
2720 !c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
2721 !c & .743,.813,.886,.947,1.138,1.377,1.896/
2722 real(kind=kind_phys) tf, tcr, tcrf
2723 parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
2724 #if HWRF==1
2725 real*8 :: gasdev,ran1 !zhang
2726 real :: rr !zhang
2727 logical,save :: pert_sas_local !zhang
2728 integer,save :: ens_random_seed_local !zhang
2729 real,save :: ens_sasamp_local !zhang
2730 data ens_random_seed_local/0/
2731 if ( ens_random_seed_local .eq. 0 ) then
2732 ens_random_seed_local=ens_random_seed
2733 pert_sas_local=pert_sas
2734 ens_sasamp_local=ens_sasamp
2735 endif
2736 #endif
2737 !
2738 !c-----------------------------------------------------------------------
2739 !
2741 km1 = km - 1
2742 !c
2743 !c initialize arrays
2744 !c
2745 do i=1,im
2746 kcnv(i)=0
2747 cnvflg(i) = .true.
2748 rn(i)=0.
2749 kbot(i)=km+1
2750 ktop(i)=0
2751 kbcon(i)=km
2752 ktcon(i)=1
2753 dtconv(i) = 3600.
2754 cldwrk(i) = 0.
2755 pdot(i) = 0.
2756 pbcdif(i)= 0.
2757 lmin(i) = 1
2758 jmin(i) = 1
2759 qlko_ktcon(i) = 0.
2760 edt(i) = 0.
2761 edto(i) = 0.
2762 edtx(i) = 0.
2763 acrt(i) = 0.
2764 acrtfct(i) = 1.
2765 aa1(i) = 0.
2766 aa2(i) = 0.
2767 xaa0(i) = 0.
2768 pwavo(i)= 0.
2769 pwevo(i)= 0.
2770 xpwav(i)= 0.
2771 xpwev(i)= 0.
2772 vshear(i) = 0.
2773 enddo
2774 ! hchuang code change
2775 ! do k = 1, km
2776 ! do i = 1, im
2777 ! ud_mf(i,k) = 0.
2778 ! dd_mf(i,k) = 0.
2779 ! dt_mf(i,k) = 0.
2780 ! enddo
2781 ! enddo
2782 !c
2783 do k = 1, 15
2784 acrit(k) = acritt(k) * (975. - pcrit(k))
2785 enddo
2786 dt2 = delt
2787 val = 1200.
2788 dtmin = max(dt2, val )
2789 val = 3600.
2790 dtmax = max(dt2, val )
2791 !c model tunable parameters are all here
2792 mbdt = 10.
2793 edtmaxl = .3
2794 edtmaxs = .3
2795 clam = .1
2796 aafac = .1
2797 ! betal = .15
2798 ! betas = .15
2799 betal = .05
2800 betas = .05
2801 !c evef = 0.07
2802 evfact = 0.3
2803 evfactl = 0.3
2804 #if ( EM_CORE == 1 )
2805 ! HAWAII TEST - ZCX
2806 BETAl = .05
2807 betas = .05
2808 evfact = 0.5
2809 evfactl = 0.5
2810 #endif
2811 !
2812 cxlamu = 1.0e-4
2813 xlamde = 1.0e-4
2814 xlamdd = 1.0e-4
2815 !
2816 fjcap = (float(jcap) / 126.) ** 2
2817 val = 1.
2818 fjcap = max(fjcap,val)
2819 fkm = (float(km) / 28.) ** 2
2820 fkm = max(fkm,val)
2821 w1l = -8.e-3
2822 w2l = -4.e-2
2823 w3l = -5.e-3
2824 w4l = -5.e-4
2825 w1s = -2.e-4
2826 w2s = -2.e-3
2827 w3s = -1.e-3
2828 w4s = -2.e-5
2829 !c
2830 !c define top layer for search of the downdraft originating layer
2831 !c and the maximum thetae for updraft
2832 !c
2833 do i=1,im
2834 kbmax(i) = km
2835 kbm(i) = km
2836 kmax(i) = km
2837 tx1(i) = 1.0 / ps(i)
2838 enddo
2839 !
2840 do k = 1, km
2841 do i=1,im
2842 IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
2843 !2011bugfix if (prsl(i,k)*tx1(i) .gt. 0.04) kmax(i) = k + 1
2844 if (prsl(i,k)*tx1(i) .gt. 0.45) kbmax(i) = k + 1
2845 if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
2846 enddo
2847 enddo
2848 do i=1,im
2849 kbmax(i) = min(kbmax(i),kmax(i))
2850 kbm(i) = min(kbm(i),kmax(i))
2851 enddo
2852 !c
2853 !c hydrostatic height assume zero terr and initially assume
2854 !c updraft entrainment rate as an inverse function of height
2855 !c
2856 do k = 1, km
2857 do i=1,im
2858 zo(i,k) = phil(i,k) / g
2859 enddo
2860 enddo
2861 do k = 1, km1
2862 do i=1,im
2863 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
2864 xlamue(i,k) = clam / zi(i,k)
2865 enddo
2866 enddo
2867 !c
2868 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2869 !c convert surface pressure to mb from cb
2870 !c
2871 do k = 1, km
2872 do i = 1, im
2873 if (k .le. kmax(i)) then
2874 pfld(i,k) = prsl(i,k) * 10.0
2875 eta(i,k) = 1.
2876 fent1(i,k)= 1.
2877 fent2(i,k)= 1.
2878 frh(i,k) = 0.
2879 hcko(i,k) = 0.
2880 qcko(i,k) = 0.
2881 ucko(i,k) = 0.
2882 vcko(i,k) = 0.
2883 etad(i,k) = 1.
2884 hcdo(i,k) = 0.
2885 qcdo(i,k) = 0.
2886 ucdo(i,k) = 0.
2887 vcdo(i,k) = 0.
2888 qrcd(i,k) = 0.
2889 qrcdo(i,k)= 0.
2890 dbyo(i,k) = 0.
2891 pwo(i,k) = 0.
2892 pwdo(i,k) = 0.
2893 dellal(i,k) = 0.
2894 to(i,k) = t1(i,k)
2895 qo(i,k) = q1(i,k)
2896 uo(i,k) = u1(i,k) * rcs(i)
2897 vo(i,k) = v1(i,k) * rcs(i)
2898 endif
2899 enddo
2900 enddo
2901 !c
2902 !c column variables
2903 !c p is pressure of the layer (mb)
2904 !c t is temperature at t-dt (k)..tn
2905 !c q is mixing ratio at t-dt (kg/kg)..qn
2906 !c to is temperature at t+dt (k)... this is after advection and turbulan
2907 !c qo is mixing ratio at t+dt (kg/kg)..q1
2908 !c
2909 do k = 1, km
2910 do i=1,im
2911 if (k .le. kmax(i)) then
2912 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
2913 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
2914 val1 = 1.e-8
2915 qeso(i,k) = max(qeso(i,k), val1)
2916 val2 = 1.e-10
2917 qo(i,k) = max(qo(i,k), val2 )
2918 ! qo(i,k) = min(qo(i,k),qeso(i,k))
2919 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
2920 endif
2921 enddo
2922 enddo
2923 !c
2924 !c compute moist static energy
2925 !c
2926 do k = 1, km
2927 do i=1,im
2928 if (k .le. kmax(i)) then
2929 ! tem = g * zo(i,k) + cp * to(i,k)
2930 tem = phil(i,k) + cp * to(i,k)
2931 heo(i,k) = tem + hvap * qo(i,k)
2932 heso(i,k) = tem + hvap * qeso(i,k)
2933 !c heo(i,k) = min(heo(i,k),heso(i,k))
2934 endif
2935 enddo
2936 enddo
2937 !c
2938 !c determine level with largest moist static energy
2939 !c this is the level where updraft starts
2940 !c
2941 do i=1,im
2942 hmax(i) = heo(i,1)
2943 kb(i) = 1
2944 enddo
2945 do k = 2, km
2946 do i=1,im
2947 if (k .le. kbm(i)) then
2948 if(heo(i,k).gt.hmax(i)) then
2949 kb(i) = k
2950 hmax(i) = heo(i,k)
2951 endif
2952 endif
2953 enddo
2954 enddo
2955 !c
2956 do k = 1, km1
2957 do i=1,im
2958 if (k .le. kmax(i)-1) then
2959 dz = .5 * (zo(i,k+1) - zo(i,k))
2960 dp = .5 * (pfld(i,k+1) - pfld(i,k))
2961 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
2962 pprime = pfld(i,k+1) + epsm1 * es
2963 qs = eps * es / pprime
2964 dqsdp = - qs / pprime
2965 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
2966 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
2967 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
2968 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
2969 dq = dqsdt * dt + dqsdp * dp
2970 to(i,k) = to(i,k+1) + dt
2971 qo(i,k) = qo(i,k+1) + dq
2972 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
2973 endif
2974 enddo
2975 enddo
2976 !
2977 do k = 1, km1
2978 do i=1,im
2979 if (k .le. kmax(i)-1) then
2980 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
2981 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
2982 val1 = 1.e-8
2983 qeso(i,k) = max(qeso(i,k), val1)
2984 val2 = 1.e-10
2985 qo(i,k) = max(qo(i,k), val2 )
2986 ! qo(i,k) = min(qo(i,k),qeso(i,k))
2987 val1 = 1.0
2988 frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), val1)
2989 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
2990 & cp * to(i,k) + hvap * qo(i,k)
2991 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
2992 & cp * to(i,k) + hvap * qeso(i,k)
2993 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
2994 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
2995 endif
2996 enddo
2997 enddo
2998 !c
2999 !c look for the level of free convection as cloud base
3000 !c
3001 do i=1,im
3002 flg(i) = .true.
3003 kbcon(i) = kmax(i)
3004 enddo
3005 do k = 1, km1
3006 do i=1,im
3007 if (flg(i).and.k.le.kbmax(i)) then
3008 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
3009 kbcon(i) = k
3010 flg(i) = .false.
3011 endif
3012 endif
3013 enddo
3014 enddo
3015 !c
3016 do i=1,im
3017 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
3018 enddo
3019 !!
3020 totflg = .true.
3021 do i=1,im
3022 totflg = totflg .and. (.not. cnvflg(i))
3023 enddo
3024 if(totflg) return
3025 !!
3026 !c
3027 !c determine critical convective inhibition
3028 !c as a function of vertical velocity at cloud base.
3029 !c
3030 do i=1,im
3031 if(cnvflg(i)) then
3032 pdot(i) = 10.* dot(i,kbcon(i))
3033 endif
3034 enddo
3035 do i=1,im
3036 if(cnvflg(i)) then
3037 if(slimsk(i).eq.1.) then
3038 w1 = w1l
3039 w2 = w2l
3040 w3 = w3l
3041 w4 = w4l
3042 else
3043 w1 = w1s
3044 w2 = w2s
3045 w3 = w3s
3046 w4 = w4s
3047 endif
3048 if(pdot(i).le.w4) then
3049 tem = (pdot(i) - w4) / (w3 - w4)
3050 elseif(pdot(i).ge.-w4) then
3051 tem = - (pdot(i) + w4) / (w4 - w3)
3052 else
3053 tem = 0.
3054 endif
3055 val1 = -1.
3056 tem = max(tem,val1)
3057 val2 = 1.
3058 tem = min(tem,val2)
3059 tem = 1. - tem
3060 tem1= .5*(cincrmax-cincrmin)
3061 cincr = cincrmax - tem * tem1
3062 pbcdif(i) = pfld(i,kb(i)) - pfld(i,kbcon(i))
3063 #if HWRF==1
3064 ! randomly perturb the convection trigger
3065 if( pert_sas_local ) then
3066 rr=2.0*ens_sasamp_local*ran1(ens_random_seed_local)-ens_sasamp_local
3067 print*, "zhang inde sas=", rr,ens_sasamp_local,ens_random_seed_local
3068 cincr=cincr+rr
3069 endif
3070 #endif
3071 if(pbcdif(i).gt.cincr) then
3072 cnvflg(i) = .false.
3073 endif
3074 endif
3075 enddo
3076 !!
3077 totflg = .true.
3078 do i=1,im
3079 totflg = totflg .and. (.not. cnvflg(i))
3080 enddo
3081 if(totflg) return
3082 !!
3083 !c
3084 !c assume that updraft entrainment rate above cloud base is
3085 !c same as that at cloud base
3086 !c
3087 do k = 2, km1
3088 do i=1,im
3089 if(cnvflg(i).and. &
3090 & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
3091 xlamue(i,k) = xlamue(i,kbcon(i))
3092 endif
3093 enddo
3094 enddo
3095 !c
3096 !c assume the detrainment rate for the updrafts to be same as
3097 !c the entrainment rate at cloud base
3098 !c
3099 do i = 1, im
3100 if(cnvflg(i)) then
3101 xlamud(i) = xlamue(i,kbcon(i))
3102 endif
3103 enddo
3104 !c
3105 !c functions rapidly decreasing with height, mimicking a cloud ensemble
3106 !c (Bechtold et al., 2008)
3107 !c
3108 do k = 2, km1
3109 do i=1,im
3110 if(cnvflg(i).and. &
3111 & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
3112 tem = qeso(i,k)/qeso(i,kbcon(i))
3113 fent1(i,k) = tem**2
3114 fent2(i,k) = tem**3
3115 endif
3116 enddo
3117 enddo
3118 !c
3119 !c final entrainment rate as the sum of turbulent part and organized entrainment
3120 !c depending on the environmental relative humidity
3121 !c (Bechtold et al., 2008)
3122 !c
3123 do k = 2, km1
3124 do i=1,im
3125 if(cnvflg(i).and. &
3126 & (k.ge.kbcon(i).and.k.lt.kmax(i))) then
3127 tem = cxlamu * frh(i,k) * fent2(i,k)
3128 xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
3129 endif
3130 enddo
3131 enddo
3132 !c
3133 !c determine updraft mass flux for the subcloud layers
3134 !c
3135 do k = km1, 1, -1
3136 do i = 1, im
3137 if (cnvflg(i)) then
3138 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
3139 dz = zi(i,k+1) - zi(i,k)
3140 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
3141 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
3142 endif
3143 endif
3144 enddo
3145 enddo
3146 !c
3147 !c compute mass flux above cloud base
3148 !c
3149 do k = 2, km1
3150 do i = 1, im
3151 if(cnvflg(i))then
3152 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
3153 dz = zi(i,k) - zi(i,k-1)
3154 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
3155 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
3156 endif
3157 endif
3158 enddo
3159 enddo
3160 !c
3161 !c compute updraft cloud properties
3162 !c
3163 do i = 1, im
3164 if(cnvflg(i)) then
3165 indx = kb(i)
3166 hcko(i,indx) = heo(i,indx)
3167 ucko(i,indx) = uo(i,indx)
3168 vcko(i,indx) = vo(i,indx)
3169 pwavo(i) = 0.
3170 endif
3171 enddo
3172 !c
3173 !c cloud property is modified by the entrainment process
3174 !c
3175 do k = 2, km1
3176 do i = 1, im
3177 if (cnvflg(i)) then
3178 if(k.gt.kb(i).and.k.lt.kmax(i)) then
3179 dz = zi(i,k) - zi(i,k-1)
3180 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3181 tem1 = 0.5 * xlamud(i) * dz
3182 factor = 1. + tem - tem1
3183 ptem = 0.5 * tem + pgcon
3184 ptem1= 0.5 * tem - pgcon
3185 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
3186 & (heo(i,k)+heo(i,k-1)))/factor
3187 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
3188 & +ptem1*uo(i,k-1))/factor
3189 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
3190 & +ptem1*vo(i,k-1))/factor
3191 dbyo(i,k) = hcko(i,k) - heso(i,k)
3192 endif
3193 endif
3194 enddo
3195 enddo
3196 !c
3197 !c taking account into convection inhibition due to existence of
3198 !c dry layers below cloud base
3199 !c
3200 do i=1,im
3201 flg(i) = cnvflg(i)
3202 kbcon1(i) = kmax(i)
3203 enddo
3204 do k = 2, km1
3205 do i=1,im
3206 if (flg(i).and.k.lt.kmax(i)) then
3207 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
3208 kbcon1(i) = k
3209 flg(i) = .false.
3210 endif
3211 endif
3212 enddo
3213 enddo
3214 do i=1,im
3215 if(cnvflg(i)) then
3216 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
3217 endif
3218 enddo
3219 do i=1,im
3220 if(cnvflg(i)) then
3221 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
3222 if(tem.gt.dthk) then
3223 cnvflg(i) = .false.
3224 endif
3225 endif
3226 enddo
3227 !!
3228 totflg = .true.
3229 do i = 1, im
3230 totflg = totflg .and. (.not. cnvflg(i))
3231 enddo
3232 if(totflg) return
3233 !!
3234 !c
3235 !c determine first guess cloud top as the level of zero buoyancy
3236 !c
3237 do i = 1, im
3238 flg(i) = cnvflg(i)
3239 ktcon(i) = 1
3240 enddo
3241 do k = 2, km1
3242 do i = 1, im
3243 if (flg(i).and.k .lt. kmax(i)) then
3244 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
3245 ktcon(i) = k
3246 flg(i) = .false.
3247 endif
3248 endif
3249 enddo
3250 enddo
3251 !c
3252 do i = 1, im
3253 if(cnvflg(i)) then
3254 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
3255 if(tem.lt.cthk) cnvflg(i) = .false.
3256 endif
3257 enddo
3258 !!
3259 totflg = .true.
3260 do i = 1, im
3261 totflg = totflg .and. (.not. cnvflg(i))
3262 enddo
3263 if(totflg) return
3264 !!
3265 !c
3266 !c search for downdraft originating level above theta-e minimum
3267 !c
3268 do i = 1, im
3269 if(cnvflg(i)) then
3270 hmin(i) = heo(i,kbcon1(i))
3271 lmin(i) = kbmax(i)
3272 jmin(i) = kbmax(i)
3273 endif
3274 enddo
3275 do k = 2, km1
3276 do i = 1, im
3277 if (cnvflg(i) .and. k .le. kbmax(i)) then
3278 if(k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
3279 lmin(i) = k + 1
3280 hmin(i) = heo(i,k)
3281 endif
3282 endif
3283 enddo
3284 enddo
3285 !c
3286 !c make sure that jmin(i) is within the cloud
3287 !c
3288 do i = 1, im
3289 if(cnvflg(i)) then
3290 jmin(i) = min(lmin(i),ktcon(i)-1)
3291 jmin(i) = max(jmin(i),kbcon1(i)+1)
3292 if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
3293 endif
3294 enddo
3295 !c
3296 !c specify upper limit of mass flux at cloud base
3297 !c
3298 do i = 1, im
3299 if(cnvflg(i)) then
3300 ! xmbmax(i) = .1
3301 !
3302 k = kbcon(i)
3303 dp = 1000. * del(i,k)
3304 xmbmax(i) = dp / (g * dt2)
3305 xmbmax(i) = min(sas_mass_flux,xmbmax(i))
3306 !
3307 ! tem = dp / (g * dt2)
3308 ! xmbmax(i) = min(tem, xmbmax(i))
3309 endif
3310 enddo
3311 !c
3312 !c compute cloud moisture property and precipitation
3313 !c
3314 do i = 1, im
3315 if (cnvflg(i)) then
3316 aa1(i) = 0.
3317 qcko(i,kb(i)) = qo(i,kb(i))
3318 ! rhbar(i) = 0.
3319 endif
3320 enddo
3321 do k = 2, km1
3322 do i = 1, im
3323 if (cnvflg(i)) then
3324 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
3325 dz = zi(i,k) - zi(i,k-1)
3326 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3327 qrch = qeso(i,k) &
3328 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3329 !cj
3330 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3331 tem1 = 0.5 * xlamud(i) * dz
3332 factor = 1. + tem - tem1
3333 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
3334 & (qo(i,k)+qo(i,k-1)))/factor
3335 !cj
3336 dq = eta(i,k) * (qcko(i,k) - qrch)
3337 !c
3338 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
3339 !c
3340 !c check if there is excess moisture to release latent heat
3341 !c
3342 if(k.ge.kbcon(i).and.dq.gt.0.) then
3343 etah = .5 * (eta(i,k) + eta(i,k-1))
3344 if(ncloud.gt.0..and.k.gt.jmin(i)) then
3345 dp = 1000. * del(i,k)
3346 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
3347 dellal(i,k) = etah * c1 * dz * qlk * g / dp
3348 else
3349 qlk = dq / (eta(i,k) + etah * c0 * dz)
3350 endif
3351 aa1(i) = aa1(i) - dz * g * qlk
3352 qcko(i,k) = qlk + qrch
3353 pwo(i,k) = etah * c0 * dz * qlk
3354 pwavo(i) = pwavo(i) + pwo(i,k)
3355 endif
3356 endif
3357 endif
3358 enddo
3359 enddo
3360 !c
3361 ! do i = 1, im
3362 ! if(cnvflg(i)) then
3363 ! indx = ktcon(i) - kb(i) - 1
3364 ! rhbar(i) = rhbar(i) / float(indx)
3365 ! endif
3366 ! enddo
3367 !c
3368 !c calculate cloud work function
3369 !c
3370 do k = 2, km1
3371 do i = 1, im
3372 if (cnvflg(i)) then
3373 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
3374 dz1 = zo(i,k+1) - zo(i,k)
3375 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3376 rfact = 1. + delta * cp * gamma &
3377 & * to(i,k) / hvap
3378 aa1(i) = aa1(i) + &
3379 & dz1 * (g / (cp * to(i,k))) &
3380 & * dbyo(i,k) / (1. + gamma) &
3381 & * rfact
3382 val = 0.
3383 aa1(i)=aa1(i)+ &
3384 & dz1 * g * delta * &
3385 & max(val,(qeso(i,k) - qo(i,k)))
3386 endif
3387 endif
3388 enddo
3389 enddo
3390 do i = 1, im
3391 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
3392 enddo
3393 !!
3394 totflg = .true.
3395 do i=1,im
3396 totflg = totflg .and. (.not. cnvflg(i))
3397 enddo
3398 if(totflg) return
3399 !!
3400 !c
3401 !c estimate the onvective overshooting as the level
3402 !c where the [aafac * cloud work function] becomes zero,
3403 !c which is the final cloud top
3404 !c
3405 do i = 1, im
3406 if (cnvflg(i)) then
3407 aa2(i) = aafac * aa1(i)
3408 endif
3409 enddo
3410 !c
3411 do i = 1, im
3412 flg(i) = cnvflg(i)
3413 ktcon1(i) = kmax(i) - 1
3414 enddo
3415 do k = 2, km1
3416 do i = 1, im
3417 if (flg(i)) then
3418 if(k.ge.ktcon(i).and.k.lt.kmax(i)) then
3419 dz1 = zo(i,k+1) - zo(i,k)
3420 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3421 rfact = 1. + delta * cp * gamma &
3422 & * to(i,k) / hvap
3423 aa2(i) = aa2(i) + &
3424 & dz1 * (g / (cp * to(i,k))) &
3425 & * dbyo(i,k) / (1. + gamma) &
3426 & * rfact
3427 if(aa2(i).lt.0.) then
3428 ktcon1(i) = k
3429 flg(i) = .false.
3430 endif
3431 endif
3432 endif
3433 enddo
3434 enddo
3435 !c
3436 !c compute cloud moisture property, detraining cloud water
3437 !c and precipitation in overshooting layers
3438 !c
3439 do k = 2, km1
3440 do i = 1, im
3441 if (cnvflg(i)) then
3442 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
3443 dz = zi(i,k) - zi(i,k-1)
3444 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3445 qrch = qeso(i,k) &
3446 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3447 !cj
3448 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3449 tem1 = 0.5 * xlamud(i) * dz
3450 factor = 1. + tem - tem1
3451 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
3452 & (qo(i,k)+qo(i,k-1)))/factor
3453 !cj
3454 dq = eta(i,k) * (qcko(i,k) - qrch)
3455 !c
3456 !c check if there is excess moisture to release latent heat
3457 !c
3458 if(dq.gt.0.) then
3459 etah = .5 * (eta(i,k) + eta(i,k-1))
3460 if(ncloud.gt.0.) then
3461 dp = 1000. * del(i,k)
3462 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
3463 dellal(i,k) = etah * c1 * dz * qlk * g / dp
3464 else
3465 qlk = dq / (eta(i,k) + etah * c0 * dz)
3466 endif
3467 qcko(i,k) = qlk + qrch
3468 pwo(i,k) = etah * c0 * dz * qlk
3469 pwavo(i) = pwavo(i) + pwo(i,k)
3470 endif
3471 endif
3472 endif
3473 enddo
3474 enddo
3475 !c
3476 !c exchange ktcon with ktcon1
3477 !c
3478 do i = 1, im
3479 if(cnvflg(i)) then
3480 kk = ktcon(i)
3481 ktcon(i) = ktcon1(i)
3482 ktcon1(i) = kk
3483 endif
3484 enddo
3485 !c
3486 !c this section is ready for cloud water
3487 !c
3488 if(ncloud.gt.0) then
3489 !c
3490 !c compute liquid and vapor separation at cloud top
3491 !c
3492 do i = 1, im
3493 if(cnvflg(i)) then
3494 k = ktcon(i) - 1
3495 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3496 qrch = qeso(i,k) &
3497 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3498 dq = qcko(i,k) - qrch
3499 !c
3500 !c check if there is excess moisture to release latent heat
3501 !c
3502 if(dq.gt.0.) then
3503 qlko_ktcon(i) = dq
3504 qcko(i,k) = qrch
3505 endif
3506 endif
3507 enddo
3508 endif
3509 !c
3510 !ccccc if(lat.eq.latd.and.lon.eq.lond.and.cnvflg(i)) then
3511 !ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
3512 !ccccc endif
3513 !c
3514 !c------- downdraft calculations
3515 !c
3516 !c--- compute precipitation efficiency in terms of windshear
3517 !c
3518 do i = 1, im
3519 if(cnvflg(i)) then
3520 vshear(i) = 0.
3521 endif
3522 enddo
3523 do k = 2, km
3524 do i = 1, im
3525 if (cnvflg(i)) then
3526 if(k.gt.kb(i).and.k.le.ktcon(i)) then
3527 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
3528 & + (vo(i,k)-vo(i,k-1)) ** 2)
3529 vshear(i) = vshear(i) + shear
3530 endif
3531 endif
3532 enddo
3533 enddo
3534 do i = 1, im
3535 if(cnvflg(i)) then
3536 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
3537 e1=1.591-.639*vshear(i) &
3538 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
3539 edt(i)=1.-e1
3540 val = .9
3541 edt(i) = min(edt(i),val)
3542 val = .0
3543 edt(i) = max(edt(i),val)
3544 edto(i)=edt(i)
3545 edtx(i)=edt(i)
3546 endif
3547 enddo
3548 !c
3549 !c determine detrainment rate between 1 and kbcon
3550 !c
3551 do i = 1, im
3552 if(cnvflg(i)) then
3553 sumx(i) = 0.
3554 endif
3555 enddo
3556 do k = 1, km1
3557 do i = 1, im
3558 if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
3559 dz = zi(i,k+1) - zi(i,k)
3560 sumx(i) = sumx(i) + dz
3561 endif
3562 enddo
3563 enddo
3564 do i = 1, im
3565 beta = betas
3566 if(slimsk(i).eq.1.) beta = betal
3567 if(cnvflg(i)) then
3568 dz = (sumx(i)+zi(i,1))/float(kbcon(i))
3569 tem = 1./float(kbcon(i))
3570 xlamd(i) = (1.-beta**tem)/dz
3571 endif
3572 enddo
3573 !c
3574 !c determine downdraft mass flux
3575 !c
3576 do k = km1, 1, -1
3577 do i = 1, im
3578 if (cnvflg(i) .and. k .le. kmax(i)-1) then
3579 if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
3580 dz = zi(i,k+1) - zi(i,k)
3581 ptem = xlamdd - xlamde
3582 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
3583 else if(k.lt.kbcon(i)) then
3584 dz = zi(i,k+1) - zi(i,k)
3585 ptem = xlamd(i) + xlamdd - xlamde
3586 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
3587 endif
3588 endif
3589 enddo
3590 enddo
3591 !c
3592 !c--- downdraft moisture properties
3593 !c
3594 do i = 1, im
3595 if(cnvflg(i)) then
3596 jmn = jmin(i)
3597 hcdo(i,jmn) = heo(i,jmn)
3598 qcdo(i,jmn) = qo(i,jmn)
3599 qrcdo(i,jmn)= qeso(i,jmn)
3600 ucdo(i,jmn) = uo(i,jmn)
3601 vcdo(i,jmn) = vo(i,jmn)
3602 pwevo(i) = 0.
3603 endif
3604 enddo
3605 !cj
3606 do k = km1, 1, -1
3607 do i = 1, im
3608 if (cnvflg(i) .and. k.lt.jmin(i)) then
3609 dz = zi(i,k+1) - zi(i,k)
3610 if(k.ge.kbcon(i)) then
3611 tem = xlamde * dz
3612 tem1 = 0.5 * xlamdd * dz
3613 else
3614 tem = xlamde * dz
3615 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
3616 endif
3617 factor = 1. + tem - tem1
3618 ptem = 0.5 * tem - pgcon
3619 ptem1= 0.5 * tem + pgcon
3620 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
3621 & (heo(i,k)+heo(i,k+1)))/factor
3622 ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1) &
3623 & +ptem1*uo(i,k))/factor
3624 vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1) &
3625 & +ptem1*vo(i,k))/factor
3626 dbyo(i,k) = hcdo(i,k) - heso(i,k)
3627 endif
3628 enddo
3629 enddo
3630 !c
3631 do k = km1, 1, -1
3632 do i = 1, im
3633 if (cnvflg(i).and.k.lt.jmin(i)) then
3634 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3635 qrcdo(i,k) = qeso(i,k)+ &
3636 & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
3637 ! detad = etad(i,k+1) - etad(i,k)
3638 !cj
3639 dz = zi(i,k+1) - zi(i,k)
3640 if(k.ge.kbcon(i)) then
3641 tem = xlamde * dz
3642 tem1 = 0.5 * xlamdd * dz
3643 else
3644 tem = xlamde * dz
3645 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
3646 endif
3647 factor = 1. + tem - tem1
3648 qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
3649 & (qo(i,k)+qo(i,k+1)))/factor
3650 !cj
3651 ! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
3652 ! & etad(i,k) * qrcdo(i,k)
3653 ! pwdo(i,k) = pwdo(i,k) - detad *
3654 ! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
3655 !cj
3656 pwdo(i,k) = etad(i,k+1) * (qcdo(i,k) - qrcdo(i,k))
3657 qcdo(i,k) = qrcdo(i,k)
3658 pwevo(i) = pwevo(i) + pwdo(i,k)
3659 endif
3660 enddo
3661 enddo
3662 !c
3663 !c--- final downdraft strength dependent on precip
3664 !c--- efficiency (edt), normalized condensate (pwav), and
3665 !c--- evaporate (pwev)
3666 !c
3667 do i = 1, im
3668 edtmax = edtmaxl
3669 if(slimsk(i).eq.0.) edtmax = edtmaxs
3670 if(cnvflg(i)) then
3671 if(pwevo(i).lt.0.) then
3672 edto(i) = -edto(i) * pwavo(i) / pwevo(i)
3673 edto(i) = min(edto(i),edtmax)
3674 else
3675 edto(i) = 0.
3676 endif
3677 endif
3678 enddo
3679 !c
3680 !c--- downdraft cloudwork functions
3681 !c
3682 do k = km1, 1, -1
3683 do i = 1, im
3684 if (cnvflg(i) .and. k .lt. jmin(i)) then
3685 gamma = el2orc * qeso(i,k) / to(i,k)**2
3686 dhh=hcdo(i,k)
3687 dt=to(i,k)
3688 dg=gamma
3689 dh=heso(i,k)
3690 dz=-1.*(zo(i,k+1)-zo(i,k))
3691 aa1(i)=aa1(i)+edto(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
3692 & *(1.+delta*cp*dg*dt/hvap)
3693 val=0.
3694 aa1(i)=aa1(i)+edto(i)* &
3695 & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
3696 endif
3697 enddo
3698 enddo
3699 do i = 1, im
3700 if(cnvflg(i).and.aa1(i).le.0.) then
3701 cnvflg(i) = .false.
3702 endif
3703 enddo
3704 !!
3705 totflg = .true.
3706 do i=1,im
3707 totflg = totflg .and. (.not. cnvflg(i))
3708 enddo
3709 if(totflg) return
3710 !!
3711 !c
3712 !c--- what would the change be, that a cloud with unit mass
3713 !c--- will do to the environment?
3714 !c
3715 do k = 1, km
3716 do i = 1, im
3717 if(cnvflg(i) .and. k .le. kmax(i)) then
3718 dellah(i,k) = 0.
3719 dellaq(i,k) = 0.
3720 dellau(i,k) = 0.
3721 dellav(i,k) = 0.
3722 endif
3723 enddo
3724 enddo
3725 do i = 1, im
3726 if(cnvflg(i)) then
3727 dp = 1000. * del(i,1)
3728 dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1) &
3729 & - heo(i,1)) * g / dp
3730 dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1) &
3731 & - qo(i,1)) * g / dp
3732 dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1) &
3733 & - uo(i,1)) * g / dp
3734 dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1) &
3735 & - vo(i,1)) * g / dp
3736 endif
3737 enddo
3738 !c
3739 !c--- changed due to subsidence and entrainment
3740 !c
3741 do k = 2, km1
3742 do i = 1, im
3743 if (cnvflg(i).and.k.lt.ktcon(i)) then
3744 aup = 1.
3745 if(k.le.kb(i)) aup = 0.
3746 adw = 1.
3747 if(k.gt.jmin(i)) adw = 0.
3748 dp = 1000. * del(i,k)
3749 dz = zi(i,k) - zi(i,k-1)
3750 !c
3751 dv1h = heo(i,k)
3752 dv2h = .5 * (heo(i,k) + heo(i,k-1))
3753 dv3h = heo(i,k-1)
3754 dv1q = qo(i,k)
3755 dv2q = .5 * (qo(i,k) + qo(i,k-1))
3756 dv3q = qo(i,k-1)
3757 dv1u = uo(i,k)
3758 dv2u = .5 * (uo(i,k) + uo(i,k-1))
3759 dv3u = uo(i,k-1)
3760 dv1v = vo(i,k)
3761 dv2v = .5 * (vo(i,k) + vo(i,k-1))
3762 dv3v = vo(i,k-1)
3763 !c
3764 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
3765 tem1 = xlamud(i)
3766 !c
3767 if(k.le.kbcon(i)) then
3768 ptem = xlamde
3769 ptem1 = xlamd(i)+xlamdd
3770 else
3771 ptem = xlamde
3772 ptem1 = xlamdd
3773 endif
3774 !cj
3775 dellah(i,k) = dellah(i,k) + &
3776 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h &
3777 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h &
3778 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz &
3779 & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
3780 & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1)) &
3781 & *dz) *g/dp
3782 !cj
3783 dellaq(i,k) = dellaq(i,k) + &
3784 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q &
3785 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q &
3786 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz &
3787 & + aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
3788 & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1)) &
3789 & *dz) *g/dp
3790 !23456789012345678901234567890123456789012345678901234567890123456789012
3791 !cj
3792 dellau(i,k) = dellau(i,k) + &
3793 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1u &
3794 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3u &
3795 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz &
3796 & + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
3797 & + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz &
3798 & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u) &
3799 & ) *g/dp
3800 !cj
3801 dellav(i,k) = dellav(i,k) + &
3802 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1v &
3803 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3v &
3804 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz &
3805 & + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
3806 & + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz &
3807 & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v) &
3808 & ) *g/dp
3809 !cj
3810 endif
3811 enddo
3812 enddo
3813 !c
3814 !c------- cloud top
3815 !c
3816 do i = 1, im
3817 if(cnvflg(i)) then
3818 indx = ktcon(i)
3819 dp = 1000. * del(i,indx)
3820 dv1h = heo(i,indx-1)
3821 dellah(i,indx) = eta(i,indx-1) * &
3822 & (hcko(i,indx-1) - dv1h) * g / dp
3823 dv1q = qo(i,indx-1)
3824 dellaq(i,indx) = eta(i,indx-1) * &
3825 & (qcko(i,indx-1) - dv1q) * g / dp
3826 dv1u = uo(i,indx-1)
3827 dellau(i,indx) = eta(i,indx-1) * &
3828 & (ucko(i,indx-1) - dv1u) * g / dp
3829 dv1v = vo(i,indx-1)
3830 dellav(i,indx) = eta(i,indx-1) * &
3831 & (vcko(i,indx-1) - dv1v) * g / dp
3832 !c
3833 !c cloud water
3834 !c
3835 dellal(i,indx) = eta(i,indx-1) * &
3836 & qlko_ktcon(i) * g / dp
3837 endif
3838 enddo
3839 !c
3840 !c------- final changed variable per unit mass flux
3841 !c
3842 do k = 1, km
3843 do i = 1, im
3844 if (cnvflg(i).and.k .le. kmax(i)) then
3845 if(k.gt.ktcon(i)) then
3846 qo(i,k) = q1(i,k)
3847 to(i,k) = t1(i,k)
3848 endif
3849 if(k.le.ktcon(i)) then
3850 qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
3851 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
3852 to(i,k) = dellat * mbdt + t1(i,k)
3853 val = 1.e-10
3854 qo(i,k) = max(qo(i,k), val )
3855 endif
3856 endif
3857 enddo
3858 enddo
3859 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3860 !c
3861 !c--- the above changed environment is now used to calulate the
3862 !c--- effect the arbitrary cloud (with unit mass flux)
3863 !c--- would have on the stability,
3864 !c--- which then is used to calculate the real mass flux,
3865 !c--- necessary to keep this change in balance with the large-scale
3866 !c--- destabilization.
3867 !c
3868 !c--- environmental conditions again, first heights
3869 !c
3870 do k = 1, km
3871 do i = 1, im
3872 if(cnvflg(i) .and. k .le. kmax(i)) then
3873 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
3874 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
3875 val = 1.e-8
3876 qeso(i,k) = max(qeso(i,k), val )
3877 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
3878 endif
3879 enddo
3880 enddo
3881 !c
3882 !c--- moist static energy
3883 !c
3884 do k = 1, km1
3885 do i = 1, im
3886 if(cnvflg(i) .and. k .le. kmax(i)-1) then
3887 dz = .5 * (zo(i,k+1) - zo(i,k))
3888 dp = .5 * (pfld(i,k+1) - pfld(i,k))
3889 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
3890 pprime = pfld(i,k+1) + epsm1 * es
3891 qs = eps * es / pprime
3892 dqsdp = - qs / pprime
3893 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
3894 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
3895 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
3896 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
3897 dq = dqsdt * dt + dqsdp * dp
3898 to(i,k) = to(i,k+1) + dt
3899 qo(i,k) = qo(i,k+1) + dq
3900 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
3901 endif
3902 enddo
3903 enddo
3904 do k = 1, km1
3905 do i = 1, im
3906 if(cnvflg(i) .and. k .le. kmax(i)-1) then
3907 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
3908 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
3909 val1 = 1.e-8
3910 qeso(i,k) = max(qeso(i,k), val1)
3911 val2 = 1.e-10
3912 qo(i,k) = max(qo(i,k), val2 )
3913 ! qo(i,k) = min(qo(i,k),qeso(i,k))
3914 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3915 & cp * to(i,k) + hvap * qo(i,k)
3916 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3917 & cp * to(i,k) + hvap * qeso(i,k)
3918 endif
3919 enddo
3920 enddo
3921 do i = 1, im
3922 if(cnvflg(i)) then
3923 k = kmax(i)
3924 heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
3925 heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
3926 !c heo(i,k) = min(heo(i,k),heso(i,k))
3927 endif
3928 enddo
3929 !c
3930 !c**************************** static control
3931 !c
3932 !c------- moisture and cloud work functions
3933 !c
3934 do i = 1, im
3935 if(cnvflg(i)) then
3936 xaa0(i) = 0.
3937 xpwav(i) = 0.
3938 endif
3939 enddo
3940 !c
3941 do i = 1, im
3942 if(cnvflg(i)) then
3943 indx = kb(i)
3944 hcko(i,indx) = heo(i,indx)
3945 qcko(i,indx) = qo(i,indx)
3946 endif
3947 enddo
3948 do k = 2, km1
3949 do i = 1, im
3950 if (cnvflg(i)) then
3951 if(k.gt.kb(i).and.k.le.ktcon(i)) then
3952 dz = zi(i,k) - zi(i,k-1)
3953 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3954 tem1 = 0.5 * xlamud(i) * dz
3955 factor = 1. + tem - tem1
3956 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
3957 & (heo(i,k)+heo(i,k-1)))/factor
3958 endif
3959 endif
3960 enddo
3961 enddo
3962 do k = 2, km1
3963 do i = 1, im
3964 if (cnvflg(i)) then
3965 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
3966 dz = zi(i,k) - zi(i,k-1)
3967 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3968 xdby = hcko(i,k) - heso(i,k)
3969 xqrch = qeso(i,k) &
3970 & + gamma * xdby / (hvap * (1. + gamma))
3971 !cj
3972 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3973 tem1 = 0.5 * xlamud(i) * dz
3974 factor = 1. + tem - tem1
3975 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
3976 & (qo(i,k)+qo(i,k-1)))/factor
3977 !cj
3978 dq = eta(i,k) * (qcko(i,k) - xqrch)
3979 !c
3980 if(k.ge.kbcon(i).and.dq.gt.0.) then
3981 etah = .5 * (eta(i,k) + eta(i,k-1))
3982 if(ncloud.gt.0..and.k.gt.jmin(i)) then
3983 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
3984 else
3985 qlk = dq / (eta(i,k) + etah * c0 * dz)
3986 endif
3987 if(k.lt.ktcon1(i)) then
3988 xaa0(i) = xaa0(i) - dz * g * qlk
3989 endif
3990 qcko(i,k) = qlk + xqrch
3991 xpw = etah * c0 * dz * qlk
3992 xpwav(i) = xpwav(i) + xpw
3993 endif
3994 endif
3995 if(k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
3996 dz1 = zo(i,k+1) - zo(i,k)
3997 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3998 rfact = 1. + delta * cp * gamma &
3999 & * to(i,k) / hvap
4000 xaa0(i) = xaa0(i) &
4001 & + dz1 * (g / (cp * to(i,k))) &
4002 & * xdby / (1. + gamma) &
4003 & * rfact
4004 val=0.
4005 xaa0(i)=xaa0(i)+ &
4006 & dz1 * g * delta * &
4007 & max(val,(qeso(i,k) - qo(i,k)))
4008 endif
4009 endif
4010 enddo
4011 enddo
4012 !c
4013 !c------- downdraft calculations
4014 !c
4015 !c--- downdraft moisture properties
4016 !c
4017 do i = 1, im
4018 if(cnvflg(i)) then
4019 jmn = jmin(i)
4020 hcdo(i,jmn) = heo(i,jmn)
4021 qcdo(i,jmn) = qo(i,jmn)
4022 qrcd(i,jmn) = qeso(i,jmn)
4023 xpwev(i) = 0.
4024 endif
4025 enddo
4026 !cj
4027 do k = km1, 1, -1
4028 do i = 1, im
4029 if (cnvflg(i) .and. k.lt.jmin(i)) then
4030 dz = zi(i,k+1) - zi(i,k)
4031 if(k.ge.kbcon(i)) then
4032 tem = xlamde * dz
4033 tem1 = 0.5 * xlamdd * dz
4034 else
4035 tem = xlamde * dz
4036 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
4037 endif
4038 factor = 1. + tem - tem1
4039 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
4040 & (heo(i,k)+heo(i,k+1)))/factor
4041 endif
4042 enddo
4043 enddo
4044 !cj
4045 do k = km1, 1, -1
4046 do i = 1, im
4047 if (cnvflg(i) .and. k .lt. jmin(i)) then
4048 dq = qeso(i,k)
4049 dt = to(i,k)
4050 gamma = el2orc * dq / dt**2
4051 dh = hcdo(i,k) - heso(i,k)
4052 qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
4053 ! detad = etad(i,k+1) - etad(i,k)
4054 !cj
4055 dz = zi(i,k+1) - zi(i,k)
4056 if(k.ge.kbcon(i)) then
4057 tem = xlamde * dz
4058 tem1 = 0.5 * xlamdd * dz
4059 else
4060 tem = xlamde * dz
4061 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
4062 endif
4063 factor = 1. + tem - tem1
4064 qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
4065 & (qo(i,k)+qo(i,k+1)))/factor
4066 !cj
4067 ! xpwd = etad(i,k+1) * qcdo(i,k+1) -
4068 ! & etad(i,k) * qrcd(i,k)
4069 ! xpwd = xpwd - detad *
4070 ! & .5 * (qrcd(i,k) + qrcd(i,k+1))
4071 !cj
4072 xpwd = etad(i,k+1) * (qcdo(i,k) - qrcd(i,k))
4073 qcdo(i,k)= qrcd(i,k)
4074 xpwev(i) = xpwev(i) + xpwd
4075 endif
4076 enddo
4077 enddo
4078 !c
4079 do i = 1, im
4080 edtmax = edtmaxl
4081 if(slimsk(i).eq.0.) edtmax = edtmaxs
4082 if(cnvflg(i)) then
4083 if(xpwev(i).ge.0.) then
4084 edtx(i) = 0.
4085 else
4086 edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
4087 edtx(i) = min(edtx(i),edtmax)
4088 endif
4089 endif
4090 enddo
4091 !c
4092 !c
4093 !c--- downdraft cloudwork functions
4094 !c
4095 !c
4096 do k = km1, 1, -1
4097 do i = 1, im
4098 if (cnvflg(i) .and. k.lt.jmin(i)) then
4099 gamma = el2orc * qeso(i,k) / to(i,k)**2
4100 dhh=hcdo(i,k)
4101 dt= to(i,k)
4102 dg= gamma
4103 dh= heso(i,k)
4104 dz=-1.*(zo(i,k+1)-zo(i,k))
4105 xaa0(i)=xaa0(i)+edtx(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
4106 & *(1.+delta*cp*dg*dt/hvap)
4107 val=0.
4108 xaa0(i)=xaa0(i)+edtx(i)* &
4109 & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
4110 endif
4111 enddo
4112 enddo
4113 !c
4114 !c calculate critical cloud work function
4115 !c
4116 do i = 1, im
4117 if(cnvflg(i)) then
4118 if(pfld(i,ktcon(i)).lt.pcrit(15))then
4119 acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i))) &
4120 & /(975.-pcrit(15))
4121 else if(pfld(i,ktcon(i)).gt.pcrit(1))then
4122 acrt(i)=acrit(1)
4123 else
4124 k = int((850. - pfld(i,ktcon(i)))/50.) + 2
4125 k = min(k,15)
4126 k = max(k,2)
4127 acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))* &
4128 & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
4129 endif
4130 endif
4131 enddo
4132 do i = 1, im
4133 if(cnvflg(i)) then
4134 if(slimsk(i).eq.1.) then
4135 w1 = w1l
4136 w2 = w2l
4137 w3 = w3l
4138 w4 = w4l
4139 else
4140 w1 = w1s
4141 w2 = w2s
4142 w3 = w3s
4143 w4 = w4s
4144 endif
4145 !c
4146 !c modify critical cloud workfunction by cloud base vertical velocity
4147 !c
4148 if(pdot(i).le.w4) then
4149 acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
4150 elseif(pdot(i).ge.-w4) then
4151 acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
4152 else
4153 acrtfct(i) = 0.
4154 endif
4155 val1 = -1.
4156 acrtfct(i) = max(acrtfct(i),val1)
4157 val2 = 1.
4158 acrtfct(i) = min(acrtfct(i),val2)
4159 acrtfct(i) = 1. - acrtfct(i)
4160 !c
4161 !c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
4162 !c
4163 !c if(rhbar(i).ge..8) then
4164 !c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
4165 !c endif
4166 !c
4167 !c modify adjustment time scale by cloud base vertical velocity
4168 !c
4169 val1=0.
4170 dtconv(i) = dt2 + max((1800. - dt2),val1) * &
4171 & (pdot(i) - w2) / (w1 - w2)
4172 !c dtconv(i) = max(dtconv(i), dt2)
4173 !c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
4174 dtconv(i) = max(dtconv(i),dtmin)
4175 dtconv(i) = min(dtconv(i),dtmax)
4176 !c
4177 endif
4178 enddo
4179 !c
4180 !c--- large scale forcing
4181 !c
4182 xmbmx1=-1.e20
4183 do i= 1, im
4184 if(cnvflg(i)) then
4185 fld(i)=(aa1(i)-acrt(i)* acrtfct(i))/dtconv(i)
4186 if(fld(i).le.0.) cnvflg(i) = .false.
4187 endif
4188 if(cnvflg(i)) then
4189 !c xaa0(i) = max(xaa0(i),0.)
4190 xk(i) = (xaa0(i) - aa1(i)) / mbdt
4191 if(xk(i).ge.0.) cnvflg(i) = .false.
4192 endif
4193 !c
4194 !c--- kernel, cloud base mass flux
4195 !c
4196 if(cnvflg(i)) then
4197 xmb(i) = -fld(i) / xk(i)
4198 xmb(i) = min(xmb(i),xmbmax(i))
4199 xmbmx1=max(xmbmx1,xmb(i))
4200 endif
4201 enddo
4202 ! if(xmbmx1.gt.0.4)print*,'qingfu test xmbmx1=',xmbmx1
4203 !!
4204 totflg = .true.
4205 do i=1,im
4206 totflg = totflg .and. (.not. cnvflg(i))
4207 enddo
4208 if(totflg) return
4209 !!
4210 !c
4211 !c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
4212 !c
4213 do k = 1, km
4214 do i = 1, im
4215 if (cnvflg(i) .and. k .le. kmax(i)) then
4216 to(i,k) = t1(i,k)
4217 qo(i,k) = q1(i,k)
4218 uo(i,k) = u1(i,k)
4219 vo(i,k) = v1(i,k)
4220 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
4221 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
4222 val = 1.e-8
4223 qeso(i,k) = max(qeso(i,k), val )
4224 endif
4225 enddo
4226 enddo
4227 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4228 !c
4229 !c--- feedback: simply the changes from the cloud with unit mass flux
4230 !c--- multiplied by the mass flux necessary to keep the
4231 !c--- equilibrium with the larger-scale.
4232 !c
4233 do i = 1, im
4234 delhbar(i) = 0.
4235 delqbar(i) = 0.
4236 deltbar(i) = 0.
4237 delubar(i) = 0.
4238 delvbar(i) = 0.
4239 qcond(i) = 0.
4240 enddo
4241 do k = 1, km
4242 do i = 1, im
4243 if (cnvflg(i) .and. k .le. kmax(i)) then
4244 if(k.le.ktcon(i)) then
4245 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
4246 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
4247 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
4248 tem = 1./rcs(i)
4249 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
4250 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
4251 dp = 1000. * del(i,k)
4252 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
4253 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
4254 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
4255 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
4256 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
4257 endif
4258 endif
4259 enddo
4260 enddo
4261 do k = 1, km
4262 do i = 1, im
4263 if (cnvflg(i) .and. k .le. kmax(i)) then
4264 if(k.le.ktcon(i)) then
4265 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
4266 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
4267 val = 1.e-8
4268 qeso(i,k) = max(qeso(i,k), val )
4269 endif
4270 endif
4271 enddo
4272 enddo
4273 !c
4274 do i = 1, im
4275 rntot(i) = 0.
4276 delqev(i) = 0.
4277 delq2(i) = 0.
4278 flg(i) = cnvflg(i)
4279 enddo
4280 do k = km, 1, -1
4281 do i = 1, im
4282 if (cnvflg(i) .and. k .le. kmax(i)) then
4283 if(k.lt.ktcon(i)) then
4284 aup = 1.
4285 if(k.le.kb(i)) aup = 0.
4286 adw = 1.
4287 if(k.ge.jmin(i)) adw = 0.
4288 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
4289 rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
4290 endif
4291 endif
4292 enddo
4293 enddo
4294 do k = km, 1, -1
4295 do i = 1, im
4296 if (k .le. kmax(i)) then
4297 deltv(i) = 0.
4298 delq(i) = 0.
4299 qevap(i) = 0.
4300 if(cnvflg(i).and.k.lt.ktcon(i)) then
4301 aup = 1.
4302 if(k.le.kb(i)) aup = 0.
4303 adw = 1.
4304 if(k.ge.jmin(i)) adw = 0.
4305 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
4306 rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
4307 endif
4308 if(flg(i).and.k.lt.ktcon(i)) then
4309 evef = edt(i) * evfact
4310 if(slimsk(i).eq.1.) evef=edt(i) * evfactl
4311 ! if(slimsk(i).eq.1.) evef=.07
4312 !c if(slimsk(i).ne.1.) evef = 0.
4313 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
4314 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
4315 dp = 1000. * del(i,k)
4316 if(rn(i).gt.0..and.qcond(i).lt.0.) then
4317 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
4318 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
4319 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
4320 endif
4321 if(rn(i).gt.0..and.qcond(i).lt.0..and. &
4322 & delq2(i).gt.rntot(i)) then
4323 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
4324 flg(i) = .false.
4325 endif
4326 if(rn(i).gt.0..and.qevap(i).gt.0.) then
4327 q1(i,k) = q1(i,k) + qevap(i)
4328 t1(i,k) = t1(i,k) - elocp * qevap(i)
4329 rn(i) = rn(i) - .001 * qevap(i) * dp / g
4330 deltv(i) = - elocp*qevap(i)/dt2
4331 delq(i) = + qevap(i)/dt2
4332 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
4333 endif
4334 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
4335 delqbar(i) = delqbar(i) + delq(i)*dp/g
4336 deltbar(i) = deltbar(i) + deltv(i)*dp/g
4337 endif
4338 endif
4339 enddo
4340 enddo
4341 !cj
4342 ! do i = 1, im
4343 ! if(me.eq.31.and.cnvflg(i)) then
4344 ! if(cnvflg(i)) then
4345 ! print *, ' deep delhbar, delqbar, deltbar = ',
4346 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
4347 ! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
4348 ! print *, ' precip =', hvap*rn(i)*1000./dt2
4349 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
4350 ! endif
4351 ! enddo
4352 !c
4353 !c precipitation rate converted to actual precip
4354 !c in unit of m instead of kg
4355 !c
4356 do i = 1, im
4357 if(cnvflg(i)) then
4358 !c
4359 !c in the event of upper level rain evaporation and lower level downdraft
4360 !c moistening, rn can become negative, in this case, we back out of the
4361 !c heating and the moistening
4362 !c
4363 if(rn(i).lt.0..and..not.flg(i)) rn(i) = 0.
4364 if(rn(i).le.0.) then
4365 rn(i) = 0.
4366 else
4367 ktop(i) = ktcon(i)
4368 kbot(i) = kbcon(i)
4369 kcnv(i) = 1
4370 cldwrk(i) = aa1(i)
4371 endif
4372 endif
4373 enddo
4374 !c
4375 !c cloud water
4376 !c
4377 if (ncloud.gt.0) then
4378 !
4379 val1=1.0
4380 val2=0.0
4381 do k = 1, km
4382 do i = 1, im
4383 if (cnvflg(i) .and. rn(i).gt.0.) then
4384 if (k.gt.kb(i).and.k.le.ktcon(i)) then
4385 tem = dellal(i,k) * xmb(i) * dt2
4386 tem1 = max(val2, min(val1, (tcr-t1(i,k))*tcrf))
4387 if (ql(i,k,2) .gt. -999.0) then
4388 ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
4389 ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
4390 else
4391 ql(i,k,1) = ql(i,k,1) + tem
4392 endif
4393 endif
4394 endif
4395 enddo
4396 enddo
4397 !
4398 endif
4399 !c
4400 do k = 1, km
4401 do i = 1, im
4402 if(cnvflg(i).and.rn(i).le.0.) then
4403 if (k .le. kmax(i)) then
4404 t1(i,k) = to(i,k)
4405 q1(i,k) = qo(i,k)
4406 u1(i,k) = uo(i,k)
4407 v1(i,k) = vo(i,k)
4408 endif
4409 endif
4410 enddo
4411 enddo
4412 !
4413 ! hchuang code change
4414 !
4415 ! do k = 1, km
4416 ! do i = 1, im
4417 ! if(cnvflg(i).and.rn(i).gt.0.) then
4418 ! if(k.ge.kb(i) .and. k.lt.ktop(i)) then
4419 ! ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
4420 ! endif
4421 ! endif
4422 ! enddo
4423 ! enddo
4424 ! do i = 1, im
4425 ! if(cnvflg(i).and.rn(i).gt.0.) then
4426 ! k = ktop(i)-1
4427 ! dt_mf(i,k) = ud_mf(i,k)
4428 ! endif
4429 ! enddo
4430 ! do k = 1, km
4431 ! do i = 1, im
4432 ! if(cnvflg(i).and.rn(i).gt.0.) then
4433 ! if(k.ge.1 .and. k.le.jmin(i)) then
4434 ! dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
4435 ! endif
4436 ! endif
4437 ! enddo
4438 ! enddo
4439 !!
4440 return
4441 end subroutine sascnvn
4443 subroutine shalcnv(im,ix,km,jcap,delt,del,prsl,ps,phil,ql, &
4444 & q1,t1,u1,v1,rcs,rn,kbot,ktop,kcnv,slimsk, &
4445 & dot,ncloud,hpbl,heat,evap,pgcon)
4446 !
4447 use MODULE_GFS_machine , only : kind_phys
4448 use MODULE_GFS_funcphys , only : fpvs
4449 use MODULE_GFS_physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
4450 &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
4451 &, rd => con_rd, cvap => con_cvap, cliq => con_cliq &
4452 &, eps => con_eps, epsm1 => con_epsm1
4453 implicit none
4454 !
4455 integer im, ix, km, jcap, ncloud, &
4456 & kbot(im), ktop(im), kcnv(im)
4457 real(kind=kind_phys) delt
4458 real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
4459 & ql(ix,km,2),q1(ix,km), t1(ix,km), &
4460 & u1(ix,km), v1(ix,km), rcs(im), &
4461 & rn(im), slimsk(im), &
4462 & dot(ix,km), phil(ix,km), hpbl(im), &
4463 & heat(im), evap(im)
4464 ! &, ud_mf(im,km),dt_mf(im,km)
4466 real ud_mf(im,km),dt_mf(im,km)
4467 !
4468 integer i,j,indx, jmn, k, kk, latd, lond, km1
4469 integer kpbl(im)
4470 !
4471 real(kind=kind_phys) c0, cpoel, dellat, delta, &
4472 & desdt, deta, detad, dg, &
4473 & dh, dhh, dlnsig, dp, &
4474 & dq, dqsdp, dqsdt, dt, &
4475 & dt2, dtmax, dtmin, dv1h, &
4476 & dv1q, dv2h, dv2q, dv1u, &
4477 & dv1v, dv2u, dv2v, dv3q, &
4478 & dv3h, dv3u, dv3v, clam, &
4479 & dz, dz1, e1, &
4480 & el2orc, elocp, aafac, cthk, &
4481 & es, etah, h1, dthk, &
4482 & evef, evfact, evfactl, fact1, &
4483 & fact2, factor, fjcap, &
4484 & g, gamma, pprime, betaw, &
4485 & qlk, qrch, qs, c1, &
4486 & rain, rfact, shear, tem1, &
4487 & tem2, terr, val, val1, &
4488 & val2, w1, w1l, w1s, &
4489 & w2, w2l, w2s, w3, &
4490 & w3l, w3s, w4, w4l, &
4491 & w4s, tem, ptem, ptem1, &
4492 & pgcon
4493 !
4494 integer kb(im), kbcon(im), kbcon1(im), &
4495 & ktcon(im), ktcon1(im), &
4496 & kbm(im), kmax(im)
4497 !
4498 real(kind=kind_phys) aa1(im), &
4499 & delhbar(im), delq(im), delq2(im), &
4500 & delqbar(im), delqev(im), deltbar(im), &
4501 & deltv(im), edt(im), &
4502 & wstar(im), sflx(im), &
4503 & pdot(im), po(im,km), &
4504 & qcond(im), qevap(im), hmax(im), &
4505 & rntot(im), vshear(im), &
4506 & xlamud(im), xmb(im), xmbmax(im), &
4507 & delubar(im), delvbar(im)
4508 !c
4509 real(kind=kind_phys) cincr, cincrmax, cincrmin
4510 !cc
4511 !c physical parameters
4512 parameter(g=grav)
4513 parameter(cpoel=cp/hvap,elocp=hvap/cp, &
4514 & el2orc=hvap*hvap/(rv*cp))
4515 parameter(terr=0.,c0=.002,c1=5.e-4,delta=fv)
4516 parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
4517 parameter(cthk=50.,cincrmax=180.,cincrmin=120.,dthk=25.)
4518 parameter(h1=0.33333333)
4519 !c local variables and arrays
4520 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), &
4521 & uo(im,km), vo(im,km), qeso(im,km)
4522 !c cloud water
4523 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), &
4524 & dbyo(im,km), zo(im,km), xlamue(im,km), &
4525 & heo(im,km), heso(im,km), &
4526 & dellah(im,km), dellaq(im,km), &
4527 & dellau(im,km), dellav(im,km), hcko(im,km), &
4528 & ucko(im,km), vcko(im,km), qcko(im,km), &
4529 & eta(im,km), zi(im,km), pwo(im,km), &
4530 & tx1(im)
4531 !
4532 logical totflg, cnvflg(im), flg(im)
4533 !
4534 real(kind=kind_phys) tf, tcr, tcrf
4535 parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
4536 !
4537 !c-----------------------------------------------------------------------
4538 !
4539 km1 = km - 1
4540 !c
4541 !c compute surface buoyancy flux
4542 !c
4543 do i=1,im
4544 sflx(i) = heat(i)+fv*t1(i,1)*evap(i)
4545 enddo
4546 !c
4547 !c initialize arrays
4548 !c
4549 do i=1,im
4550 cnvflg(i) = .true.
4551 if(kcnv(i).eq.1) cnvflg(i) = .false.
4552 if(sflx(i).le.0.) cnvflg(i) = .false.
4553 if(cnvflg(i)) then
4554 kbot(i)=km+1
4555 ktop(i)=0
4556 endif
4557 rn(i)=0.
4558 kbcon(i)=km
4559 ktcon(i)=1
4560 kb(i)=km
4561 pdot(i) = 0.
4562 qlko_ktcon(i) = 0.
4563 edt(i) = 0.
4564 aa1(i) = 0.
4565 vshear(i) = 0.
4566 enddo
4567 ! hchuang code change
4568 do k = 1, km
4569 do i = 1, im
4570 ud_mf(i,k) = 0.
4571 dt_mf(i,k) = 0.
4572 enddo
4573 enddo
4574 !!
4575 totflg = .true.
4576 do i=1,im
4577 totflg = totflg .and. (.not. cnvflg(i))
4578 enddo
4579 if(totflg) return
4580 !!
4581 !c
4582 dt2 = delt
4583 val = 1200.
4584 dtmin = max(dt2, val )
4585 val = 3600.
4586 dtmax = max(dt2, val )
4587 !c model tunable parameters are all here
4588 clam = .3
4589 aafac = .1
4590 betaw = .03
4591 !c evef = 0.07
4592 evfact = 0.3
4593 evfactl = 0.3
4594 !
4595 fjcap = (float(jcap) / 126.) ** 2
4596 val = 1.
4597 fjcap = max(fjcap,val)
4598 w1l = -8.e-3
4599 w2l = -4.e-2
4600 w3l = -5.e-3
4601 w4l = -5.e-4
4602 w1s = -2.e-4
4603 w2s = -2.e-3
4604 w3s = -1.e-3
4605 w4s = -2.e-5
4606 !c
4607 !c define top layer for search of the downdraft originating layer
4608 !c and the maximum thetae for updraft
4609 !c
4610 do i=1,im
4611 kbm(i) = km
4612 kmax(i) = km
4613 tx1(i) = 1.0 / ps(i)
4614 enddo
4615 !
4616 do k = 1, km
4617 do i=1,im
4618 if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
4619 if (prsl(i,k)*tx1(i) .gt. 0.60) kmax(i) = k + 1
4620 enddo
4621 enddo
4622 do i=1,im
4623 kbm(i) = min(kbm(i),kmax(i))
4624 enddo
4625 !c
4626 !!c hydrostatic height assume zero terr and compute
4627 !c updraft entrainment rate as an inverse function of height
4628 !c
4629 do k = 1, km
4630 do i=1,im
4631 zo(i,k) = phil(i,k) / g
4632 enddo
4633 enddo
4634 do k = 1, km1
4635 do i=1,im
4636 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
4637 xlamue(i,k) = clam / zi(i,k)
4638 enddo
4639 enddo
4640 do i=1,im
4641 xlamue(i,km) = xlamue(i,km1)
4642 enddo
4643 !c
4644 !c pbl height
4645 !c
4646 do i=1,im
4647 flg(i) = cnvflg(i)
4648 kpbl(i)= 1
4649 enddo
4650 do k = 2, km1
4651 do i=1,im
4652 if (flg(i).and.zo(i,k).le.hpbl(i)) then
4653 kpbl(i) = k
4654 else
4655 flg(i) = .false.
4656 endif
4657 enddo
4658 enddo
4659 do i=1,im
4660 kpbl(i)= min(kpbl(i),kbm(i))
4661 enddo
4662 !c
4663 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4664 !c convert surface pressure to mb from cb
4665 !c
4666 do k = 1, km
4667 do i = 1, im
4668 if (cnvflg(i) .and. k .le. kmax(i)) then
4669 pfld(i,k) = prsl(i,k) * 10.0
4670 eta(i,k) = 1.
4671 hcko(i,k) = 0.
4672 qcko(i,k) = 0.
4673 ucko(i,k) = 0.
4674 vcko(i,k) = 0.
4675 dbyo(i,k) = 0.
4676 pwo(i,k) = 0.
4677 dellal(i,k) = 0.
4678 to(i,k) = t1(i,k)
4679 qo(i,k) = q1(i,k)
4680 uo(i,k) = u1(i,k) * rcs(i)
4681 vo(i,k) = v1(i,k) * rcs(i)
4682 endif
4683 enddo
4684 enddo
4685 !c
4686 !c column variables
4687 !c p is pressure of the layer (mb)
4688 !c t is temperature at t-dt (k)..tn
4689 !c q is mixing ratio at t-dt (kg/kg)..qn
4690 !c to is temperature at t+dt (k)... this is after advection and turbulan
4691 !c qo is mixing ratio at t+dt (kg/kg)..q1
4692 !c
4693 do k = 1, km
4694 do i=1,im
4695 if (cnvflg(i) .and. k .le. kmax(i)) then
4696 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
4697 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
4698 val1 = 1.e-8
4699 qeso(i,k) = max(qeso(i,k), val1)
4700 val2 = 1.e-10
4701 qo(i,k) = max(qo(i,k), val2 )
4702 ! qo(i,k) = min(qo(i,k),qeso(i,k))
4703 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
4704 endif
4705 enddo
4706 enddo
4707 !c
4708 !c compute moist static energy
4709 !c
4710 do k = 1, km
4711 do i=1,im
4712 if (cnvflg(i) .and. k .le. kmax(i)) then
4713 ! tem = g * zo(i,k) + cp * to(i,k)
4714 tem = phil(i,k) + cp * to(i,k)
4715 heo(i,k) = tem + hvap * qo(i,k)
4716 heso(i,k) = tem + hvap * qeso(i,k)
4717 !c heo(i,k) = min(heo(i,k),heso(i,k))
4718 endif
4719 enddo
4720 enddo
4721 !c
4722 !c determine level with largest moist static energy within pbl
4723 !c this is the level where updraft starts
4724 !c
4725 do i=1,im
4726 if (cnvflg(i)) then
4727 hmax(i) = heo(i,1)
4728 kb(i) = 1
4729 endif
4730 enddo
4731 do k = 2, km
4732 do i=1,im
4733 if (cnvflg(i).and.k.le.kpbl(i)) then
4734 if(heo(i,k).gt.hmax(i)) then
4735 kb(i) = k
4736 hmax(i) = heo(i,k)
4737 endif
4738 endif
4739 enddo
4740 enddo
4741 !c
4742 do k = 1, km1
4743 do i=1,im
4744 if (cnvflg(i) .and. k .le. kmax(i)-1) then
4745 dz = .5 * (zo(i,k+1) - zo(i,k))
4746 dp = .5 * (pfld(i,k+1) - pfld(i,k))
4747 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
4748 pprime = pfld(i,k+1) + epsm1 * es
4749 qs = eps * es / pprime
4750 dqsdp = - qs / pprime
4751 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
4752 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
4753 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
4754 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
4755 dq = dqsdt * dt + dqsdp * dp
4756 to(i,k) = to(i,k+1) + dt
4757 qo(i,k) = qo(i,k+1) + dq
4758 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
4759 endif
4760 enddo
4761 enddo
4762 !
4763 do k = 1, km1
4764 do i=1,im
4765 if (cnvflg(i) .and. k .le. kmax(i)-1) then
4766 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
4767 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
4768 val1 = 1.e-8
4769 qeso(i,k) = max(qeso(i,k), val1)
4770 val2 = 1.e-10
4771 qo(i,k) = max(qo(i,k), val2 )
4772 ! qo(i,k) = min(qo(i,k),qeso(i,k))
4773 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
4774 & cp * to(i,k) + hvap * qo(i,k)
4775 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
4776 & cp * to(i,k) + hvap * qeso(i,k)
4777 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
4778 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
4779 endif
4780 enddo
4781 enddo
4782 !c
4783 !c look for the level of free convection as cloud base
4784 !c
4785 do i=1,im
4786 flg(i) = cnvflg(i)
4787 if(flg(i)) kbcon(i) = kmax(i)
4788 enddo
4789 do k = 2, km1
4790 do i=1,im
4791 if (flg(i).and.k.lt.kbm(i)) then
4792 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
4793 kbcon(i) = k
4794 flg(i) = .false.
4795 endif
4796 endif
4797 enddo
4798 enddo
4799 !c
4800 do i=1,im
4801 if(cnvflg(i)) then
4802 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
4803 endif
4804 enddo
4805 !!
4806 totflg = .true.
4807 do i=1,im
4808 totflg = totflg .and. (.not. cnvflg(i))
4809 enddo
4810 if(totflg) return
4811 !!
4812 !c
4813 !c determine critical convective inhibition
4814 !c as a function of vertical velocity at cloud base.
4815 !c
4816 do i=1,im
4817 if(cnvflg(i)) then
4818 pdot(i) = 10.* dot(i,kbcon(i))
4819 endif
4820 enddo
4821 do i=1,im
4822 if(cnvflg(i)) then
4823 if(slimsk(i).eq.1.) then
4824 w1 = w1l
4825 w2 = w2l
4826 w3 = w3l
4827 w4 = w4l
4828 else
4829 w1 = w1s
4830 w2 = w2s
4831 w3 = w3s
4832 w4 = w4s
4833 endif
4834 if(pdot(i).le.w4) then
4835 ptem = (pdot(i) - w4) / (w3 - w4)
4836 elseif(pdot(i).ge.-w4) then
4837 ptem = - (pdot(i) + w4) / (w4 - w3)
4838 else
4839 ptem = 0.
4840 endif
4841 val1 = -1.
4842 ptem = max(ptem,val1)
4843 val2 = 1.
4844 ptem = min(ptem,val2)
4845 ptem = 1. - ptem
4846 ptem1= .5*(cincrmax-cincrmin)
4847 cincr = cincrmax - ptem * ptem1
4848 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
4849 if(tem1.gt.cincr) then
4850 cnvflg(i) = .false.
4851 endif
4852 endif
4853 enddo
4854 !!
4855 totflg = .true.
4856 do i=1,im
4857 totflg = totflg .and. (.not. cnvflg(i))
4858 enddo
4859 if(totflg) return
4860 !!
4861 !c
4862 !c assume the detrainment rate for the updrafts to be same as
4863 !c the entrainment rate at cloud base
4864 !c
4865 do i = 1, im
4866 if(cnvflg(i)) then
4867 xlamud(i) = xlamue(i,kbcon(i))
4868 endif
4869 enddo
4870 !c
4871 !c determine updraft mass flux for the subcloud layers
4872 !c
4873 do k = km1, 1, -1
4874 do i = 1, im
4875 if (cnvflg(i)) then
4876 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
4877 dz = zi(i,k+1) - zi(i,k)
4878 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
4879 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
4880 endif
4881 endif
4882 enddo
4883 enddo
4884 !c
4885 !c compute mass flux above cloud base
4886 !c
4887 do k = 2, km1
4888 do i = 1, im
4889 if(cnvflg(i))then
4890 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
4891 dz = zi(i,k) - zi(i,k-1)
4892 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
4893 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
4894 endif
4895 endif
4896 enddo
4897 enddo
4898 !c
4899 !c compute updraft cloud property
4900 !c
4901 do i = 1, im
4902 if(cnvflg(i)) then
4903 indx = kb(i)
4904 hcko(i,indx) = heo(i,indx)
4905 ucko(i,indx) = uo(i,indx)
4906 vcko(i,indx) = vo(i,indx)
4907 endif
4908 enddo
4909 !c
4910 do k = 2, km1
4911 do i = 1, im
4912 if (cnvflg(i)) then
4913 if(k.gt.kb(i).and.k.lt.kmax(i)) then
4914 dz = zi(i,k) - zi(i,k-1)
4915 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
4916 tem1 = 0.5 * xlamud(i) * dz
4917 factor = 1. + tem - tem1
4918 ptem = 0.5 * tem + pgcon
4919 ptem1= 0.5 * tem - pgcon
4920 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
4921 & (heo(i,k)+heo(i,k-1)))/factor
4922 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
4923 & +ptem1*uo(i,k-1))/factor
4924 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
4925 & +ptem1*vo(i,k-1))/factor
4926 dbyo(i,k) = hcko(i,k) - heso(i,k)
4927 endif
4928 endif
4929 enddo
4930 enddo
4931 !c
4932 !c taking account into convection inhibition due to existence of
4933 !c dry layers below cloud base
4934 !c
4935 do i=1,im
4936 flg(i) = cnvflg(i)
4937 kbcon1(i) = kmax(i)
4938 enddo
4939 do k = 2, km1
4940 do i=1,im
4941 if (flg(i).and.k.lt.kbm(i)) then
4942 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
4943 kbcon1(i) = k
4944 flg(i) = .false.
4945 endif
4946 endif
4947 enddo
4948 enddo
4949 do i=1,im
4950 if(cnvflg(i)) then
4951 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
4952 endif
4953 enddo
4954 do i=1,im
4955 if(cnvflg(i)) then
4956 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
4957 if(tem.gt.dthk) then
4958 cnvflg(i) = .false.
4959 endif
4960 endif
4961 enddo
4962 !!
4963 totflg = .true.
4964 do i = 1, im
4965 totflg = totflg .and. (.not. cnvflg(i))
4966 enddo
4967 if(totflg) return
4968 !!
4969 !c
4970 !c determine first guess cloud top as the level of zero buoyancy
4971 !c limited to the level of sigma=0.7
4972 !c
4973 do i = 1, im
4974 flg(i) = cnvflg(i)
4975 if(flg(i)) ktcon(i) = kbm(i)
4976 enddo
4977 do k = 2, km1
4978 do i=1,im
4979 if (flg(i).and.k .lt. kbm(i)) then
4980 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
4981 ktcon(i) = k
4982 flg(i) = .false.
4983 endif
4984 endif
4985 enddo
4986 enddo
4987 !c
4988 !c turn off shallow convection if cloud top is less than pbl top
4989 !c
4990 do i=1,im
4991 if(cnvflg(i)) then
4992 kk = kpbl(i)+1
4993 if(ktcon(i).le.kk) cnvflg(i) = .false.
4994 endif
4995 enddo
4996 ! c
4997 ! c turn off shallow convection if cloud depth is less than
4998 ! c a threshold value (cthk)
4999 ! c
5000 do i = 1, im
5001 if(cnvflg(i)) then
5002 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
5003 if(tem.lt.cthk) cnvflg(i) = .false.
5004 endif
5005 enddo
5006 !!
5007 totflg = .true.
5008 do i = 1, im
5009 totflg = totflg .and. (.not. cnvflg(i))
5010 enddo
5011 if(totflg) return
5012 !!
5013 !c
5014 !c specify upper limit of mass flux at cloud base
5015 !c
5016 do i = 1, im
5017 if(cnvflg(i)) then
5018 ! xmbmax(i) = .1
5019 !
5020 k = kbcon(i)
5021 dp = 1000. * del(i,k)
5022 xmbmax(i) = dp / (g * dt2)
5023 !
5024 ! tem = dp / (g * dt2)
5025 ! xmbmax(i) = min(tem, xmbmax(i))
5026 endif
5027 enddo
5028 !c
5029 !c compute cloud moisture property and precipitation
5030 !c
5031 do i = 1, im
5032 if (cnvflg(i)) then
5033 aa1(i) = 0.
5034 qcko(i,kb(i)) = qo(i,kb(i))
5035 endif
5036 enddo
5037 do k = 2, km1
5038 do i = 1, im
5039 if (cnvflg(i)) then
5040 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
5041 dz = zi(i,k) - zi(i,k-1)
5042 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5043 qrch = qeso(i,k) &
5044 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
5045 !cj
5046 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
5047 tem1 = 0.5 * xlamud(i) * dz
5048 factor = 1. + tem - tem1
5049 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
5050 & (qo(i,k)+qo(i,k-1)))/factor
5051 !cj
5052 dq = eta(i,k) * (qcko(i,k) - qrch)
5053 !c
5054 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
5055 !c
5056 !c below lfc check if there is excess moisture to release latent heat
5057 !c
5058 if(k.ge.kbcon(i).and.dq.gt.0.) then
5059 etah = .5 * (eta(i,k) + eta(i,k-1))
5060 if(ncloud.gt.0.) then
5061 dp = 1000. * del(i,k)
5062 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
5063 dellal(i,k) = etah * c1 * dz * qlk * g / dp
5064 else
5065 qlk = dq / (eta(i,k) + etah * c0 * dz)
5066 endif
5067 aa1(i) = aa1(i) - dz * g * qlk
5068 qcko(i,k)= qlk + qrch
5069 pwo(i,k) = etah * c0 * dz * qlk
5070 endif
5071 endif
5072 endif
5073 enddo
5074 enddo
5075 !c
5076 !c calculate cloud work function
5077 !c
5078 do k = 2, km1
5079 do i = 1, im
5080 if (cnvflg(i)) then
5081 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
5082 dz1 = zo(i,k+1) - zo(i,k)
5083 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5084 rfact = 1. + delta * cp * gamma &
5085 & * to(i,k) / hvap
5086 aa1(i) = aa1(i) + &
5087 & dz1 * (g / (cp * to(i,k))) &
5088 & * dbyo(i,k) / (1. + gamma) &
5089 & * rfact
5090 val = 0.
5091 aa1(i)=aa1(i)+ &
5092 & dz1 * g * delta * &
5093 & max(val,(qeso(i,k) - qo(i,k)))
5094 endif
5095 endif
5096 enddo
5097 enddo
5098 do i = 1, im
5099 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
5100 enddo
5101 !!
5102 totflg = .true.
5103 do i=1,im
5104 totflg = totflg .and. (.not. cnvflg(i))
5105 enddo
5106 if(totflg) return
5107 !!
5108 !c
5109 !c estimate the onvective overshooting as the level
5110 !c where the [aafac * cloud work function] becomes zero,
5111 !c which is the final cloud top
5112 !c limited to the level of sigma=0.7
5113 !c
5114 do i = 1, im
5115 if (cnvflg(i)) then
5116 aa1(i) = aafac * aa1(i)
5117 endif
5118 enddo
5119 !c
5120 do i = 1, im
5121 flg(i) = cnvflg(i)
5122 ktcon1(i) = kbm(i)
5123 enddo
5124 do k = 2, km1
5125 do i = 1, im
5126 if (flg(i)) then
5127 if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
5128 dz1 = zo(i,k+1) - zo(i,k)
5129 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5130 rfact = 1. + delta * cp * gamma &
5131 & * to(i,k) / hvap
5132 aa1(i) = aa1(i) + &
5133 & dz1 * (g / (cp * to(i,k))) &
5134 & * dbyo(i,k) / (1. + gamma) &
5135 & * rfact
5136 if(aa1(i).lt.0.) then
5137 ktcon1(i) = k
5138 flg(i) = .false.
5139 endif
5140 endif
5141 endif
5142 enddo
5143 enddo
5144 !c
5145 !c compute cloud moisture property, detraining cloud water
5146 !c and precipitation in overshooting layers
5147 !c
5148 do k = 2, km1
5149 do i = 1, im
5150 if (cnvflg(i)) then
5151 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
5152 dz = zi(i,k) - zi(i,k-1)
5153 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5154 qrch = qeso(i,k) &
5155 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
5156 !cj
5157 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
5158 tem1 = 0.5 * xlamud(i) * dz
5159 factor = 1. + tem - tem1
5160 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
5161 & (qo(i,k)+qo(i,k-1)))/factor
5162 !cj
5163 dq = eta(i,k) * (qcko(i,k) - qrch)
5164 !c
5165 !c check if there is excess moisture to release latent heat
5166 !c
5167 if(dq.gt.0.) then
5168 etah = .5 * (eta(i,k) + eta(i,k-1))
5169 if(ncloud.gt.0.) then
5170 dp = 1000. * del(i,k)
5171 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
5172 dellal(i,k) = etah * c1 * dz * qlk * g / dp
5173 else
5174 qlk = dq / (eta(i,k) + etah * c0 * dz)
5175 endif
5176 qcko(i,k) = qlk + qrch
5177 pwo(i,k) = etah * c0 * dz * qlk
5178 endif
5179 endif
5180 endif
5181 enddo
5182 enddo
5183 !c
5184 !c exchange ktcon with ktcon1
5185 !c
5186 do i = 1, im
5187 if(cnvflg(i)) then
5188 kk = ktcon(i)
5189 ktcon(i) = ktcon1(i)
5190 ktcon1(i) = kk
5191 endif
5192 enddo
5193 !c
5194 !c this section is ready for cloud water
5195 !c
5196 if(ncloud.gt.0) then
5197 !c
5198 !c compute liquid and vapor separation at cloud top
5199 !c
5200 do i = 1, im
5201 if(cnvflg(i)) then
5202 k = ktcon(i) - 1
5203 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5204 qrch = qeso(i,k) &
5205 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
5206 dq = qcko(i,k) - qrch
5207 !c
5208 !c check if there is excess moisture to release latent heat
5209 !c
5210 if(dq.gt.0.) then
5211 qlko_ktcon(i) = dq
5212 qcko(i,k) = qrch
5213 endif
5214 endif
5215 enddo
5216 endif
5217 !!c
5218 !c--- compute precipitation efficiency in terms of windshear
5219 !c
5220 do i = 1, im
5221 if(cnvflg(i)) then
5222 vshear(i) = 0.
5223 endif
5224 enddo
5225 do k = 2, km
5226 do i = 1, im
5227 if (cnvflg(i)) then
5228 if(k.gt.kb(i).and.k.le.ktcon(i)) then
5229 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
5230 & + (vo(i,k)-vo(i,k-1)) ** 2)
5231 vshear(i) = vshear(i) + shear
5232 endif
5233 endif
5234 enddo
5235 enddo
5236 do i = 1, im
5237 if(cnvflg(i)) then
5238 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
5239 e1=1.591-.639*vshear(i) &
5240 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
5241 edt(i)=1.-e1
5242 val = .9
5243 edt(i) = min(edt(i),val)
5244 val = .0
5245 edt(i) = max(edt(i),val)
5246 endif
5247 enddo
5248 !c
5249 !c--- what would the change be, that a cloud with unit mass
5250 !c--- will do to the environment?
5251 !c
5252 do k = 1, km
5253 do i = 1, im
5254 if(cnvflg(i) .and. k .le. kmax(i)) then
5255 dellah(i,k) = 0.
5256 dellaq(i,k) = 0.
5257 dellau(i,k) = 0.
5258 dellav(i,k) = 0.
5259 endif
5260 enddo
5261 enddo
5262 !c
5263 !c--- changed due to subsidence and entrainment
5264 !c
5265 do k = 2, km1
5266 do i = 1, im
5267 if (cnvflg(i)) then
5268 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
5269 dp = 1000. * del(i,k)
5270 dz = zi(i,k) - zi(i,k-1)
5271 !c
5272 dv1h = heo(i,k)
5273 dv2h = .5 * (heo(i,k) + heo(i,k-1))
5274 dv3h = heo(i,k-1)
5275 dv1q = qo(i,k)
5276 dv2q = .5 * (qo(i,k) + qo(i,k-1))
5277 dv3q = qo(i,k-1)
5278 dv1u = uo(i,k)
5279 dv2u = .5 * (uo(i,k) + uo(i,k-1))
5280 dv3u = uo(i,k-1)
5281 dv1v = vo(i,k)
5282 dv2v = .5 * (vo(i,k) + vo(i,k-1))
5283 dv3v = vo(i,k-1)
5284 !c
5285 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
5286 tem1 = xlamud(i)
5287 !cj
5288 dellah(i,k) = dellah(i,k) + &
5289 & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h &
5290 & - tem*eta(i,k-1)*dv2h*dz &
5291 & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
5292 & ) *g/dp
5293 !cj
5294 dellaq(i,k) = dellaq(i,k) + &
5295 & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q &
5296 & - tem*eta(i,k-1)*dv2q*dz &
5297 & + tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
5298 & ) *g/dp
5299 !cj
5300 dellau(i,k) = dellau(i,k) + &
5301 & ( eta(i,k)*dv1u - eta(i,k-1)*dv3u &
5302 & - tem*eta(i,k-1)*dv2u*dz &
5303 & + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
5304 & - pgcon*eta(i,k-1)*(dv1u-dv3u) &
5305 & ) *g/dp
5306 !cj
5307 dellav(i,k) = dellav(i,k) + &
5308 & ( eta(i,k)*dv1v - eta(i,k-1)*dv3v &
5309 & - tem*eta(i,k-1)*dv2v*dz &
5310 & + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
5311 & - pgcon*eta(i,k-1)*(dv1v-dv3v) &
5312 & ) *g/dp
5313 !cj
5314 endif
5315 endif
5316 enddo
5317 enddo
5318 !c
5319 !c------- cloud top
5320 !c
5321 do i = 1, im
5322 if(cnvflg(i)) then
5323 indx = ktcon(i)
5324 dp = 1000. * del(i,indx)
5325 dv1h = heo(i,indx-1)
5326 dellah(i,indx) = eta(i,indx-1) * &
5327 & (hcko(i,indx-1) - dv1h) * g / dp
5328 dv1q = qo(i,indx-1)
5329 dellaq(i,indx) = eta(i,indx-1) * &
5330 & (qcko(i,indx-1) - dv1q) * g / dp
5331 dv1u = uo(i,indx-1)
5332 dellau(i,indx) = eta(i,indx-1) * &
5333 & (ucko(i,indx-1) - dv1u) * g / dp
5334 dv1v = vo(i,indx-1)
5335 dellav(i,indx) = eta(i,indx-1) * &
5336 & (vcko(i,indx-1) - dv1v) * g / dp
5337 !c
5338 !c cloud water
5339 !c
5340 dellal(i,indx) = eta(i,indx-1) * &
5341 & qlko_ktcon(i) * g / dp
5342 endif
5343 enddo
5344 !c
5345 !c mass flux at cloud base for shallow convection
5346 !c (Grant, 2001)
5347 !c
5348 do i= 1, im
5349 if(cnvflg(i)) then
5350 k = kbcon(i)
5351 ! ptem = g*sflx(i)*zi(i,k)/t1(i,1)
5352 ptem = g*sflx(i)*hpbl(i)/t1(i,1)
5353 wstar(i) = ptem**h1
5354 tem = po(i,k)*100. / (rd*t1(i,k))
5355 xmb(i) = betaw*tem*wstar(i)
5356 xmb(i) = min(xmb(i),xmbmax(i))
5357 endif
5358 enddo
5359 !c
5360 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5361 !c
5362 do k = 1, km
5363 do i = 1, im
5364 if (cnvflg(i) .and. k .le. kmax(i)) then
5365 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
5366 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
5367 val = 1.e-8
5368 qeso(i,k) = max(qeso(i,k), val )
5369 endif
5370 enddo
5371 enddo
5372 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5373 !c
5374 do i = 1, im
5375 delhbar(i) = 0.
5376 delqbar(i) = 0.
5377 deltbar(i) = 0.
5378 delubar(i) = 0.
5379 delvbar(i) = 0.
5380 qcond(i) = 0.
5381 enddo
5382 do k = 1, km
5383 do i = 1, im
5384 if (cnvflg(i)) then
5385 if(k.gt.kb(i).and.k.le.ktcon(i)) then
5386 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
5387 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
5388 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
5389 tem = 1./rcs(i)
5390 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
5391 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
5392 dp = 1000. * del(i,k)
5393 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
5394 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
5395 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
5396 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
5397 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
5398 endif
5399 endif
5400 enddo
5401 enddo
5402 do k = 1, km
5403 do i = 1, im
5404 if (cnvflg(i)) then
5405 if(k.gt.kb(i).and.k.le.ktcon(i)) then
5406 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
5407 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
5408 val = 1.e-8
5409 qeso(i,k) = max(qeso(i,k), val )
5410 endif
5411 endif
5412 enddo
5413 enddo
5414 !c
5415 do i = 1, im
5416 rntot(i) = 0.
5417 delqev(i) = 0.
5418 delq2(i) = 0.
5419 flg(i) = cnvflg(i)
5420 enddo
5421 do k = km, 1, -1
5422 do i = 1, im
5423 if (cnvflg(i)) then
5424 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
5425 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
5426 endif
5427 endif
5428 enddo
5429 enddo
5430 !c
5431 !c evaporating rain
5432 !c
5433 do k = km, 1, -1
5434 do i = 1, im
5435 if (k .le. kmax(i)) then
5436 deltv(i) = 0.
5437 delq(i) = 0.
5438 qevap(i) = 0.
5439 if(cnvflg(i)) then
5440 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
5441 rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
5442 endif
5443 endif
5444 if(flg(i).and.k.lt.ktcon(i)) then
5445 evef = edt(i) * evfact
5446 if(slimsk(i).eq.1.) evef=edt(i) * evfactl
5447 ! if(slimsk(i).eq.1.) evef=.07
5448 !c if(slimsk(i).ne.1.) evef = 0.
5449 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
5450 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
5451 dp = 1000. * del(i,k)
5452 if(rn(i).gt.0..and.qcond(i).lt.0.) then
5453 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
5454 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
5455 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
5456 endif
5457 if(rn(i).gt.0..and.qcond(i).lt.0..and. &
5458 & delq2(i).gt.rntot(i)) then
5459 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
5460 flg(i) = .false.
5461 endif
5462 if(rn(i).gt.0..and.qevap(i).gt.0.) then
5463 tem = .001 * dp / g
5464 tem1 = qevap(i) * tem
5465 if(tem1.gt.rn(i)) then
5466 qevap(i) = rn(i) / tem
5467 rn(i) = 0.
5468 else
5469 rn(i) = rn(i) - tem1
5470 endif
5471 q1(i,k) = q1(i,k) + qevap(i)
5472 t1(i,k) = t1(i,k) - elocp * qevap(i)
5473 deltv(i) = - elocp*qevap(i)/dt2
5474 delq(i) = + qevap(i)/dt2
5475 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
5476 endif
5477 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
5478 delqbar(i) = delqbar(i) + delq(i)*dp/g
5479 deltbar(i) = deltbar(i) + deltv(i)*dp/g
5480 endif
5481 endif
5482 enddo
5483 enddo
5484 !cj
5485 ! do i = 1, im
5486 ! if(me.eq.31.and.cnvflg(i)) then
5487 ! if(cnvflg(i)) then
5488 ! print *, ' shallow delhbar, delqbar, deltbar = ',
5489 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
5490 ! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
5491 ! print *, ' precip =', hvap*rn(i)*1000./dt2
5492 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
5493 ! endif
5494 ! enddo
5495 !cj
5496 do i = 1, im
5497 if(cnvflg(i)) then
5498 if(rn(i).lt.0..or..not.flg(i)) rn(i) = 0.
5499 ktop(i) = ktcon(i)
5500 kbot(i) = kbcon(i)
5501 kcnv(i) = 0
5502 endif
5503 enddo
5504 !c
5505 !c cloud water
5506 !c
5507 if (ncloud.gt.0) then
5508 !
5509 val1 = 1.0
5510 val2 = 0.
5511 do k = 1, km1
5512 do i = 1, im
5513 if (cnvflg(i)) then
5514 if (k.gt.kb(i).and.k.le.ktcon(i)) then
5515 tem = dellal(i,k) * xmb(i) * dt2
5516 ! tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
5517 tem1 = max(val2, min(val1, (tcr-t1(i,k))*tcrf))
5518 if (ql(i,k,2) .gt. -999.0) then
5519 ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
5520 ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
5521 else
5522 ql(i,k,1) = ql(i,k,1) + tem
5523 endif
5524 endif
5525 endif
5526 enddo
5527 enddo
5528 !
5529 endif
5530 !
5531 ! hchuang code change
5532 !
5533 do k = 1, km
5534 do i = 1, im
5535 if(cnvflg(i)) then
5536 if(k.ge.kb(i) .and. k.lt.ktop(i)) then
5537 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
5538 endif
5539 endif
5540 enddo
5541 enddo
5542 do i = 1, im
5543 if(cnvflg(i)) then
5544 k = ktop(i)-1
5545 dt_mf(i,k) = ud_mf(i,k)
5546 endif
5547 enddo
5548 !!
5549 return
5550 end subroutine shalcnv
5551 END MODULE module_cu_sas
5553 FUNCTION ran1(idum)
5554 implicit none
5555 integer idum,ia,im,iq,ir,ntab,ndiv
5556 real*8 am,eps,rnmx
5557 real*8 ran1
5558 parameter (ia=16807,im=2147483647,am=1./im,iq=127773,ir=2836, &
5559 & ntab=32,ndiv=1+(im-1)/ntab,eps=1.2e-7,rnmx=1.-eps)
5560 integer j,k,iv(32),iy,junk
5561 common /random/ junk,iv,iy
5562 data iv /ntab*0/, iy /0/
5563 if (idum.le.0.or.iy.eq.0) then
5564 idum=max(-idum,1)
5565 do j=ntab+8,1,-1
5566 k=idum/iq
5567 idum=ia*(idum-k*iq)-ir*k
5568 if (idum.lt.0) idum=idum+im
5569 if (j.le.ntab) iv(j)=idum
5570 enddo
5571 iy=iv(1)
5572 endif
5573 k=idum/iq
5574 idum=ia*(idum-k*iq)-ir*k
5575 if (idum.lt.0) idum=idum+im
5576 j=1+iy/ndiv
5577 iy=iv(j)
5578 iv(j)=idum
5579 ran1=min(am*iy,rnmx)
5580 return
5581 END FUNCTION ran1
5583 FUNCTION gasdev(idum)
5584 INTEGER idum
5585 REAL*8 gasdev
5586 !CU USES ran1
5587 INTEGER iset
5588 REAL*8 fac,gset,rsq,v1,v2,ran1
5589 SAVE iset,gset
5590 DATA iset/0/
5591 if (iset.eq.0) then
5592 1 v1=2.*ran1(idum)-1.
5593 v2=2.*ran1(idum)-1.
5594 rsq=v1**2+v2**2
5595 if(rsq.ge.1..or.rsq.eq.0.)goto 1
5596 fac=sqrt(-2.*log(rsq)/rsq)
5597 gset=v1*fac
5598 gasdev=v2*fac
5599 iset=1
5600 else
5601 gasdev=gset
5602 iset=0
5603 endif
5604 return
5605 END FUNCTION gasdev