| blob | ec11287d73f179e26a1d0e1af76d546c6f601503 |
1 !!
2 !! LOG
3 !! 2015-12-10 Weiguo Wang added gfs scale-aware SAS scheme to HWRF
4 !! A description of this updated Scale-aware SAS can be found in the HWRF Scientific Documentation:
5 !! https://dtcenter.org/HurrWRF/users/docs/scientific_documents/HWRFv3.9a_ScientificDoc.pdf
8 MODULE module_cu_scalesas
10 CONTAINS
12 !-----------------------------------------------------------------
13 SUBROUTINE CU_SCALESAS(DT,ITIMESTEP,STEPCU, &
14 RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
15 RUCUTEN,RVCUTEN, &
16 RAINCV,PRATEC,HTOP,HBOT, &
17 U3D,V3D,W,T3D,QV3D,QC3D,QI3D,PI3D,RHO3D, &
18 DZ8W,PCPS,P8W,XLAND,CU_ACT_FLAG, &
19 P_QC, &
20 MOMMIX, & ! gopal's doing
21 PGCON,sas_mass_flux, &
22 pert_sas, ens_random_seed, ens_sasamp, &
23 shalconv,shal_pgcon, &
24 HPBL2D,EVAP2D,HEAT2D, & !Kwon for shallow convection
25 P_QI,P_FIRST_SCALAR, &
26 DX2D, DY, & ! Wang W for scale-aware cnv
27 SCALEFUN, SCALEFUN1, & !cnv scale functions
28 SIGMU,SIGMU1, &
29 ids,ide, jds,jde, kds,kde, &
30 ims,ime, jms,jme, kms,kme, &
31 its,ite, jts,jte, kts,kte )
33 !-------------------------------------------------------------------
34 USE MODULE_GFS_MACHINE , ONLY : kind_phys
35 USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
36 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
37 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
38 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
39 &, EPS => con_eps, EPSM1 => con_epsm1 &
40 &, ROVCP => con_rocp, RD => con_rd
41 !-------------------------------------------------------------------
42 IMPLICIT NONE
43 !-------------------------------------------------------------------
44 !-- U3D 3D u-velocity interpolated to theta points (m/s)
45 !-- V3D 3D v-velocity interpolated to theta points (m/s)
46 !-- TH3D 3D potential temperature (K)
47 !-- T3D temperature (K)
48 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
49 !-- QC3D 3D cloud mixing ratio (Kg/Kg)
50 !-- QI3D 3D ice mixing ratio (Kg/Kg)
51 !-- P8w 3D pressure at full levels (Pa)
52 !-- Pcps 3D pressure (Pa)
53 !-- PI3D 3D exner function (dimensionless)
54 !-- rr3D 3D dry air density (kg/m^3)
55 !-- RUBLTEN U tendency due to
56 ! PBL parameterization (m/s^2)
57 !-- RVBLTEN V tendency due to
58 ! PBL parameterization (m/s^2)
59 !-- RTHBLTEN Theta tendency due to
60 ! PBL parameterization (K/s)
61 !-- RQVBLTEN Qv tendency due to
62 ! PBL parameterization (kg/kg/s)
63 !-- RQCBLTEN Qc tendency due to
64 ! PBL parameterization (kg/kg/s)
65 !-- RQIBLTEN Qi tendency due to
66 ! PBL parameterization (kg/kg/s)
67 !
68 !-- MOMMIX MOMENTUM MIXING COEFFICIENT (can be set in the namelist)
69 !-- RUCUTEN U tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
70 !-- RVCUTEN V tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
71 !
72 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
73 !-- GRAV acceleration due to gravity (m/s^2)
74 !-- ROVCP R/CP
75 !-- RD gas constant for dry air (J/kg/K)
76 !-- ROVG R/G
77 !-- P_QI species index for cloud ice
78 !-- dz8w dz between full levels (m)
79 !-- z height above sea level (m)
80 !-- PSFC pressure at the surface (Pa)
81 !-- UST u* in similarity theory (m/s)
82 !-- PBL PBL height (m)
83 !-- PSIM similarity stability function for momentum
84 !-- PSIH similarity stability function for heat
85 !-- HFX upward heat flux at the surface (W/m^2)
86 !-- QFX upward moisture flux at the surface (kg/m^2/s)
87 !-- TSK surface temperature (K)
88 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
89 !-- WSPD wind speed at lowest model level (m/s)
90 !-- BR bulk Richardson number in surface layer
91 !-- DT time step (s)
92 !-- rvovrd R_v divided by R_d (dimensionless)
93 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
94 !-- KARMAN Von Karman constant
95 !-- ids start index for i in domain
96 !-- ide end index for i in domain
97 !-- jds start index for j in domain
98 !-- jde end index for j in domain
99 !-- kds start index for k in domain
100 !-- kde end index for k in domain
101 !-- ims start index for i in memory
102 !-- ime end index for i in memory
103 !-- jms start index for j in memory
104 !-- jme end index for j in memory
105 !-- kms start index for k in memory
106 !-- kme end index for k in memory
107 !-- its start index for i in tile
108 !-- ite end index for i in tile
109 !-- jts start index for j in tile
110 !-- jte end index for j in tile
111 !-- kts start index for k in tile
112 !-- kte end index for k in tile
113 !-------------------------------------------------------------------
115 INTEGER :: ICLDCK
117 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
118 ims,ime, jms,jme, kms,kme, &
119 its,ite, jts,jte, kts,kte, &
120 ITIMESTEP, & !NSTD
121 P_FIRST_SCALAR, &
122 P_QC, &
123 P_QI, &
124 STEPCU
126 REAL, INTENT(IN) :: &
127 DT
129 REAL, INTENT(IN) :: DY
130 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DX2D
131 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: SCALEFUN,SCALEFUN1, & !updaft area fraction for deep and shallow cnv
132 SIGMU, SIGMU1 !updaft area fraction for deep and shallow cnv
135 REAL, OPTIONAL, INTENT(IN) :: PGCON,sas_mass_flux,shal_pgcon
136 INTEGER, OPTIONAL, INTENT(IN) :: shalconv
137 REAL(kind=kind_phys) :: PGCON_USE,SHAL_PGCON_USE,massf
138 logical,optional,intent(in) :: pert_sas
139 integer,optional,intent(in) :: ens_random_seed
140 real,optional,intent(in) :: ens_sasamp
141 INTEGER :: shalconv_use
142 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
143 RQCCUTEN, &
144 RQICUTEN, &
145 RQVCUTEN, &
146 RTHCUTEN
147 REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(INOUT) :: &
148 RUCUTEN, &
149 RVCUTEN
150 REAL, OPTIONAL, INTENT(IN) :: MOMMIX
152 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
153 INTENT(IN) :: HPBL2D,EVAP2D,HEAT2D !Kwon for sha
155 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
156 XLAND
158 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
159 RAINCV, PRATEC
161 REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
162 HBOT, &
163 HTOP
165 LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) :: &
166 CU_ACT_FLAG
169 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
170 DZ8W, &
171 P8w, &
172 Pcps, &
173 PI3D, &
174 QC3D, &
175 QI3D, &
176 QV3D, &
177 RHO3D, &
178 T3D, &
179 U3D, &
180 V3D, &
181 W
183 !--------------------------- LOCAL VARS ------------------------------
185 REAL, DIMENSION(ims:ime, jms:jme) :: &
186 PSFC
188 REAL, DIMENSION(its:ite, jts:jte) :: &
189 RAINCV1, PRATEC1
190 REAL, DIMENSION(its:ite, jts:jte) :: &
191 RAINCV2, PRATEC2
193 REAL (kind=kind_phys) :: &
194 DELT, &
195 DPSHC, &
196 RDELT, &
197 RSEED
199 REAL (kind=kind_phys), DIMENSION(its:ite) :: &
200 CLDWRK, &
201 PS, &
202 RCS, &
203 RN, &
204 SLIMSK, &
205 HPBL,EVAP,HEAT ! for shallow convection
207 REAL (kind=kind_phys), DIMENSION(its:ite) :: garea ! grid box area in m^2,
208 ! scale-aware cnv
210 REAL (kind=kind_phys), DIMENSION(its:ite) :: SCALEFUN_out,SCALEFUN1_out,& !updraft area fraction for deep and shallow cnv
211 SIGMU_out,SIGMU1_out !updraft area fraction for deep and shallow cnv
214 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
215 PRSI
217 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
218 DEL, &
219 DOT, &
220 PHIL, &
221 PRSL, &
222 PRSLK, &
223 Q1, &
224 T1, &
225 U1, &
226 V1, &
227 ZI, &
228 ZL
230 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
231 cnvw,cnvc,ud_mf,dd_mf,dt_mf ! *_mf not useful for HWRF. it is for transport
233 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) :: &
234 QL
236 INTEGER, DIMENSION(its:ite) :: &
237 KBOT, &
238 KTOP, &
239 KCNV
241 INTEGER :: &
242 I, &
243 IGPVS, &
244 IM, &
245 J, &
246 JCAP, &
247 K, &
248 KM, &
249 KP, &
250 KX, &
251 NCLOUD
253 DATA IGPVS/0/
255 !-----------------------------------------------------------------------
256 !
257 ! write(0,*)' in scale-aware sas'
259 if(present(shalconv)) then
260 shalconv_use=shalconv
261 else
262 #if (NMM_CORE==1)
263 shalconv_use=0
264 #else
265 #if (EM_CORE==1)
266 shalconv_use=1
267 #else
268 shalconv_use=0
269 #endif
270 #endif
271 endif
273 if(present(pgcon)) then
274 pgcon_use = pgcon
275 else
276 ! pgcon_use = 0.7 ! Gregory et al. (1997, QJRMS)
277 pgcon_use = 0.55 ! Zhang & Wu (2003,JAS), used in GFS (25km res spectral)
278 ! pgcon_use = 0.2 ! HWRF, for model tuning purposes
279 ! pgcon_use = 0.3 ! GFDL, or so I am told
281 ! For those attempting to tune pgcon:
283 ! The value of 0.55 comes from an observational study of
284 ! synoptic-scale deep convection and 0.7 came from an
285 ! incorrect fit to the same data. That value is likely
286 ! correct for deep convection at gridscales near that of GFS,
287 ! but is questionable in shallow convection, or for scales
288 ! much finer than synoptic scales.
290 ! Then again, the assumptions of SAS break down when the
291 ! gridscale is near the convection scale anyway. In a large
292 ! storm such as a hurricane, there is often no environment to
293 ! detrain into since adjancent gridsquares are also undergoing
294 ! active convection. Each gridsquare will no longer have many
295 ! updrafts and downdrafts. At sub-convective timescales, you
296 ! will find unstable columns for many (say, 5 second length)
297 ! timesteps in a real atmosphere during a convection cell's
298 ! lifetime, so forcing it to be neutrally stable is unphysical.
300 ! Hence, in scales near the convection scale (cells have
301 ! ~0.5-4km diameter in hurricanes), this parameter is more of a
302 ! tuning parameter to get a scheme that is inappropriate for
303 ! that resolution to do a reasonable job.
305 ! Your mileage might vary.
307 ! - Sam Trahan
308 endif
310 if(present(sas_mass_flux)) then
311 massf=sas_mass_flux
312 ! Use this to reduce the fluxes added by SAS to prevent
313 ! computational instability as a result of large fluxes.
314 else
315 massf=9e9 ! large number to disable check
316 endif
318 if(present(shal_pgcon)) then
319 if(shal_pgcon>=0) then
320 shal_pgcon_use = shal_pgcon
321 else
322 ! shal_pgcon<0 means use deep pgcon
323 shal_pgcon_use = pgcon_use
324 endif
325 else
326 ! Default: Same as deep convection pgcon
327 shal_pgcon_use = pgcon_use
328 ! Read the warning above though. It may be advisable for
329 ! these to be different.
330 endif
332 DO J=JTS,JTE
333 DO I=ITS,ITE
334 CU_ACT_FLAG(I,J)=.TRUE.
335 ENDDO
336 ENDDO
338 IM=ITE-ITS+1
339 KX=KTE-KTS+1
340 JCAP=126
341 DPSHC=30_kind_phys
342 DELT=DT*STEPCU
343 RDELT=1./DELT
344 NCLOUD=1
347 DO J=jms,jme
348 DO I=ims,ime
349 PSFC(i,j)=P8w(i,kms,j)
350 ENDDO
351 ENDDO
353 if(igpvs.eq.0) CALL GFUNCPHYS
354 igpvs=1
356 !------------- J LOOP (OUTER) --------------------------------------------------
358 big_outer_j_loop: DO J=jts,jte
360 ! --------------- compute zi and zl -----------------------------------------
361 DO i=its,ite
362 ZI(I,KTS)=0.0
363 ENDDO
365 DO k=kts+1,kte
366 KM=K-1
367 DO i=its,ite
368 ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
369 ENDDO
370 ENDDO
372 DO k=kts+1,kte
373 KM=K-1
374 DO i=its,ite
375 ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
376 ENDDO
377 ENDDO
379 DO i=its,ite
380 ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
381 ENDDO
383 ! --------------- end compute zi and zl -------------------------------------
385 DO i=its,ite
386 !!PS(i)=PSFC(i,j)*.001
387 PS(i)=PSFC(i,j) ! wang, in Pa
388 RCS(i)=1.
389 SLIMSK(i)=ABS(XLAND(i,j)-2.)
391 garea(I)=DX2D(i,j)*DY*2.0 ! grid box area in m^2
392 KCNV(I)=0 ! initialize KCNV
393 ENDDO
395 #if (NMM_CORE == 1)
396 if(shalconv_use==1) then
397 DO i=its,ite
398 HPBL(I) = HPBL2D(I,J) ! for shallow convection
399 EVAP(I) = EVAP2D(I,J) ! for shallow convection
400 HEAT(I) = HEAT2D(I,J) ! for shallow convection
401 ENDDO
402 endif
403 #endif
405 DO i=its,ite
406 PRSI(i,kts)=PS(i)
407 ENDDO
409 DO k=kts,kte
410 kp=k+1
411 DO i=its,ite
412 !PRSL(I,K)=Pcps(i,k,j)*.001
413 PRSL(I,K)=Pcps(i,k,j) ! in Pa
414 PHIL(I,K)=ZL(I,K)*GRAV
415 !DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
416 DOT(i,k)=-0.5*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j)) ! in Pa
417 ENDDO
418 ENDDO
420 DO k=kts,kte
421 DO i=its,ite
422 DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
423 U1(i,k)=U3D(i,k,j)
424 V1(i,k)=V3D(i,k,j)
425 Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
426 T1(i,k)=T3D(i,k,j)
427 QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
428 QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
429 !PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
430 PRSLK(I,K)=(PRSL(i,k)*1.0e-5)**ROVCP ! prsl in Pa
431 ENDDO
432 ENDDO
434 DO k=kts+1,kte+1
435 km=k-1
436 DO i=its,ite
437 PRSI(i,k)=PRSI(i,km)-del(i,km)
438 ENDDO
439 ENDDO
441 ! write(0,*)'dx2d=',dx2d(its,jts:jts+5)
442 ! write(0,*)'dy=',dy
443 ! write(0,*)'ps=',ps(its)
444 ! write(0,*)'del=',del(its,kts:kts+5)
445 ! write(0,*)'prsl=',prsl(its,kts:kts+5)
446 ! write(0,*)'dot=',dot(its,kts:kts+5)
448 ! 2016-03-22 near final version of scale-aware convection
449 call mfdeepcnv(im,im,kx,delt,del,prsl,ps,phil,ql, &
450 & q1,t1,u1,v1,cldwrk,rn,kbot,ktop,kcnv,nint(slimsk),garea, &
451 & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc, SIGMU_out,SCALEFUN_out, &
452 & pert_sas, ens_random_seed, ens_sasamp)
454 do i=its,ite
455 RAINCV1(I,J)=RN(I)*1000./STEPCU
456 PRATEC1(I,J)=RN(I)*1000./(STEPCU * DT)
457 enddo
458 !
459 do i=its,ite
460 RAINCV2(I,J)=0.
461 PRATEC2(I,J)=0.
462 enddo
463 !
464 !
465 if_shallow_conv: if(shalconv_use==1) then
466 #if (NMM_CORE == 1)
467 ! NMM calls the new shallow convection developed by J Han
468 ! (Added to WRF by Y.Kwon)
470 ! 2016-03-22 near final version of scale-aware convection
471 call mfshalcnv(im,im,kx,delt,del,prsl,ps,phil,ql, &
472 & q1,t1,u1,v1,rn,kbot,ktop,kcnv,nint(slimsk),garea, &
473 & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc,SIGMU1_out,SCALEFUN1_out, &
474 & pert_sas, ens_random_seed, ens_sasamp)
478 DO I=ITS,ITE
479 RAINCV2(I,J)=RN(I)*1000./STEPCU
480 PRATEC2(I,J)=RN(I)*1000./(STEPCU * DT)
481 ENDDO
482 !
483 #else
484 #if (EM_CORE == 1)
485 ! NOTE: ARW should be able to call the new shalcnv here, but
486 ! they need to add the three new variables, so I'm leaving the
487 ! old shallow convection call here - Sam Trahan
489 !! removed CALL OLD_ARW_SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KCNV,Q1,T1,DPSHC)
490 #else
491 ! Shallow convection is untested for other cores.
492 #endif
493 #endif
494 endif if_shallow_conv
496 !! output updraft area fractions for deep and shallow cnv
497 do i=its,ite
498 SCALEFUN(I,J)=SCALEFUN_OUT(i)
499 SCALEFUN1(I,J)=SCALEFUN1_OUT(i)
500 SIGMU(I,J)=SIGMU_OUT(i)
501 SIGMU1(I,J)=SIGMU1_OUT(i)
502 enddo
505 DO I=ITS,ITE
506 RAINCV(I,J)= RAINCV1(I,J) + RAINCV2(I,J)
507 PRATEC(I,J)= PRATEC1(I,J) + PRATEC2(I,J)
508 HBOT(I,J)=KBOT(I)
509 HTOP(I,J)=KTOP(I)
510 ENDDO
512 DO K=KTS,KTE
513 DO I=ITS,ITE
514 RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
515 RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
516 ENDDO
517 ENDDO
519 !===============================================================================
520 ! ADD MOMENTUM MIXING TERM AS TENDENCIES. This is gopal's doing for SAS
521 ! MOMMIX is the reduction factor set to 0.7 by default. Because NMM has
522 ! divergence damping term, a reducion factor for cumulum mixing may be
523 ! required otherwise storms were too weak.
524 !===============================================================================
525 !
526 #if (NMM_CORE == 1)
527 DO K=KTS,KTE
528 DO I=ITS,ITE
529 ! RUCUTEN(I,J,K)=MOMMIX*(U1(I,K)-U3D(I,K,J))*RDELT
530 ! RVCUTEN(I,J,K)=MOMMIX*(V1(I,K)-V3D(I,K,J))*RDELT
531 RUCUTEN(I,J,K)=(U1(I,K)-U3D(I,K,J))*RDELT
532 RVCUTEN(I,J,K)=(V1(I,K)-V3D(I,K,J))*RDELT
533 ENDDO
534 ENDDO
535 #endif
538 IF(P_QC .ge. P_FIRST_SCALAR)THEN
539 DO K=KTS,KTE
540 DO I=ITS,ITE
541 RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
542 ENDDO
543 ENDDO
544 ENDIF
546 IF(P_QI .ge. P_FIRST_SCALAR)THEN
547 DO K=KTS,KTE
548 DO I=ITS,ITE
549 RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
550 ENDDO
551 ENDDO
552 ENDIF
554 ENDDO big_outer_j_loop ! Outer most J loop
556 END SUBROUTINE CU_SCALESAS
558 !====================================================================
559 SUBROUTINE scalesasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
560 RUCUTEN,RVCUTEN, &
561 RESTART,P_QC,P_QI,P_FIRST_SCALAR, &
562 allowed_to_read, &
563 ids, ide, jds, jde, kds, kde, &
564 ims, ime, jms, jme, kms, kme, &
565 its, ite, jts, jte, kts, kte )
566 !--------------------------------------------------------------------
567 IMPLICIT NONE
568 !--------------------------------------------------------------------
569 LOGICAL , INTENT(IN) :: allowed_to_read,restart
570 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
571 ims, ime, jms, jme, kms, kme, &
572 its, ite, jts, jte, kts, kte
573 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
575 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
576 RTHCUTEN, &
577 RQVCUTEN, &
578 RQCCUTEN, &
579 RQICUTEN
580 REAL, DIMENSION( ims:ime , jms:jme , kms:kme ) , INTENT(OUT) :: &
581 RUCUTEN, & ! gopal's doing for SAS
582 RVCUTEN
584 INTEGER :: i, j, k, itf, jtf, ktf
586 jtf=min0(jte,jde-1)
587 ktf=min0(kte,kde-1)
588 itf=min0(ite,ide-1)
590 #ifdef HWRF
591 !zhang's doing
592 IF(.not.restart .or. .not.allowed_to_read)THEN
593 !end of zhang's doing
594 #else
595 IF(.not.restart)THEN
596 #endif
597 DO j=jts,jtf
598 DO k=kts,ktf
599 DO i=its,itf
600 RTHCUTEN(i,k,j)=0.
601 RQVCUTEN(i,k,j)=0.
602 RUCUTEN(i,j,k)=0.
603 RVCUTEN(i,j,k)=0.
604 ENDDO
605 ENDDO
606 ENDDO
608 IF (P_QC .ge. P_FIRST_SCALAR) THEN
609 DO j=jts,jtf
610 DO k=kts,ktf
611 DO i=its,itf
612 RQCCUTEN(i,k,j)=0.
613 ENDDO
614 ENDDO
615 ENDDO
616 ENDIF
618 IF (P_QI .ge. P_FIRST_SCALAR) THEN
619 DO j=jts,jtf
620 DO k=kts,ktf
621 DO i=its,itf
622 RQICUTEN(i,k,j)=0.
623 ENDDO
624 ENDDO
625 ENDDO
626 ENDIF
627 ENDIF
629 END SUBROUTINE scalesasinit
631 !-----------------------------------------------------------------------
632 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
633 !!! scale aware SAS
634 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
636 subroutine mfdeepcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, &
637 & q1,t1,u1,v1,cldwrk,rn,kbot,ktop,kcnv,islimsk,garea, &
638 & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc, sigmuout,scaldfunc, &
639 & pert_sas, ens_random_seed, ens_sasamp)
641 !
642 ! use machine , only : kind_phys
643 ! use funcphys , only : fpvs
644 ! use physcons, grav => con_g, cp => con_cp, hvap => con_hvap
645 USE MODULE_GFS_MACHINE, ONLY : kind_phys
646 USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
647 USE MODULE_GFS_PHYSCONS, grav => con_g, cp => con_cp &
648 & , hvap => con_hvap &
649 &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
650 &, rd => con_rd, cvap => con_cvap, cliq => con_cliq &
651 &, eps => con_eps, epsm1 => con_epsm1
655 implicit none
656 !
657 integer im, ix, km, ncloud, &
658 & kbot(im), ktop(im), kcnv(im)
659 ! &, me
660 real(kind=kind_phys) delt
661 real(kind=kind_phys) psp(im), delp(ix,km), prslp(ix,km)
662 real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
663 & ql(ix,km,2),q1(ix,km), t1(ix,km), &
664 & u1(ix,km), v1(ix,km), &
665 ! & u1(ix,km), v1(ix,km), rcs(im),
666 & cldwrk(im), rn(im), garea(im), &
667 & dot(ix,km), phil(ix,km), &
668 & cnvw(ix,km),cnvc(ix,km), &
669 ! hchuang code change mass flux output
670 & ud_mf(im,km),dd_mf(im,km),dt_mf(im,km)
671 !
672 integer i, indx, jmn, k, kk, km1, n
673 integer, dimension(im), intent(in) :: islimsk
674 ! integer latd,lond
675 !
676 real(kind=kind_phys) clam, cxlamu, cxlamd, &
677 & xlamde, xlamdd, &
678 & crtlamu, crtlamd
679 !
680 ! real(kind=kind_phys) detad
681 real(kind=kind_phys) adw, aup, aafac, &
682 & beta, betal, betas, &
683 & c0l, c0s, d0, &
684 & c1, asolfac, &
685 & dellat, delta, desdt, dg, &
686 & dh, dhh, dp, &
687 & dq, dqsdp, dqsdt, dt, &
688 & dt2, dtmax, dtmin, &
689 & dxcrtas, dxcrtuf, &
690 & dv1h, dv2h, dv3h, &
691 & dv1q, dv2q, dv3q, &
692 & dz, dz1, e1, edtmax, &
693 & edtmaxl, edtmaxs, el2orc, elocp, &
694 & es, etah, &
695 & cthk, dthk, &
696 & evef, evfact, evfactl, fact1, &
697 & fact2, factor, &
698 & g, gamma, pprime, cm, &
699 & qlk, qrch, qs, &
700 & rain, rfact, shear, tfac, &
701 & val, val1, val2, &
702 & w1, w1l, w1s, w2, &
703 & w2l, w2s, w3, w3l, &
704 & w3s, w4, w4l, w4s, &
705 & rho, betaw, &
706 & xdby, xpw, xpwd, &
707 ! & xqrch, mbdt, tem, &
708 & xqrch, tem, tem1, tem2, &
709 & ptem, ptem1, ptem2, &
710 & pgcon
711 !
712 logical,optional,intent(in) :: pert_sas
713 integer,optional,intent(in) :: ens_random_seed
714 real,optional,intent(in) :: ens_sasamp
716 integer kb(im), kbcon(im), kbcon1(im), &
717 & ktcon(im), ktcon1(im), ktconn(im), &
718 & jmin(im), lmin(im), kbmax(im), &
719 & kbm(im), kmax(im)
720 !
721 ! real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im), &
722 real(kind=kind_phys) aa1(im), &
723 & umean(im), tauadv(im), gdx(im), &
724 & delhbar(im), delq(im), delq2(im), &
725 & delqbar(im), delqev(im), deltbar(im), &
726 & deltv(im), dtconv(im), edt(im), &
727 & edto(im), edtx(im), fld(im), &
728 & hcdo(im,km), hmax(im), hmin(im), &
729 & ucdo(im,km), vcdo(im,km),aa2(im), &
730 & pdot(im), po(im,km), &
731 & pwavo(im), pwevo(im), mbdt(im), &
732 & qcdo(im,km), qcond(im), qevap(im), &
733 & rntot(im), vshear(im), xaa0(im), &
734 & xk(im), xlamd(im), cina(im), &
735 & xmb(im), xmbmax(im), xpwav(im), &
736 & xpwev(im), xlamx(im), &
737 & delubar(im),delvbar(im)
738 !
739 real(kind=kind_phys) c0(im)
740 !cj
741 real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn, &
742 & cinacr, cinacrmx, cinacrmn
743 !cj
744 !
745 ! parameters for updraft velocity calculation
746 real(kind=kind_phys) bet1, cd1, f1, gam1, &
747 & bb1, bb2, wucb
748 !
749 !c physical parameters
750 parameter(g=grav,asolfac=0.89)
751 parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
752 parameter(c0s=.002,c1=.002,d0=.01)
753 parameter(c0l=c0s*asolfac)
754 !
755 ! asolfac: aerosol-aware parameter based on Lim & Hong (2012)
756 ! asolfac= cx / c0s(=.002)
757 ! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s)
758 ! Nccn: CCN number concentration in cm^(-3)
759 ! Until a realistic Nccn is provided, typical Nccns are assumed
760 ! as Nccn=100 for sea and Nccn=7000 for land
761 !
762 parameter(cm=1.0,delta=fv)
763 parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
764 parameter(cthk=200.,dthk=25.)
765 parameter(cinpcrmx=180.,cinpcrmn=120.)
766 parameter(cinacrmx=-120.,cinacrmn=-120.)
767 parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
768 parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3)
769 !
770 ! local variables and arrays
771 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), &
772 & uo(im,km), vo(im,km), qeso(im,km)
773 ! for updraft velocity calculation
774 real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km)
775 real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im)
777 real(kind=kind_phys) sigmuout(im)
778 !
779 !c cloud water
780 ! real(kind=kind_phys) tvo(im,km)
781 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km), &
782 & dbyo(im,km), zo(im,km), &
783 & xlamue(im,km), xlamud(im,km), &
784 & fent1(im,km), fent2(im,km), frh(im,km), &
785 & heo(im,km), heso(im,km), &
786 & qrcd(im,km), dellah(im,km), dellaq(im,km), &
787 & dellau(im,km), dellav(im,km), hcko(im,km), &
788 & ucko(im,km), vcko(im,km), qcko(im,km), &
789 & eta(im,km), etad(im,km), zi(im,km), &
790 & qrcko(im,km), qrcdo(im,km), &
791 & pwo(im,km), pwdo(im,km), c0t(im,km), &
792 & tx1(im), sumx(im), cnvwt(im,km)
793 ! &, rhbar(im)
794 !
795 logical totflg, cnvflg(im), asqecflg(im), flg(im)
796 !
797 ! asqecflg: flag for the quasi-equilibrium assumption of Arakawa-Schubert
798 !
799 ! real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
800 !! save pcrit, acritt
801 ! data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,
802 ! & 350.,300.,250.,200.,150./
803 ! data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,
804 ! & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
805 !c gdas derived acrit
806 !c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
807 !c & .743,.813,.886,.947,1.138,1.377,1.896/
808 real(kind=kind_phys) tf, tcr, tcrf
809 parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
812 #if HWRF==1
813 real*8 :: gasdev,ran1 !zhang
814 real :: rr !zhang
815 logical,save :: pert_sas_local !zhang
816 integer,save :: ens_random_seed_local,env_pp_local !zhang
817 integer :: ensda_physics_pert !zhang
818 real,save :: ens_sasamp_local !zhang
819 data ens_random_seed_local/0/
820 data env_pp_local/0/
821 CHARACTER(len=3) :: env_memb,env_pp
822 if ( ens_random_seed_local .eq. 0 ) then
823 CALL nl_get_ensda_physics_pert(1,ensda_physics_pert)
824 ens_random_seed_local=ens_random_seed
825 env_pp_local=ensda_physics_pert
826 pert_sas_local=.false.
827 ens_sasamp_local=0.0
828 ! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99
829 if ( env_pp_local .eq. 1 ) then
830 if ( ens_random_seed .ne. 99 ) then
831 pert_sas_local=.true.
832 ens_sasamp_local=ens_sasamp
833 else
834 ! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero
835 ens_random_seed_local=ens_random_seed
836 pert_sas_local=pert_sas
837 ens_sasamp_local=ens_sasamp
838 endif
839 else
840 ens_random_seed_local=ens_random_seed
841 pert_sas_local=pert_sas
842 ens_sasamp_local=ens_sasamp
843 endif
844 print*, "DESAS ==", ens_random_seed_local,pert_sas_local,ens_sasamp_local,ensda_physics_pert
845 endif
846 #endif
847 !
848 !c-----------------------------------------------------------------------
849 !
850 !************************************************************************
851 ! convert input Pa terms to Cb terms -- Moorthi
852 ps = psp * 0.001
853 prsl = prslp * 0.001
854 del = delp * 0.001
855 !************************************************************************
856 !
857 !
858 km1 = km - 1
859 !c
860 !c initialize arrays
861 !c
862 do i=1,im
863 cnvflg(i) = .true.
864 rn(i)=0.
865 mbdt(i)=10.
866 kbot(i)=km+1
867 ktop(i)=0
868 kbcon(i)=km
869 ktcon(i)=1
870 ktconn(i)=1
871 dtconv(i) = 3600.
872 cldwrk(i) = 0.
873 pdot(i) = 0.
874 lmin(i) = 1
875 jmin(i) = 1
876 qlko_ktcon(i) = 0.
877 edt(i) = 0.
878 edto(i) = 0.
879 edtx(i) = 0.
880 ! acrt(i) = 0.
881 ! acrtfct(i) = 1.
882 aa1(i) = 0.
883 aa2(i) = 0.
884 xaa0(i) = 0.
885 cina(i) = 0.
886 pwavo(i)= 0.
887 pwevo(i)= 0.
888 xpwav(i)= 0.
889 xpwev(i)= 0.
890 vshear(i) = 0.
891 gdx(i) = sqrt(garea(i))
893 scaldfunc(i)=-1.0 ! initialized wang
894 sigmagfm(i)=-1.0
895 sigmuout(i)=-1.0
897 enddo
898 !
899 do i=1,im
900 if(islimsk(i) == 1) then
901 c0(i) = c0l
902 else
903 c0(i) = c0s
904 endif
905 enddo
906 do k = 1, km
907 do i = 1, im
908 if(t1(i,k) > 273.16) then
909 c0t(i,k) = c0(i)
910 else
911 tem = d0 * (t1(i,k) - 273.16)
912 tem1 = exp(tem)
913 c0t(i,k) = c0(i) * tem1
914 endif
915 enddo
916 enddo
917 !
918 do k = 1, km
919 do i = 1, im
920 cnvw(i,k) = 0.
921 cnvc(i,k) = 0.
922 enddo
923 enddo
924 ! hchuang code change
925 do k = 1, km
926 do i = 1, im
927 ud_mf(i,k) = 0.
928 dd_mf(i,k) = 0.
929 dt_mf(i,k) = 0.
930 enddo
931 enddo
932 !c
933 ! do k = 1, 15
934 ! acrit(k) = acritt(k) * (975. - pcrit(k))
935 ! enddo
936 !
937 dt2 = delt
938 ! val = 1200.
939 val = 600.
940 dtmin = max(dt2, val )
941 ! val = 5400.
942 val = 10800.
943 dtmax = max(dt2, val )
944 !c model tunable parameters are all here
945 edtmaxl = .3
946 edtmaxs = .3
947 clam = .1
948 aafac = .1
949 ! betal = .15
950 ! betas = .15
951 betal = .05
952 betas = .05
953 !c evef = 0.07
954 evfact = 0.3
955 evfactl = 0.3
956 !
957 crtlamu = 1.0e-4
958 crtlamd = 1.0e-4
959 !
960 cxlamu = 1.0e-3
961 cxlamd = 1.0e-4
962 xlamde = 1.0e-4
963 xlamdd = 1.0e-4
964 !
965 ! pgcon = 0.7 ! Gregory et al. (1997, QJRMS)
966 pgcon = 0.55 ! Zhang & Wu (2003,JAS)
967 !
968 w1l = -8.e-3
969 w2l = -4.e-2
970 w3l = -5.e-3
971 w4l = -5.e-4
972 w1s = -2.e-4
973 w2s = -2.e-3
974 w3s = -1.e-3
975 w4s = -2.e-5
976 !c
977 !c define top layer for search of the downdraft originating layer
978 !c and the maximum thetae for updraft
979 !c
980 do i=1,im
981 kbmax(i) = km
982 kbm(i) = km
983 kmax(i) = km
984 tx1(i) = 1.0 / ps(i)
985 enddo
986 !
987 do k = 1, km
988 do i=1,im
989 if (prsl(i,k)*tx1(i) > 0.04) kmax(i) = k + 1
990 if (prsl(i,k)*tx1(i) > 0.45) kbmax(i) = k + 1
991 if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1
992 enddo
993 enddo
994 do i=1,im
995 kmax(i) = min(km,kmax(i))
996 kbmax(i) = min(kbmax(i),kmax(i))
997 kbm(i) = min(kbm(i),kmax(i))
998 enddo
999 !c
1000 !c hydrostatic height assume zero terr and initially assume
1001 !c updraft entrainment rate as an inverse function of height
1002 !c
1003 do k = 1, km
1004 do i=1,im
1005 zo(i,k) = phil(i,k) / g
1006 enddo
1007 enddo
1008 do k = 1, km1
1009 do i=1,im
1010 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
1011 xlamue(i,k) = clam / zi(i,k)
1012 ! xlamue(i,k) = max(xlamue(i,k), crtlamu)
1013 enddo
1014 enddo
1015 !c
1016 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1017 !c convert surface pressure to mb from cb
1018 !c
1019 do k = 1, km
1020 do i = 1, im
1021 if (k <= kmax(i)) then
1022 pfld(i,k) = prsl(i,k) * 10.0
1023 eta(i,k) = 1.
1024 fent1(i,k)= 1.
1025 fent2(i,k)= 1.
1026 frh(i,k) = 0.
1027 hcko(i,k) = 0.
1028 qcko(i,k) = 0.
1029 qrcko(i,k)= 0.
1030 ucko(i,k) = 0.
1031 vcko(i,k) = 0.
1032 etad(i,k) = 1.
1033 hcdo(i,k) = 0.
1034 qcdo(i,k) = 0.
1035 ucdo(i,k) = 0.
1036 vcdo(i,k) = 0.
1037 qrcd(i,k) = 0.
1038 qrcdo(i,k)= 0.
1039 dbyo(i,k) = 0.
1040 pwo(i,k) = 0.
1041 pwdo(i,k) = 0.
1042 dellal(i,k) = 0.
1043 to(i,k) = t1(i,k)
1044 qo(i,k) = q1(i,k)
1045 uo(i,k) = u1(i,k)
1046 vo(i,k) = v1(i,k)
1047 ! uo(i,k) = u1(i,k) * rcs(i)
1048 ! vo(i,k) = v1(i,k) * rcs(i)
1049 wu2(i,k) = 0.
1050 buo(i,k) = 0.
1051 drag(i,k) = 0.
1052 cnvwt(i,k)= 0.
1053 endif
1054 enddo
1055 enddo
1056 !c
1057 !c column variables
1058 !c p is pressure of the layer (mb)
1059 !c t is temperature at t-dt (k)..tn
1060 !c q is mixing ratio at t-dt (kg/kg)..qn
1061 !c to is temperature at t+dt (k)... this is after advection and turbulan
1062 !c qo is mixing ratio at t+dt (kg/kg)..q1
1063 !c
1064 do k = 1, km
1065 do i=1,im
1066 if (k <= kmax(i)) then
1067 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1068 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
1069 val1 = 1.e-8
1070 qeso(i,k) = max(qeso(i,k), val1)
1071 val2 = 1.e-10
1072 qo(i,k) = max(qo(i,k), val2 )
1073 ! qo(i,k) = min(qo(i,k),qeso(i,k))
1074 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
1075 endif
1076 enddo
1077 enddo
1078 !c
1079 !c compute moist static energy
1080 !c
1081 do k = 1, km
1082 do i=1,im
1083 if (k <= kmax(i)) then
1084 ! tem = g * zo(i,k) + cp * to(i,k)
1085 tem = phil(i,k) + cp * to(i,k)
1086 heo(i,k) = tem + hvap * qo(i,k)
1087 heso(i,k) = tem + hvap * qeso(i,k)
1088 !c heo(i,k) = min(heo(i,k),heso(i,k))
1089 endif
1090 enddo
1091 enddo
1092 !c
1093 !c determine level with largest moist static energy
1094 !c this is the level where updraft starts
1095 !c
1096 do i=1,im
1097 hmax(i) = heo(i,1)
1098 kb(i) = 1
1099 enddo
1100 do k = 2, km
1101 do i=1,im
1102 if (k <= kbm(i)) then
1103 if(heo(i,k) > hmax(i)) then
1104 kb(i) = k
1105 hmax(i) = heo(i,k)
1106 endif
1107 endif
1108 enddo
1109 enddo
1110 !c
1111 do k = 1, km1
1112 do i=1,im
1113 if (k <= kmax(i)-1) then
1114 dz = .5 * (zo(i,k+1) - zo(i,k))
1115 dp = .5 * (pfld(i,k+1) - pfld(i,k))
1116 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
1117 pprime = pfld(i,k+1) + epsm1 * es
1118 qs = eps * es / pprime
1119 dqsdp = - qs / pprime
1120 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1121 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
1122 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1123 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
1124 dq = dqsdt * dt + dqsdp * dp
1125 to(i,k) = to(i,k+1) + dt
1126 qo(i,k) = qo(i,k+1) + dq
1127 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
1128 endif
1129 enddo
1130 enddo
1131 !
1132 do k = 1, km1
1133 do i=1,im
1134 if (k <= kmax(i)-1) then
1135 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
1136 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
1137 val1 = 1.e-8
1138 qeso(i,k) = max(qeso(i,k), val1)
1139 val2 = 1.e-10
1140 qo(i,k) = max(qo(i,k), val2 )
1141 ! qo(i,k) = min(qo(i,k),qeso(i,k))
1142 frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), 1._kind_phys)
1143 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
1144 & cp * to(i,k) + hvap * qo(i,k)
1145 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
1146 & cp * to(i,k) + hvap * qeso(i,k)
1147 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
1148 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
1149 endif
1150 enddo
1151 enddo
1152 !c
1153 !c look for the level of free convection as cloud base
1154 !c
1155 do i=1,im
1156 flg(i) = .true.
1157 kbcon(i) = kmax(i)
1158 enddo
1159 do k = 1, km1
1160 do i=1,im
1161 if (flg(i) .and. k <= kbmax(i)) then
1162 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
1163 kbcon(i) = k
1164 flg(i) = .false.
1165 endif
1166 endif
1167 enddo
1168 enddo
1169 !c
1170 do i=1,im
1171 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
1172 enddo
1173 !!
1174 totflg = .true.
1175 do i=1,im
1176 totflg = totflg .and. (.not. cnvflg(i))
1177 enddo
1178 if(totflg) return
1179 !!
1180 do i=1,im
1181 if(cnvflg(i)) then
1182 ! pdot(i) = 10.* dot(i,kbcon(i))
1183 pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
1184 endif
1185 enddo
1186 !c
1187 !c turn off convection if pressure depth between parcel source level
1188 !c and cloud base is larger than a critical value, cinpcr
1189 !c
1190 do i=1,im
1191 if(cnvflg(i)) then
1192 if(islimsk(i) == 1) then
1193 w1 = w1l
1194 w2 = w2l
1195 w3 = w3l
1196 w4 = w4l
1197 else
1198 w1 = w1s
1199 w2 = w2s
1200 w3 = w3s
1201 w4 = w4s
1202 endif
1203 if(pdot(i) <= w4) then
1204 tem = (pdot(i) - w4) / (w3 - w4)
1205 elseif(pdot(i) >= -w4) then
1206 tem = - (pdot(i) + w4) / (w4 - w3)
1207 else
1208 tem = 0.
1209 endif
1210 val1 = -1.
1211 tem = max(tem,val1)
1212 val2 = 1.
1213 tem = min(tem,val2)
1214 ptem = 1. - tem
1215 ptem1= .5*(cinpcrmx-cinpcrmn)
1216 cinpcr = cinpcrmx - ptem * ptem1
1217 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
1218 #if HWRF==1
1219 ! randomly perturb the convection trigger
1220 !zz if( pert_sas_local .and. ens_random_seed_local .gt. 0 ) then
1221 if( pert_sas_local ) then
1222 !zz print*,"ens_random_seed==",ens_random_seed,ens_random_seed_local
1223 ens_random_seed_local=ran1(-ens_random_seed_local)*1000
1224 rr=2.0*ens_sasamp_local*ran1(-ens_random_seed_local)-ens_sasamp_local
1225 !zz print*, "zhang inde desas=a", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr
1226 cinpcr=cinpcr+rr
1227 !zz print*, "zhang inde desas=b", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr
1228 endif
1229 #endif
1230 if(tem1 > cinpcr) then
1231 cnvflg(i) = .false.
1232 endif
1233 endif
1234 enddo
1235 !!
1236 totflg = .true.
1237 do i=1,im
1238 totflg = totflg .and. (.not. cnvflg(i))
1239 enddo
1240 if(totflg) return
1241 !!
1242 !c
1243 !c assume that updraft entrainment rate above cloud base is
1244 !c same as that at cloud base
1245 !c
1246 do i=1,im
1247 if(cnvflg(i)) then
1248 xlamx(i) = xlamue(i,kbcon(i))
1249 endif
1250 enddo
1251 do k = 2, km1
1252 do i=1,im
1253 if(cnvflg(i).and. &
1254 & (k > kbcon(i) .and. k < kmax(i))) then
1255 xlamue(i,k) = xlamx(i)
1256 endif
1257 enddo
1258 enddo
1259 !c
1260 !c specify a background (turbulent) detrainment rate for the updrafts
1261 !c
1262 do k = 1, km1
1263 do i=1,im
1264 if(cnvflg(i) .and. k < kmax(i)) then
1265 xlamud(i,k) = xlamx(i)
1266 ! xlamud(i,k) = crtlamd
1267 endif
1268 enddo
1269 enddo
1270 !c
1271 !c functions rapidly decreasing with height, mimicking a cloud ensemble
1272 !c (Bechtold et al., 2008)
1273 !c
1274 do k = 2, km1
1275 do i=1,im
1276 if(cnvflg(i).and. &
1277 & (k > kbcon(i) .and. k < kmax(i))) then
1278 tem = qeso(i,k)/qeso(i,kbcon(i))
1279 fent1(i,k) = tem**2
1280 fent2(i,k) = tem**3
1281 endif
1282 enddo
1283 enddo
1284 !c
1285 !c final entrainment and detrainment rates as the sum of turbulent part and
1286 !c organized entrainment depending on the environmental relative humidity
1287 !c (Bechtold et al., 2008)
1288 !c
1289 do k = 2, km1
1290 do i=1,im
1291 if(cnvflg(i) .and. &
1292 & (k > kbcon(i) .and. k < kmax(i))) then
1293 tem = cxlamu * frh(i,k) * fent2(i,k)
1294 xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
1295 ! tem1 = cxlamd * frh(i,k)
1296 ! xlamud(i,k) = xlamud(i,k) + tem1
1297 endif
1298 enddo
1299 enddo
1300 !
1301 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1302 !c
1303 !c determine updraft mass flux for the subcloud layers
1304 !c
1305 do k = km1, 1, -1
1306 do i = 1, im
1307 if (cnvflg(i)) then
1308 if(k < kbcon(i) .and. k >= kb(i)) then
1309 dz = zi(i,k+1) - zi(i,k)
1310 tem = 0.5*(xlamud(i,k)+xlamud(i,k+1))
1311 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-tem
1312 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
1313 endif
1314 endif
1315 enddo
1316 enddo
1317 !c
1318 !c compute mass flux above cloud base
1319 !c
1320 do i = 1, im
1321 flg(i) = cnvflg(i)
1322 enddo
1323 do k = 2, km1
1324 do i = 1, im
1325 if(flg(i))then
1326 if(k > kbcon(i) .and. k < kmax(i)) then
1327 dz = zi(i,k) - zi(i,k-1)
1328 tem = 0.5*(xlamud(i,k)+xlamud(i,k-1))
1329 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-tem
1330 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
1331 if(eta(i,k) <= 0.) then
1332 kmax(i) = k
1333 ktconn(i) = k
1334 flg(i) = .false.
1335 endif
1336 endif
1337 endif
1338 enddo
1339 enddo
1340 !c
1341 !c compute updraft cloud properties
1342 !c
1343 do i = 1, im
1344 if(cnvflg(i)) then
1345 indx = kb(i)
1346 hcko(i,indx) = heo(i,indx)
1347 ucko(i,indx) = uo(i,indx)
1348 vcko(i,indx) = vo(i,indx)
1349 pwavo(i) = 0.
1350 endif
1351 enddo
1352 !c
1353 !c cloud property is modified by the entrainment process
1354 !c
1355 ! cm is an enhancement factor in entrainment rates for momentum
1356 !
1357 do k = 2, km1
1358 do i = 1, im
1359 if (cnvflg(i)) then
1360 if(k > kb(i) .and. k < kmax(i)) then
1361 dz = zi(i,k) - zi(i,k-1)
1362 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1363 tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
1364 factor = 1. + tem - tem1
1365 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
1366 & (heo(i,k)+heo(i,k-1)))/factor
1367 dbyo(i,k) = hcko(i,k) - heso(i,k)
1368 !
1369 tem = 0.5 * cm * tem
1370 factor = 1. + tem
1371 ptem = tem + pgcon
1372 ptem1= tem - pgcon
1373 ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*uo(i,k) &
1374 & +ptem1*uo(i,k-1))/factor
1375 vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*vo(i,k) &
1376 & +ptem1*vo(i,k-1))/factor
1377 endif
1378 endif
1379 enddo
1380 enddo
1381 !c
1382 !c taking account into convection inhibition due to existence of
1383 !c dry layers below cloud base
1384 !c
1385 do i=1,im
1386 flg(i) = cnvflg(i)
1387 kbcon1(i) = kmax(i)
1388 enddo
1389 do k = 2, km1
1390 do i=1,im
1391 if (flg(i) .and. k < kmax(i)) then
1392 if(k >= kbcon(i) .and. dbyo(i,k) > 0.) then
1393 kbcon1(i) = k
1394 flg(i) = .false.
1395 endif
1396 endif
1397 enddo
1398 enddo
1399 do i=1,im
1400 if(cnvflg(i)) then
1401 if(kbcon1(i) == kmax(i)) cnvflg(i) = .false.
1402 endif
1403 enddo
1404 do i=1,im
1405 if(cnvflg(i)) then
1406 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
1407 if(tem > dthk) then
1408 cnvflg(i) = .false.
1409 endif
1410 endif
1411 enddo
1412 !!
1413 totflg = .true.
1414 do i = 1, im
1415 totflg = totflg .and. (.not. cnvflg(i))
1416 enddo
1417 if(totflg) return
1418 !!
1419 !c
1420 !c calculate convective inhibition
1421 !c
1422 do k = 2, km1
1423 do i = 1, im
1424 if (cnvflg(i)) then
1425 if(k > kb(i) .and. k < kbcon1(i)) then
1426 dz1 = zo(i,k+1) - zo(i,k)
1427 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1428 rfact = 1. + delta * cp * gamma &
1429 & * to(i,k) / hvap
1430 cina(i) = cina(i) + &
1431 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
1432 & dz1 * (g / (cp * to(i,k))) &
1433 & * dbyo(i,k) / (1. + gamma) &
1434 & * rfact
1435 val = 0.
1436 cina(i) = cina(i) + &
1437 ! & dz1 * eta(i,k) * g * delta *
1438 & dz1 * g * delta * &
1439 & max(val,(qeso(i,k) - qo(i,k)))
1440 endif
1441 endif
1442 enddo
1443 enddo
1444 do i = 1, im
1445 if(cnvflg(i)) then
1446 !
1447 ! if(islimsk(i) == 1) then
1448 ! w1 = w1l
1449 ! w2 = w2l
1450 ! w3 = w3l
1451 ! w4 = w4l
1452 ! else
1453 ! w1 = w1s
1454 ! w2 = w2s
1455 ! w3 = w3s
1456 ! w4 = w4s
1457 ! endif
1458 ! if(pdot(i) <= w4) then
1459 ! tem = (pdot(i) - w4) / (w3 - w4)
1460 ! elseif(pdot(i) >= -w4) then
1461 ! tem = - (pdot(i) + w4) / (w4 - w3)
1462 ! else
1463 ! tem = 0.
1464 ! endif
1465 !
1466 ! val1 = -1.
1467 ! tem = max(tem,val1)
1468 ! val2 = 1.
1469 ! tem = min(tem,val2)
1470 ! tem = 1. - tem
1471 ! tem1= .5*(cinacrmx-cinacrmn)
1472 ! cinacr = cinacrmx - tem * tem1
1473 !
1474 cinacr = cinacrmx
1475 if(cina(i) < cinacr) cnvflg(i) = .false.
1476 endif
1477 enddo
1478 !!
1479 totflg = .true.
1480 do i=1,im
1481 totflg = totflg .and. (.not. cnvflg(i))
1482 enddo
1483 if(totflg) return
1484 !!
1485 !c
1486 !c determine first guess cloud top as the level of zero buoyancy
1487 !c
1488 do i = 1, im
1489 flg(i) = cnvflg(i)
1490 ktcon(i) = 1
1491 enddo
1492 do k = 2, km1
1493 do i = 1, im
1494 if (flg(i) .and. k < kmax(i)) then
1495 if(k > kbcon1(i) .and. dbyo(i,k) < 0.) then
1496 ktcon(i) = k
1497 flg(i) = .false.
1498 endif
1499 endif
1500 enddo
1501 enddo
1502 !c
1503 do i = 1, im
1504 if(cnvflg(i)) then
1505 if(ktcon(i) == 1 .and. ktconn(i) > 1) then
1506 ktcon(i) = ktconn(i)
1507 endif
1508 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
1509 if(tem < cthk) cnvflg(i) = .false.
1510 endif
1511 enddo
1512 !!
1513 totflg = .true.
1514 do i = 1, im
1515 totflg = totflg .and. (.not. cnvflg(i))
1516 enddo
1517 if(totflg) return
1518 !!
1519 !c
1520 !c search for downdraft originating level above theta-e minimum
1521 !c
1522 do i = 1, im
1523 if(cnvflg(i)) then
1524 hmin(i) = heo(i,kbcon1(i))
1525 lmin(i) = kbmax(i)
1526 jmin(i) = kbmax(i)
1527 endif
1528 enddo
1529 do k = 2, km1
1530 do i = 1, im
1531 if (cnvflg(i) .and. k <= kbmax(i)) then
1532 if(k > kbcon1(i) .and. heo(i,k) < hmin(i)) then
1533 lmin(i) = k + 1
1534 hmin(i) = heo(i,k)
1535 endif
1536 endif
1537 enddo
1538 enddo
1539 !c
1540 !c make sure that jmin(i) is within the cloud
1541 !c
1542 do i = 1, im
1543 if(cnvflg(i)) then
1544 jmin(i) = min(lmin(i),ktcon(i)-1)
1545 jmin(i) = max(jmin(i),kbcon1(i)+1)
1546 if(jmin(i) >= ktcon(i)) cnvflg(i) = .false.
1547 endif
1548 enddo
1549 !c
1550 !c specify upper limit of mass flux at cloud base
1551 !c
1552 do i = 1, im
1553 if(cnvflg(i)) then
1554 ! xmbmax(i) = .1
1555 !
1556 k = kbcon(i)
1557 dp = 1000. * del(i,k)
1558 xmbmax(i) = dp / (g * dt2)
1559 ! mbdt(i) = 0.1 * dp / g
1560 !
1561 ! tem = dp / (g * dt2)
1562 ! xmbmax(i) = min(tem, xmbmax(i))
1563 endif
1564 enddo
1565 !c
1566 !c compute cloud moisture property and precipitation
1567 !c
1568 do i = 1, im
1569 if (cnvflg(i)) then
1570 ! aa1(i) = 0.
1571 qcko(i,kb(i)) = qo(i,kb(i))
1572 qrcko(i,kb(i)) = qo(i,kb(i))
1573 ! rhbar(i) = 0.
1574 endif
1575 enddo
1576 do k = 2, km1
1577 do i = 1, im
1578 if (cnvflg(i)) then
1579 if(k > kb(i) .and. k < ktcon(i)) then
1580 dz = zi(i,k) - zi(i,k-1)
1581 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1582 qrch = qeso(i,k) &
1583 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1584 !cj
1585 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1586 tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
1587 factor = 1. + tem - tem1
1588 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
1589 & (qo(i,k)+qo(i,k-1)))/factor
1590 qrcko(i,k) = qcko(i,k)
1591 !cj
1592 dq = eta(i,k) * (qcko(i,k) - qrch)
1593 !c
1594 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
1595 !c
1596 !c check if there is excess moisture to release latent heat
1597 !c
1598 if(k >= kbcon(i) .and. dq > 0.) then
1599 etah = .5 * (eta(i,k) + eta(i,k-1))
1600 if(ncloud > 0 .and. k > jmin(i)) then
1601 dp = 1000. * del(i,k)
1602 ptem = c0t(i,k) + c1
1603 qlk = dq / (eta(i,k) + etah * ptem * dz)
1604 dellal(i,k) = etah * c1 * dz * qlk * g / dp
1605 else
1606 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1607 endif
1608 ! aa1(i) = aa1(i) - dz * g * qlk * etah
1609 ! aa1(i) = aa1(i) - dz * g * qlk
1610 buo(i,k) = buo(i,k) - g * qlk
1611 qcko(i,k) = qlk + qrch
1612 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1613 pwavo(i) = pwavo(i) + pwo(i,k)
1614 ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
1615 cnvwt(i,k) = etah * qlk * g / dp
1616 endif
1617 !
1618 ! compute buoyancy and drag for updraft velocity
1619 !
1620 if(k >= kbcon(i)) then
1621 rfact = 1. + delta * cp * gamma &
1622 & * to(i,k) / hvap
1623 buo(i,k) = buo(i,k) + (g / (cp * to(i,k))) &
1624 & * dbyo(i,k) / (1. + gamma) &
1625 & * rfact
1626 val = 0.
1627 buo(i,k) = buo(i,k) + g * delta * &
1628 & max(val,(qeso(i,k) - qo(i,k)))
1629 drag(i,k) = max(xlamue(i,k),xlamud(i,k))
1630 endif
1631 !
1632 endif
1633 endif
1634 enddo
1635 enddo
1636 !c
1637 ! do i = 1, im
1638 ! if(cnvflg(i)) then
1639 ! indx = ktcon(i) - kb(i) - 1
1640 ! rhbar(i) = rhbar(i) / float(indx)
1641 ! endif
1642 ! enddo
1643 !c
1644 !c calculate cloud work function
1645 !c
1646 ! do k = 2, km1
1647 ! do i = 1, im
1648 ! if (cnvflg(i)) then
1649 ! if(k >= kbcon(i) .and. k < ktcon(i)) then
1650 ! dz1 = zo(i,k+1) - zo(i,k)
1651 ! gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1652 ! rfact = 1. + delta * cp * gamma
1653 ! & * to(i,k) / hvap
1654 ! aa1(i) = aa1(i) +
1655 !! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
1656 ! & dz1 * (g / (cp * to(i,k)))
1657 ! & * dbyo(i,k) / (1. + gamma)
1658 ! & * rfact
1659 ! val = 0.
1660 ! aa1(i) = aa1(i) +
1661 !! & dz1 * eta(i,k) * g * delta *
1662 ! & dz1 * g * delta *
1663 ! & max(val,(qeso(i,k) - qo(i,k)))
1664 ! endif
1665 ! endif
1666 ! enddo
1667 ! enddo
1668 !
1669 ! calculate cloud work function
1670 !
1671 do i = 1, im
1672 if (cnvflg(i)) then
1673 aa1(i) = 0.
1674 endif
1675 enddo
1676 do k = 2, km1
1677 do i = 1, im
1678 if (cnvflg(i)) then
1679 if(k >= kbcon(i) .and. k < ktcon(i)) then
1680 dz1 = zo(i,k+1) - zo(i,k)
1681 ! aa1(i) = aa1(i) + buo(i,k) * dz1 * eta(i,k)
1682 aa1(i) = aa1(i) + buo(i,k) * dz1
1683 endif
1684 endif
1685 enddo
1686 enddo
1687 !
1688 do i = 1, im
1689 if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
1690 enddo
1691 !!
1692 totflg = .true.
1693 do i=1,im
1694 totflg = totflg .and. (.not. cnvflg(i))
1695 enddo
1696 if(totflg) return
1697 !!
1698 !c
1699 !c estimate the onvective overshooting as the level
1700 !c where the [aafac * cloud work function] becomes zero,
1701 !c which is the final cloud top
1702 !c
1703 do i = 1, im
1704 if (cnvflg(i)) then
1705 aa2(i) = aafac * aa1(i)
1706 endif
1707 enddo
1708 !c
1709 do i = 1, im
1710 flg(i) = cnvflg(i)
1711 ktcon1(i) = kmax(i)
1712 enddo
1713 do k = 2, km1
1714 do i = 1, im
1715 if (flg(i)) then
1716 if(k >= ktcon(i) .and. k < kmax(i)) then
1717 dz1 = zo(i,k+1) - zo(i,k)
1718 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1719 rfact = 1. + delta * cp * gamma &
1720 & * to(i,k) / hvap
1721 aa2(i) = aa2(i) + &
1722 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
1723 & dz1 * (g / (cp * to(i,k))) &
1724 & * dbyo(i,k) / (1. + gamma) &
1725 & * rfact
1726 ! val = 0.
1727 ! aa2(i) = aa2(i) +
1728 !! & dz1 * eta(i,k) * g * delta *
1729 ! & dz1 * g * delta *
1730 ! & max(val,(qeso(i,k) - qo(i,k)))
1731 if(aa2(i) < 0.) then
1732 ktcon1(i) = k
1733 flg(i) = .false.
1734 endif
1735 endif
1736 endif
1737 enddo
1738 enddo
1739 !c
1740 !c compute cloud moisture property, detraining cloud water
1741 !c and precipitation in overshooting layers
1742 !c
1743 do k = 2, km1
1744 do i = 1, im
1745 if (cnvflg(i)) then
1746 if(k >= ktcon(i) .and. k < ktcon1(i)) then
1747 dz = zi(i,k) - zi(i,k-1)
1748 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1749 qrch = qeso(i,k) &
1750 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1751 !cj
1752 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1753 tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
1754 factor = 1. + tem - tem1
1755 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
1756 & (qo(i,k)+qo(i,k-1)))/factor
1757 qrcko(i,k) = qcko(i,k)
1758 !cj
1759 dq = eta(i,k) * (qcko(i,k) - qrch)
1760 !c
1761 !c check if there is excess moisture to release latent heat
1762 !c
1763 if(dq > 0.) then
1764 etah = .5 * (eta(i,k) + eta(i,k-1))
1765 if(ncloud > 0) then
1766 dp = 1000. * del(i,k)
1767 ptem = c0t(i,k) + c1
1768 qlk = dq / (eta(i,k) + etah * ptem * dz)
1769 dellal(i,k) = etah * c1 * dz * qlk * g / dp
1770 else
1771 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1772 endif
1773 qcko(i,k) = qlk + qrch
1774 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1775 pwavo(i) = pwavo(i) + pwo(i,k)
1776 ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * g / dp
1777 cnvwt(i,k) = etah * qlk * g / dp
1778 endif
1779 endif
1780 endif
1781 enddo
1782 enddo
1783 !
1784 ! compute updraft velocity square(wu2)
1785 !
1786 ! bb1 = 2. * (1.+bet1*cd1)
1787 ! bb2 = 2. / (f1*(1.+gam1))
1788 !
1789 ! bb1 = 3.9
1790 ! bb2 = 0.67
1791 !
1792 ! bb1 = 2.0
1793 ! bb2 = 4.0
1794 !
1795 bb1 = 4.0
1796 bb2 = 0.8
1797 !
1798 do i = 1, im
1799 if (cnvflg(i)) then
1800 k = kbcon1(i)
1801 tem = po(i,k) / (rd * to(i,k))
1802 wucb = -0.01 * dot(i,k) / (tem * g)
1803 if(wucb > 0.) then
1804 wu2(i,k) = wucb * wucb
1805 else
1806 wu2(i,k) = 0.
1807 endif
1808 endif
1809 enddo
1810 do k = 2, km1
1811 do i = 1, im
1812 if (cnvflg(i)) then
1813 if(k > kbcon1(i) .and. k < ktcon(i)) then
1814 dz = zi(i,k) - zi(i,k-1)
1815 tem = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz
1816 tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz
1817 ptem = (1. - tem) * wu2(i,k-1)
1818 ptem1 = 1. + tem
1819 wu2(i,k) = (ptem + tem1) / ptem1
1820 wu2(i,k) = max(wu2(i,k), 0._kind_phys)
1821 endif
1822 endif
1823 enddo
1824 enddo
1825 !
1826 ! compute updraft velocity average over the whole cumulus
1827 !
1828 do i = 1, im
1829 wc(i) = 0.
1830 sumx(i) = 0.
1831 enddo
1832 do k = 2, km1
1833 do i = 1, im
1834 if (cnvflg(i)) then
1835 if(k > kbcon1(i) .and. k < ktcon(i)) then
1836 dz = zi(i,k) - zi(i,k-1)
1837 tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1838 wc(i) = wc(i) + tem * dz
1839 sumx(i) = sumx(i) + dz
1840 endif
1841 endif
1842 enddo
1843 enddo
1844 do i = 1, im
1845 if(cnvflg(i)) then
1846 if(sumx(i) == 0.) then
1847 cnvflg(i)=.false.
1848 else
1849 wc(i) = wc(i) / sumx(i)
1850 endif
1851 val = 1.e-4
1852 if (wc(i) < val) cnvflg(i)=.false.
1853 endif
1854 enddo
1855 !c
1856 !c exchange ktcon with ktcon1
1857 !c
1858 do i = 1, im
1859 if(cnvflg(i)) then
1860 kk = ktcon(i)
1861 ktcon(i) = ktcon1(i)
1862 ktcon1(i) = kk
1863 endif
1864 enddo
1865 !c
1866 !c this section is ready for cloud water
1867 !c
1868 if(ncloud > 0) then
1869 !c
1870 !c compute liquid and vapor separation at cloud top
1871 !c
1872 do i = 1, im
1873 if(cnvflg(i)) then
1874 k = ktcon(i) - 1
1875 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1876 qrch = qeso(i,k) &
1877 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1878 dq = qcko(i,k) - qrch
1879 !c
1880 !c check if there is excess moisture to release latent heat
1881 !c
1882 if(dq > 0.) then
1883 qlko_ktcon(i) = dq
1884 qcko(i,k) = qrch
1885 endif
1886 endif
1887 enddo
1888 endif
1889 !c
1890 !ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then
1891 !ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
1892 !ccccc endif
1893 !c
1894 !c------- downdraft calculations
1895 !c
1896 !c--- compute precipitation efficiency in terms of windshear
1897 !c
1898 do i = 1, im
1899 if(cnvflg(i)) then
1900 vshear(i) = 0.
1901 endif
1902 enddo
1903 do k = 2, km
1904 do i = 1, im
1905 if (cnvflg(i)) then
1906 if(k > kb(i) .and. k <= ktcon(i)) then
1907 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
1908 & + (vo(i,k)-vo(i,k-1)) ** 2)
1909 vshear(i) = vshear(i) + shear
1910 endif
1911 endif
1912 enddo
1913 enddo
1914 do i = 1, im
1915 if(cnvflg(i)) then
1916 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
1917 e1=1.591-.639*vshear(i) &
1918 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1919 edt(i)=1.-e1
1920 val = .9
1921 edt(i) = min(edt(i),val)
1922 val = .0
1923 edt(i) = max(edt(i),val)
1924 edto(i)=edt(i)
1925 edtx(i)=edt(i)
1926 endif
1927 enddo
1928 !c
1929 !c determine detrainment rate between 1 and kbcon
1930 !c
1931 do i = 1, im
1932 if(cnvflg(i)) then
1933 sumx(i) = 0.
1934 endif
1935 enddo
1936 do k = 1, km1
1937 do i = 1, im
1938 if(cnvflg(i)) then
1939 if(k >= 1 .and. k < kbcon(i)) then
1940 dz = zi(i,k+1) - zi(i,k)
1941 sumx(i) = sumx(i) + dz
1942 endif
1943 endif
1944 enddo
1945 enddo
1946 do i = 1, im
1947 beta = betas
1948 if(islimsk(i) == 1) beta = betal
1949 if(cnvflg(i)) then
1950 dz = (sumx(i)+zi(i,1))/float(kbcon(i))
1951 tem = 1./float(kbcon(i))
1952 xlamd(i) = (1.-beta**tem)/dz
1953 endif
1954 enddo
1955 !c
1956 !c determine downdraft mass flux
1957 !c
1958 do k = km1, 1, -1
1959 do i = 1, im
1960 if (cnvflg(i) .and. k <= kmax(i)-1) then
1961 if(k < jmin(i) .and. k >= kbcon(i)) then
1962 dz = zi(i,k+1) - zi(i,k)
1963 ptem = xlamdd - xlamde
1964 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1965 else if(k < kbcon(i)) then
1966 dz = zi(i,k+1) - zi(i,k)
1967 ptem = xlamd(i) + xlamdd - xlamde
1968 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
1969 endif
1970 endif
1971 enddo
1972 enddo
1973 !c
1974 !c--- downdraft moisture properties
1975 !c
1976 do i = 1, im
1977 if(cnvflg(i)) then
1978 jmn = jmin(i)
1979 hcdo(i,jmn) = heo(i,jmn)
1980 qcdo(i,jmn) = qo(i,jmn)
1981 qrcdo(i,jmn)= qo(i,jmn)
1982 ucdo(i,jmn) = uo(i,jmn)
1983 vcdo(i,jmn) = vo(i,jmn)
1984 pwevo(i) = 0.
1985 endif
1986 enddo
1987 !cj
1988 do k = km1, 1, -1
1989 do i = 1, im
1990 if (cnvflg(i) .and. k < jmin(i)) then
1991 dz = zi(i,k+1) - zi(i,k)
1992 if(k >= kbcon(i)) then
1993 tem = xlamde * dz
1994 tem1 = 0.5 * xlamdd * dz
1995 else
1996 tem = xlamde * dz
1997 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1998 endif
1999 factor = 1. + tem - tem1
2000 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
2001 & (heo(i,k)+heo(i,k+1)))/factor
2002 dbyo(i,k) = hcdo(i,k) - heso(i,k)
2003 !
2004 tem = 0.5 * cm * tem
2005 factor = 1. + tem
2006 ptem = tem - pgcon
2007 ptem1= tem + pgcon
2008 ucdo(i,k) = ((1.-tem)*ucdo(i,k+1)+ptem*uo(i,k+1) &
2009 & +ptem1*uo(i,k))/factor
2010 vcdo(i,k) = ((1.-tem)*vcdo(i,k+1)+ptem*vo(i,k+1) &
2011 & +ptem1*vo(i,k))/factor
2012 endif
2013 enddo
2014 enddo
2015 !c
2016 do k = km1, 1, -1
2017 do i = 1, im
2018 if (cnvflg(i) .and. k < jmin(i)) then
2019 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2020 qrcdo(i,k) = qeso(i,k)+ &
2021 & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
2022 ! detad = etad(i,k+1) - etad(i,k)
2023 !cj
2024 dz = zi(i,k+1) - zi(i,k)
2025 if(k >= kbcon(i)) then
2026 tem = xlamde * dz
2027 tem1 = 0.5 * xlamdd * dz
2028 else
2029 tem = xlamde * dz
2030 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
2031 endif
2032 factor = 1. + tem - tem1
2033 qcdo(i,k) = ((1.-tem1)*qrcdo(i,k+1)+tem*0.5* &
2034 & (qo(i,k)+qo(i,k+1)))/factor
2035 !cj
2036 ! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
2037 ! & etad(i,k) * qrcdo(i,k)
2038 ! pwdo(i,k) = pwdo(i,k) - detad *
2039 ! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
2040 !cj
2041 pwdo(i,k) = etad(i,k) * (qcdo(i,k) - qrcdo(i,k))
2042 pwevo(i) = pwevo(i) + pwdo(i,k)
2043 endif
2044 enddo
2045 enddo
2046 !c
2047 !c--- final downdraft strength dependent on precip
2048 !c--- efficiency (edt), normalized condensate (pwav), and
2049 !c--- evaporate (pwev)
2050 !c
2051 do i = 1, im
2052 edtmax = edtmaxl
2053 if(islimsk(i) == 0) edtmax = edtmaxs
2054 if(cnvflg(i)) then
2055 if(pwevo(i) < 0.) then
2056 edto(i) = -edto(i) * pwavo(i) / pwevo(i)
2057 edto(i) = min(edto(i),edtmax)
2058 else
2059 edto(i) = 0.
2060 endif
2061 endif
2062 enddo
2063 !c
2064 !c--- downdraft cloudwork functions
2065 !c
2066 do k = km1, 1, -1
2067 do i = 1, im
2068 if (cnvflg(i) .and. k < jmin(i)) then
2069 gamma = el2orc * qeso(i,k) / to(i,k)**2
2070 dhh=hcdo(i,k)
2071 dt=to(i,k)
2072 dg=gamma
2073 dh=heso(i,k)
2074 dz=-1.*(zo(i,k+1)-zo(i,k))
2075 ! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k)
2076 aa1(i)=aa1(i)+edto(i)*dz &
2077 & *(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
2078 & *(1.+delta*cp*dg*dt/hvap)
2079 val=0.
2080 ! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k)
2081 aa1(i)=aa1(i)+edto(i)*dz &
2082 & *g*delta*max(val,(qeso(i,k)-qo(i,k)))
2083 endif
2084 enddo
2085 enddo
2086 do i = 1, im
2087 if(cnvflg(i) .and. aa1(i) <= 0.) then
2088 cnvflg(i) = .false.
2089 endif
2090 enddo
2091 !!
2092 totflg = .true.
2093 do i=1,im
2094 totflg = totflg .and. (.not. cnvflg(i))
2095 enddo
2096 if(totflg) return
2097 !!
2098 !c
2099 !c--- what would the change be, that a cloud with unit mass
2100 !c--- will do to the environment?
2101 !c
2102 do k = 1, km
2103 do i = 1, im
2104 if(cnvflg(i) .and. k <= kmax(i)) then
2105 dellah(i,k) = 0.
2106 dellaq(i,k) = 0.
2107 dellau(i,k) = 0.
2108 dellav(i,k) = 0.
2109 endif
2110 enddo
2111 enddo
2112 do i = 1, im
2113 if(cnvflg(i)) then
2114 dp = 1000. * del(i,1)
2115 dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1) &
2116 & - heo(i,1)) * g / dp
2117 dellaq(i,1) = edto(i) * etad(i,1) * (qrcdo(i,1) &
2118 & - qo(i,1)) * g / dp
2119 dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1) &
2120 & - uo(i,1)) * g / dp
2121 dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1) &
2122 & - vo(i,1)) * g / dp
2123 endif
2124 enddo
2125 !c
2126 !c--- changed due to subsidence and entrainment
2127 !c
2128 do k = 2, km1
2129 do i = 1, im
2130 if (cnvflg(i) .and. k < ktcon(i)) then
2131 aup = 1.
2132 if(k <= kb(i)) aup = 0.
2133 adw = 1.
2134 if(k > jmin(i)) adw = 0.
2135 dp = 1000. * del(i,k)
2136 dz = zi(i,k) - zi(i,k-1)
2137 !c
2138 dv1h = heo(i,k)
2139 dv2h = .5 * (heo(i,k) + heo(i,k-1))
2140 dv3h = heo(i,k-1)
2141 dv1q = qo(i,k)
2142 dv2q = .5 * (qo(i,k) + qo(i,k-1))
2143 dv3q = qo(i,k-1)
2144 !c
2145 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
2146 tem1 = 0.5 * (xlamud(i,k)+xlamud(i,k-1))
2147 !c
2148 if(k <= kbcon(i)) then
2149 ptem = xlamde
2150 ptem1 = xlamd(i)+xlamdd
2151 else
2152 ptem = xlamde
2153 ptem1 = xlamdd
2154 endif
2155 !cj
2156 dellah(i,k) = dellah(i,k) + &
2157 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h &
2158 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h &
2159 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz &
2160 & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
2161 & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz &
2162 & ) *g/dp
2163 !cj
2164 dellaq(i,k) = dellaq(i,k) + &
2165 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q &
2166 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q &
2167 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz &
2168 & + aup*tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz &
2169 & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qcdo(i,k-1))*dz &
2170 & ) *g/dp
2171 !cj
2172 tem1=eta(i,k)*(uo(i,k)-ucko(i,k))
2173 tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1))
2174 ptem1=etad(i,k)*(uo(i,k)-ucdo(i,k))
2175 ptem2=etad(i,k-1)*(uo(i,k-1)-ucdo(i,k-1))
2176 dellau(i,k) = dellau(i,k) + &
2177 & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*g/dp
2178 !cj
2179 tem1=eta(i,k)*(vo(i,k)-vcko(i,k))
2180 tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1))
2181 ptem1=etad(i,k)*(vo(i,k)-vcdo(i,k))
2182 ptem2=etad(i,k-1)*(vo(i,k-1)-vcdo(i,k-1))
2183 dellav(i,k) = dellav(i,k) + &
2184 & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*g/dp
2185 !cj
2186 endif
2187 enddo
2188 enddo
2189 !c
2190 !c------- cloud top
2191 !c
2192 do i = 1, im
2193 if(cnvflg(i)) then
2194 indx = ktcon(i)
2195 dp = 1000. * del(i,indx)
2196 dv1h = heo(i,indx-1)
2197 dellah(i,indx) = eta(i,indx-1) * &
2198 & (hcko(i,indx-1) - dv1h) * g / dp
2199 dv1q = qo(i,indx-1)
2200 dellaq(i,indx) = eta(i,indx-1) * &
2201 & (qcko(i,indx-1) - dv1q) * g / dp
2202 dellau(i,indx) = eta(i,indx-1) * &
2203 & (ucko(i,indx-1) - uo(i,indx-1)) * g / dp
2204 dellav(i,indx) = eta(i,indx-1) * &
2205 & (vcko(i,indx-1) - vo(i,indx-1)) * g / dp
2206 !c
2207 !c cloud water
2208 !c
2209 dellal(i,indx) = eta(i,indx-1) * &
2210 & qlko_ktcon(i) * g / dp
2211 endif
2212 enddo
2213 !c
2214 !c------- final changed variable per unit mass flux
2215 !c
2216 ! if grid size is less than a threshold value (dxcrtas),
2217 ! the quasi-equilibrium assumption of Arakawa-Schubert is not
2218 ! used any longer.
2219 !
2220 do i = 1, im
2221 asqecflg(i) = cnvflg(i)
2222 if(asqecflg(i) .and. gdx(i) < dxcrtas) then
2223 asqecflg(i) = .false.
2224 endif
2225 enddo
2226 !
2227 do k = 1, km
2228 do i = 1, im
2229 if (asqecflg(i) .and. k <= kmax(i)) then
2230 if(k > ktcon(i)) then
2231 qo(i,k) = q1(i,k)
2232 to(i,k) = t1(i,k)
2233 endif
2234 if(k <= ktcon(i)) then
2235 qo(i,k) = dellaq(i,k) * mbdt(i) + q1(i,k)
2236 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
2237 to(i,k) = dellat * mbdt(i) + t1(i,k)
2238 val = 1.e-10
2239 qo(i,k) = max(qo(i,k), val )
2240 endif
2241 endif
2242 enddo
2243 enddo
2244 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2245 !c
2246 !c--- the above changed environment is now used to calulate the
2247 !c--- effect the arbitrary cloud (with unit mass flux)
2248 !c--- would have on the stability,
2249 !c--- which then is used to calculate the real mass flux,
2250 !c--- necessary to keep this change in balance with the large-scale
2251 !c--- destabilization.
2252 !c
2253 !c--- environmental conditions again, first heights
2254 !c
2255 do k = 1, km
2256 do i = 1, im
2257 if(asqecflg(i) .and. k <= kmax(i)) then
2258 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
2259 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
2260 val = 1.e-8
2261 qeso(i,k) = max(qeso(i,k), val )
2262 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
2263 endif
2264 enddo
2265 enddo
2266 !c
2267 !c--- moist static energy
2268 !c
2269 do k = 1, km1
2270 do i = 1, im
2271 if(asqecflg(i) .and. k <= kmax(i)-1) then
2272 dz = .5 * (zo(i,k+1) - zo(i,k))
2273 dp = .5 * (pfld(i,k+1) - pfld(i,k))
2274 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
2275 pprime = pfld(i,k+1) + epsm1 * es
2276 qs = eps * es / pprime
2277 dqsdp = - qs / pprime
2278 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
2279 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
2280 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
2281 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
2282 dq = dqsdt * dt + dqsdp * dp
2283 to(i,k) = to(i,k+1) + dt
2284 qo(i,k) = qo(i,k+1) + dq
2285 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
2286 endif
2287 enddo
2288 enddo
2289 do k = 1, km1
2290 do i = 1, im
2291 if(asqecflg(i) .and. k <= kmax(i)-1) then
2292 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
2293 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
2294 val1 = 1.e-8
2295 qeso(i,k) = max(qeso(i,k), val1)
2296 val2 = 1.e-10
2297 qo(i,k) = max(qo(i,k), val2 )
2298 ! qo(i,k) = min(qo(i,k),qeso(i,k))
2299 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
2300 & cp * to(i,k) + hvap * qo(i,k)
2301 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
2302 & cp * to(i,k) + hvap * qeso(i,k)
2303 endif
2304 enddo
2305 enddo
2306 do i = 1, im
2307 if(asqecflg(i)) then
2308 k = kmax(i)
2309 heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
2310 heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
2311 !c heo(i,k) = min(heo(i,k),heso(i,k))
2312 endif
2313 enddo
2314 !c
2315 !c**************************** static control
2316 !c
2317 !c------- moisture and cloud work functions
2318 !c
2319 do i = 1, im
2320 if(asqecflg(i)) then
2321 xaa0(i) = 0.
2322 xpwav(i) = 0.
2323 endif
2324 enddo
2325 !c
2326 do i = 1, im
2327 if(asqecflg(i)) then
2328 indx = kb(i)
2329 hcko(i,indx) = heo(i,indx)
2330 qcko(i,indx) = qo(i,indx)
2331 endif
2332 enddo
2333 do k = 2, km1
2334 do i = 1, im
2335 if (asqecflg(i)) then
2336 if(k > kb(i) .and. k <= ktcon(i)) then
2337 dz = zi(i,k) - zi(i,k-1)
2338 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2339 tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
2340 factor = 1. + tem - tem1
2341 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
2342 & (heo(i,k)+heo(i,k-1)))/factor
2343 endif
2344 endif
2345 enddo
2346 enddo
2347 do k = 2, km1
2348 do i = 1, im
2349 if (asqecflg(i)) then
2350 if(k > kb(i) .and. k < ktcon(i)) then
2351 dz = zi(i,k) - zi(i,k-1)
2352 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2353 xdby = hcko(i,k) - heso(i,k)
2354 xqrch = qeso(i,k) &
2355 & + gamma * xdby / (hvap * (1. + gamma))
2356 !cj
2357 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2358 tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz
2359 factor = 1. + tem - tem1
2360 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
2361 & (qo(i,k)+qo(i,k-1)))/factor
2362 !cj
2363 dq = eta(i,k) * (qcko(i,k) - xqrch)
2364 !c
2365 if(k >= kbcon(i) .and. dq > 0.) then
2366 etah = .5 * (eta(i,k) + eta(i,k-1))
2367 if(ncloud > 0 .and. k > jmin(i)) then
2368 ptem = c0t(i,k) + c1
2369 qlk = dq / (eta(i,k) + etah * ptem * dz)
2370 else
2371 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
2372 endif
2373 if(k < ktcon1(i)) then
2374 ! xaa0(i) = xaa0(i) - dz * g * qlk * etah
2375 xaa0(i) = xaa0(i) - dz * g * qlk
2376 endif
2377 qcko(i,k) = qlk + xqrch
2378 xpw = etah * c0t(i,k) * dz * qlk
2379 xpwav(i) = xpwav(i) + xpw
2380 endif
2381 endif
2382 if(k >= kbcon(i) .and. k < ktcon1(i)) then
2383 dz1 = zo(i,k+1) - zo(i,k)
2384 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2385 rfact = 1. + delta * cp * gamma &
2386 & * to(i,k) / hvap
2387 xaa0(i) = xaa0(i) &
2388 ! & + dz1 * eta(i,k) * (g / (cp * to(i,k)))
2389 & + dz1 * (g / (cp * to(i,k))) &
2390 & * xdby / (1. + gamma) &
2391 & * rfact
2392 val=0.
2393 xaa0(i) = xaa0(i) + &
2394 ! & dz1 * eta(i,k) * g * delta *
2395 & dz1 * g * delta * &
2396 & max(val,(qeso(i,k) - qo(i,k)))
2397 endif
2398 endif
2399 enddo
2400 enddo
2401 !c
2402 !c------- downdraft calculations
2403 !c
2404 !c--- downdraft moisture properties
2405 !c
2406 do i = 1, im
2407 if(asqecflg(i)) then
2408 jmn = jmin(i)
2409 hcdo(i,jmn) = heo(i,jmn)
2410 qcdo(i,jmn) = qo(i,jmn)
2411 qrcd(i,jmn) = qo(i,jmn)
2412 xpwev(i) = 0.
2413 endif
2414 enddo
2415 !cj
2416 do k = km1, 1, -1
2417 do i = 1, im
2418 if (asqecflg(i) .and. k < jmin(i)) then
2419 dz = zi(i,k+1) - zi(i,k)
2420 if(k >= kbcon(i)) then
2421 tem = xlamde * dz
2422 tem1 = 0.5 * xlamdd * dz
2423 else
2424 tem = xlamde * dz
2425 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
2426 endif
2427 factor = 1. + tem - tem1
2428 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
2429 & (heo(i,k)+heo(i,k+1)))/factor
2430 endif
2431 enddo
2432 enddo
2433 !cj
2434 do k = km1, 1, -1
2435 do i = 1, im
2436 if (asqecflg(i) .and. k < jmin(i)) then
2437 dq = qeso(i,k)
2438 dt = to(i,k)
2439 gamma = el2orc * dq / dt**2
2440 dh = hcdo(i,k) - heso(i,k)
2441 qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
2442 ! detad = etad(i,k+1) - etad(i,k)
2443 !cj
2444 dz = zi(i,k+1) - zi(i,k)
2445 if(k >= kbcon(i)) then
2446 tem = xlamde * dz
2447 tem1 = 0.5 * xlamdd * dz
2448 else
2449 tem = xlamde * dz
2450 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
2451 endif
2452 factor = 1. + tem - tem1
2453 qcdo(i,k) = ((1.-tem1)*qrcd(i,k+1)+tem*0.5* &
2454 & (qo(i,k)+qo(i,k+1)))/factor
2455 !cj
2456 ! xpwd = etad(i,k+1) * qcdo(i,k+1) -
2457 ! & etad(i,k) * qrcd(i,k)
2458 ! xpwd = xpwd - detad *
2459 ! & .5 * (qrcd(i,k) + qrcd(i,k+1))
2460 !cj
2461 xpwd = etad(i,k) * (qcdo(i,k) - qrcd(i,k))
2462 xpwev(i) = xpwev(i) + xpwd
2463 endif
2464 enddo
2465 enddo
2466 !c
2467 do i = 1, im
2468 edtmax = edtmaxl
2469 if(islimsk(i) == 0) edtmax = edtmaxs
2470 if(asqecflg(i)) then
2471 if(xpwev(i) >= 0.) then
2472 edtx(i) = 0.
2473 else
2474 edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
2475 edtx(i) = min(edtx(i),edtmax)
2476 endif
2477 endif
2478 enddo
2479 !c
2480 !c
2481 !c--- downdraft cloudwork functions
2482 !c
2483 !c
2484 do k = km1, 1, -1
2485 do i = 1, im
2486 if (asqecflg(i) .and. k < jmin(i)) then
2487 gamma = el2orc * qeso(i,k) / to(i,k)**2
2488 dhh=hcdo(i,k)
2489 dt= to(i,k)
2490 dg= gamma
2491 dh= heso(i,k)
2492 dz=-1.*(zo(i,k+1)-zo(i,k))
2493 ! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k)
2494 xaa0(i)=xaa0(i)+edtx(i)*dz &
2495 & *(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
2496 & *(1.+delta*cp*dg*dt/hvap)
2497 val=0.
2498 ! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k)
2499 xaa0(i)=xaa0(i)+edtx(i)*dz &
2500 & *g*delta*max(val,(qeso(i,k)-qo(i,k)))
2501 endif
2502 enddo
2503 enddo
2504 !c
2505 !c calculate critical cloud work function
2506 !c
2507 ! do i = 1, im
2508 ! if(cnvflg(i)) then
2509 ! if(pfld(i,ktcon(i)) < pcrit(15))then
2510 ! acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i)))
2511 ! & /(975.-pcrit(15))
2512 ! else if(pfld(i,ktcon(i)) > pcrit(1))then
2513 ! acrt(i)=acrit(1)
2514 ! else
2515 ! k = int((850. - pfld(i,ktcon(i)))/50.) + 2
2516 ! k = min(k,15)
2517 ! k = max(k,2)
2518 ! acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*
2519 ! & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
2520 ! endif
2521 ! endif
2522 ! enddo
2523 ! do i = 1, im
2524 ! if(cnvflg(i)) then
2525 ! if(islimsk(i) == 1) then
2526 ! w1 = w1l
2527 ! w2 = w2l
2528 ! w3 = w3l
2529 ! w4 = w4l
2530 ! else
2531 ! w1 = w1s
2532 ! w2 = w2s
2533 ! w3 = w3s
2534 ! w4 = w4s
2535 ! endif
2536 !c
2537 !c modify critical cloud workfunction by cloud base vertical velocity
2538 !c
2539 ! if(pdot(i) <= w4) then
2540 ! acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
2541 ! elseif(pdot(i) >= -w4) then
2542 ! acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
2543 ! else
2544 ! acrtfct(i) = 0.
2545 ! endif
2546 ! val1 = -1.
2547 ! acrtfct(i) = max(acrtfct(i),val1)
2548 ! val2 = 1.
2549 ! acrtfct(i) = min(acrtfct(i),val2)
2550 ! acrtfct(i) = 1. - acrtfct(i)
2551 !c
2552 !c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
2553 !c
2554 !c if(rhbar(i) >= .8) then
2555 !c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
2556 !c endif
2557 !c
2558 !c modify adjustment time scale by cloud base vertical velocity
2559 !c
2560 ! dtconv(i) = dt2 + max((1800. - dt2),0.) *
2561 ! & (pdot(i) - w2) / (w1 - w2)
2562 !c dtconv(i) = max(dtconv(i), dt2)
2563 !c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
2564 !
2565 ! dtconv(i) = max(dtconv(i),dtmin)
2566 ! dtconv(i) = min(dtconv(i),dtmax)
2567 !c
2568 ! endif
2569 ! enddo
2570 !
2571 ! compute convective turn-over time
2572 !
2573 do i= 1, im
2574 if(cnvflg(i)) then
2575 tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
2576 dtconv(i) = tem / wc(i)
2577 dtconv(i) = max(dtconv(i),dtmin)
2578 dtconv(i) = min(dtconv(i),dtmax)
2579 endif
2580 enddo
2581 !
2582 ! compute advective time scale using a mean cloud layer wind speed
2583 !
2584 do i= 1, im
2585 if(cnvflg(i)) then
2586 sumx(i) = 0.
2587 umean(i) = 0.
2588 endif
2589 enddo
2590 do k = 2, km1
2591 do i = 1, im
2592 if(cnvflg(i)) then
2593 if(k >= kbcon1(i) .and. k < ktcon1(i)) then
2594 dz = zi(i,k) - zi(i,k-1)
2595 tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
2596 umean(i) = umean(i) + tem * dz
2597 sumx(i) = sumx(i) + dz
2598 endif
2599 endif
2600 enddo
2601 enddo
2602 do i= 1, im
2603 if(cnvflg(i)) then
2604 umean(i) = umean(i) / sumx(i)
2605 umean(i) = max(umean(i), 1._kind_phys)
2606 tauadv(i) = gdx(i) / umean(i)
2607 endif
2608 enddo
2609 !c
2610 !c compute cloud base mass flux as a function of the mean
2611 !c updraft velcoity for the grid sizes where
2612 !c the quasi-equilibrium assumption of Arakawa-Schubert is not
2613 !c valid any longer.
2614 !c
2615 do i= 1, im
2616 if(cnvflg(i) .and. .not.asqecflg(i)) then
2617 k = kbcon(i)
2618 rho = po(i,k)*100. / (rd*to(i,k))
2619 tfac = tauadv(i) / dtconv(i)
2620 tfac = min(tfac, 1._kind_phys)
2621 xmb(i) = tfac*betaw*rho*wc(i)
2622 endif
2623 enddo
2624 !c
2625 !c compute cloud base mass flux using
2626 !c the quasi-equilibrium assumption of Arakawa-Schubert
2627 !c
2628 do i= 1, im
2629 if(asqecflg(i)) then
2630 ! fld(i)=(aa1(i)-acrt(i)*acrtfct(i))/dtconv(i)
2631 fld(i)=aa1(i)/dtconv(i)
2632 if(fld(i) <= 0.) then
2633 asqecflg(i) = .false.
2634 cnvflg(i) = .false.
2635 endif
2636 endif
2637 if(asqecflg(i)) then
2638 !c xaa0(i) = max(xaa0(i),0.)
2639 xk(i) = (xaa0(i) - aa1(i)) / mbdt(i)
2640 if(xk(i) >= 0.) then
2641 asqecflg(i) = .false.
2642 cnvflg(i) = .false.
2643 endif
2644 endif
2645 !c
2646 !c--- kernel, cloud base mass flux
2647 !c
2648 if(asqecflg(i)) then
2649 tfac = tauadv(i) / dtconv(i)
2650 tfac = min(tfac, 1._kind_phys)
2651 xmb(i) = -tfac * fld(i) / xk(i)
2652 ! xmb(i) = min(xmb(i),xmbmax(i))
2653 endif
2654 enddo
2655 !!
2656 totflg = .true.
2657 do i=1,im
2658 totflg = totflg .and. (.not. cnvflg(i))
2659 enddo
2660 if(totflg) return
2661 !!
2662 !
2663 !--- modified Grell & Freitas' (2014) updraft fraction which uses
2664 ! actual entrainment rate at cloud base
2665 !
2666 do i = 1, im
2667 if(cnvflg(i)) then
2668 tem = min(max(xlamx(i), 7.e-5_kind_phys), 3.e-4_kind_phys)
2669 tem = 0.2 / tem
2670 tem1 = 3.14 * tem * tem
2671 sigmagfm(i) = tem1 / garea(i)
2672 sigmagfm(i) = max(sigmagfm(i), 0.001_kind_phys)
2673 sigmagfm(i) = min(sigmagfm(i), 0.999_kind_phys)
2674 endif
2675 enddo
2676 !
2677 !--- compute scale-aware function based on Arakawa & Wu (2013)
2678 !
2679 do i = 1, im
2680 if(cnvflg(i)) then
2681 if (gdx(i) < dxcrtuf) then
2682 scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
2683 scaldfunc(i) = max(min(scaldfunc(i), 1.0_kind_phys), 0._kind_phys)
2684 sigmuout(i)=sigmagfm(i)
2685 else
2686 scaldfunc(i) = 1.0
2687 endif
2688 xmb(i) = xmb(i) * scaldfunc(i)
2689 xmb(i) = min(xmb(i),xmbmax(i))
2690 endif
2691 enddo
2692 !c
2693 !c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
2694 !c
2695 do k = 1, km
2696 do i = 1, im
2697 if (cnvflg(i) .and. k <= kmax(i)) then
2698 to(i,k) = t1(i,k)
2699 qo(i,k) = q1(i,k)
2700 uo(i,k) = u1(i,k)
2701 vo(i,k) = v1(i,k)
2702 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
2703 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
2704 val = 1.e-8
2705 qeso(i,k) = max(qeso(i,k), val )
2706 endif
2707 enddo
2708 enddo
2709 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2710 !c
2711 !c--- feedback: simply the changes from the cloud with unit mass flux
2712 !c--- multiplied by the mass flux necessary to keep the
2713 !c--- equilibrium with the larger-scale.
2714 !c
2715 do i = 1, im
2716 delhbar(i) = 0.
2717 delqbar(i) = 0.
2718 deltbar(i) = 0.
2719 delubar(i) = 0.
2720 delvbar(i) = 0.
2721 qcond(i) = 0.
2722 enddo
2723 do k = 1, km
2724 do i = 1, im
2725 if (cnvflg(i) .and. k <= kmax(i)) then
2726 if(k <= ktcon(i)) then
2727 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
2728 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
2729 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
2730 ! tem = 1./rcs(i)
2731 ! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
2732 ! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
2733 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
2734 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
2735 dp = 1000. * del(i,k)
2736 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
2737 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
2738 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
2739 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
2740 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
2741 endif
2742 endif
2743 enddo
2744 enddo
2745 do k = 1, km
2746 do i = 1, im
2747 if (cnvflg(i) .and. k <= kmax(i)) then
2748 if(k <= ktcon(i)) then
2749 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
2750 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
2751 val = 1.e-8
2752 qeso(i,k) = max(qeso(i,k), val )
2753 endif
2754 endif
2755 enddo
2756 enddo
2757 !c
2758 do i = 1, im
2759 rntot(i) = 0.
2760 delqev(i) = 0.
2761 delq2(i) = 0.
2762 flg(i) = cnvflg(i)
2763 enddo
2764 do k = km, 1, -1
2765 do i = 1, im
2766 if (cnvflg(i) .and. k <= kmax(i)) then
2767 if(k < ktcon(i)) then
2768 aup = 1.
2769 if(k <= kb(i)) aup = 0.
2770 adw = 1.
2771 if(k >= jmin(i)) adw = 0.
2772 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
2773 rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
2774 endif
2775 endif
2776 enddo
2777 enddo
2778 do k = km, 1, -1
2779 do i = 1, im
2780 if (k <= kmax(i)) then
2781 deltv(i) = 0.
2782 delq(i) = 0.
2783 qevap(i) = 0.
2784 if(cnvflg(i) .and. k < ktcon(i)) then
2785 aup = 1.
2786 if(k <= kb(i)) aup = 0.
2787 adw = 1.
2788 if(k >= jmin(i)) adw = 0.
2789 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
2790 rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
2791 endif
2792 if(flg(i) .and. k < ktcon(i)) then
2793 evef = edt(i) * evfact
2794 if(islimsk(i) == 1) evef=edt(i) * evfactl
2795 ! if(islimsk(i) == 1) evef=.07
2796 !c if(islimsk(i) == 1) evef = 0.
2797 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
2798 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
2799 dp = 1000. * del(i,k)
2800 if(rn(i) > 0. .and. qcond(i) < 0.) then
2801 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
2802 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
2803 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
2804 endif
2805 if(rn(i) > 0. .and. qcond(i) < 0. .and. &
2806 & delq2(i) > rntot(i)) then
2807 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
2808 flg(i) = .false.
2809 endif
2810 if(rn(i) > 0. .and. qevap(i) > 0.) then
2811 q1(i,k) = q1(i,k) + qevap(i)
2812 t1(i,k) = t1(i,k) - elocp * qevap(i)
2813 rn(i) = rn(i) - .001 * qevap(i) * dp / g
2814 deltv(i) = - elocp*qevap(i)/dt2
2815 delq(i) = + qevap(i)/dt2
2816 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
2817 endif
2818 delqbar(i) = delqbar(i) + delq(i)*dp/g
2819 deltbar(i) = deltbar(i) + deltv(i)*dp/g
2820 endif
2821 endif
2822 enddo
2823 enddo
2824 !cj
2825 ! do i = 1, im
2826 ! if(me == 31 .and. cnvflg(i)) then
2827 ! if(cnvflg(i)) then
2828 ! print *, ' deep delhbar, delqbar, deltbar = ',
2829 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
2830 ! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
2831 ! print *, ' precip =', hvap*rn(i)*1000./dt2
2832 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
2833 ! endif
2834 ! enddo
2835 !c
2836 !c precipitation rate converted to actual precip
2837 !c in unit of m instead of kg
2838 !c
2839 do i = 1, im
2840 if(cnvflg(i)) then
2841 !c
2842 !c in the event of upper level rain evaporation and lower level downdraft
2843 !c moistening, rn can become negative, in this case, we back out of the
2844 !c heating and the moistening
2845 !c
2846 if(rn(i) < 0. .and. .not.flg(i)) rn(i) = 0.
2847 if(rn(i) <= 0.) then
2848 rn(i) = 0.
2849 else
2850 ktop(i) = ktcon(i)
2851 kbot(i) = kbcon(i)
2852 kcnv(i) = 1
2853 cldwrk(i) = aa1(i)
2854 endif
2855 endif
2856 enddo
2857 !c
2858 !c convective cloud water
2859 !c
2860 do k = 1, km
2861 do i = 1, im
2862 if (cnvflg(i) .and. rn(i) > 0.) then
2863 if (k >= kbcon(i) .and. k < ktcon(i)) then
2864 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
2865 endif
2866 endif
2867 enddo
2868 enddo
2869 !c
2870 !c convective cloud cover
2871 !c
2872 do k = 1, km
2873 do i = 1, im
2874 if (cnvflg(i) .and. rn(i) > 0.) then
2875 if (k >= kbcon(i) .and. k < ktcon(i)) then
2876 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
2877 cnvc(i,k) = min(cnvc(i,k), 0.6_kind_phys)
2878 cnvc(i,k) = max(cnvc(i,k), 0.0_kind_phys)
2879 endif
2880 endif
2881 enddo
2882 enddo
2884 !c
2885 !c cloud water
2886 !c
2887 if (ncloud > 0) then
2888 !
2889 do k = 1, km
2890 do i = 1, im
2891 if (cnvflg(i) .and. rn(i) > 0.) then
2892 ! if (k > kb(i) .and. k <= ktcon(i)) then
2893 if (k >= kbcon(i) .and. k <= ktcon(i)) then
2894 tem = dellal(i,k) * xmb(i) * dt2
2895 tem1 = max(0.0_kind_phys, min(1.0_kind_phys, (tcr-t1(i,k))*tcrf))
2896 if (ql(i,k,2) > -999.0) then
2897 ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
2898 ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
2899 else
2900 ql(i,k,1) = ql(i,k,1) + tem
2901 endif
2902 endif
2903 endif
2904 enddo
2905 enddo
2906 !
2907 endif
2908 !c
2909 do k = 1, km
2910 do i = 1, im
2911 if(cnvflg(i) .and. rn(i) <= 0.) then
2912 if (k <= kmax(i)) then
2913 t1(i,k) = to(i,k)
2914 q1(i,k) = qo(i,k)
2915 u1(i,k) = uo(i,k)
2916 v1(i,k) = vo(i,k)
2917 endif
2918 endif
2919 enddo
2920 enddo
2921 !
2922 ! hchuang code change
2923 !
2924 do k = 1, km
2925 do i = 1, im
2926 if(cnvflg(i) .and. rn(i) > 0.) then
2927 if(k >= kb(i) .and. k < ktop(i)) then
2928 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
2929 endif
2930 endif
2931 enddo
2932 enddo
2933 do i = 1, im
2934 if(cnvflg(i) .and. rn(i) > 0.) then
2935 k = ktop(i)-1
2936 dt_mf(i,k) = ud_mf(i,k)
2937 endif
2938 enddo
2939 do k = 1, km
2940 do i = 1, im
2941 if(cnvflg(i) .and. rn(i) > 0.) then
2942 if(k >= 1 .and. k <= jmin(i)) then
2943 dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
2944 endif
2945 endif
2946 enddo
2947 enddo
2948 !!
2949 return
2950 end subroutine mfdeepcnv
2951 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2953 subroutine mfshalcnv(im,ix,km,delt,delp,prslp,psp,phil,ql, &
2954 & q1,t1,u1,v1,rn,kbot,ktop,kcnv,islimsk,garea, &
2955 & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc,sigmagfm,scaldfunc, &
2956 & pert_sas, ens_random_seed, ens_sasamp)
2958 !
2959 ! use machine , only : kind_phys
2960 ! use funcphys , only : fpvs
2961 ! use physcons, grav => con_g, cp => con_cp, hvap => con_hvap
2962 USE MODULE_GFS_MACHINE, ONLY : kind_phys
2963 USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
2964 USE MODULE_GFS_PHYSCONS, grav => con_g, cp => con_cp &
2965 & , hvap => con_hvap &
2966 &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
2967 &, rd => con_rd, cvap => con_cvap, cliq => con_cliq &
2968 &, eps => con_eps, epsm1 => con_epsm1
2972 implicit none
2973 !
2974 integer im, ix, km, ncloud, &
2975 & kbot(im), ktop(im), kcnv(im)
2976 ! &, me
2977 real(kind=kind_phys) delt
2978 real(kind=kind_phys) psp(im), delp(ix,km), prslp(ix,km)
2979 real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
2980 & ql(ix,km,2),q1(ix,km), t1(ix,km), &
2981 & u1(ix,km), v1(ix,km), &
2982 ! & u1(ix,km), v1(ix,km), rcs(im),
2983 & rn(im), garea(im), &
2984 & dot(ix,km), phil(ix,km), hpbl(im), &
2985 & cnvw(ix,km),cnvc(ix,km) &
2986 ! hchuang code change mass flux output
2987 &, ud_mf(im,km),dt_mf(im,km)
2988 !
2989 integer i,j,indx, k, kk, km1, n
2990 integer kpbl(im)
2991 integer, dimension(im), intent(in) :: islimsk
2992 !
2993 real(kind=kind_phys) dellat, delta, &
2994 & c0l, c0s, d0, &
2995 & c1, asolfac, &
2996 & desdt, dp, &
2997 & dq, dqsdp, dqsdt, dt, &
2998 & dt2, dtmax, dtmin, dxcrt, &
2999 & dv1h, dv2h, dv3h, &
3000 & dv1q, dv2q, dv3q, &
3001 & dz, dz1, e1, clam, &
3002 & el2orc, elocp, aafac, cm, &
3003 & es, etah, h1, &
3004 & evef, evfact, evfactl, fact1, &
3005 & fact2, factor, dthk, &
3006 & g, gamma, pprime, betaw, &
3007 & qlk, qrch, qs, &
3008 & rfact, shear, tfac, &
3009 & val, val1, val2, &
3010 & w1, w1l, w1s, w2, &
3011 & w2l, w2s, w3, w3l, &
3012 & w3s, w4, w4l, w4s, &
3013 & rho, tem, tem1, tem2, &
3014 & ptem, ptem1, &
3015 & pgcon
3016 !
3017 logical,optional,intent(in) :: pert_sas
3018 integer,optional,intent(in) :: ens_random_seed
3019 real,optional,intent(in) :: ens_sasamp
3021 integer kb(im), kbcon(im), kbcon1(im), &
3022 & ktcon(im), ktcon1(im), ktconn(im), &
3023 & kbm(im), kmax(im)
3024 !
3025 real(kind=kind_phys) aa1(im), cina(im), &
3026 & umean(im), tauadv(im), gdx(im), &
3027 & delhbar(im), delq(im), delq2(im), &
3028 & delqbar(im), delqev(im), deltbar(im), &
3029 & deltv(im), dtconv(im), edt(im), &
3030 & pdot(im), po(im,km), &
3031 & qcond(im), qevap(im), hmax(im), &
3032 & rntot(im), vshear(im), &
3033 & xlamud(im), xmb(im), xmbmax(im), &
3034 & delubar(im), delvbar(im)
3035 !
3036 real(kind=kind_phys) c0(im)
3037 !c
3038 real(kind=kind_phys) crtlamd
3039 !
3040 real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn, &
3041 & cinacr, cinacrmx, cinacrmn
3042 !
3043 ! parameters for updraft velocity calculation
3044 real(kind=kind_phys) bet1, cd1, f1, gam1, &
3045 & bb1, bb2, wucb
3046 !cc
3047 !c physical parameters
3048 parameter(g=grav,asolfac=0.89)
3049 parameter(elocp=hvap/cp, &
3050 & el2orc=hvap*hvap/(rv*cp))
3051 parameter(c0s=0.002,c1=5.e-4,d0=.01)
3052 parameter(c0l=c0s*asolfac)
3053 !
3054 ! asolfac: aerosol-aware parameter based on Lim & Hong (2012)
3055 ! asolfac= cx / c0s(=.002)
3056 ! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s)
3057 ! Nccn: CCN number concentration in cm^(-3)
3058 ! Until a realistic Nccn is provided, typical Nccns are assumed
3059 ! as Nccn=100 for sea and Nccn=7000 for land
3060 !
3061 parameter(cm=1.0,delta=fv)
3062 parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
3063 parameter(dthk=25.)
3064 parameter(cinpcrmx=180.,cinpcrmn=120.)
3065 parameter(cinacrmx=-120.,cinacrmn=-120.)
3066 parameter(crtlamd=3.e-4)
3067 parameter(dtmax=10800.,dtmin=600.)
3068 parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
3069 parameter(betaw=.03,dxcrt=15.e3)
3070 parameter(h1=0.33333333)
3071 !c local variables and arrays
3072 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), &
3073 & uo(im,km), vo(im,km), qeso(im,km)
3074 ! for updraft velocity calculation
3075 real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km)
3076 real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im)
3077 !
3078 !c cloud water
3079 ! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
3080 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), &
3081 & dbyo(im,km), zo(im,km), xlamue(im,km), &
3082 & heo(im,km), heso(im,km), &
3083 & dellah(im,km), dellaq(im,km), &
3084 & dellau(im,km), dellav(im,km), hcko(im,km), &
3085 & ucko(im,km), vcko(im,km), qcko(im,km), &
3086 & qrcko(im,km), eta(im,km), &
3087 & zi(im,km), pwo(im,km), c0t(im,km), &
3088 & sumx(im), tx1(im), cnvwt(im,km)
3089 !
3090 logical totflg, cnvflg(im), flg(im)
3091 !
3092 real(kind=kind_phys) tf, tcr, tcrf
3093 parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
3095 #if HWRF==1
3096 real*8 :: gasdev,ran1 !zhang
3097 real :: rr !zhang
3098 logical,save :: pert_sas_local !zhang
3099 integer,save :: ens_random_seed_local,env_pp_local !zhang
3100 integer :: ensda_physics_pert !zhang
3101 real,save :: ens_sasamp_local !zhang
3102 data ens_random_seed_local/0/
3103 data env_pp_local/0/
3104 CHARACTER(len=3) :: env_memb,env_pp
3105 if ( ens_random_seed_local .eq. 0 ) then
3106 CALL nl_get_ensda_physics_pert(1,ensda_physics_pert)
3107 ens_random_seed_local=ens_random_seed
3108 env_pp_local=ensda_physics_pert
3109 pert_sas_local=.false.
3110 ens_sasamp_local=0.0
3111 ! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99
3112 if ( env_pp_local .eq. 1 ) then
3113 if ( ens_random_seed .ne. 99 ) then
3114 pert_sas_local=.true.
3115 ens_sasamp_local=ens_sasamp
3116 else
3117 ! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero
3118 ens_random_seed_local=ens_random_seed
3119 pert_sas_local=pert_sas
3120 ens_sasamp_local=ens_sasamp
3121 endif
3122 else
3123 ens_random_seed_local=ens_random_seed
3124 pert_sas_local=pert_sas
3125 ens_sasamp_local=ens_sasamp
3126 endif
3128 print*, "SHSAS ==", ens_random_seed_local,pert_sas_local,ens_sasamp_local,ensda_physics_pert
3129 endif
3130 #endif
3132 !
3133 !c-----------------------------------------------------------------------
3134 !
3135 !************************************************************************
3136 ! convert input Pa terms to Cb terms -- Moorthi
3137 ps = psp * 0.001
3138 prsl = prslp * 0.001
3139 del = delp * 0.001
3140 !************************************************************************
3141 !
3142 km1 = km - 1
3143 !c
3144 !c initialize arrays
3145 !c
3146 do i=1,im
3147 cnvflg(i) = .true.
3148 if(kcnv(i) == 1) cnvflg(i) = .false.
3149 if(cnvflg(i)) then
3150 kbot(i)=km+1
3151 ktop(i)=0
3152 endif
3153 rn(i)=0.
3154 kbcon(i)=km
3155 ktcon(i)=1
3156 ktconn(i)=1
3157 kb(i)=km
3158 pdot(i) = 0.
3159 qlko_ktcon(i) = 0.
3160 edt(i) = 0.
3161 aa1(i) = 0.
3162 cina(i) = 0.
3163 vshear(i) = 0.
3164 gdx(i) = sqrt(garea(i))
3165 scaldfunc(i)=-1.0 ! wang initialized
3166 sigmagfm(i)=-1.0
3167 enddo
3168 !!
3169 totflg = .true.
3170 do i=1,im
3171 totflg = totflg .and. (.not. cnvflg(i))
3172 enddo
3173 if(totflg) return
3174 !!
3175 do i=1,im
3176 if(islimsk(i) == 1) then
3177 c0(i) = c0l
3178 else
3179 c0(i) = c0s
3180 endif
3181 enddo
3182 !
3183 do k = 1, km
3184 do i = 1, im
3185 if(t1(i,k) > 273.16) then
3186 c0t(i,k) = c0(i)
3187 else
3188 tem = d0 * (t1(i,k) - 273.16)
3189 tem1 = exp(tem)
3190 c0t(i,k) = c0(i) * tem1
3191 endif
3192 enddo
3193 enddo
3194 !
3195 do k = 1, km
3196 do i = 1, im
3197 cnvw(i,k) = 0.
3198 cnvc(i,k) = 0.
3199 enddo
3200 enddo
3201 ! hchuang code change
3202 do k = 1, km
3203 do i = 1, im
3204 ud_mf(i,k) = 0.
3205 dt_mf(i,k) = 0.
3206 enddo
3207 enddo
3208 !c
3209 dt2 = delt
3210 !
3211 !c model tunable parameters are all here
3212 clam = .3
3213 aafac = .1
3214 !c evef = 0.07
3215 evfact = 0.3
3216 evfactl = 0.3
3217 !
3218 ! pgcon = 0.7 ! Gregory et al. (1997, QJRMS)
3219 pgcon = 0.55 ! Zhang & Wu (2003,JAS)
3220 w1l = -8.e-3
3221 w2l = -4.e-2
3222 w3l = -5.e-3
3223 w4l = -5.e-4
3224 w1s = -2.e-4
3225 w2s = -2.e-3
3226 w3s = -1.e-3
3227 w4s = -2.e-5
3228 !c
3229 !c define top layer for search of the downdraft originating layer
3230 !c and the maximum thetae for updraft
3231 !c
3232 do i=1,im
3233 kbm(i) = km
3234 kmax(i) = km
3235 tx1(i) = 1.0 / ps(i)
3236 enddo
3237 !
3238 do k = 1, km
3239 do i=1,im
3240 if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1
3241 if (prsl(i,k)*tx1(i) > 0.60) kmax(i) = k + 1
3242 enddo
3243 enddo
3244 do i=1,im
3245 kbm(i) = min(kbm(i),kmax(i))
3246 enddo
3247 !c
3248 !c hydrostatic height assume zero terr and compute
3249 !c updraft entrainment rate as an inverse function of height
3250 !c
3251 do k = 1, km
3252 do i=1,im
3253 zo(i,k) = phil(i,k) / g
3254 enddo
3255 enddo
3256 do k = 1, km1
3257 do i=1,im
3258 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
3259 xlamue(i,k) = clam / zi(i,k)
3260 enddo
3261 enddo
3262 do i=1,im
3263 xlamue(i,km) = xlamue(i,km1)
3264 enddo
3265 !c
3266 !c pbl height
3267 !c
3268 do i=1,im
3269 flg(i) = cnvflg(i)
3270 kpbl(i)= 1
3271 enddo
3272 do k = 2, km1
3273 do i=1,im
3274 if (flg(i) .and. zo(i,k) <= hpbl(i)) then
3275 kpbl(i) = k
3276 else
3277 flg(i) = .false.
3278 endif
3279 enddo
3280 enddo
3281 do i=1,im
3282 kpbl(i)= min(kpbl(i),kbm(i))
3283 enddo
3284 !c
3285 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3286 !c convert surface pressure to mb from cb
3287 !c
3288 do k = 1, km
3289 do i = 1, im
3290 if (cnvflg(i) .and. k <= kmax(i)) then
3291 pfld(i,k) = prsl(i,k) * 10.0
3292 eta(i,k) = 1.
3293 hcko(i,k) = 0.
3294 qcko(i,k) = 0.
3295 qrcko(i,k)= 0.
3296 ucko(i,k) = 0.
3297 vcko(i,k) = 0.
3298 dbyo(i,k) = 0.
3299 pwo(i,k) = 0.
3300 dellal(i,k) = 0.
3301 to(i,k) = t1(i,k)
3302 qo(i,k) = q1(i,k)
3303 uo(i,k) = u1(i,k)
3304 vo(i,k) = v1(i,k)
3305 ! uo(i,k) = u1(i,k) * rcs(i)
3306 ! vo(i,k) = v1(i,k) * rcs(i)
3307 wu2(i,k) = 0.
3308 buo(i,k) = 0.
3309 drag(i,k) = 0.
3310 cnvwt(i,k) = 0.
3311 endif
3312 enddo
3313 enddo
3314 !c
3315 !c column variables
3316 !c p is pressure of the layer (mb)
3317 !c t is temperature at t-dt (k)..tn
3318 !c q is mixing ratio at t-dt (kg/kg)..qn
3319 !c to is temperature at t+dt (k)... this is after advection and turbulan
3320 !c qo is mixing ratio at t+dt (kg/kg)..q1
3321 !c
3322 do k = 1, km
3323 do i=1,im
3324 if (cnvflg(i) .and. k <= kmax(i)) then
3325 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
3326 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
3327 val1 = 1.e-8
3328 qeso(i,k) = max(qeso(i,k), val1)
3329 val2 = 1.e-10
3330 qo(i,k) = max(qo(i,k), val2 )
3331 ! qo(i,k) = min(qo(i,k),qeso(i,k))
3332 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
3333 endif
3334 enddo
3335 enddo
3336 !c
3337 !c compute moist static energy
3338 !c
3339 do k = 1, km
3340 do i=1,im
3341 if (cnvflg(i) .and. k <= kmax(i)) then
3342 ! tem = g * zo(i,k) + cp * to(i,k)
3343 tem = phil(i,k) + cp * to(i,k)
3344 heo(i,k) = tem + hvap * qo(i,k)
3345 heso(i,k) = tem + hvap * qeso(i,k)
3346 !c heo(i,k) = min(heo(i,k),heso(i,k))
3347 endif
3348 enddo
3349 enddo
3350 !c
3351 !c determine level with largest moist static energy within pbl
3352 !c this is the level where updraft starts
3353 !c
3354 do i=1,im
3355 if (cnvflg(i)) then
3356 hmax(i) = heo(i,1)
3357 kb(i) = 1
3358 endif
3359 enddo
3360 do k = 2, km
3361 do i=1,im
3362 if (cnvflg(i) .and. k <= kpbl(i)) then
3363 if(heo(i,k) > hmax(i)) then
3364 kb(i) = k
3365 hmax(i) = heo(i,k)
3366 endif
3367 endif
3368 enddo
3369 enddo
3370 !c
3371 do k = 1, km1
3372 do i=1,im
3373 if (cnvflg(i) .and. k <= kmax(i)-1) then
3374 dz = .5 * (zo(i,k+1) - zo(i,k))
3375 dp = .5 * (pfld(i,k+1) - pfld(i,k))
3376 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
3377 pprime = pfld(i,k+1) + epsm1 * es
3378 qs = eps * es / pprime
3379 dqsdp = - qs / pprime
3380 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
3381 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
3382 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
3383 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
3384 dq = dqsdt * dt + dqsdp * dp
3385 to(i,k) = to(i,k+1) + dt
3386 qo(i,k) = qo(i,k+1) + dq
3387 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
3388 endif
3389 enddo
3390 enddo
3391 !
3392 do k = 1, km1
3393 do i=1,im
3394 if (cnvflg(i) .and. k <= kmax(i)-1) then
3395 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
3396 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
3397 val1 = 1.e-8
3398 qeso(i,k) = max(qeso(i,k), val1)
3399 val2 = 1.e-10
3400 qo(i,k) = max(qo(i,k), val2 )
3401 ! qo(i,k) = min(qo(i,k),qeso(i,k))
3402 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3403 & cp * to(i,k) + hvap * qo(i,k)
3404 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3405 & cp * to(i,k) + hvap * qeso(i,k)
3406 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
3407 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
3408 endif
3409 enddo
3410 enddo
3411 !c
3412 !c look for the level of free convection as cloud base
3413 !c
3414 do i=1,im
3415 flg(i) = cnvflg(i)
3416 if(flg(i)) kbcon(i) = kmax(i)
3417 enddo
3418 do k = 2, km1
3419 do i=1,im
3420 if (flg(i) .and. k < kbm(i)) then
3421 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k)) then
3422 kbcon(i) = k
3423 flg(i) = .false.
3424 endif
3425 endif
3426 enddo
3427 enddo
3428 !c
3429 do i=1,im
3430 if(cnvflg(i)) then
3431 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
3432 endif
3433 enddo
3434 !!
3435 totflg = .true.
3436 do i=1,im
3437 totflg = totflg .and. (.not. cnvflg(i))
3438 enddo
3439 if(totflg) return
3440 !!
3441 do i=1,im
3442 if(cnvflg(i)) then
3443 ! pdot(i) = 10.* dot(i,kbcon(i))
3444 pdot(i) = 0.01 * dot(i,kbcon(i)) ! Now dot is in Pa/s
3445 endif
3446 enddo
3447 !c
3448 !c turn off convection if pressure depth between parcel source level
3449 !c and cloud base is larger than a critical value, cinpcr
3450 !c
3451 do i=1,im
3452 if(cnvflg(i)) then
3453 if(islimsk(i) == 1) then
3454 w1 = w1l
3455 w2 = w2l
3456 w3 = w3l
3457 w4 = w4l
3458 else
3459 w1 = w1s
3460 w2 = w2s
3461 w3 = w3s
3462 w4 = w4s
3463 endif
3464 if(pdot(i) <= w4) then
3465 tem = (pdot(i) - w4) / (w3 - w4)
3466 elseif(pdot(i) >= -w4) then
3467 tem = - (pdot(i) + w4) / (w4 - w3)
3468 else
3469 tem = 0.
3470 endif
3471 val1 = -1.
3472 tem = max(tem,val1)
3473 val2 = 1.
3474 tem = min(tem,val2)
3475 ptem = 1. - tem
3476 ptem1= .5*(cinpcrmx-cinpcrmn)
3477 cinpcr = cinpcrmx - ptem * ptem1
3478 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
3479 #if HWRF==1
3480 ! randomly perturb the convection trigger
3481 !zzz if( pert_sas_local .and. ens_random_seed_local .gt. 0 ) then
3482 if( pert_sas_local ) then
3483 !zz print*, "zhang inde ens_random_seed=", ens_random_seed,ens_random_seed_local
3484 ens_random_seed_local=ran1(-ens_random_seed_local)*1000
3485 rr=2.0*ens_sasamp_local*ran1(-ens_random_seed_local)-ens_sasamp_local
3486 !zz print*, "zhang inde shsas=a", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr
3487 cinpcr=cinpcr+rr
3488 !zz print*, "zhang inde shsas=b", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr
3489 endif
3490 #endif
3491 if(tem1 > cinpcr) then
3492 cnvflg(i) = .false.
3493 endif
3494 endif
3495 enddo
3496 !!
3497 totflg = .true.
3498 do i=1,im
3499 totflg = totflg .and. (.not. cnvflg(i))
3500 enddo
3501 if(totflg) return
3502 !!
3503 !c
3504 !c specify the detrainment rate for the updrafts
3505 !c
3506 do i = 1, im
3507 if(cnvflg(i)) then
3508 xlamud(i) = xlamue(i,kbcon(i))
3509 ! xlamud(i) = crtlamd
3510 endif
3511 enddo
3512 !c
3513 !c determine updraft mass flux for the subcloud layers
3514 !c
3515 do k = km1, 1, -1
3516 do i = 1, im
3517 if (cnvflg(i)) then
3518 if(k < kbcon(i) .and. k >= kb(i)) then
3519 dz = zi(i,k+1) - zi(i,k)
3520 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
3521 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
3522 endif
3523 endif
3524 enddo
3525 enddo
3526 !c
3527 !c compute mass flux above cloud base
3528 !c
3529 do i = 1, im
3530 flg(i) = cnvflg(i)
3531 enddo
3532 do k = 2, km1
3533 do i = 1, im
3534 if(flg(i))then
3535 if(k > kbcon(i) .and. k < kmax(i)) then
3536 dz = zi(i,k) - zi(i,k-1)
3537 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
3538 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
3539 if(eta(i,k) <= 0.) then
3540 kmax(i) = k
3541 ktconn(i) = k
3542 kbm(i) = min(kbm(i),kmax(i))
3543 flg(i) = .false.
3544 endif
3545 endif
3546 endif
3547 enddo
3548 enddo
3549 !c
3550 !c compute updraft cloud property
3551 !c
3552 do i = 1, im
3553 if(cnvflg(i)) then
3554 indx = kb(i)
3555 hcko(i,indx) = heo(i,indx)
3556 ucko(i,indx) = uo(i,indx)
3557 vcko(i,indx) = vo(i,indx)
3558 endif
3559 enddo
3560 !c
3561 ! cm is an enhancement factor in entrainment rates for momentum
3562 !
3563 do k = 2, km1
3564 do i = 1, im
3565 if (cnvflg(i)) then
3566 if(k > kb(i) .and. k < kmax(i)) then
3567 dz = zi(i,k) - zi(i,k-1)
3568 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3569 tem1 = 0.5 * xlamud(i) * dz
3570 factor = 1. + tem - tem1
3571 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
3572 & (heo(i,k)+heo(i,k-1)))/factor
3573 dbyo(i,k) = hcko(i,k) - heso(i,k)
3574 !
3575 tem = 0.5 * cm * tem
3576 factor = 1. + tem
3577 ptem = tem + pgcon
3578 ptem1= tem - pgcon
3579 ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*uo(i,k) &
3580 & +ptem1*uo(i,k-1))/factor
3581 vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*vo(i,k) &
3582 & +ptem1*vo(i,k-1))/factor
3583 endif
3584 endif
3585 enddo
3586 enddo
3587 !c
3588 !c taking account into convection inhibition due to existence of
3589 !c dry layers below cloud base
3590 !c
3591 do i=1,im
3592 flg(i) = cnvflg(i)
3593 kbcon1(i) = kmax(i)
3594 enddo
3595 do k = 2, km1
3596 do i=1,im
3597 if (flg(i) .and. k < kbm(i)) then
3598 if(k >= kbcon(i) .and. dbyo(i,k) > 0.) then
3599 kbcon1(i) = k
3600 flg(i) = .false.
3601 endif
3602 endif
3603 enddo
3604 enddo
3605 do i=1,im
3606 if(cnvflg(i)) then
3607 if(kbcon1(i) == kmax(i)) cnvflg(i) = .false.
3608 endif
3609 enddo
3610 do i=1,im
3611 if(cnvflg(i)) then
3612 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
3613 if(tem > dthk) then
3614 cnvflg(i) = .false.
3615 endif
3616 endif
3617 enddo
3618 !!
3619 totflg = .true.
3620 do i = 1, im
3621 totflg = totflg .and. (.not. cnvflg(i))
3622 enddo
3623 if(totflg) return
3624 !!
3625 !c
3626 !c calculate convective inhibition
3627 !c
3628 do k = 2, km1
3629 do i = 1, im
3630 if (cnvflg(i)) then
3631 if(k > kb(i) .and. k < kbcon1(i)) then
3632 dz1 = zo(i,k+1) - zo(i,k)
3633 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3634 rfact = 1. + delta * cp * gamma &
3635 & * to(i,k) / hvap
3636 cina(i) = cina(i) + &
3637 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
3638 & dz1 * (g / (cp * to(i,k))) &
3639 & * dbyo(i,k) / (1. + gamma) &
3640 & * rfact
3641 val = 0.
3642 cina(i) = cina(i) + &
3643 ! & dz1 * eta(i,k) * g * delta *
3644 & dz1 * g * delta * &
3645 & max(val,(qeso(i,k) - qo(i,k)))
3646 endif
3647 endif
3648 enddo
3649 enddo
3650 do i = 1, im
3651 if(cnvflg(i)) then
3652 !
3653 ! if(islimsk(i) == 1) then
3654 ! w1 = w1l
3655 ! w2 = w2l
3656 ! w3 = w3l
3657 ! w4 = w4l
3658 ! else
3659 ! w1 = w1s
3660 ! w2 = w2s
3661 ! w3 = w3s
3662 ! w4 = w4s
3663 ! endif
3664 ! if(pdot(i) <= w4) then
3665 ! tem = (pdot(i) - w4) / (w3 - w4)
3666 ! elseif(pdot(i) >= -w4) then
3667 ! tem = - (pdot(i) + w4) / (w4 - w3)
3668 ! else
3669 ! tem = 0.
3670 ! endif
3671 !
3672 ! val1 = -1.
3673 ! tem = max(tem,val1)
3674 ! val2 = 1.
3675 ! tem = min(tem,val2)
3676 ! tem = 1. - tem
3677 ! tem1= .5*(cinacrmx-cinacrmn)
3678 ! cinacr = cinacrmx - tem * tem1
3679 !
3680 cinacr = cinacrmx
3681 if(cina(i) < cinacr) cnvflg(i) = .false.
3682 endif
3683 enddo
3684 !!
3685 totflg = .true.
3686 do i=1,im
3687 totflg = totflg .and. (.not. cnvflg(i))
3688 enddo
3689 if(totflg) return
3690 !!
3691 !c
3692 !c determine first guess cloud top as the level of zero buoyancy
3693 !c limited to the level of P/Ps=0.7
3694 !c
3695 do i = 1, im
3696 flg(i) = cnvflg(i)
3697 if(flg(i)) ktcon(i) = kbm(i)
3698 enddo
3699 do k = 2, km1
3700 do i=1,im
3701 if (flg(i) .and. k < kbm(i)) then
3702 if(k > kbcon1(i) .and. dbyo(i,k) < 0.) then
3703 ktcon(i) = k
3704 flg(i) = .false.
3705 endif
3706 endif
3707 enddo
3708 enddo
3709 !c
3710 !c specify upper limit of mass flux at cloud base
3711 !c
3712 do i = 1, im
3713 if(cnvflg(i)) then
3714 ! xmbmax(i) = .1
3715 !
3716 k = kbcon(i)
3717 dp = 1000. * del(i,k)
3718 xmbmax(i) = dp / (g * dt2)
3719 !
3720 ! tem = dp / (g * dt2)
3721 ! xmbmax(i) = min(tem, xmbmax(i))
3722 endif
3723 enddo
3724 !c
3725 !c compute cloud moisture property and precipitation
3726 !c
3727 do i = 1, im
3728 if (cnvflg(i)) then
3729 aa1(i) = 0.
3730 qcko(i,kb(i)) = qo(i,kb(i))
3731 qrcko(i,kb(i)) = qo(i,kb(i))
3732 endif
3733 enddo
3734 do k = 2, km1
3735 do i = 1, im
3736 if (cnvflg(i)) then
3737 if(k > kb(i) .and. k < ktcon(i)) then
3738 dz = zi(i,k) - zi(i,k-1)
3739 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3740 qrch = qeso(i,k) &
3741 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3742 !cj
3743 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3744 tem1 = 0.5 * xlamud(i) * dz
3745 factor = 1. + tem - tem1
3746 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
3747 & (qo(i,k)+qo(i,k-1)))/factor
3748 qrcko(i,k) = qcko(i,k)
3749 !cj
3750 dq = eta(i,k) * (qcko(i,k) - qrch)
3751 !c
3752 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
3753 !c
3754 !c below lfc check if there is excess moisture to release latent heat
3755 !c
3756 if(k >= kbcon(i) .and. dq > 0.) then
3757 etah = .5 * (eta(i,k) + eta(i,k-1))
3758 if(ncloud > 0) then
3759 dp = 1000. * del(i,k)
3760 ptem = c0t(i,k) + c1
3761 qlk = dq / (eta(i,k) + etah * ptem * dz)
3762 dellal(i,k) = etah * c1 * dz * qlk * g / dp
3763 else
3764 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
3765 endif
3766 buo(i,k) = buo(i,k) - g * qlk
3767 qcko(i,k)= qlk + qrch
3768 pwo(i,k) = etah * c0t(i,k) * dz * qlk
3769 cnvwt(i,k) = etah * qlk * g / dp
3770 endif
3771 !
3772 ! compute buoyancy and drag for updraft velocity
3773 !
3774 if(k >= kbcon(i)) then
3775 rfact = 1. + delta * cp * gamma &
3776 & * to(i,k) / hvap
3777 buo(i,k) = buo(i,k) + (g / (cp * to(i,k))) &
3778 & * dbyo(i,k) / (1. + gamma) &
3779 & * rfact
3780 val = 0.
3781 buo(i,k) = buo(i,k) + g * delta * &
3782 & max(val,(qeso(i,k) - qo(i,k)))
3783 drag(i,k) = max(xlamue(i,k),xlamud(i))
3784 endif
3785 !
3786 endif
3787 endif
3788 enddo
3789 enddo
3790 !c
3791 !c calculate cloud work function
3792 !c
3793 ! do k = 2, km1
3794 ! do i = 1, im
3795 ! if (cnvflg(i)) then
3796 ! if(k >= kbcon(i) .and. k < ktcon(i)) then
3797 ! dz1 = zo(i,k+1) - zo(i,k)
3798 ! gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3799 ! rfact = 1. + delta * cp * gamma
3800 ! & * to(i,k) / hvap
3801 ! aa1(i) = aa1(i) +
3802 !! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
3803 ! & dz1 * (g / (cp * to(i,k)))
3804 ! & * dbyo(i,k) / (1. + gamma)
3805 ! & * rfact
3806 ! val = 0.
3807 ! aa1(i) = aa1(i) +
3808 !! & dz1 * eta(i,k) * g * delta *
3809 ! & dz1 * g * delta *
3810 ! & max(val,(qeso(i,k) - qo(i,k)))
3811 ! endif
3812 ! endif
3813 ! enddo
3814 ! enddo
3815 ! do i = 1, im
3816 ! if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
3817 ! enddo
3818 !
3819 ! calculate cloud work function
3820 !
3821 do i = 1, im
3822 if (cnvflg(i)) then
3823 aa1(i) = 0.
3824 endif
3825 enddo
3826 do k = 2, km1
3827 do i = 1, im
3828 if (cnvflg(i)) then
3829 if(k >= kbcon(i) .and. k < ktcon(i)) then
3830 dz1 = zo(i,k+1) - zo(i,k)
3831 aa1(i) = aa1(i) + buo(i,k) * dz1
3832 endif
3833 endif
3834 enddo
3835 enddo
3836 do i = 1, im
3837 if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
3838 enddo
3839 !!
3840 totflg = .true.
3841 do i=1,im
3842 totflg = totflg .and. (.not. cnvflg(i))
3843 enddo
3844 if(totflg) return
3845 !!
3846 !c
3847 !c estimate the onvective overshooting as the level
3848 !c where the [aafac * cloud work function] becomes zero,
3849 !c which is the final cloud top
3850 !c limited to the level of P/Ps=0.7
3851 !c
3852 do i = 1, im
3853 if (cnvflg(i)) then
3854 aa1(i) = aafac * aa1(i)
3855 endif
3856 enddo
3857 !c
3858 do i = 1, im
3859 flg(i) = cnvflg(i)
3860 ktcon1(i) = kbm(i)
3861 enddo
3862 do k = 2, km1
3863 do i = 1, im
3864 if (flg(i)) then
3865 if(k >= ktcon(i) .and. k < kbm(i)) then
3866 dz1 = zo(i,k+1) - zo(i,k)
3867 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3868 rfact = 1. + delta * cp * gamma &
3869 & * to(i,k) / hvap
3870 aa1(i) = aa1(i) + &
3871 ! & dz1 * eta(i,k) * (g / (cp * to(i,k)))
3872 & dz1 * (g / (cp * to(i,k))) &
3873 & * dbyo(i,k) / (1. + gamma) &
3874 & * rfact
3875 ! val = 0.
3876 ! aa1(i) = aa1(i) +
3877 !! & dz1 * eta(i,k) * g * delta *
3878 ! & dz1 * g * delta *
3879 ! & max(val,(qeso(i,k) - qo(i,k)))
3880 if(aa1(i) < 0.) then
3881 ktcon1(i) = k
3882 flg(i) = .false.
3883 endif
3884 endif
3885 endif
3886 enddo
3887 enddo
3888 !c
3889 !c compute cloud moisture property, detraining cloud water
3890 !c and precipitation in overshooting layers
3891 !c
3892 do k = 2, km1
3893 do i = 1, im
3894 if (cnvflg(i)) then
3895 if(k >= ktcon(i) .and. k < ktcon1(i)) then
3896 dz = zi(i,k) - zi(i,k-1)
3897 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3898 qrch = qeso(i,k) &
3899 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3900 !cj
3901 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3902 tem1 = 0.5 * xlamud(i) * dz
3903 factor = 1. + tem - tem1
3904 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
3905 & (qo(i,k)+qo(i,k-1)))/factor
3906 qrcko(i,k) = qcko(i,k)
3907 !cj
3908 dq = eta(i,k) * (qcko(i,k) - qrch)
3909 !c
3910 !c check if there is excess moisture to release latent heat
3911 !c
3912 if(dq > 0.) then
3913 etah = .5 * (eta(i,k) + eta(i,k-1))
3914 if(ncloud > 0) then
3915 dp = 1000. * del(i,k)
3916 ptem = c0t(i,k) + c1
3917 qlk = dq / (eta(i,k) + etah * ptem * dz)
3918 dellal(i,k) = etah * c1 * dz * qlk * g / dp
3919 else
3920 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
3921 endif
3922 qcko(i,k) = qlk + qrch
3923 pwo(i,k) = etah * c0t(i,k) * dz * qlk
3924 cnvwt(i,k) = etah * qlk * g / dp
3925 endif
3926 endif
3927 endif
3928 enddo
3929 enddo
3930 !
3931 ! compute updraft velocity square(wu2)
3932 !
3933 ! bb1 = 2. * (1.+bet1*cd1)
3934 ! bb2 = 2. / (f1*(1.+gam1))
3935 !
3936 ! bb1 = 3.9
3937 ! bb2 = 0.67
3938 !
3939 ! bb1 = 2.0
3940 ! bb2 = 4.0
3941 !
3942 bb1 = 4.0
3943 bb2 = 0.8
3944 !
3945 do i = 1, im
3946 if (cnvflg(i)) then
3947 k = kbcon1(i)
3948 tem = po(i,k) / (rd * to(i,k))
3949 wucb = -0.01 * dot(i,k) / (tem * g)
3950 if(wucb > 0.) then
3951 wu2(i,k) = wucb * wucb
3952 else
3953 wu2(i,k) = 0.
3954 endif
3955 endif
3956 enddo
3957 do k = 2, km1
3958 do i = 1, im
3959 if (cnvflg(i)) then
3960 if(k > kbcon1(i) .and. k < ktcon(i)) then
3961 dz = zi(i,k) - zi(i,k-1)
3962 tem = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz
3963 tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz
3964 ptem = (1. - tem) * wu2(i,k-1)
3965 ptem1 = 1. + tem
3966 wu2(i,k) = (ptem + tem1) / ptem1
3967 wu2(i,k) = max(wu2(i,k), 0._kind_phys)
3968 endif
3969 endif
3970 enddo
3971 enddo
3972 !
3973 ! compute updraft velocity averaged over the whole cumulus
3974 !
3975 do i = 1, im
3976 wc(i) = 0.
3977 sumx(i) = 0.
3978 enddo
3979 do k = 2, km1
3980 do i = 1, im
3981 if (cnvflg(i)) then
3982 if(k > kbcon1(i) .and. k < ktcon(i)) then
3983 dz = zi(i,k) - zi(i,k-1)
3984 tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
3985 wc(i) = wc(i) + tem * dz
3986 sumx(i) = sumx(i) + dz
3987 endif
3988 endif
3989 enddo
3990 enddo
3991 do i = 1, im
3992 if(cnvflg(i)) then
3993 if(sumx(i) == 0.) then
3994 cnvflg(i)=.false.
3995 else
3996 wc(i) = wc(i) / sumx(i)
3997 endif
3998 val = 1.e-4
3999 if (wc(i) < val) cnvflg(i)=.false.
4000 endif
4001 enddo
4002 !c
4003 !c exchange ktcon with ktcon1
4004 !c
4005 do i = 1, im
4006 if(cnvflg(i)) then
4007 kk = ktcon(i)
4008 ktcon(i) = ktcon1(i)
4009 ktcon1(i) = kk
4010 endif
4011 enddo
4012 !c
4013 !c this section is ready for cloud water
4014 !c
4015 if(ncloud > 0) then
4016 !c
4017 !c compute liquid and vapor separation at cloud top
4018 !c
4019 do i = 1, im
4020 if(cnvflg(i)) then
4021 k = ktcon(i) - 1
4022 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
4023 qrch = qeso(i,k) &
4024 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
4025 dq = qcko(i,k) - qrch
4026 !c
4027 !c check if there is excess moisture to release latent heat
4028 !c
4029 if(dq > 0.) then
4030 qlko_ktcon(i) = dq
4031 qcko(i,k) = qrch
4032 endif
4033 endif
4034 enddo
4035 endif
4036 !c
4037 !c--- compute precipitation efficiency in terms of windshear
4038 !c
4039 do i = 1, im
4040 if(cnvflg(i)) then
4041 vshear(i) = 0.
4042 endif
4043 enddo
4044 do k = 2, km
4045 do i = 1, im
4046 if (cnvflg(i)) then
4047 if(k > kb(i) .and. k <= ktcon(i)) then
4048 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
4049 & + (vo(i,k)-vo(i,k-1)) ** 2)
4050 vshear(i) = vshear(i) + shear
4051 endif
4052 endif
4053 enddo
4054 enddo
4055 do i = 1, im
4056 if(cnvflg(i)) then
4057 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
4058 e1=1.591-.639*vshear(i) &
4059 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
4060 edt(i)=1.-e1
4061 val = .9
4062 edt(i) = min(edt(i),val)
4063 val = .0
4064 edt(i) = max(edt(i),val)
4065 endif
4066 enddo
4067 !c
4068 !c--- what would the change be, that a cloud with unit mass
4069 !c--- will do to the environment?
4070 !c
4071 do k = 1, km
4072 do i = 1, im
4073 if(cnvflg(i) .and. k <= kmax(i)) then
4074 dellah(i,k) = 0.
4075 dellaq(i,k) = 0.
4076 dellau(i,k) = 0.
4077 dellav(i,k) = 0.
4078 endif
4079 enddo
4080 enddo
4081 !c
4082 !c--- changed due to subsidence and entrainment
4083 !c
4084 do k = 2, km1
4085 do i = 1, im
4086 if (cnvflg(i)) then
4087 if(k > kb(i) .and. k < ktcon(i)) then
4088 dp = 1000. * del(i,k)
4089 dz = zi(i,k) - zi(i,k-1)
4090 !c
4091 dv1h = heo(i,k)
4092 dv2h = .5 * (heo(i,k) + heo(i,k-1))
4093 dv3h = heo(i,k-1)
4094 dv1q = qo(i,k)
4095 dv2q = .5 * (qo(i,k) + qo(i,k-1))
4096 dv3q = qo(i,k-1)
4097 !c
4098 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
4099 tem1 = xlamud(i)
4100 !cj
4101 dellah(i,k) = dellah(i,k) + &
4102 & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h &
4103 & - tem*eta(i,k-1)*dv2h*dz &
4104 & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
4105 & ) *g/dp
4106 !cj
4107 dellaq(i,k) = dellaq(i,k) + &
4108 & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q &
4109 & - tem*eta(i,k-1)*dv2q*dz &
4110 & + tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz &
4111 & ) *g/dp
4112 !cj
4113 tem1=eta(i,k)*(uo(i,k)-ucko(i,k))
4114 tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1))
4115 dellau(i,k) = dellau(i,k) + (tem1-tem2) * g/dp
4116 !cj
4117 tem1=eta(i,k)*(vo(i,k)-vcko(i,k))
4118 tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1))
4119 dellav(i,k) = dellav(i,k) + (tem1-tem2) * g/dp
4120 !cj
4121 endif
4122 endif
4123 enddo
4124 enddo
4125 !c
4126 !c------- cloud top
4127 !c
4128 do i = 1, im
4129 if(cnvflg(i)) then
4130 indx = ktcon(i)
4131 dp = 1000. * del(i,indx)
4132 dv1h = heo(i,indx-1)
4133 dellah(i,indx) = eta(i,indx-1) * &
4134 & (hcko(i,indx-1) - dv1h) * g / dp
4135 dv1q = qo(i,indx-1)
4136 dellaq(i,indx) = eta(i,indx-1) * &
4137 & (qcko(i,indx-1) - dv1q) * g / dp
4138 dellau(i,indx) = eta(i,indx-1) * &
4139 & (ucko(i,indx-1) - uo(i,indx-1)) * g / dp
4140 dellav(i,indx) = eta(i,indx-1) * &
4141 & (vcko(i,indx-1) - vo(i,indx-1)) * g / dp
4142 !c
4143 !c cloud water
4144 !c
4145 dellal(i,indx) = eta(i,indx-1) * &
4146 & qlko_ktcon(i) * g / dp
4147 endif
4148 enddo
4149 !
4150 !
4151 ! compute convective turn-over time
4152 !
4153 do i= 1, im
4154 if(cnvflg(i)) then
4155 tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
4156 dtconv(i) = tem / wc(i)
4157 dtconv(i) = max(dtconv(i),dtmin)
4158 dtconv(i) = max(dtconv(i),dt2)
4159 dtconv(i) = min(dtconv(i),dtmax)
4160 endif
4161 enddo
4162 !
4163 ! compute advective time scale using a mean cloud layer wind speed
4164 !
4165 do i= 1, im
4166 if(cnvflg(i)) then
4167 sumx(i) = 0.
4168 umean(i) = 0.
4169 endif
4170 enddo
4171 do k = 2, km1
4172 do i = 1, im
4173 if(cnvflg(i)) then
4174 if(k >= kbcon1(i) .and. k < ktcon1(i)) then
4175 dz = zi(i,k) - zi(i,k-1)
4176 tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
4177 umean(i) = umean(i) + tem * dz
4178 sumx(i) = sumx(i) + dz
4179 endif
4180 endif
4181 enddo
4182 enddo
4183 do i= 1, im
4184 if(cnvflg(i)) then
4185 umean(i) = umean(i) / sumx(i)
4186 umean(i) = max(umean(i), 1._kind_phys)
4187 tauadv(i) = gdx(i) / umean(i)
4188 endif
4189 enddo
4190 !c
4191 !c compute cloud base mass flux as a function of the mean
4192 !c updraft velcoity
4193 !c
4194 do i= 1, im
4195 if(cnvflg(i)) then
4196 k = kbcon(i)
4197 rho = po(i,k)*100. / (rd*to(i,k))
4198 tfac = tauadv(i) / dtconv(i)
4199 tfac = min(tfac, 1._kind_phys)
4200 xmb(i) = tfac*betaw*rho*wc(i)
4201 endif
4202 enddo
4203 !
4204 !--- modified Grell & Freitas' (2014) updraft fraction which uses
4205 ! actual entrainment rate at cloud base
4206 !
4207 do i = 1, im
4208 if(cnvflg(i)) then
4209 tem = min(max(xlamue(i,kbcon(i)), 2.e-4_kind_phys), 6.e-4_kind_phys)
4210 tem = 0.2 / tem
4211 tem1 = 3.14 * tem * tem
4212 sigmagfm(i) = tem1 / garea(i)
4213 sigmagfm(i) = max(sigmagfm(i), 0.001_kind_phys)
4214 sigmagfm(i) = min(sigmagfm(i), 0.999_kind_phys)
4215 endif
4216 enddo
4217 !
4218 !--- compute scale-aware function based on Arakawa & Wu (2013)
4219 !
4220 !
4221 do i = 1, im
4222 if(cnvflg(i)) then
4223 if (gdx(i) < dxcrt) then
4224 scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
4225 scaldfunc(i) = max(min(scaldfunc(i), 1._kind_phys), 0._kind_phys)
4226 else
4227 scaldfunc(i) = 1.0
4228 endif
4229 xmb(i) = xmb(i) * scaldfunc(i)
4230 xmb(i) = min(xmb(i),xmbmax(i))
4231 endif
4232 enddo
4233 !c
4234 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4235 !c
4236 do k = 1, km
4237 do i = 1, im
4238 if (cnvflg(i) .and. k <= kmax(i)) then
4239 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
4240 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
4241 val = 1.e-8
4242 qeso(i,k) = max(qeso(i,k), val )
4243 endif
4244 enddo
4245 enddo
4246 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4247 !c
4248 do i = 1, im
4249 delhbar(i) = 0.
4250 delqbar(i) = 0.
4251 deltbar(i) = 0.
4252 delubar(i) = 0.
4253 delvbar(i) = 0.
4254 qcond(i) = 0.
4255 enddo
4256 do k = 1, km
4257 do i = 1, im
4258 if (cnvflg(i)) then
4259 if(k > kb(i) .and. k <= ktcon(i)) then
4260 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
4261 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
4262 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
4263 ! tem = 1./rcs(i)
4264 ! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
4265 ! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
4266 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
4267 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
4268 dp = 1000. * del(i,k)
4269 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
4270 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
4271 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
4272 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
4273 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
4274 endif
4275 endif
4276 enddo
4277 enddo
4278 !
4279 do k = 1, km
4280 do i = 1, im
4281 if (cnvflg(i)) then
4282 if(k > kb(i) .and. k <= ktcon(i)) then
4283 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
4284 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
4285 val = 1.e-8
4286 qeso(i,k) = max(qeso(i,k), val )
4287 endif
4288 endif
4289 enddo
4290 enddo
4291 !c
4292 do i = 1, im
4293 rntot(i) = 0.
4294 delqev(i) = 0.
4295 delq2(i) = 0.
4296 flg(i) = cnvflg(i)
4297 enddo
4298 do k = km, 1, -1
4299 do i = 1, im
4300 if (cnvflg(i)) then
4301 if(k < ktcon(i) .and. k > kb(i)) then
4302 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
4303 endif
4304 endif
4305 enddo
4306 enddo
4307 !c
4308 !c evaporating rain
4309 !c
4310 do k = km, 1, -1
4311 do i = 1, im
4312 if (k <= kmax(i)) then
4313 deltv(i) = 0.
4314 delq(i) = 0.
4315 qevap(i) = 0.
4316 if(cnvflg(i)) then
4317 if(k < ktcon(i) .and. k > kb(i)) then
4318 rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
4319 endif
4320 endif
4321 if(flg(i) .and. k < ktcon(i)) then
4322 evef = edt(i) * evfact
4323 if(islimsk(i) == 1) evef=edt(i) * evfactl
4324 ! if(islimsk(i) == 1) evef=.07
4325 !c if(islimsk(i) == 1) evef = 0.
4326 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
4327 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
4328 dp = 1000. * del(i,k)
4329 if(rn(i) > 0. .and. qcond(i) < 0.) then
4330 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
4331 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
4332 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
4333 endif
4334 if(rn(i) > 0. .and. qcond(i) < 0. .and. &
4335 & delq2(i) > rntot(i)) then
4336 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
4337 flg(i) = .false.
4338 endif
4339 if(rn(i) > 0. .and. qevap(i) > 0.) then
4340 tem = .001 * dp / g
4341 tem1 = qevap(i) * tem
4342 if(tem1 > rn(i)) then
4343 qevap(i) = rn(i) / tem
4344 rn(i) = 0.
4345 else
4346 rn(i) = rn(i) - tem1
4347 endif
4348 q1(i,k) = q1(i,k) + qevap(i)
4349 t1(i,k) = t1(i,k) - elocp * qevap(i)
4350 deltv(i) = - elocp*qevap(i)/dt2
4351 delq(i) = + qevap(i)/dt2
4352 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
4353 endif
4354 delqbar(i) = delqbar(i) + delq(i)*dp/g
4355 deltbar(i) = deltbar(i) + deltv(i)*dp/g
4356 endif
4357 endif
4358 enddo
4359 enddo
4360 !cj
4361 ! do i = 1, im
4362 ! if(me == 31 .and. cnvflg(i)) then
4363 ! if(cnvflg(i)) then
4364 ! print *, ' shallow delhbar, delqbar, deltbar = ',
4365 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
4366 ! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
4367 ! print *, ' precip =', hvap*rn(i)*1000./dt2
4368 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
4369 ! endif
4370 ! enddo
4371 !cj
4372 do i = 1, im
4373 if(cnvflg(i)) then
4374 if(rn(i) < 0. .or. .not.flg(i)) rn(i) = 0.
4375 ktop(i) = ktcon(i)
4376 kbot(i) = kbcon(i)
4377 kcnv(i) = 2
4378 endif
4379 enddo
4380 !c
4381 !c convective cloud water
4382 !c
4383 do k = 1, km
4384 do i = 1, im
4385 if (cnvflg(i)) then
4386 if (k >= kbcon(i) .and. k < ktcon(i)) then
4387 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
4388 endif
4389 endif
4390 enddo
4391 enddo
4393 !c
4394 !c convective cloud cover
4395 !c
4396 do k = 1, km
4397 do i = 1, im
4398 if (cnvflg(i)) then
4399 if (k >= kbcon(i) .and. k < ktcon(i)) then
4400 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
4401 cnvc(i,k) = min(cnvc(i,k), 0.2_kind_phys)
4402 cnvc(i,k) = max(cnvc(i,k), 0.0_kind_phys)
4403 endif
4404 endif
4405 enddo
4406 enddo
4408 !c
4409 !c cloud water
4410 !c
4411 if (ncloud > 0) then
4412 !
4413 do k = 1, km1
4414 do i = 1, im
4415 if (cnvflg(i)) then
4416 ! if (k > kb(i) .and. k <= ktcon(i)) then
4417 if (k >= kbcon(i) .and. k <= ktcon(i)) then
4418 tem = dellal(i,k) * xmb(i) * dt2
4419 tem1 = max(0.0_kind_phys, min(1.0_kind_phys, (tcr-t1(i,k))*tcrf))
4420 if (ql(i,k,2) > -999.0) then
4421 ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
4422 ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
4423 else
4424 ql(i,k,1) = ql(i,k,1) + tem
4425 endif
4426 endif
4427 endif
4428 enddo
4429 enddo
4430 !
4431 endif
4432 !
4433 ! hchuang code change
4434 !
4435 do k = 1, km
4436 do i = 1, im
4437 if(cnvflg(i)) then
4438 if(k >= kb(i) .and. k < ktop(i)) then
4439 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
4440 endif
4441 endif
4442 enddo
4443 enddo
4444 do i = 1, im
4445 if(cnvflg(i)) then
4446 k = ktop(i)-1
4447 dt_mf(i,k) = ud_mf(i,k)
4448 endif
4449 enddo
4450 !!
4451 return
4452 end subroutine mfshalcnv
4453 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4454 END MODULE module_cu_scalesas