<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_SFC_WTQ'><A href='../../html_code/physics/da_sfc_wtq.inc.html#DA_SFC_WTQ' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_sfc_wtq (psfc, tg, ps, ts, qs, us, vs, & 3,3
hs, roughness, xland, dx, u10, v10, t2, q2, regime, &
has_lsm, regime_wrf, qsfc_wrf, znt_wrf, ust_wrf, mol_wrf, hfx, qfx, pblh)
!---------------------------------------------------------------------------
! Purpose: Calculate the 10m wind, 2m temperature and moisture based on the
! similarity theory/
!
! The unit for pressure : psfc, ps is Pa.
! The unit for temperature: tg, ts, t2 is K.
! The unit for moisture : qs, q2 is kg/kg.
! The unit for wind : us, vs, u10, v10 is m/s.
! The unit for height : hs, roughness is m.
! xland and regime are dimensionless.
!
! History: Nov 2010 - improve calculation consistency with WRF model (Eric Chiang)
! Jul 2015 - further improvement on consistency
!
! Reference:
! ---------
!
! input Variables:
!
! psfc, tg : surface pressure and ground temperature
! ps, ts, qs, us, vs, hs : model variable at lowlest half sigma level
! dx (m) : horizontal resolution
!
!
! Constants:
!
! hs : height at the lowest half sigma level
! roughness : roughness
! xland : land-water-mask (=2 water, =1 land)
!
! output Variables:
!
! regime : PBL regime
! u10, v10 : 10-m high observed wind components
! t2 , q2 : 2-m high observed temperature and mixing ratio
!
!---------------------------------------------------------------------------
!
! psim : mechanical psi at lowlest sigma level
! psim2 : mechanical psi at 2m
! psimz : mechanical psi at 10m
!
!---------------------------------------------------------------------------
implicit none
real, intent (in) :: ps , ts , qs , us, vs
real, intent (in) :: psfc, tg
real, intent (in) :: hs, roughness, xland
real, intent (out) :: regime
real, intent (out) :: u10, v10, t2, q2
logical, intent(in), optional :: has_lsm
real, intent(in), optional :: regime_wrf, qsfc_wrf, znt_wrf, ust_wrf, mol_wrf
real, intent(in), optional :: hfx, qfx, pblh
logical :: use_table = .true.
logical :: use_ust_wrf = .false.
logical :: vconv_wrf
integer :: nn, nz, n2
real :: rr, rz, r2
real :: cqs2, chs2, rho, rhox, fluxc, visc, restar, z0t, z0q
! h10 is the height of 10m where the wind observed
! h2 is the height of 2m where the temperature and
! moisture observed.
real, parameter :: h10 = 10.0, h2 = 2.0
! Default roughness over the land
real, parameter :: zint0 = 0.01
! Von Karman constant
real, parameter :: k_kar = 0.4
! Working variables
real :: Vc2, Va2, V2, vc, wspd
real :: rib, rcp, xx, yy, cc
real :: psiw, psiz, mol, ust, hol, holz, hol2
real :: psim, psimz, psim2, psih, psihz, psih2
real :: psit, psit2, psiq, psiq2
real :: gzsoz0, gz10oz0, gz2oz0
real :: eg, qg, tvg, tvs, tvs2
real :: ths, thg, thvs, thvg, thvs2, vsgd, vsgd2, dx
real :: zq0, z0
real, parameter :: ka = 2.4E-5
if (trace_use_dull) call da_trace_entry
("da_sfc_wtq")
rcp = gas_constant/cp
! 1 Compute the roughness length based upon season and land use
! 1.1 Define the roughness length
z0 = roughness
if (z0 < 0.0001) z0 = 0.0001
if ( present(znt_wrf) ) then
if ( znt_wrf > 0.0 ) then
z0 = znt_wrf
end if
end if
! 1.2 Define the rouhgness length for moisture
if (xland .ge. 1.5) then
zq0 = z0
else
zq0 = zint0
end if
! 1.3 Define the some constant variable for psi
gzsoz0 = log(hs/z0)
gz10oz0 = log(h10/z0)
gz2oz0 = log(h2/z0)
! 2. Calculate the virtual temperature
! 2.1 Compute Virtual temperature on the lowest half sigma level
tvs = ts * (1.0 + 0.608 * qs)
! 2.2 Convert ground virtual temperature assuming it's saturated
call da_tp_to_qs(tg, psfc, eg, qg) !output qg is specific humidity
qg = qg*(1.0-qg) !hcl convert to mixing ratio
if ( present(qsfc_wrf) ) then
if ( qsfc_wrf > 0.0 ) then
qg = qsfc_wrf
end if
endif
tvg = tg * (1.0 + 0.608 * qg)
! 3. Compute the potential temperature
! 3.1 Potential temperature on the lowest half sigma level
ths = ts * (1000.0 / (ps/100.0)) ** rcp
! 3.2 Potential temperature at the ground
thg = tg * (1000.0 / (psfc/100.0)) ** rcp
! 4. Virtual potential temperature
! 4.1 Virtual potential temperature on the lowest half sigma level
thvs = tvs * (1000.0 / (ps/100.0)) ** rcp
! 4.2 Virtual potential temperature at ground
thvg = tvg * (1000.0 / (psfc/100.0)) ** rcp
! 5. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
! 5.1 Velocity
! Wind speed:
Va2 = us*us + vs*vs
!
! Convective velocity:
vconv_wrf = .false.
if ( present(hfx) .and. present(qfx) .and. present(pblh) ) then
! calculating vconv over land following wrf method
if ( pblh > 0.0 ) then
vconv_wrf = .true.
end if
end if
if (thvg >= thvs) then
! prior to V3.7, Vc2 = 4.0 * (thvg - thvs)
Vc2 = thvg - thvs
else
Vc2 = 0.0
end if
if ( xland < 1.5 ) then !land
if ( vconv_wrf ) then
! following the calculation as in module_sf_sfclay.F
rhox = psfc/(gas_constant*tvg)
fluxc = max(hfx/rhox/cp+0.608*tvg*qfx/rhox, 0.0)
vc = (gravity/tg*pblh*fluxc)**0.33
vc2 = vc*vc
end if
end if
! Calculate Mahrt and Sun low-res correction ! Add by Eric Chiang ( July 2010 )
vsgd = 0.32 * (max(dx/5000.-1.,0.))**0.33 ! Add by Eric Chiang ( July 2010 )
vsgd2 = vsgd * vsgd ! Add by Eric Chiang ( July 2010 )
V2 = Va2 + Vc2 + vsgd2 ! Add by Eric Chiang ( July 2010 )
wspd = sqrt(v2)
wspd = max(wspd,0.1)
v2 = wspd*wspd
! 5.2 Bulk richardson number
rib = (gravity * hs / ths) * (thvs - thvg) / V2
! if previously unstable, do not let into regime 1 and 2
if ( present(mol_wrf) ) then
if ( mol_wrf < 0.0 ) rib = min(rib, 0.0)
end if
! Calculate ust, m/L (mol), h/L (hol)
psim = 0.0
psih = 0.0
! Friction speed
if ( present(ust_wrf) ) then
if ( ust_wrf > 0.0 ) then
use_ust_wrf = .true.
ust = ust_wrf
end if
end if
if ( .not. use_ust_wrf ) then
!ust = 0.0001 !init value as in phys/module_physics_init.F
ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
end if
! Heat flux factor
if ( present(mol_wrf) ) then
mol = mol_wrf
else
mol = k_kar * (ths - thg)/(gzsoz0 - psih)
!mol = 0.0
end if
! set regimes based on rib
if (rib .GE. 0.2) then
! Stable conditions (REGIME 1)
regime = 1.1
else if ((rib .LT. 0.2) .AND. (rib .GT. 0.0)) then
! Mechanically driven turbulence (REGIME 2)
regime = 2.1
else if (rib .EQ. 0.0) then
! Unstable Forced convection (REGIME 3)
regime = 3.1
else
! Free convection (REGIME 4)
regime = 4.1
end if
if ( present(regime_wrf) ) then
if ( regime_wrf > 0.0 ) then
regime = regime_wrf
end if
end if
! 6. CALCULATE PSI BASED UPON REGIME
!if (rib .GE. 0.2) then
if ( nint(regime) == 1 ) then
! 6.1 Stable conditions (REGIME 1)
! ---------------------------
regime = 1.1
psim = -10.0*gzsoz0
psim = max(psim,-10.0)
psimz = h10/hs * psim
psimz = max(psimz,-10.0)
psim2 = h2 /hs * psim
psim2 = max(psim2,-10.0)
psih = psim
psihz = psimz
psih2 = psim2
!else if ((rib .LT. 0.2) .AND. (rib .GT. 0.0)) then
else if ( nint(regime) == 2 ) then
! 6.2 Mechanically driven turbulence (REGIME 2)
regime = 2.1
psim = (-5.0 * rib) * gzsoz0 / (1.1 - 5.0*rib)
psim = max(psim,-10.0)
psimz = h10/hs * psim
psimz = max(psimz,-10.0)
psim2 = h2 /hs * psim
psim2 = max(psim2,-10.0)
psih = psim
psihz = psimz
psih2 = psim2
!else if (rib .EQ. 0.0) then
else if ( nint(regime) == 3 ) then
! 6.3 Unstable Forced convection (REGIME 3)
regime = 3.1
psim = 0.0
psimz = 0.0
psim2 = 0.0
psih = psim
psihz = psimz
psih2 = psim2
else
! 6.4 Free convection (REGIME 4)
regime = 4.1
cc = 2.0 * atan(1.0)
! Ratio of PBL height to Monin-Obukhov length
if (ust .LT. 0.01) then
hol = rib * gzsoz0
else
hol = k_kar * gravity * hs * mol / (ths * ust * ust)
end if
! 6.4.2 Calculate n, nz, R, Rz
holz = (h10 / hs) * hol
hol2 = (h2 / hs) * hol
hol = min(hol,0.0)
hol = max(hol,-9.9999)
holz = min(holz,0.0)
holz = max(holz,-9.9999)
hol2 = min(hol2,0.0)
hol2 = max(hol2,-9.9999)
! 6.4.3 Calculate Psim & psih
if ( use_table ) then
! Using the look-up table:
nn = int(-100.0 * hol)
rr = (-100.0 * hol) - nn
psim = psimtb(nn) + rr * (psimtb(nn+1) - psimtb(nn))
psih = psihtb(nn) + rr * (psihtb(nn+1) - psihtb(nn))
else
! Using the continuous function:
xx = (1.0 - 16.0 * hol) ** 0.25
yy = log((1.0+xx*xx)/2.0)
psim = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
psih = 2.0 * yy
end if
if ( use_table ) then
! Using the look-up table:
nz = int(-100.0 * holz)
rz = (-100.0 * holz) - nz
psimz = psimtb(nz) + rz * (psimtb(nz+1) - psimtb(nz))
psihz = psihtb(nz) + rz * (psihtb(nz+1) - psihtb(nz))
else
! Using the continuous function:
xx = (1.0 - 16.0 * holz) ** 0.25
yy = log((1.0+xx*xx)/2.0)
psimz = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
psihz = 2.0 * yy
end if
if ( use_table ) then
! Using the look-up table:
n2 = int(-100.0 * hol2)
r2 = (-100.0 * hol2) - n2
psim2 = psimtb(n2) + r2 * (psimtb(n2+1) - psimtb(n2))
psih2 = psihtb(n2) + r2 * (psihtb(n2+1) - psihtb(n2))
else
! Using the continuous function:
xx = (1.0 - 16.0 * hol2) ** 0.25
yy = log((1.0+xx*xx)/2.0)
psim2 = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
psih2 = 2.0 * yy
end if
! 6.4.4 Define the limit value for psim & psih
psim = min(psim,0.9*gzsoz0)
psimz = min(psimz,0.9*gz10oz0)
psim2 = min(psim2,0.9*gz2oz0)
psih = min(psih,0.9*gzsoz0)
psihz = min(psihz,0.9*gz10oz0)
psih2 = min(psih2,0.9*gz2oz0)
end if ! Regime
! 7. Calculate psi for wind, temperature and moisture
psiw = gzsoz0 - psim
psiz = gz10oz0 - psimz
psit = max(gzsoz0-psih, 2.0)
psit2 = gz2oz0 - psih2
if ( .not. use_ust_wrf ) then
! re-calculate ust since psim is now available
ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
end if
psiq = log(k_kar*ust*hs/ka + hs / zq0) - psih
psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2
!V3.7, as in module_sf_sfclay.F
if ( xland >= 1.5 ) then !water
visc = (1.32+0.009*(ts-273.15))*1.e-5
restar = ust*z0/visc
z0t = (5.5e-5)*(restar**(-0.60))
z0t = min(z0t,1.0e-4)
z0t = max(z0t,2.0e-9)
z0q = z0t
psiq = max(log((hs+z0q)/z0q)-psih, 2.)
psit = max(log((hs+z0t)/z0t)-psih, 2.)
psiq2 = max(log((2.+z0q)/z0q)-psih2, 2.)
psit2 = max(log((2.+z0t)/z0t)-psih2, 2.)
end if
! 8. Calculate 10m wind, 2m temperature and moisture
u10 = us * psiz / psiw
v10 = vs * psiz / psiw
t2 = (thg + (ths - thg)*psit2/psit)*((psfc/100.0)/1000.0)**rcp
q2 = qg + (qs - qg)*psiq2/psiq
if ( present(has_lsm) ) then
if ( has_lsm ) then
!cqs2: 2m surface exchange coefficient for moisture
!chs2: 2m surface exchange coefficient for heat
cqs2 = ust*k_kar/psiq2
if (xland .ge. 1.5) then
!water
chs2 = ust*k_kar/psit2
else
!land
chs2 = cqs2 !as in subroutine lsm in phys/module_sf_noahdrv.F
end if
!re-calculate T2/Q2 as in module_sf_sfcdiags.F
rho = psfc/(gas_constant*tg)
if ( cqs2 < 1.e-5 ) then
q2 = qg
else
if ( present(qfx) ) then
q2 = qg - qfx/(rho*cqs2)
end if
end if
if ( chs2 < 1.e-5 ) then
t2 = tg
else
if ( present(hfx) ) then
t2 = tg - hfx/(rho*cp*chs2)
end if
end if
end if
end if
if (trace_use_dull) call da_trace_exit("da_sfc_wtq")
end subroutine da_sfc_wtq