! translated from NN f77 to F90 and put into WRF by Mariusz Pagowski
! NOAA/GSD & CIRA/CSU, Feb 2008
! changes to original code:
! 1. code is 1d (in z)
! 2. no advection of TKE, covariances and variances
! 3. Cranck-Nicholson replaced with the implicit scheme
! 4. removed terrain dependent grid since input in WRF in actual
! distances in z[m]
! 5. cosmetic changes to adhere to WRF standard (remove common blocks,
! intent etc)
MODULE module_bl_mynn 
(docs) 3
USE module_model_constants
, only: &
&karman, g, p1000mb, &
&cp, r_d, rcp, xlv, &
&svp1, svp2, svp3, svpt0, ep_1, ep_2
IMPLICIT NONE
! The parameters below depend on stability functions of module_sf_mynn.
REAL, PARAMETER :: cphm_st=5.0, cphm_unst=16.0, &
cphh_st=5.0, cphh_unst=16.0
REAL, PARAMETER :: xlvcp=xlv/cp, ev=xlv, rd=r_d, rk=cp/rd, &
&svp11=svp1*1.e3, p608=ep_1, ep_3=1.-ep_2
REAL, PARAMETER :: tref=300.0 ! reference temperature (K)
REAL, PARAMETER :: tv0=p608*tref, tv1=(1.+p608)*tref, gtr=g/tref
! Closure constants
REAL, PARAMETER :: &
&vk = karman, &
&pr = 0.74, &
&g1 = 0.235, &
&b1 = 24.0, &
&b2 = 15.0, &
&c2 = 0.75, &
&c3 = 0.352, &
&c4 = 0.0, &
&c5 = 0.2, &
&a1 = b1*( 1.0-3.0*g1 )/6.0, &
! &c1 = g1 -1.0/( 3.0*a1*b1**(1.0/3.0) ), &
&c1 = g1 -1.0/( 3.0*a1*2.88449914061481660), &
&a2 = a1*( g1-c1 )/( g1*pr ), &
&g2 = b2/b1*( 1.0-c3 ) +2.0*a1/b1*( 3.0-2.0*c2 )
REAL, PARAMETER :: &
&cc2 = 1.0-c2, &
&cc3 = 1.0-c3, &
&e1c = 3.0*a2*b2*cc3, &
&e2c = 9.0*a1*a2*cc2, &
&e3c = 9.0*a2*a2*cc2*( 1.0-c5 ), &
&e4c = 12.0*a1*a2*cc2, &
&e5c = 6.0*a1*a1
! Constants for length scale
REAL, PARAMETER :: qmin=0.0, zmax=1.0, cns=2.7, &
&alp1=0.23, alp2=1.0, alp3=5.0, alp4=100.0
! Constants for gravitational settling
! REAL, PARAMETER :: gno=1.e6/(1.e8)**(2./3.), gpw=5./3., qcgmin=1.e-8
REAL, PARAMETER :: gno=4.64158883361278196
REAL, PARAMETER :: gpw=5./3., qcgmin=1.e-8,qkemin=1.e-12
REAL, PARAMETER :: rr2=0.7071068, rrp=0.3989423
INTEGER :: mynn_level
CONTAINS
! **********************************************************************
! * An improved Mellor-Yamada turbulence closure model *
! * *
! * Aug/2005 M. Nakanishi (N.D.A) *
! * Modified: Dec/2005 M. Nakanishi (N.D.A) *
! * naka@nda.ac.jp *
! * *
! * Contents: *
! * 1. mym_initialize (to be called once initially) *
! * gives the closure constants and initializes the turbulent *
! * quantities. *
! * (2) mym_level2 (called in the other subroutines) *
! * calculates the stability functions at Level 2. *
! * (3) mym_length (called in the other subroutines) *
! * calculates the master length scale. *
! * 4. mym_turbulence *
! * calculates the vertical diffusivity coefficients and the *
! * production terms for the turbulent quantities. *
! * 5. mym_predict *
! * predicts the turbulent quantities at the next step. *
! * 6. mym_condensation *
! * determines the liquid water content and the cloud fraction *
! * diagnostically. *
! * *
! * call mym_initialize *
! * | *
! * |<----------------+ *
! * | | *
! * call mym_condensation | *
! * call mym_turbulence | *
! * call mym_predict | *
! * | | *
! * |-----------------+ *
! * | *
! * end *
! * *
! * Variables worthy of special mention: *
! * tref : Reference temperature *
! * thl : Liquid water potential temperature *
! * qw : Total water (water vapor+liquid water) content *
! * ql : Liquid water content *
! * vt, vq : Functions for computing the buoyancy flux *
! * *
! * If the water contents are unnecessary, e.g., in the case of *
! * ocean models, thl is the potential temperature and qw, ql, vt *
! * and vq are all zero. *
! * *
! * Grid arrangement: *
! * k+1 +---------+ *
! * | | i = 1 - nx *
! * (k) | * | j = 1 - ny *
! * | | k = 1 - nz *
! * k +---------+ *
! * i (i) i+1 *
! * *
! * All the predicted variables are defined at the center (*) of *
! * the grid boxes. The diffusivity coefficients are, however, *
! * defined on the walls of the grid boxes. *
! * # Upper boundary values are given at k=nz. *
! * *
! * References: *
! * 1. Nakanishi, M., 2001: *
! * Boundary-Layer Meteor., 99, 349-378. *
! * 2. Nakanishi, M. and H. Niino, 2004: *
! * Boundary-Layer Meteor., 112, 1-31. *
! * 3. Nakanishi, M. and H. Niino, 2006: *
! * Boundary-Layer Meteor., (in press). *
! **********************************************************************
!
!
! SUBROUTINE mym_initialize:
!
! Input variables:
! iniflag : <>0; turbulent quantities will be initialized
! = 0; turbulent quantities have been already
! given, i.e., they will not be initialized
! mx, my : Maximum numbers of grid boxes
! in the x and y directions, respectively
! nx, ny, nz : Numbers of the actual grid boxes
! in the x, y and z directions, respectively
! tref : Reference temperature (K)
! dz(nz) : Vertical grid spacings (m)
! # dz(nz)=dz(nz-1)
! zw(nz+1) : Heights of the walls of the grid boxes (m)
! # zw(1)=0.0 and zw(k)=zw(k-1)+dz(k-1)
! h(mx,ny) : G^(1/2) in the terrain-following coordinate
! # h=1-zg/zt, where zg is the height of the
! terrain and zt the top of the model domain
! pi0(mx,my,nz) : Exner function at zw*h+zg (J/kg K)
! defined by c_p*( p_basic/1000hPa )^kappa
! This is usually computed by integrating
! d(pi0)/dz = -h*g/tref.
! rmo(mx,ny) : Inverse of the Obukhov length (m^(-1))
! flt, flq(mx,ny) : Turbulent fluxes of sensible and latent heat,
! respectively, e.g., flt=-u_*Theta_* (K m/s)
!! flt - liquid water potential temperature surface flux
!! flq - total water flux surface flux
! ust(mx,ny) : Friction velocity (m/s)
! pmz(mx,ny) : phi_m-zeta at z1*h+z0, where z1 (=0.5*dz(1))
! is the first grid point above the surafce, z0
! the roughness length and zeta=(z1*h+z0)*rmo
! phh(mx,ny) : phi_h at z1*h+z0
! u, v(mx,my,nz): Components of the horizontal wind (m/s)
! thl(mx,my,nz) : Liquid water potential temperature
! (K)
! qw(mx,my,nz) : Total water content Q_w (kg/kg)
!
! Output variables:
! ql(mx,my,nz) : Liquid water content (kg/kg)
! v?(mx,my,nz) : Functions for computing the buoyancy flux
! qke(mx,my,nz) : Twice the turbulent kineti! energy q^2
! (m^2/s^2)
! tsq(mx,my,nz) : Variance of Theta_l (K^2)
! qsq(mx,my,nz) : Variance of Q_w
! cov(mx,my,nz) : Covariance of Theta_l and Q_w (K)
! el(mx,my,nz) : Master length scale L (m)
! defined on the walls of the grid boxes
! bsh : no longer used
! via common : Closure constants
!
! Work arrays: see subroutine mym_level2
! pd?(mx,my,nz) : Half of the production terms at Level 2
! defined on the walls of the grid boxes
! qkw(mx,my,nz) : q on the walls of the grid boxes (m/s)
!
! # As to dtl, ...gh, see subroutine mym_turbulence.
!
SUBROUTINE mym_initialize (docs) ( kts,kte,& 1
& dz, zw, &
& u, v, thl, qw, &
! & ust, rmo, pmz, phh, flt, flq,&
& ust, rmo, &
& Qke, Tsq, Qsq, Cov)
!
INTEGER, INTENT(IN) :: kts,kte
! REAL, INTENT(IN) :: ust, rmo, pmz, phh, flt, flq
REAL, INTENT(IN) :: ust, rmo
REAL, DIMENSION(kts:kte), INTENT(in) :: dz
REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw
REAL, DIMENSION(kts:kte), INTENT(out) :: qke,tsq,qsq,cov
REAL, DIMENSION(kts:kte) :: &
&ql,el,pdk,pdt,pdq,pdc,dtl,dqw,dtv,&
&gm,gh,sm,sh,qkw,vt,vq
INTEGER :: k,l,lmax
REAL :: phm,vkz,elq,elv,b1l,b2l,pmz=1.,phh=1.,flt=0.,flq=0.,tmpq
! ** At first ql, vt and vq are set to zero. **
DO k = kts,kte
ql(k) = 0.0
vt(k) = 0.0
vq(k) = 0.0
END DO
!
CALL mym_level2 ( kts,kte,&
& dz, &
& u, v, thl, qw, &
& ql, vt, vq, &
& dtl, dqw, dtv, gm, gh, sm, sh )
!
! ** Preliminary setting **
el (kts) = 0.0
qke(kts) = ust**2 * ( b1*pmz )**(2.0/3.0)
!
phm = phh*b2 / ( b1*pmz )**(1.0/3.0)
tsq(kts) = phm*( flt/ust )**2
qsq(kts) = phm*( flq/ust )**2
cov(kts) = phm*( flt/ust )*( flq/ust )
!
DO k = kts+1,kte
vkz = vk*zw(k)
el (k) = vkz/( 1.0 + vkz/100.0 )
qke(k) = 0.0
!
tsq(k) = 0.0
qsq(k) = 0.0
cov(k) = 0.0
END DO
!
! ** Initialization with an iterative manner **
! ** lmax is the iteration count. This is arbitrary. **
lmax = 5 !!kte +3
!
DO l = 1,lmax
!
CALL mym_length ( kts,kte,&
& dz, zw, &
& rmo, flt, flq, &
& vt, vq, &
& qke, &
& dtv, &
& el, &
& qkw)
!
DO k = kts+1,kte
elq = el(k)*qkw(k)
pdk(k) = elq*( sm(k)*gm (k)+&
&sh(k)*gh (k) )
pdt(k) = elq* sh(k)*dtl(k)**2
pdq(k) = elq* sh(k)*dqw(k)**2
pdc(k) = elq* sh(k)*dtl(k)*dqw(k)
END DO
!
! ** Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 ) **
vkz = vk*0.5*dz(kts)
!
elv = 0.5*( el(kts+1)+el(kts) ) / vkz
qke(kts) = ust**2 * ( b1*pmz*elv )**(2.0/3.0)
!
phm = phh*b2 / ( b1*pmz/elv**2 )**(1.0/3.0)
tsq(kts) = phm*( flt/ust )**2
qsq(kts) = phm*( flq/ust )**2
cov(kts) = phm*( flt/ust )*( flq/ust )
!
DO k = kts+1,kte-1
b1l = b1*0.25*( el(k+1)+el(k) )
tmpq=MAX(b1l*( pdk(k+1)+pdk(k) ),qkemin)
! PRINT *,'tmpqqqqq',tmpq,pdk(k+1),pdk(k)
qke(k) = tmpq**(2.0/3.0)
!
IF ( qke(k) .LE. 0.0 ) THEN
b2l = 0.0
ELSE
b2l = b2*( b1l/b1 ) / SQRT( qke(k) )
END IF
!
tsq(k) = b2l*( pdt(k+1)+pdt(k) )
qsq(k) = b2l*( pdq(k+1)+pdq(k) )
cov(k) = b2l*( pdc(k+1)+pdc(k) )
END DO
!
END DO
!! qke(kts)=qke(kts+1)
!! tsq(kts)=tsq(kts+1)
!! qsq(kts)=qsq(kts+1)
!! cov(kts)=cov(kts+1)
qke(kte)=qke(kte-1)
tsq(kte)=tsq(kte-1)
qsq(kte)=qsq(kte-1)
cov(kte)=cov(kte-1)
!
! RETURN
END SUBROUTINE mym_initialize
!
! SUBROUTINE mym_level2:
!
! Input variables: see subroutine mym_initialize
!
! Output variables:
! dtl(mx,my,nz) : Vertical gradient of Theta_l (K/m)
! dqw(mx,my,nz) : Vertical gradient of Q_w
! dtv(mx,my,nz) : Vertical gradient of Theta_V (K/m)
! gm (mx,my,nz) : G_M divided by L^2/q^2 (s^(-2))
! gh (mx,my,nz) : G_H divided by L^2/q^2 (s^(-2))
! sm (mx,my,nz) : Stability function for momentum, at Level 2
! sh (mx,my,nz) : Stability function for heat, at Level 2
!
! These are defined on the walls of the grid boxes.
!
SUBROUTINE mym_level2 (docs) (kts,kte,& 2
& dz, &
& u, v, thl, qw, &
& ql, vt, vq, &
& dtl, dqw, dtv, gm, gh, sm, sh )
!
INTEGER, INTENT(IN) :: kts,kte
REAL, DIMENSION(kts:kte), INTENT(in) :: dz
REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,ql,vt,vq
REAL, DIMENSION(kts:kte), INTENT(out) :: &
&dtl,dqw,dtv,gm,gh,sm,sh
INTEGER :: k
REAL :: rfc,f1,f2,rf1,rf2,smc,shc,&
&ri1,ri2,ri3,ri4,duz,dtz,dqz,vtt,vqq,dtq,dzk,afk,abk,ri,rf
! ev = 2.5e6
! tv0 = 0.61*tref
! tv1 = 1.61*tref
! gtr = 9.81/tref
!
rfc = g1/( g1+g2 )
f1 = b1*( g1-c1 ) +3.0*a2*( 1.0 -c2 )*( 1.0-c5 ) &
& +2.0*a1*( 3.0-2.0*c2 )
f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 )
rf1 = b1*( g1-c1 )/f1
rf2 = b1* g1 /f2
smc = a1 /a2* f1/f2
shc = 3.0*a2*( g1+g2 )
!
ri1 = 0.5/smc
ri2 = rf1*smc
ri3 = 4.0*rf2*smc -2.0*ri2
ri4 = ri2**2
!
DO k = kts+1,kte
dzk = 0.5 *( dz(k)+dz(k-1) )
afk = dz(k)/( dz(k)+dz(k-1) )
abk = 1.0 -afk
duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2
duz = duz /dzk**2
dtz = ( thl(k)-thl(k-1) )/( dzk )
dqz = ( qw(k)-qw(k-1) )/( dzk )
!
vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
vqq = tv0 +vq(k)*abk +vq(k-1)*afk
dtq = vtt*dtz +vqq*dqz
!
dtl(k) = dtz
dqw(k) = dqz
dtv(k) = dtq
!? dtv(i,j,k) = dtz +tv0*dqz
!? : +( ev/pi0(i,j,k)-tv1 )
!? : *( ql(i,j,k)-ql(i,j,k-1) )/( dzk*h(i,j) )
!
gm (k) = duz
gh (k) = -dtq*gtr
!
! ** Gradient Richardson number **
ri = -gh(k)/MAX( duz, 1.0e-10 )
! ** Flux Richardson number **
rf = MIN( ri1*( ri+ri2-SQRT(ri**2-ri3*ri+ri4) ), rfc )
!
sh (k) = shc*( rfc-rf )/( 1.0-rf )
sm (k) = smc*( rf1-rf )/( rf2-rf ) * sh(k)
END DO
!
RETURN
END SUBROUTINE mym_level2
!
! SUBROUTINE mym_length:
!
! Input variables: see subroutine mym_initialize
!
! Output variables: see subroutine mym_initialize
!
! Work arrays:
! elt(mx,ny) : Length scale depending on the PBL depth (m)
! vsc(mx,ny) : Velocity scale q_c (m/s)
! at first, used for computing elt
!
SUBROUTINE mym_length (docs) ( kts,kte,& 2
& dz, zw, &
& rmo, flt, flq, &
& vt, vq, &
& qke, &
& dtv, &
& el, &
& qkw)
INTEGER, INTENT(IN) :: kts,kte
REAL, DIMENSION(kts:kte), INTENT(in) :: dz
REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
REAL, INTENT(in) :: rmo,flt,flq
REAL, DIMENSION(kts:kte), INTENT(IN) :: qke,vt,vq
REAL, DIMENSION(kts:kte), INTENT(out) :: &
&qkw, el
REAL, DIMENSION(kts:kte), INTENT(in) :: dtv
REAL :: elt,vsc
INTEGER :: i,j,k
REAL :: afk,abk,zwk,dzk,qdz,vflx,bv,elb,els,elf
! tv0 = 0.61*tref
! gtr = 9.81/tref
!
DO k = kts+1,kte
afk = dz(k)/( dz(k)+dz(k-1) )
abk = 1.0 -afk
qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*&
&afk,1.0e-10))
END DO
!
elt = 1.0e-5
vsc = 1.0e-5
!
! ** Strictly, zwk*h(i,j) -> ( zwk*h(i,j)+z0 ) **
DO k = kts+1,kte-1
zwk = zw(k)
dzk = 0.5*( dz(k)+dz(k-1) )
qdz = MAX( qkw(k)-qmin, 0.0 )*dzk
elt = elt +qdz*zwk
vsc = vsc +qdz
END DO
!
vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq
elt = alp1*elt/vsc
vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**(1.0/3.0)
!
! ** Strictly, el(i,j,1) is not zero. **
el(kts) = 0.0
!
DO k = kts+1,kte
zwk = zw(k)
! ** Length scale limited by the buoyancy effect **
IF ( dtv(k) .GT. 0.0 ) THEN
bv = SQRT( gtr*dtv(k) )
elb = alp2*qkw(k) / bv &
& *( 1.0 + alp3/alp2*&
&SQRT( vsc/( bv*elt ) ) )
elf=elb
elf = alp2 * qkw(k)/bv
ELSE
elb = 1.0e10
elf =elb
END IF
!
! ** Length scale in the surface layer **
IF ( rmo .GT. 0.0 ) THEN
els = vk*zwk &
& /( 1.0 + cns*MIN( zwk*rmo, zmax ) )
ELSE
els = vk*zwk &
& *( 1.0 - alp4* zwk*rmo )**0.2
END IF
!
el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf)
!! el(k) = elb/( elb/elt+elb/els+1.0 )
END DO
!
RETURN
END SUBROUTINE mym_length
!
! SUBROUTINE mym_turbulence:
!
! Input variables: see subroutine mym_initialize
! levflag : <>3; Level 2.5
! = 3; Level 3
!
! # ql, vt, vq, qke, tsq, qsq and cov are changed to input variables.
!
! Output variables: see subroutine mym_initialize
! dfm(mx,my,nz) : Diffusivity coefficient for momentum,
! divided by dz (not dz*h(i,j)) (m/s)
! dfh(mx,my,nz) : Diffusivity coefficient for heat,
! divided by dz (not dz*h(i,j)) (m/s)
! dfq(mx,my,nz) : Diffusivity coefficient for q^2,
! divided by dz (not dz*h(i,j)) (m/s)
! tcd(mx,my,nz) : Countergradient diffusion term for Theta_l
! (K/s)
! qcd(mx,my,nz) : Countergradient diffusion term for Q_w
! (kg/kg s)
! pd?(mx,my,nz) : Half of the production terms
!
! Only tcd and qcd are defined at the center of the grid boxes
!
! # DO NOT forget that tcd and qcd are added on the right-hand side
! of the equations for Theta_l and Q_w, respectively.
!
! Work arrays: see subroutine mym_initialize and level2
!
! # dtl, dqw, dtv, gm and gh are allowed to share storage units with
! dfm, dfh, dfq, tcd and qcd, respectively, for saving memory.
!
SUBROUTINE mym_turbulence (docs) ( kts,kte,& 1
& levflag, &
& dz, zw, &
& u, v, thl, ql, qw, &
& qke, tsq, qsq, cov, &
& vt, vq,&
& rmo, flt, flq, &
& El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc)
!
INTEGER, INTENT(IN) :: kts,kte
INTEGER, INTENT(IN) :: levflag
REAL, DIMENSION(kts:kte), INTENT(in) :: dz
REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
REAL, INTENT(in) :: rmo,flt,flq
REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,&
&ql,vt,vq,qke,tsq,qsq,cov
REAL, DIMENSION(kts:kte), INTENT(out) :: dfm,dfh,dfq,&
&pdk,pdt,pdq,pdc,tcd,qcd,el
REAL, DIMENSION(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh
INTEGER :: k
! REAL :: cc2,cc3,e1c,e2c,e3c,e4c,e5c
REAL :: e6c,dzk,afk,abk,vtt,vqq,&
&cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh
DOUBLE PRECISION q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel
DOUBLE PRECISION q3sq, t3sq, r3sq, c3sq, dlsq, qdiv
DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden
!
! tv0 = 0.61*tref
! gtr = 9.81/tref
!
! cc2 = 1.0-c2
! cc3 = 1.0-c3
! e1c = 3.0*a2*b2*cc3
! e2c = 9.0*a1*a2*cc2
! e3c = 9.0*a2*a2*cc2*( 1.0-c5 )
! e4c = 12.0*a1*a2*cc2
! e5c = 6.0*a1*a1
!
CALL mym_level2 (kts,kte,&
& dz, &
& u, v, thl, qw, &
& ql, vt, vq, &
& dtl, dqw, dtv, gm, gh, sm, sh )
!
CALL mym_length (kts,kte, &
& dz, zw, &
& rmo, flt, flq, &
& vt, vq, &
& qke, &
& dtv, &
& el, &
& qkw)
!
DO k = kts+1,kte
dzk = 0.5 *( dz(k)+dz(k-1) )
afk = dz(k)/( dz(k)+dz(k-1) )
abk = 1.0 -afk
elsq = el (k)**2
q2sq = b1*elsq*( sm(k)*gm(k)+sh(k)*gh(k) )
q3sq = qkw(k)**2
!
! Modified: Dec/22/2005, from here, (dlsq -> elsq)
gmel = gm (k)*elsq
ghel = gh (k)*elsq
! Modified: Dec/22/2005, up to here
!
! ** Since qkw is set to more than 0.0, q3sq > 0.0. **
IF ( q3sq .LT. q2sq ) THEN
qdiv = SQRT( q3sq/q2sq )
sm(k) = sm(k) * qdiv
sh(k) = sh(k) * qdiv
!
e1 = q3sq - e1c*ghel * qdiv**2
e2 = q3sq - e2c*ghel * qdiv**2
e3 = e1 + e3c*ghel * qdiv**2
e4 = e1 - e4c*ghel * qdiv**2
eden = e2*e4 + e3*e5c*gmel * qdiv**2
eden = MAX( eden, 1.0d-20 )
ELSE
e1 = q3sq - e1c*ghel
e2 = q3sq - e2c*ghel
e3 = e1 + e3c*ghel
e4 = e1 - e4c*ghel
eden = e2*e4 + e3*e5c*gmel
eden = MAX( eden, 1.0d-20 )
!
qdiv = 1.0
sm(k) = q3sq*a1*( e3-3.0*c1*e4 )/eden
sh(k) = q3sq*a2*( e2+3.0*c1*e5c*gmel )/eden
END IF
!
! ** Level 3 : start **
IF ( levflag .EQ. 3 ) THEN
t2sq = qdiv*b2*elsq*sh(k)*dtl(k)**2
r2sq = qdiv*b2*elsq*sh(k)*dqw(k)**2
c2sq = qdiv*b2*elsq*sh(k)*dtl(k)*dqw(k)
t3sq = MAX( tsq(k)*abk+tsq(k-1)*afk, 0.0 )
r3sq = MAX( qsq(k)*abk+qsq(k-1)*afk, 0.0 )
c3sq = cov(k)*abk+cov(k-1)*afk
!
! Modified: Dec/22/2005, from here
c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
!
vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
vqq = tv0 +vq(k)*abk +vq(k-1)*afk
t2sq = vtt*t2sq +vqq*c2sq
r2sq = vtt*c2sq +vqq*r2sq
c2sq = MAX( vtt*t2sq+vqq*r2sq, 0.0d0 )
t3sq = vtt*t3sq +vqq*c3sq
r3sq = vtt*c3sq +vqq*r3sq
c3sq = MAX( vtt*t3sq+vqq*r3sq, 0.0d0 )
!
cw25 = e1*( e2 + 3.0*c1*e5c*gmel*qdiv**2 )/( 3.0*eden )
!
! ** Limitation on q, instead of L/q **
dlsq = elsq
IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
!
! ** Limitation on c3sq (0.12 =< cw =< 0.76) **
e2 = q3sq - e2c*ghel * qdiv**2
e3 = q3sq + e3c*ghel * qdiv**2
e4 = q3sq - e4c*ghel * qdiv**2
eden = e2*e4 + e3 *e5c*gmel * qdiv**2
!
wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2 &
& *( e2*e4c - e3c*e5c*gmel * qdiv**2 )
!
IF ( wden .NE. 0.0 ) THEN
clow = q3sq*( 0.12-cw25 )*eden/wden
cupp = q3sq*( 0.76-cw25 )*eden/wden
!
IF ( wden .GT. 0.0 ) THEN
c3sq = MIN( MAX( c3sq, c2sq+clow ), c2sq+cupp )
ELSE
c3sq = MAX( MIN( c3sq, c2sq+clow ), c2sq+cupp )
END IF
END IF
!
e1 = e2 + e5c*gmel * qdiv**2
eden = MAX( eden, 1.0d-20 )
! Modified: Dec/22/2005, up to here
!
e6c = 3.0*a2*cc3*gtr * dlsq/elsq
!
! ** for Gamma_theta **
!! enum = qdiv*e6c*( t3sq-t2sq )
IF ( t2sq .GE. 0.0 ) THEN
enum = MAX( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
ELSE
enum = MIN( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
ENDIF
gamt =-e1 *enum /eden
!
! ** for Gamma_q **
!! enum = qdiv*e6c*( r3sq-r2sq )
IF ( r2sq .GE. 0.0 ) THEN
enum = MAX( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
ELSE
enum = MIN( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
ENDIF
gamq =-e1 *enum /eden
!
! ** for Sm' and Sh'd(Theta_V)/dz **
!! enum = qdiv*e6c*( c3sq-c2sq )
enum = MAX( qdiv*e6c*( c3sq-c2sq ), 0.0d0)
smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c+e4c)*a1/a2
gamv = e1 *enum*gtr/eden
!
sm(k) = sm(k) +smd
!
! ** For elh (see below), qdiv at Level 3 is reset to 1.0. **
qdiv = 1.0
! ** Level 3 : end **
!
ELSE
! ** At Level 2.5, qdiv is not reset. **
gamt = 0.0
gamq = 0.0
gamv = 0.0
END IF
!
elq = el(k)*qkw(k)
elh = elq*qdiv
!
pdk(k) = elq*( sm(k)*gm (k) &
& +sh(k)*gh (k)+gamv )
pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k)
pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k)
pdc(k) = elh*( sh(k)*dtl(k)+gamt )&
&*dqw(k)*0.5 &
&+elh*( sh(k)*dqw(k)+gamq )*dtl(k)*0.5
!
tcd(k) = elq*gamt
qcd(k) = elq*gamq
!
dfm(k) = elq*sm (k) / dzk
dfh(k) = elq*sh (k) / dzk
! Modified: Dec/22/2005, from here
! ** In sub.mym_predict, dfq for the TKE and scalar variance **
! ** are set to 3.0*dfm and 1.0*dfm, respectively. **
dfq(k) = dfm(k)
! Modified: Dec/22/2005, up to here
END DO
!
dfm(kts) = 0.0
dfh(kts) = 0.0
dfq(kts) = 0.0
tcd(kts) = 0.0
qcd(kts) = 0.0
tcd(kte) = 0.0
qcd(kte) = 0.0
!
DO k = kts,kte-1
dzk = dz(k)
tcd(k) = ( tcd(k+1)-tcd(k) )/( dzk )
qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk )
END DO
!
RETURN
END SUBROUTINE mym_turbulence
!
! SUBROUTINE mym_predict:
!
! Input variables: see subroutine mym_initialize and turbulence
! qke(mx,my,nz) : qke at (n)th time level
! tsq, ...cov : ditto
!
! Output variables:
! qke(mx,my,nz) : qke at (n+1)th time level
! tsq, ...cov : ditto
!
! Work arrays:
! qkw(mx,my,nz) : q at the center of the grid boxes (m/s)
! bp (mx,my,nz) : = 1/2*F, see below
! rp (mx,my,nz) : = P-1/2*F*Q, see below
!
! # The equation for a turbulent quantity Q can be expressed as
! dQ/dt + Ah + Av = Dh + Dv + P - F*Q, (1)
! where A is the advection, D the diffusion, P the production,
! F*Q the dissipation and h and v denote horizontal and vertical,
! respectively. If Q is q^2, F is 2q/B_1L.
! Using the Crank-Nicholson scheme for Av, Dv and F*Q, a finite
! difference equation is written as
! Q{n+1} - Q{n} = dt *( Dh{n} - Ah{n} + P{n} )
! + dt/2*( Dv{n} - Av{n} - F*Q{n} )
! + dt/2*( Dv{n+1} - Av{n+1} - F*Q{n+1} ), (2)
! where n denotes the time level.
! When the advection and diffusion terms are discretized as
! dt/2*( Dv - Av ) = a(k)Q(k+1) - b(k)Q(k) + c(k)Q(k-1), (3)
! Eq.(2) can be rewritten as
! - a(k)Q(k+1) + [ 1 + b(k) + dt/2*F ]Q(k) - c(k)Q(k-1)
! = Q{n} + dt *( Dh{n} - Ah{n} + P{n} )
! + dt/2*( Dv{n} - Av{n} - F*Q{n} ), (4)
! where Q on the left-hand side is at (n+1)th time level.
!
! In this subroutine, a(k), b(k) and c(k) are obtained from
! subprogram coefvu and are passed to subprogram tinteg via
! common. 1/2*F and P-1/2*F*Q are stored in bp and rp,
! respectively. Subprogram tinteg solves Eq.(4).
!
! Modify this subroutine according to your numerical integration
! scheme (program).
!
SUBROUTINE mym_predict (docs) (kts,kte,& 1
& levflag, &
& delt,&
& dz, &
& ust, flt, flq, pmz, phh, &
& el, dfq, &
& pdk, pdt, pdq, pdc,&
& qke, tsq, qsq, cov)
INTEGER, INTENT(IN) :: kts,kte
INTEGER, INTENT(IN) :: levflag
REAL, INTENT(IN) :: delt
REAL, DIMENSION(kts:kte), INTENT(IN) :: dz, dfq,el
REAL, DIMENSION(kts:kte), INTENT(INOUT) :: pdk, pdt, pdq, pdc
REAL, INTENT(IN) :: flt, flq, ust, pmz, phh
REAL, DIMENSION(kts:kte), INTENT(INOUT) :: qke,tsq, qsq, cov
INTEGER :: k,nz
REAL, DIMENSION(kts:kte) :: qkw, bp, rp, df3q
REAL :: vkz,pdk1,phm,pdt1,pdq1,pdc1,b1l,b2l
REAL, DIMENSION(kts:kte) :: dtz
REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
nz=kte-kts+1
! ** Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 ) **
vkz = vk*0.5*dz(kts)
!
! Modified: Dec/22/2005, from here
! ** dfq for the TKE is 3.0*dfm. **
! CALL coefvu ( dfq, 3.0 ) ! make change here
! Modified: Dec/22/2005, up to here
!
DO k = kts,kte
!! qke(k) = MAX(qke(k), 0.0)
qkw(k) = SQRT( MAX( qke(k), 0.0 ) )
df3q(k)=3.*dfq(k)
dtz(k)=delt/dz(k)
END DO
!
pdk1 = 2.0*ust**3*pmz/( vkz )
phm = 2.0/ust *phh/( vkz )
pdt1 = phm*flt**2
pdq1 = phm*flq**2
pdc1 = phm*flt*flq
!
! ** pdk(i,j,1)+pdk(i,j,2) corresponds to pdk1. **
pdk(kts) = pdk1 -pdk(kts+1)
!! pdt(kts) = pdt1 -pdt(kts+1)
!! pdq(kts) = pdq1 -pdq(kts+1)
!! pdc(kts) = pdc1 -pdc(kts+1)
pdt(kts) = pdt(kts+1)
pdq(kts) = pdq(kts+1)
pdc(kts) = pdc(kts+1)
!
! ** Prediction of twice the turbulent kinetic energy **
!! DO k = kts+1,kte-1
DO k = kts,kte-1
b1l = b1*0.5*( el(k+1)+el(k) )
bp(k) = 2.*qkw(k) / b1l
rp(k) = pdk(k+1) + pdk(k)
END DO
!! a(1)=0.
!! b(1)=1.
!! c(1)=-1.
!! d(1)=0.
! Since df3q(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*df3q(k+1)+bp(k)*delt.
DO k=kts,kte-1
a(k-kts+1)=-dtz(k)*df3q(k)
b(k-kts+1)=1.+dtz(k)*(df3q(k)+df3q(k+1))+bp(k)*delt
c(k-kts+1)=-dtz(k)*df3q(k+1)
d(k-kts+1)=rp(k)*delt + qke(k)
ENDDO
!! DO k=kts+1,kte-1
!! a(k-kts+1)=-dtz(k)*df3q(k)
!! b(k-kts+1)=1.+dtz(k)*(df3q(k)+df3q(k+1))
!! c(k-kts+1)=-dtz(k)*df3q(k+1)
!! d(k-kts+1)=rp(k)*delt + qke(k) - qke(k)*bp(k)*delt
!! ENDDO
a(nz)=-1. !0.
b(nz)=1.
c(nz)=0.
d(nz)=0.
CALL tridiag(nz,a,b,c,d)
DO k=kts,kte
qke(k)=d(k-kts+1)
ENDDO
IF ( levflag .EQ. 3 ) THEN
!
! Modified: Dec/22/2005, from here
! ** dfq for the scalar variance is 1.0*dfm. **
! CALL coefvu ( dfq, 1.0 ) make change here
! Modified: Dec/22/2005, up to here
!
! ** Prediction of the temperature variance **
!! DO k = kts+1,kte-1
DO k = kts,kte-1
b2l = b2*0.5*( el(k+1)+el(k) )
bp(k) = 2.*qkw(k) / b2l
rp(k) = pdt(k+1) + pdt(k)
END DO
!zero gradient for tsq at bottom and top
!! a(1)=0.
!! b(1)=1.
!! c(1)=-1.
!! d(1)=0.
! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
DO k=kts,kte-1
a(k-kts+1)=-dtz(k)*dfq(k)
b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
c(k-kts+1)=-dtz(k)*dfq(k+1)
d(k-kts+1)=rp(k)*delt + tsq(k)
ENDDO
!! DO k=kts+1,kte-1
!! a(k-kts+1)=-dtz(k)*dfq(k)
!! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
!! c(k-kts+1)=-dtz(k)*dfq(k+1)
!! d(k-kts+1)=rp(k)*delt + tsq(k) - tsq(k)*bp(k)*delt
!! ENDDO
a(nz)=-1. !0.
b(nz)=1.
c(nz)=0.
d(nz)=0.
CALL tridiag(nz,a,b,c,d)
DO k=kts,kte
tsq(k)=d(k-kts+1)
ENDDO
! ** Prediction of the moisture variance **
!! DO k = kts+1,kte-1
DO k = kts,kte-1
b2l = b2*0.5*( el(k+1)+el(k) )
bp(k) = 2.*qkw(k) / b2l
rp(k) = pdq(k+1) +pdq(k)
END DO
!zero gradient for qsq at bottom and top
!! a(1)=0.
!! b(1)=1.
!! c(1)=-1.
!! d(1)=0.
! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
DO k=kts,kte-1
a(k-kts+1)=-dtz(k)*dfq(k)
b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
c(k-kts+1)=-dtz(k)*dfq(k+1)
d(k-kts+1)=rp(k)*delt + qsq(k)
ENDDO
!! DO k=kts+1,kte-1
!! a(k-kts+1)=-dtz(k)*dfq(k)
!! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
!! c(k-kts+1)=-dtz(k)*dfq(k+1)
!! d(k-kts+1)=rp(k)*delt + qsq(k) -qsq(k)*bp(k)*delt
!! ENDDO
a(nz)=-1. !0.
b(nz)=1.
c(nz)=0.
d(nz)=0.
CALL tridiag(nz,a,b,c,d)
DO k=kts,kte
qsq(k)=d(k-kts+1)
ENDDO
! ** Prediction of the temperature-moisture covariance **
!! DO k = kts+1,kte-1
DO k = kts,kte-1
b2l = b2*0.5*( el(k+1)+el(k) )
bp(k) = 2.*qkw(k) / b2l
rp(k) = pdc(k+1) + pdc(k)
END DO
!zero gradient for tqcov at bottom and top
!! a(1)=0.
!! b(1)=1.
!! c(1)=-1.
!! d(1)=0.
! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
DO k=kts,kte-1
a(k-kts+1)=-dtz(k)*dfq(k)
b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
c(k-kts+1)=-dtz(k)*dfq(k+1)
d(k-kts+1)=rp(k)*delt + cov(k)
ENDDO
!! DO k=kts+1,kte-1
!! a(k-kts+1)=-dtz(k)*dfq(k)
!! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
!! c(k-kts+1)=-dtz(k)*dfq(k+1)
!! d(k-kts+1)=rp(k)*delt + cov(k) - cov(k)*bp(k)*delt
!! ENDDO
a(nz)=-1. !0.
b(nz)=1.
c(nz)=0.
d(nz)=0.
CALL tridiag(nz,a,b,c,d)
DO k=kts,kte
cov(k)=d(k-kts+1)
ENDDO
ELSE
!! DO k = kts+1,kte-1
DO k = kts,kte-1
IF ( qkw(k) .LE. 0.0 ) THEN
b2l = 0.0
ELSE
b2l = b2*0.25*( el(k+1)+el(k) )/qkw(k)
END IF
!
tsq(k) = b2l*( pdt(k+1)+pdt(k) )
qsq(k) = b2l*( pdq(k+1)+pdq(k) )
cov(k) = b2l*( pdc(k+1)+pdc(k) )
END DO
!! tsq(kts)=tsq(kts+1)
!! qsq(kts)=qsq(kts+1)
!! cov(kts)=cov(kts+1)
tsq(kte)=tsq(kte-1)
qsq(kte)=qsq(kte-1)
cov(kte)=cov(kte-1)
END IF
END SUBROUTINE mym_predict
! SUBROUTINE mym_condensation:
!
! Input variables: see subroutine mym_initialize and turbulence
! pi (mx,my,nz) : Perturbation of the Exner function (J/kg K)
! defined on the walls of the grid boxes
! This is usually computed by integrating
! d(pi)/dz = h*g*tv/tref**2
! from the upper boundary, where tv is the
! virtual potential temperature minus tref.
!
! Output variables: see subroutine mym_initialize
! cld(mx,my,nz) : Cloud fraction
!
! Work arrays:
! qmq(mx,my,nz) : Q_w-Q_{sl}, where Q_{sl} is the saturation
! specific humidity at T=Tl
! alp(mx,my,nz) : Functions in the condensation process
! bet(mx,my,nz) : ditto
! sgm(mx,my,nz) : Combined standard deviation sigma_s
! multiplied by 2/alp
!
! # qmq, alp, bet and sgm are allowed to share storage units with
! any four of other work arrays for saving memory.
!
! # Results are sensitive particularly to values of cp and rd.
! Set these values to those adopted by you.
!
SUBROUTINE mym_condensation (docs) (kts,kte, & 2
& levflag, &
& dz, &
& thl, qw, &
& p,exner, &
& tsq, qsq, cov, &
& Vt, Vq)
INTEGER, INTENT(IN) :: kts,kte
INTEGER, INTENT(in) :: levflag
REAL, DIMENSION(kts:kte), INTENT(IN) :: dz
REAL, DIMENSION(kts:kte), INTENT(IN) :: p,exner, thl, qw, &
&tsq, qsq, cov
REAL, DIMENSION(kts:kte), INTENT(OUT) :: vt,vq
REAL, DIMENSION(kts:kte) :: qmq,alp,bet,sgm,ql,cld
DOUBLE PRECISION :: t3sq, r3sq, c3sq
!
REAL :: p2a,t,esl,qsl,dqsl,q1,cld0,eq1,qll,&
&q2p,pt,rac,qt
INTEGER :: i,j,k
REAL :: erf
! Note: kte needs to be larger than kts, i.e., kte >= kts+1.
DO k = kts,kte-1
p2a = exner(k)
t = thl(k)*p2a
!x if ( ct .gt. 0.0 ) then
! a = 17.27
! b = 237.3
!x else
!x a = 21.87
!x b = 265.5
!x end if
!
! ** 3.8 = 0.622*6.11 (hPa) **
esl=svp11*EXP(svp2*(t-svpt0)/(t-svp3))
qsl=ep_2*esl/(p(k)-ep_3*esl)
! qsl = 3.8*EXP( a*ct/( b+ct ) ) / ( 1000.0*p2a**rk )
dqsl = qsl*ep_2*ev/( rd*t**2 )
!
qmq(k) = qw(k) -qsl
alp(k) = 1.0/( 1.0+dqsl*xlvcp )
bet(k) = dqsl*p2a
!
t3sq = MAX( tsq(k), 0.0 )
r3sq = MAX( qsq(k), 0.0 )
c3sq = cov(k)
c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
!
r3sq = r3sq +bet(k)**2*t3sq -2.0*bet(k)*c3sq
sgm(k) = SQRT( MAX( r3sq, 1.0d-10 ) )
END DO
!
DO k = kts,kte-1
q1 = qmq(k) / sgm(k)
cld0 = 0.5*( 1.0+erf( q1*rr2 ) )
! q1=0.
! cld0=0.
eq1 = rrp*EXP( -0.5*q1*q1 )
qll = MAX( cld0*q1 + eq1, 0.0 )
cld(k) = cld0
ql (k) = alp(k)*sgm(k)*qll
!
q2p = xlvcp/exner( k )
pt = thl(k) +q2p*ql(k)
qt = 1.0 +p608*qw(k) -(1.+p608)*ql(k)
rac = alp(k)*( cld0-qll*eq1 )*( q2p*qt-(1.+p608)*pt )
!
vt (k) = qt-1.0 -rac*bet(k)
vq (k) = p608*pt-tv0 +rac
END DO
!
cld(kte) = cld(kte-1)
ql(kte) = ql(kte-1)
vt(kte) = vt(kte-1)
vq(kte) = vq(kte-1)
RETURN
END SUBROUTINE mym_condensation
SUBROUTINE mynn_tendencies (docs) (kts,kte,& 1,6
&levflag,grav_settling,&
&delt,&
&dz,&
&u,v,th,qv,qc,p,exner,&
&thl,sqv,sqc,sqw,&
&ust,flt,flq,wspd,qcg,&
&tsq,qsq,cov,&
&tcd,qcd,&
&dfm,dfh,dfq,&
&Du,Dv,Dth,Dqv,Dqc)
INTEGER, INTENT(in) :: kts,kte
INTEGER, INTENT(in) :: grav_settling,levflag
!! grav_settling = 1 for gravitational settling of droplets
!! grav_settling = 0 otherwise
! thl - liquid water potential temperature
! qw - total water
! dfm,dfh,dfq - as above
! flt - surface flux of thl
! flq - surface flux of qw
REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,th,qv,qc,p,exner,&
&dfm,dfh,dfq,dz,tsq,qsq,cov,tcd,qcd
REAL, DIMENSION(kts:kte), INTENT(inout) :: thl,sqw,sqv,sqc
REAL, DIMENSION(kts:kte), INTENT(out) :: du,dv,dth,dqv,dqc
REAL, INTENT(IN) :: delt,ust,flt,flq,wspd,qcg
! REAL, INTENT(IN) :: delt,ust,flt,flq,qcg,&
! &gradu_top,gradv_top,gradth_top,gradqv_top
!local vars
REAL, DIMENSION(kts:kte) :: dtz,vt,vq
REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
REAL :: rhs,gfluxm,gfluxp,dztop
INTEGER :: k,kk,nz
nz=kte-kts+1
dztop=.5*(dz(kte)+dz(kte-1))
DO k=kts,kte
dtz(k)=delt/dz(k)
ENDDO
!! u
k=kts
a(1)=0.
b(1)=1.+dtz(k)*(dfm(k+1)+ust**2/wspd)
c(1)=-dtz(k)*dfm(k+1)
d(1)=u(k)
!! a(1)=0.
!! b(1)=1.+dtz(k)*dfm(k+1)
!! c(1)=-dtz(k)*dfm(k+1)
!! d(1)=u(k)*(1.-ust**2/wspd*dtz(k))
DO k=kts+1,kte-1
kk=k-kts+1
a(kk)=-dtz(k)*dfm(k)
b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
c(kk)=-dtz(k)*dfm(k+1)
d(kk)=u(k)
ENDDO
!! no flux at the top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=0.
!! specified gradient at the top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=gradu_top*dztop
!! prescribed value
a(nz)=0
b(nz)=1.
c(nz)=0.
d(nz)=u(kte)
CALL tridiag(nz,a,b,c,d)
DO k=kts,kte
du(k)=(d(k-kts+1)-u(k))/delt
ENDDO
!! v
k=kts
a(1)=0.
b(1)=1.+dtz(k)*(dfm(k+1)+ust**2/wspd)
c(1)=-dtz(k)*dfm(k+1)
d(1)=v(k)
!! a(1)=0.
!! b(1)=1.+dtz(k)*dfm(k+1)
!! c(1)=-dtz(k)*dfm(k+1)
!! d(1)=v(k)*(1.-ust**2/wspd*dtz(k))
DO k=kts+1,kte-1
kk=k-kts+1
a(kk)=-dtz(k)*dfm(k)
b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
c(kk)=-dtz(k)*dfm(k+1)
d(kk)=v(k)
ENDDO
!! no flux at the top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=0.
!! specified gradient at the top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=gradv_top*dztop
!! prescribed value
a(nz)=0
b(nz)=1.
c(nz)=0.
d(nz)=v(kte)
CALL tridiag(nz,a,b,c,d)
DO k=kts,kte
dv(k)=(d(k-kts+1)-v(k))/delt
ENDDO
!! thl
k=kts
a(1)=0.
b(1)=1.+dtz(k)*dfh(k+1)
c(1)=-dtz(k)*dfh(k+1)
! if qcg not used then assume constant flux in the surface layer
IF (qcg < qcgmin) THEN
IF (sqc(k) > qcgmin) THEN
gfluxm=grav_settling*gno*sqc(k)**gpw
ELSE
gfluxm=0.
ENDIF
ELSE
gfluxm=grav_settling*gno*(qcg/(1.+qcg))**gpw
ENDIF
IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
ELSE
gfluxp=0.
ENDIF
rhs=-xlvcp/exner(k)&
&*( &
&(gfluxp - gfluxm)/dz(k)&
& ) + tcd(k)
d(1)=thl(k)+dtz(k)*flt+rhs*delt
DO k=kts+1,kte-1
kk=k-kts+1
a(kk)=-dtz(k)*dfh(k)
b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
c(kk)=-dtz(k)*dfh(k+1)
IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
ELSE
gfluxp=0.
ENDIF
IF (.5*(sqc(k-1)+sqc(k)) > qcgmin) THEN
gfluxm=grav_settling*gno*(.5*(sqc(k-1)+sqc(k)))**gpw
ELSE
gfluxm=0.
ENDIF
rhs=-xlvcp/exner(k)&
&*( &
&(gfluxp - gfluxm)/dz(k)&
& ) + tcd(k)
d(kk)=thl(k)+rhs*delt
ENDDO
!! no flux at the top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=0.
!! specified gradient at the top
!assume gradthl_top=gradth_top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=gradth_top*dztop
!! prescribed value
a(nz)=0.
b(nz)=1.
c(nz)=0.
d(nz)=thl(kte)
CALL tridiag(nz,a,b,c,d)
DO k=kts,kte
thl(k)=d(k-kts+1)
ENDDO
!! total water
k=kts
a(1)=0.
b(1)=1.+dtz(k)*dfh(k+1)
c(1)=-dtz(k)*dfh(k+1)
IF (qcg < qcgmin) THEN
IF (sqc(k) > qcgmin) THEN
gfluxm=grav_settling*gno*sqc(k)**gpw
ELSE
gfluxm=0.
ENDIF
ELSE
gfluxm=grav_settling*gno*(qcg/(1.+qcg))**gpw
ENDIF
IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
ELSE
gfluxp=0.
ENDIF
rhs=&
&( &
&(gfluxp - gfluxm)/dz(k)&
& ) + qcd(k)
d(1)=sqw(k)+dtz(k)*flq+rhs*delt
DO k=kts+1,kte-1
kk=k-kts+1
a(kk)=-dtz(k)*dfh(k)
b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
c(kk)=-dtz(k)*dfh(k+1)
IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
ELSE
gfluxp=0.
ENDIF
IF (.5*(sqc(k-1)+sqc(k)) > qcgmin) THEN
gfluxm=grav_settling*gno*(.5*(sqc(k-1)+sqc(k)))**gpw
ELSE
gfluxm=0.
ENDIF
rhs=&
&( &
&(gfluxp - gfluxm)/dz(k)&
& ) + qcd(k)
d(kk)=sqw(k) + rhs*delt
ENDDO
!! no flux at the top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=0.
!! specified gradient at the top
!assume gradqw_top=gradqv_top
! a(nz)=-1.
! b(nz)=1.
! c(nz)=0.
! d(nz)=gradqv_top*dztop
!! prescribed value
a(nz)=0.
b(nz)=1.
c(nz)=0.
d(nz)=sqw(kte)
CALL tridiag(nz,a,b,c,d)
!convert to mixing ratios for wrf
DO k=kts,kte
sqw(k)=d(k-kts+1)
sqv(k)=sqw(k)-sqc(k)
Dqv(k)=(sqv(k)/(1.-sqv(k))-qv(k))/delt
! Dqc(k)=(sqc(k)/(1.-sqc(k))-qc(k))/delt
Dqc(k)=0.
Dth(k)=(thl(k)+xlvcp/exner(k)*sqc(k)-th(k))/delt
ENDDO
END SUBROUTINE mynn_tendencies
SUBROUTINE retrieve_exchange_coeffs (docs) (kts,kte,& 1
&dfm,dfh,dfq,dz,&
&K_m,K_h,K_q)
INTEGER , INTENT(in) :: kts,kte
REAL, DIMENSION(KtS:KtE), INTENT(in) :: dz,dfm,dfh,dfq
REAL, DIMENSION(KtS:KtE), INTENT(out) :: &
&K_m, K_h, K_q
INTEGER :: k
REAL :: dzk
K_m(kts)=0.
K_h(kts)=0.
K_q(kts)=0.
DO k=kts+1,kte
dzk = 0.5 *( dz(k)+dz(k-1) )
K_m(k)=dfm(k)*dzk
K_h(k)=dfh(k)*dzk
K_q(k)=dfq(k)*dzk
ENDDO
END SUBROUTINE retrieve_exchange_coeffs
SUBROUTINE tridiag (docs) (n,a,b,c,d) 8,2
!! to solve system of linear eqs on tridiagonal matrix n times n
!! after Peaceman and Rachford, 1955
!! a,b,c,d - are vectors of order n
!! a,b,c - are coefficients on the LHS
!! d - is initially RHS on the output becomes a solution vector
INTEGER, INTENT(in):: n
REAL, DIMENSION(n), INTENT(in) :: a,b
REAL, DIMENSION(n), INTENT(inout) :: c,d
INTEGER :: i
REAL :: p
REAL, DIMENSION(n) :: q
c(n)=0.
q(1)=-c(1)/b(1)
d(1)=d(1)/b(1)
DO i=2,n
p=1./(b(i)+a(i)*q(i-1))
q(i)=-c(i)*p
d(i)=(d(i)-a(i)*d(i-1))*p
ENDDO
DO i=n-1,1,-1
d(i)=d(i)+q(i)*d(i+1)
ENDDO
END SUBROUTINE tridiag
SUBROUTINE mynn_bl_driver (docs) (& 1,6
&initflag,&
&grav_settling,&
&delt,&
&dz,&
&u,v,th,qv,qc,&
&p,exner,rho,&
&xland,ts,qsfc,qcg,ps,&
&ust,ch,hfx,qfx,rmol,wspd,&
&Qke,Tsq,Qsq,Cov,&
&Du,Dv,Dth,&
&Dqv,Dqc,&
! &K_m,K_h,K_q&
&K_h,k_m&
&,IDS,IDE,JDS,JDE,KDS,KDE &
&,IMS,IME,JMS,JME,KMS,KME &
&,ITS,ITE,JTS,JTE,KTS,KTE)
INTEGER, INTENT(in) :: initflag
INTEGER, INTENT(in) :: grav_settling
INTEGER,INTENT(IN) :: &
& IDS,IDE,JDS,JDE,KDS,KDE &
&,IMS,IME,JMS,JME,KMS,KME &
&,ITS,ITE,JTS,JTE,KTS,KTE
! initflag > 0 for TRUE
! else for FALSE
! levflag : <>3; Level 2.5
! = 3; Level 3
! grav_settling = 1 when gravitational settling accounted for
! grav_settling = 0 when gravitational settling NOT accounted for
REAL, INTENT(in) :: delt
REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(in) :: dz,&
&u,v,th,qv,qc,p,exner,rho
REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(in) :: xland,ust,&
&ch,rmol,ts,qsfc,qcg,ps,hfx,qfx, wspd
!
REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: &
&Qke,Tsq,Qsq,Cov
REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
&Du,Dv,Dth,Dqv,Dqc
REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
&K_h,K_m
!local vars
INTEGER :: ITF,JTF,KTF
INTEGER :: i,j,k
REAL, DIMENSION(KMS:KME) :: thl,sqv,sqc,sqw,&
&El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, Vt, Vq
REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: K_q
REAL, DIMENSION(KMS:KME+1) :: zw
REAL :: cpm,sqcg,flt,flq,pmz,phh,exnerg,zet
INTEGER, SAVE :: levflag
JTF=MIN0(JTE,JDE-1)
ITF=MIN0(ITE,IDE-1)
KTF=MIN0(KTE,KDE-1)
IF (initflag > 0) THEN
levflag=mynn_level
DO j=JTS,JTF
DO i=ITS,ITF
DO k=KTS,KTF
sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
thl(k)=th(i,k,j)
IF (k==kts) THEN
zw(k)=0.
ELSE
zw(k)=zw(k-1)+dz(i,k-1,j)
ENDIF
k_m(i,k,j)=0.
k_h(i,k,j)=0.
k_q(i,k,j)=0.
ENDDO
zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
CALL mym_initialize ( kts,kte,&
&dz(i,kts:kte,j), zw, &
&u(i,kts:kte,j), v(i,kts:kte,j), &
&thl, sqv,&
&ust(i,j), rmol(i,j),&
&Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
&Qsq(i,kts:kte,j), Cov(i,kts:kte,j))
ENDDO
ENDDO
ENDIF
DO j=JTS,JTF
DO i=ITS,ITF
DO k=KTS,KTF
sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
sqc(k)=qc(i,k,j)/(1.+qc(i,k,j))
sqw(k)=sqv(k)+sqc(k)
thl(k)=th(i,k,j)-xlvcp/exner(i,k,j)*sqc(k)
IF (k==kts) THEN
zw(k)=0.
ELSE
zw(k)=zw(k-1)+dz(i,k-1,j)
ENDIF
ENDDO
zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
sqcg=qcg(i,j)/(1.+qcg(i,j))
cpm=cp*(1.+0.8*qv(i,kts,j))
! The exchange coefficient for cloud water is assumed to be the same as
! that for heat. CH is multiplied by WSPD. See module_sf_mynn.F
exnerg=(ps(i,j)/p1000mb)**rcp
flt = hfx(i,j)/( rho(i,kts,j)*cpm ) &
+xlvcp*ch(i,j)*(sqc(kts)/exner(i,kts,j)-sqcg/exnerg)
flq = qfx(i,j)/ rho(i,kts,j) &
-ch(i,j)*(sqc(kts) -sqcg )
!!!!!
zet = 0.5*dz(i,kts,j)*rmol(i,j)
if ( zet >= 0.0 ) then
pmz = 1.0 + (cphm_st-1.0) * zet
phh = 1.0 + cphh_st * zet
else
pmz = 1.0/ (1.0-cphm_unst*zet)**0.25 - zet
phh = 1.0/SQRT(1.0-cphh_unst*zet)
end if
!!!!!
CALL mym_condensation ( kts,kte,&
&levflag, &
&dz(i,kts:kte,j), &
&thl, sqw, &
&p(i,kts:kte,j),exner(i,kts:kte,j), &
&tsq(i,kts:kte,j), qsq(i,kts:kte,j), cov(i,kts:kte,j), &
&Vt, Vq)
CALL mym_turbulence ( kts,kte,&
&levflag, &
&dz(i,kts:kte,j), zw, &
&u(i,kts:kte,j), v(i,kts:kte,j), thl, sqc, sqw, &
&qke(i,kts:kte,j), tsq(i,kts:kte,j), &
&qsq(i,kts:kte,j), cov(i,kts:kte,j), &
&vt, vq,&
&rmol(i,j), flt, flq, &
&El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc)
CALL mym_predict (kts,kte,&
&levflag, &
&delt,&
&dz(i,kts:kte,j), &
&ust(i,j), flt, flq, pmz, phh, &
&el, dfq, &
&pdk, pdt, pdq, pdc,&
&Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
&Qsq(i,kts:kte,j), Cov(i,kts:kte,j))
CALL mynn_tendencies(kts,kte,&
&levflag,grav_settling,&
&delt,&
&dz(i,kts:kte,j),&
&u(i,kts:kte,j),v(i,kts:kte,j),&
&th(i,kts:kte,j),qv(i,kts:kte,j),qc(i,kts:kte,j),&
&p(i,kts:kte,j),exner(i,kts:kte,j),&
&thl,sqv,sqc,sqw,&
&ust(i,j),flt,flq,wspd(i,j),qcg(i,j),&
&tsq(i,kts:kte,j),qsq(i,kts:kte,j),cov(i,kts:kte,j),&
&tcd,qcd,&
&dfm,dfh,dfq,&
&Du(i,kts:kte,j),Dv(i,kts:kte,j),Dth(i,kts:kte,j),&
&Dqv(i,kts:kte,j),Dqc(i,kts:kte,j))
CALL retrieve_exchange_coeffs(kts,kte,&
&dfm,dfh,dfq,dz(i,kts:kte,j),&
&K_m(i,kts:kte,j),K_h(i,kts:kte,j),K_q(i,kts:kte,j))
ENDDO
ENDDO
END SUBROUTINE mynn_bl_driver
SUBROUTINE mynn_bl_init_driver (docs) (& 1
&Du,Dv,Dth,&
&Dqv,Dqc&
&,RESTART,ALLOWED_TO_READ,LEVEL&
&,IDS,IDE,JDS,JDE,KDS,KDE &
&,IMS,IME,JMS,JME,KMS,KME &
&,ITS,ITE,JTS,JTE,KTS,KTE)
LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
INTEGER,INTENT(IN) :: LEVEL
INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
& IMS,IME,JMS,JME,KMS,KME, &
& ITS,ITE,JTS,JTE,KTS,KTE
REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: &
&Du,Dv,Dth,Dqv,Dqc
INTEGER :: I,J,K,ITF,JTF,KTF
JTF=MIN0(JTE,JDE-1)
KTF=MIN0(KTE,KDE-1)
ITF=MIN0(ITE,IDE-1)
IF(.NOT.RESTART)THEN
DO J=JTS,JTF
DO K=KTS,KTF
DO I=ITS,ITF
Du(i,k,j)=0.
Dv(i,k,j)=0.
Dth(i,k,j)=0.
Dqv(i,k,j)=0.
Dqc(i,k,j)=0.
ENDDO
ENDDO
ENDDO
ENDIF
mynn_level=level
END SUBROUTINE mynn_bl_init_driver
END MODULE module_bl_mynn