1 !WRF:MODEL_LAYER:PHYSICS
5 ! USE module_model_constants
7 REAL, PARAMETER :: RIC = 0.25 ! critical Richardson number
8 REAL, PARAMETER :: CRANKP = 0.5 ! CRANK-NIC PARAMETER
12 !-----------------------------------------------------------------------
13 !-----------------------------------------------------------------------
14 SUBROUTINE ACMPBL(XTIME, DTPBL, ZNW, SIGMAH, &
15 U3D, V3D, PP3D, DZ8W, TH3D, T3D, &
16 QV3D, QC3D, QI3D, RR3D, &
20 PBLH, KPBL2D, REGIME, &
21 GZ1OZ0, WSPD, PSIM, MUT, &
22 RUBLTEN, RVBLTEN, RTHBLTEN, &
23 RQVBLTEN, RQCBLTEN, RQIBLTEN, &
24 ids,ide, jds,jde, kds,kde, &
25 ims,ime, jms,jme, kms,kme, &
26 its,ite, jts,jte, kts,kte)
27 !-----------------------------------------------------------------------
28 !-----------------------------------------------------------------------
30 ! THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO
31 ! THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2 (ACM2), WHICH IS A COMBINED
32 ! LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992)
35 ! Pleim (2007) A combined local and non-local closure model for the atmospheric
36 ! boundary layer. Part1: Model description and testing.
38 ! Pleim (2007) A combined local and non-local closure model for the atmospheric
39 ! boundary layer. Part2: Application and evaluation in a mesoscale
40 ! meteorology model. JAMC, 46, 1396-1409
43 ! AX 3/2005 - developed WRF version based on the MM5 PX LSM
44 ! RG and JP 7/2006 - Finished WRF adaptation
45 ! JP 12/2011 12/2011 - ACM2 modified so it's not dependent on first layer thickness.
47 !**********************************************************************
51 !-- XTIME Time since simulation start (min)
52 !-- DTPBL PBL time step
53 !-- ZNW Sigma at full layer
54 !-- SIGMAH Sigma at half layer
55 !-- U3D 3D u-velocity interpolated to theta points (m/s)
56 !-- V3D 3D v-velocity interpolated to theta points (m/s)
57 !-- PP3D Pressure at half levels (Pa)
58 !-- DZ8W dz between full levels (m)
59 !-- TH3D Potential Temperature (K)
60 !-- T3D Temperature (K)
61 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
62 !-- QC3D 3D cloud mixing ratio (Kg/Kg)
63 !-- QI3D 3D ice mixing ratio (Kg/Kg)
64 !-- RR3D 3D dry air density (kg/m^3)
65 !-- UST Friction Velocity (m/s)
66 !-- HFX Upward heat flux at the surface (w/m^2)
67 !-- QFX Upward moisture flux at the surface (Kg/m^2/s)
68 !-- TSK Surface temperature (K)
69 !-- PSFC Pressure at the surface (Pa)
70 !-- EP1 Constant for virtual temperature (r_v/r_d-1) (dimensionless)
73 !-- RD gas constant for dry air (j/kg/k)
74 !-- CPD heat capacity at constant pressure for dry air (j/kg/k)
75 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
76 !-- WSPD wind speed at lowest model level (m/s)
77 !-- PSIM similarity stability function for momentum
78 !-- MUT Total Mu : Psfc - Ptop
80 !-- ids start index for i in domain
81 !-- ide end index for i in domain
82 !-- jds start index for j in domain
83 !-- jde end index for j in domain
84 !-- kds start index for k in domain
85 !-- kde end index for k in domain
86 !-- ims start index for i in memory
87 !-- ime end index for i in memory
88 !-- jms start index for j in memory
89 !-- jme end index for j in memory
90 !-- kms start index for k in memory
91 !-- kme end index for k in memory
92 !-- jts start index for j in tile
93 !-- jte end index for j in tile
94 !-- kts start index for k in tile
95 !-- kte end index for k in tile
98 !-- PBLH PBL height (m)
99 !-- KPBL2D K index for PBL layer
100 !-- REGIME Flag indicating PBL regime (stable, unstable, etc.)
101 !-- RUBLTEN U tendency due to PBL parameterization (m/s^2)
102 !-- RVBLTEN V tendency due to PBL parameterization (m/s^2)
103 !-- RTHBLTEN Theta tendency due to PBL parameterization (K/s)
104 !-- RQVBLTEN Qv tendency due to PBL parameterization (kg/kg/s)
105 !-- RQCBLTEN Qc tendency due to PBL parameterization (kg/kg/s)
106 !-- RQIBLTEN Qi tendency due to PBL parameterization (kg/kg/s)
107 !-----------------------------------------------------------------------
108 !-----------------------------------------------------------------------
112 ! DECLARATIONS - INTEGER
113 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
114 ims,ime, jms,jme, kms,kme, &
115 its,ite, jts,jte, kts,kte, XTIME
117 ! DECLARATIONS - REAL
118 REAL, INTENT(IN) :: DTPBL, EP1, &
121 REAL, DIMENSION( kms:kme ), INTENT(IN) :: ZNW, SIGMAH
123 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
124 INTENT(IN) :: U3D, V3D, &
129 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: PSIM, GZ1OZ0, &
133 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PBLH, REGIME, UST
135 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
136 INTENT(INOUT) :: RUBLTEN, RVBLTEN, &
137 RTHBLTEN, RQVBLTEN, &
140 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: KPBL2D
147 REAL, DIMENSION( kts:kte ) :: DSIGH, DSIGHI, DSIGFI
148 REAL, DIMENSION( 0:kte ) :: SIGMAF
150 REAL, PARAMETER :: KARMAN = 0.4
161 DSIGH(K) = SIGMAF(K) - SIGMAF(K-1)
162 DSIGHI(K) = 1.0 / DSIGH(K)
166 DSIGFI(K) = 1.0 / (SIGMAH(K+1) - SIGMAH(K))
169 DSIGFI(kte) = DSIGFI(kte-1)
172 CALL ACM2D(j=J,xtime=XTIME, dtpbl=DTPBL, sigmaf=SIGMAF, sigmah=SIGMAH &
173 ,dsigfi=DSIGFI,dsighi=DSIGHI,dsigh=DSIGH &
174 ,us=u3d(ims,kms,j),vs=v3d(ims,kms,j) &
175 ,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j) &
176 ,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j) &
177 ,qis=qi3d(ims,kms,j),dzf=DZ8W(ims,kms,j) &
178 ,densx=RR3D(ims,kms,j) &
179 ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) &
180 ,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j) &
181 ,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
182 ,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt &
183 ,psfcpa=psfc(ims,j),ust=ust(ims,j) &
185 ,regime=regime(ims,j),psim=psim(ims,j) &
186 ,hfx=hfx(ims,j),qfx=qfx(ims,j) &
187 ,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j) &
188 ,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j) &
190 ,ep1=ep1,karman=karman &
191 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
192 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
193 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
196 END SUBROUTINE ACMPBL
197 !-----------------------------------------------------------------------
198 !-----------------------------------------------------------------------
201 !-----------------------------------------------------------------------
202 !-----------------------------------------------------------------------
203 SUBROUTINE ACM2D(j,XTIME, DTPBL, sigmaf, sigmah &
204 ,dsigfi,dsighi,dsigh &
205 ,us,vs,theta,tt,qvs,qcs,qis &
206 ,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp &
207 ,cpd,g,rovcp,rd,rdt,psfcpa,ust &
209 ,hfx,qfx,tg,gz1oz0,wspd ,klpbl &
212 ,ids,ide, jds,jde, kds,kde &
213 ,ims,ime, jms,jme, kms,kme &
214 ,its,ite, jts,jte, kts,kte )
216 !-----------------------------------------------------------------------
217 !-----------------------------------------------------------------------
223 REAL, DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF
224 REAL, DIMENSION( kms:kme ), INTENT(IN) :: SIGMAH
225 REAL, DIMENSION( kts:kte ), INTENT(IN) :: DSIGH, DSIGHI, DSIGFI
226 REAL , INTENT(IN) :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT
227 REAL , DIMENSION( ims:ime ), INTENT(INOUT) :: PBL, UST
229 REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, TT, &
231 REAL, DIMENSION( ims:ime, kms:kme ), intent(in) :: DZF
232 REAL, DIMENSION( ims:ime, kms:kme ), intent(inout) :: utnp, &
238 real, dimension( ims:ime ), intent(in ) :: psfcpa
239 real, dimension( ims:ime ), intent(in ) :: tg
240 real, dimension( ims:ime ), intent(inout) :: regime
241 real, dimension( ims:ime ), intent(in) :: wspd, psim, gz1oz0
242 real, dimension( ims:ime ), intent(in) :: hfx, qfx
243 real, dimension( ims:ime ), intent(in) :: mut
246 INTEGER, DIMENSION( ims:ime ), INTENT(OUT):: KLPBL
247 INTEGER, INTENT(IN) :: XTIME
248 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
249 ims,ime, jms,jme, kms,kme, &
250 its,ite, jts,jte, kts,kte, j
251 !--------------------------------------------------------------------
255 INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV
258 REAL :: TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ
260 REAL, DIMENSION( its:ite ) :: FINT, PSTAR, CPAIR
261 REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB, &
262 EDDYZ, UX, VX, THETAX, &
264 REAL, DIMENSION( its:ite, 0:kte ) :: ZF
265 REAL, DIMENSION( its:ite) :: WST, TST, QST, USTM, TSTV
266 REAL, DIMENSION( its:ite ) :: PBLSIG, MOL
267 REAL :: FINTT, ZMIX, UMIX, VMIX
268 REAL :: TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH
269 REAL :: A,TST12,RL,ZFUNC
270 ! REAL, PARAMETER :: KARMAN = 0.4
273 INTEGER :: KL, jtf, ktf, itf, KMIX, KSRC
275 character*256 :: message
276 !-----initialize vertical tendencies and
299 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
300 !!! Compute Micromet Scaling variables, not availiable in WRF for ACM
302 CPAIR(I) = CPD * (1.0 + 0.84 * QVS(I,1)) ! J/(K KG)
303 TMPFX = HFX(I) / (cpair(i) * DENSX(I,1))
304 TMPVTCON = 1.0 + EP1 * QVS(I,1) ! COnversion factor for virtual temperature
305 WS1 = SQRT(US(I,1)**2 + VS(I,1)**2) ! Level 1 wind speed
306 TST(I) = -TMPFX / UST(I)
307 QST(I) = -QFX(I) / (UST(I)*DENSX(I,1))
308 USTM(I) = UST(I) * WS1 / wspd(i)
309 THV1 = TMPVTCON * THETA(I,1)
310 TSTV(I) = TST(I)*TMPVTCON + THV1*EP1*QST(I)
311 IF(ABS(TSTV(I)).LT.1.0E-6) THEN
312 TSTV(I) = SIGN(1.0E-6,TSTV(I))
314 MOL(I) = THV1 * UST(i)**2/(KARMAN*G*TSTV(I))
315 WST(I) = UST(I) * (PBL(I)/(KARMAN*ABS(MOL(I)))) ** 0.333333
316 PSTAR(I) = MUT(I)/1000. ! P* in cb
318 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
320 !... Compute PBL height
322 !... compute the height of full- and half-sigma level above ground level
330 ZF(I,K) = DZF(I,K) + ZF(I,K-1)
331 ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1))
337 TVCON = 1.0 + EP1 * QVS(I,K)
338 THETAV(I,K) = THETA(I,K) * TVCON
343 !... COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990
347 IF (SIGMAF(K).lT.0.9955) GO TO 69
352 TH1 = TH1 + THETAV(I,K)
355 IF(MOL(I).LT.0.0 .AND. XTIME.GT.1) then
356 WSS = (UST(I) ** 3 + 0.6 * WST(I) ** 3) ** 0.33333
357 TCONV = -8.5 * UST(I) * TSTV(I) / WSS
363 DTMP = THETAV(I,K) - TH1
364 IF (DTMP.LT.0.0) KMIX = K
367 FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1) &
369 ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX)
370 UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX)
371 VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX)
378 DTMP = THETAV(I,K) - TH1
379 TOG = 0.5 * (THETAV(I,K) + TH1) / G
380 WSSQ = (US(I,K)-UMIX)**2 &
382 IF (KMIX == 1) WSSQ = WSSQ + 100.*UST(I)*UST(I)
383 WSSQ = MAX( WSSQ, 0.1 )
384 RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ)
385 IF (RIB(I,K) .GE. RIC) GO TO 201
388 write (message, *)' RIBX never exceeds RIC, RIB(i,kte) = ',rib(i,5), &
389 ' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i), &
390 ' TCONV = ',TCONV,' WST = ',WST(I), &
391 ' KMIX = ',kmix,' UST = ',UST(I), &
392 ' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1), &
394 CALL wrf_error_fatal ( message )
402 IF (KPBLH(I) .NE. 1) THEN
403 !---------INTERPOLATE BETWEEN LEVELS -- jp 7/93
404 FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) - &
406 IF (FINT(I) .GT. 0.5) THEN
408 FINT(I) = FINT(I) - 0.5
410 KPBLHT = KPBLH(I) - 1
411 FINT(I) = FINT(I) + 0.5
413 PBL(I) = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) + &
416 PBLSIG(I) = FINT(I) * DSIGH(KPBLHT) + SIGMAF(KPBLHT-1) ! sigma at PBL height
420 PBLSIG(I) = SIGMAF(1)
428 ! Check for CBL and identify conv. vs. non-conv cells
429 IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3 &
430 .AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN
432 REGIME(I) = 4.0 ! FREE CONVECTIVE - ACM
437 CALL EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &
438 US, VS, TT, THETAV, DENSX, PSTAR, &
439 QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, &
440 EDDYZ, its,ite, kts,kte,ims,ime, kms,kme)
442 CALL ACM(DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, J, &
443 KLPBL, PBL, PBLSIG, MOL, UST, &
444 TST, QST, USTM, EDDYZ, DENSX, &
445 US, VS, THETA, QVS, QCS, QIS, &
446 UX, VX, THETAX, QVX, QCX, QIX, &
447 ids,ide, jds,jde, kds,kde, &
448 ims,ime, jms,jme, kms,kme, &
449 its,ite, jts,jte, kts,kte)
451 !... Calculate tendency due to PBL parameterization
455 UTNP(I,K) = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT
456 VTNP(I,K) = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT
457 TTNP(I,K) = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT
458 QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT
459 QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT
460 QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT
465 !-----------------------------------------------------------------------
466 !-----------------------------------------------------------------------
468 !-----------------------------------------------------------------------
469 !-----------------------------------------------------------------------
470 SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
471 RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
472 restart, allowed_to_read , &
473 ids, ide, jds, jde, kds, kde, &
474 ims, ime, jms, jme, kms, kme, &
475 its, ite, jts, jte, kts, kte )
476 !-----------------------------------------------------------------------
478 ! This subroutine is for preparing ACM PBL variables.
479 ! Called from module_physics_init.F
482 ! AX 3/2005 - Originally developed
483 !-----------------------------------------------------------------------
486 !-----------------------------------------------------------------------
487 !-----------------------------------------------------------------------
490 LOGICAL , INTENT(IN) :: restart , allowed_to_read
492 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
493 ims, ime, jms, jme, kms, kme, &
494 its, ite, jts, jte, kts, kte
496 INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR
498 ! REAL , DIMENSION( kms:kme ), INTENT(IN) :: SHALF
499 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
508 INTEGER :: i, j, k, itf, jtf, ktf
529 IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
540 END SUBROUTINE acminit
541 !-----------------------------------------------------------------------
542 !-----------------------------------------------------------------------
544 !-----------------------------------------------------------------------
545 !-------------------------------------------------------------------
546 SUBROUTINE EDDYX(DTPBL, ZF, ZA, MOL, PBL, UST, &
547 US, VS, TT, THETAV, DENSX, PSTAR, &
548 QVS, QCS, QIS, DSIGFI, G, RD, CPAIR, &
549 EDDYZ, its,ite, kts,kte,ims,ime,kms,kme )
552 !**********************************************************************
553 ! Two methods for computing Kz:
554 ! 1. Boundary scaling similar to Holtslag and Boville (1993)
555 ! 2. Local Kz computed as function of local Richardson # and vertical
556 ! wind shear, similar to LIU & CARROLL (1996)
558 !**********************************************************************
560 !-- DTPBL time step of the minor loop for the land-surface/pbl model
561 !-- ZF height of full sigma level
562 !-- ZA height of half sigma level
563 !-- MOL Monin-Obukhov length in 1D form
564 !-- PBL PBL height in 1D form
565 !-- UST friction velocity U* in 1D form (m/s)
569 !-- THETAV potential virtual temperature
570 !-- DENSX dry air density (kg/m^3)
571 !-- PSTAR P*=Psfc-Ptop
572 !-- QVS water vapor mixing ratio (Kg/Kg)
573 !-- QCS cloud mixing ratio (Kg/Kg)
574 !-- QIS ice mixing ratio (Kg/Kg)
575 !-- DSIGFI inverse of sigma layer delta
577 !-- RD gas constant for dry air (j/kg/k)
578 !-- CPAIR specific heat of moist air (M^2 S^-2 K^-1)
579 !-- EDDYZ eddy diffusivity KZ
580 !-----------------------------------------------------------------------
581 !-----------------------------------------------------------------------
588 INTEGER, INTENT(IN) :: its,ite, kts,kte,ims,ime,kms,kme
590 REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST
591 REAL , INTENT(IN) :: DTPBL, G, RD
592 REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGFI
593 REAL , DIMENSION( its:ite ), INTENT(IN) :: MOL, PSTAR, CPAIR
595 REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, TT, &
597 REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV
598 REAL, DIMENSION( its:ite, 0:kte ) , INTENT(IN) :: ZF
600 REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ
602 !.......Local variables
605 INTEGER :: ILX, KL, KLM, K, I
608 REAL :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ
609 REAL :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF, KZO
612 REAL, PARAMETER :: RV = 461.5
613 REAL, PARAMETER :: RC = 0.25
614 REAL, PARAMETER :: RLAM = 80.0
615 REAL, PARAMETER :: GAMH = 16.0 !15.0 ! Holtslag and Boville (1993)
616 REAL, PARAMETER :: BETAH = 5.0 ! Holtslag and Boville (1993)
617 REAL, PARAMETER :: KARMAN = 0.4
618 REAL, PARAMETER :: EDYZ0 = 0.01 ! New Min Kz
619 ! REAL, PARAMETER :: EDYZ0 = 0.1
620 !-- IMVDIF imvdif=1 for moist adiabat vertical diffusion
621 INTEGER, PARAMETER :: imvdif = 1
631 DZF = ZA(I,K+1) - ZA(I,K)
633 !--------------------------------------------------------------------------
634 IF (ZF(I,K) .LT. PBL(I)) THEN
635 ZOVL = ZF(I,K) / MOL(I)
636 IF (ZOVL .LT. 0.0) THEN
637 IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN
638 PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL)
641 ZSOL = 0.1 * PBL(I) / MOL(I)
642 PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL)
645 ELSE IF (ZOVL .LT. 1.0) THEN
646 PHIH = 1.0 + BETAH * ZOVL
652 ZFUNC = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** 2
653 EDYZ = KARMAN * WT * ZFUNC
655 !--------------------------------------------------------------------------
656 SS = ((US(I,K+1) - US(I,K)) ** 2 + (VS(I,K+1) - VS(I,K)) ** 2) &
657 / (DZF * DZF) + 1.0E-9
658 GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K))
659 RI = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS)
660 !--------------------------------------------------------------------------
661 ! Adjustment to vert diff in Moist air
663 IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+ &
664 QIS(I,K+1)) .GT. 0.01E-3) THEN
665 QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1))
666 TMEAN = 0.5 * (TT(I,K) + TT(I,K+1))
667 XLV = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6
668 ALPH = XLV * QMEAN / RD / TMEAN
669 CHI = XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN
670 RI = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) * &
671 ((CHI - ALPH) / (1.0 + CHI)))
674 !--------------------------------------------------------------------------
677 SQL = (ZK * RLAM / (RLAM + ZK)) ** 2
679 IF (RI .GE. 0.0) THEN
680 IF (ZF(I,K).LT.PBL(I).AND.ZOVL.GT.0.0) THEN
681 FH = MAX((1.-ZF(I,K)/PBL(I))**2,0.01) * PHIH **(-2)
684 FH = (MAX(1.-RI/RC,0.01))**2
686 EDDYZ(I,K) = KZO + SQRT(SS) * FH * SQL
688 EDDYZ(I,K) = KZO + SQRT(SS * (1.0 - 25.0 * RI)) * SQL
691 IF(EDYZ.GT.EDDYZ(I,K)) THEN
695 EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K))
696 EDDYZ(I,K) = MAX(KZO,EDDYZ(I,K))
698 DENSF = 0.5 * (DENSX(I,K+1) + DENSX(I,K))
700 EDDYZ(I,K) = EDDYZ(I,K) * (DENSF * G / PSTAR(I)) ** 2 * &
701 DTPBL * DSIGFI(K)*1E-6
707 EDDYZ(I,KL) = 0.0 ! EDDYZ(I,KLM) -- changed jp 3/08
711 !-----------------------------------------------------------------------
712 !-----------------------------------------------------------------------
714 !-----------------------------------------------------------------------
715 !-------------------------------------------------------------------
716 SUBROUTINE ACM (DTPBL, PSTAR, NOCONV, SIGMAF, DSIGH, DSIGHI, JX, &
717 KLPBL, PBL, PBLSIG, MOL, UST, &
718 TST, QST, USTM, EDDYZ, DENSX, &
719 US, VS, THETA, QVS, QCS, QIS, &
720 UX, VX, THETAX, QVX, QCX, QIX, &
721 ids,ide, jds,jde, kds,kde, &
722 ims,ime, jms,jme, kms,kme, &
723 its,ite, jts,jte, kts,kte)
724 !**********************************************************************
725 ! PBL model called the Asymmetric Convective Model, Version 2 (ACM2)
726 ! -- See top of module for summary and references
728 !---- REVISION HISTORY:
729 ! AX 3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM
730 ! JP and RG 8/2006 - updates
732 !**********************************************************************
734 !-- DTPBL PBL time step
735 !-- PSTAR Psurf - Ptop in cb
736 !-- NOCONV If free convection =0, no; =1, yes
737 !-- SIGMAF Sigma for full layer
738 !-- DSIGH Sigma thickness
739 !-- DSIGHI Inverse of sigma thickness
741 !-- KLPBL PBL level at K index
742 !-- PBL PBL height in m
743 !-- PBLSIG Sigma level for PBL
744 !-- MOL Monin-Obukhov length in 1D form
745 !-- UST U* in 1D form
746 !-- TST Theta* in 1D form
747 !-- QST Q* in 1D form
748 !-- USTM U* for computation of momemtum flux
749 !-- EDDYZ eddy diffusivity KZ
750 !-- DENSX dry air density (kg/m^3)
753 !-- THETA potential temperature
754 !-- QVS water vapor mixing ratio (Kg/Kg)
755 !-- QCS cloud mixing ratio (Kg/Kg)
756 !-- QIS ice mixing ratio (Kg/Kg)
759 !-- THETAX new potential temperature
760 !-- QVX new water vapor mixing ratio (Kg/Kg)
761 !-- QCX new cloud mixing ratio (Kg/Kg)
762 !-- QIX new ice mixing ratio (Kg/Kg)
763 !-----------------------------------------------------------------------
764 !-----------------------------------------------------------------------
771 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
772 ims,ime, jms,jme, kms,kme, &
773 its,ite, jts,jte, kts,kte, JX
774 INTEGER, DIMENSION( its:ite ), INTENT(IN) :: NOCONV
775 INTEGER, DIMENSION( ims:ime ), INTENT(IN) :: KLPBL
778 REAL , DIMENSION( ims:ime ), INTENT(IN) :: PBL, UST
779 REAL , INTENT(IN) :: DTPBL
780 REAL , DIMENSION( its:ite ), INTENT(IN) :: PSTAR, PBLSIG, &
783 REAL , DIMENSION( kts:kte ), INTENT(IN) :: DSIGHI, DSIGH
784 REAL , DIMENSION( 0:kte ), INTENT(IN) :: SIGMAF
785 REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT) :: EDDYZ
786 REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN) :: US,VS, THETA, &
788 REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX, THETAX, &
790 !.......Local variables
793 INTEGER, PARAMETER :: NSP = 6
795 !......ACM2 Parameters
796 ! INTEGER, PARAMETER :: IFACM = 0
798 REAL, PARAMETER :: G1000 = 9.8 * 1.0E-3
799 REAL, PARAMETER :: XX = 0.5 ! FACTOR APPLIED TO CONV MIXING TIME STEP
800 REAL, PARAMETER :: KARMAN = 0.4
803 INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L
805 INTEGER, DIMENSION( its:ite ) :: KCBL
808 REAL :: G1000I, MBMAX, HOVL, MEDDY, MBAR
809 REAL :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1
810 REAL, DIMENSION( its:ite ) :: PSTARI, FSACM, DTLIM
811 REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN
812 REAL, DIMENSION( 1:NSP, its:ite ) :: FS, BCBOTN
813 REAL, DIMENSION( kts:kte ) :: XPLUS, XMINUS
815 REAL, DIMENSION( 1:NSP,its:ite,kts:kte ) :: VCI
817 REAL, DIMENSION( kts:kte ) :: AI, BI, CI, EI !, Y
818 REAL, DIMENSION( 1:NSP,kts:kte ) :: DI, UI
820 !--Start Exicutable ----
830 !---COMPUTE ACM MIXING RATE
833 PSTARI(I) = 1.0 / PSTAR(I)
837 IF (NOCONV(I) .EQ. 1) THEN
840 !-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE
841 !--New couple ACM & EDDY-------------------------------------------------------------
842 HOVL = -PBL(I) / MOL(I)
843 FSACM(I) = 1./(1.+((KARMAN/(HOVL))**0.3333)/(0.72*KARMAN))
844 MEDDY = EDDYZ(I,1) / (DTPBL * (PBLSIG(I) - SIGMAF(1)))
845 MBAR = MEDDY * FSACM(I)
847 EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I))
850 MBMAX = AMAX1(MBMAX,MBAR)
853 MDWN(K,I) = MBAR * (PBLSIG(I) - SIGMAF(K-1)) * DSIGHI(K)
856 MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)
857 MDWN(KCBL(I)+1,I) = 0.0
859 ENDDO ! end of I loop
863 EKZ = EDDYZ(I,K) / DTPBL * DSIGHI(K)
864 DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))
869 IF (NOCONV(I) .EQ. 1) THEN
870 KCBLMX = AMAX0(KLPBL(I),KCBLMX)
871 RZ = (SIGMAF(KCBL(I)) - SIGMAF(1)) * DSIGHI(1)
872 DTLIM(I) = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))
878 VCI(1,I,K) = THETA(I,K)
879 VCI(2,I,K) = QVS(I,K)
882 ! -- Also mix cloud water and ice IF necessary
883 ! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN !!! Check other PBL models
884 VCI(5,I,K) = QCS(I,K)
885 VCI(6,I,K) = QIS(I,K)
892 FS(1,I) = -UST(I) * TST(I) * DENSX(I,1) * PSTARI(I)
893 FS(2,I) = -UST(I) * QST(I) * DENSX(I,1) * PSTARI(I)
894 FM = -USTM(I) * USTM(I) * DENSX(I,1) * PSTARI(I)
895 WSPD = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9
896 FS(3,I) = FM * US(I,1) / WSPD
897 FS(4,I) = FM * VS(I,1) / WSPD
899 FS(6,I) = 0.0 ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0
902 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
905 NLP = INT(DTPBL / DTLIM(I) + 1.0)
908 DO NL = 1,NLP ! LOOP OVER SUB TIME LOOP
910 !-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES
920 EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DSIGH(K) * DSIGHI(K-1)
921 BI(K) = 1.0 + CRANKP * MDWN(K,I) * DTS
922 AI(K) = -CRANKP * MBARKS(K,I) * DTS
925 EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DSIGHI(1 )* DTRAT
926 AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DSIGHI(2) * DTRAT
933 XPLUS(K) = EDDYZ(I,K) * DSIGHI(K) * DTRAT
934 XMINUS(K) = EDDYZ(I,K-1) * DSIGHI(K) * DTRAT
935 CI(K) = - XMINUS(K) * CRANKP
936 EI(K) = EI(K) - XPLUS(K) * CRANKP
937 BI(K) = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP
940 IF (NOCONV(I) .EQ. 1) THEN
941 BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBLSIG(I) - SIGMAF(1)) * &
942 DTS * DSIGHI(1) + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
944 BI(1) = 1.0 + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
954 !** COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION
957 DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) * &
958 VCI(L,I,K) + DSIGH(K+1) * DSIGHI(K) * &
959 MDWN(K+1,I) * VCI(L,I,K+1))
960 DI(L,K) = VCI(L,I,K) + (1.0 - CRANKP) * DELC
973 DI(L,K) = DI(L,K) - (1.0 - CRANKP) * XMINUS(K) * &
974 (VCI(L,I,K) - VCI(L,I,K-1))
978 DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) * &
979 (VCI(L,I,K+1) - VCI(L,I,K)) - &
980 (1.0 - CRANKP) * XMINUS(K) * &
981 (VCI(L,I,K) - VCI(L,I,K-1))
986 IF (NOCONV(I) .EQ. 1) THEN
988 F1 = -G1000I * (MBARKS(1,I) * &
989 (PBLSIG(I) - SIGMAF(1)) * VCI(L,I,1) - &
990 MDWN(2,I) * VCI(L,I,2) * DSIGH(2))
992 DI(L,1) = VCI(L,I,1) - G1000 * (FS(L,I) - (1.0 - CRANKP) &
993 * F1) * DSIGHI(1) * DTS
997 DI(L,1) = VCI(L,I,1) - G1000 * FS(L,I) * DSIGHI(1) * DTS
1001 DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DSIGHI(1) &
1002 * DTRAT * (VCI(L,I,2) - VCI(L,I,1))
1004 IF ( NOCONV(I) .EQ. 1 ) THEN
1005 CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)
1007 CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)
1010 !-- COMPUTE NEW THETAV AND Q
1013 VCI(L,I,K) = UI(L,K)
1018 ENDDO ! END SUB TIME LOOP
1019 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1024 THETAX(I,K) = VCI(1,I,K)
1025 QVX(I,K) = VCI(2,I,K)
1026 UX(I,K) = VCI(3,I,K)
1027 VX(I,K) = VCI(4,I,K)
1033 QCX(I,K) = VCI(5,I,K)
1034 QIX(I,K) = VCI(6,I,K)
1039 !-----------------------------------------------------------------------
1040 !-----------------------------------------------------------------------
1042 !-----------------------------------------------------------------------
1043 !-----------------------------------------------------------------------
1044 SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP)
1046 !-----------------------------------------------------------------------
1047 !-----------------------------------------------------------------------
1050 !-- Bordered band diagonal matrix solver for ACM2
1052 !-- ACM2 Matrix is in this form:
1074 ! L61 L62 L63 L64 L65 1
1075 !---------------------------------------------------------
1077 INTEGER, INTENT(IN) :: KL
1078 INTEGER, INTENT(IN) :: NSP
1079 REAL A(KL),B(KL),E(KL)
1080 REAL C(KL),D(NSP,KL),X(NSP,KL)
1083 REAL Y(NSP,KL),AIJ,SUM
1084 REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL)
1087 !-- Define Upper and Lower matrices
1102 L(I,J) = (AIJ-L(I,J-1)*E(J-1))/ &
1103 (B(J)-L(J,J-1)*E(J-1))
1109 UII(I) = B(I)-L(I,I-1)*E(I-1)
1113 !-- Forward sub for Ly=d
1119 SUM = SUM - L(I,J)*Y(V,J)
1125 !-- Back sub for Ux=y
1128 X(V,KL) = Y(V,KL)*RUII(KL)
1132 X(V,I) = (Y(V,I)-UIIP1(I)*X(V,I+1))*RUII(I)
1136 END SUBROUTINE MATRIX
1139 !-----------------------------------------------------------------------
1140 !-----------------------------------------------------------------------
1141 SUBROUTINE TRI ( L, D, U, B, X,KL,NSP)
1142 !-----------------------------------------------------------------------
1143 !-----------------------------------------------------------------------
1146 ! Solves tridiagonal system by Thomas algorithm.
1147 ! The associated tri-diagonal system is stored in 3 arrays
1150 ! U : super-diagonal
1151 ! B : right hand side function
1152 ! X : return solution from tridiagonal solver
1154 ! [ D(1) U(1) 0 0 0 ... 0 ]
1155 ! [ L(2) D(2) U(2) 0 0 ... . ]
1156 ! [ 0 L(3) D(3) U(3) 0 ... . ]
1157 ! [ . . . . . ] X(i) = B(i)
1162 !-----------------------------------------------------------------------
1168 INTEGER, INTENT(IN) :: KL
1169 INTEGER, INTENT(IN) :: NSP
1171 REAL L( KL ) ! subdiagonal
1172 REAL D(KL) ! diagonal
1173 REAL U( KL ) ! superdiagonal
1174 REAL B(NSP,KL ) ! R.H. side
1175 REAL X( NSP,KL ) ! solution
1183 ! Decomposition and forward substitution:
1186 X( V,1 ) = BET * B(V,1 )
1190 GAM(K ) = BET * U( K-1 )
1191 BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) )
1193 X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) )
1197 ! Back-substitution:
1199 DO K = KL - 1, 1, -1
1201 X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 )
1206 !-----------------------------------------------------------------------
1207 !-----------------------------------------------------------------------
1209 END MODULE module_bl_acm