subroutine da_sfc_wtq_lin(psfc, tg, ps, ts, qs, us, vs, regime, & 1,5
psfc_prime,tg_prime,ps_prime,ts_prime,qs_prime, &
us_prime,vs_prime, hs,roughness,xland,dx, & ! Modified by Eric Chiang (JULY 2010)
u10_prime,v10_prime,t2_prime,q2_prime)
!---------------------------------------------------------------------------
! Purpose: Calculate the 10m wind, 2m temperature and moisture based on the
! similarity theory/
!
! Reference:
! ---------
!
! input Variables(basic state):
!
! psfc, tg : surface pressure and ground temperature
! ps, ts, qs, us, vs, hs : model variable at lowlest half sigma level
! regime : PBL regime
!
! input Variables(pertubation):
!
! psfc_prime, tg_prime : Surface pressure and ground temperature
! ps_prime, ts_prime, : Model variables at the lowest half sigma
! qs_prime, us_prime, : level
! vs_prime :
!
! Constants:
!
! hs : height at the lowest half sigma level
! roughness : roughness
! xland : land-water-mask (=2 water, =1 land)
!
! output Variables(pertubation):
!
! u10_prime, v10_prime : 10-m high observed wind components
! t2_prime , q2_prime : 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) :: regime
real, intent (in) :: ps , ts , qs , us, vs, psfc, tg
real, intent (in) :: ps_prime, ts_prime, qs_prime, us_prime, vs_prime, psfc_prime, tg_prime
real, intent (in) :: hs, roughness, xland
real, intent (out) :: u10_prime, v10_prime, t2_prime, q2_prime
! Maximum number of iterations in computing psim, psih
integer, parameter :: k_iteration = 10
! integer, parameter :: k_iteration = 1
! 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
real :: rib, rcp, xx, yy, cc, Pi
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
real :: ths, thg, thvs, thvg, vsgd, vsgd2, dx
real :: zq0, z0
real :: Vc2_prime, Va2_prime, V2_prime
real :: rib_prime, xx_prime, yy_prime
real :: psiw_prime, psiz_prime, mol_prime, ust_prime, &
hol_prime, holz_prime, hol2_prime
real :: psim_prime, psimz_prime, psim2_prime, &
psih_prime, psihz_prime, psih2_prime
real :: psit_prime, psit2_prime, &
psiq_prime, psiq2_prime
real :: qg_prime, tvg_prime, tvs_prime
real :: ths_prime, thg_prime, thvs_prime, thvg_prime
real, parameter :: ka = 2.4E-5
integer :: iregime
if (trace_use) call da_trace_entry
("da_sfc_wtq_lin")
rcp = gas_constant/cp
! 1 Compute the roughness length based upon season and land use
! =====================================
! 1.1 Define the rouhness length
! -----------------
z0 = roughness
if (z0 < 0.0001) z0 = 0.0001
! 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_prime = ts_prime * (1.0 + 0.608 * qs) + 0.608 * ts * qs_prime
tvs = ts * (1.0 + 0.608 * qs)
! 2.2 Compute the ground saturated mixing ratio and the ground virtual
! temperature
! ----------------------------------------------------------------
call da_tp_to_qs(tg, psfc, eg, qg)
call da_tp_to_qs_lin1(tg, psfc, eg, tg_prime, psfc_prime, qg_prime)
qg_prime = qg_prime * qg
tvg_prime = tg_prime * (1.0 + 0.608 * qg) + 0.608 * tg * qg_prime
tvg = tg * (1.0 + 0.608 * qg)
! 3. Compute the potential temperature and virtual potential temperature
! =======================================================================
! 3.1 Potential temperature on the lowest half sigma level
! ----------------------------------------------------
Pi = (100000.0 / ps) ** rcp
ths_prime = (ts_prime - ps_prime * rcp * ts/ps) * Pi
ths = ts * Pi
! 3.2 Virtual potential temperature on the lowest half sigma level
! ------------------------------------------------------------
thvs_prime = (tvs_prime - ps_prime * rcp * tvs/ps) * Pi
thvs = tvs * Pi
! 3.3 Potential temperature at the ground
! -----------------------------------
Pi = (100000.0 / psfc) ** rcp
thg_prime = (tg_prime - psfc_prime * rcp * tg/psfc) * Pi
thg = tg * Pi
! 3.4 Virtual potential temperature at ground
! ---------------------------------------
thvg_prime = (tvg_prime - psfc_prime * rcp * tvg/psfc) * Pi
thvg = tvg * Pi
! 4. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
! =================================================
! 4.1 Velocity
! --------
! Wind speed:
Va2_prime = 2.0*us*us_prime + 2.0*vs*vs_prime
Va2 = us*us + vs*vs
! Convective velocity:
if (thvg >= thvs) then
Vc2_prime = 4.0 * (thvg_prime - thvs_prime)
Vc2 = 4.0 * (thvg - thvs)
else
Vc2_prime = 0.0
Vc2 = 0.0
end if
! Calculate Mahrt and Sun low-res correction ! Add by Eric Chiang ( July 2010 )
! dx is a constant, so vsgd is also a constant,
! the perturnations of vsgd and vsgd2 are zero. (YRG, 09/15/2011)
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_prime should be computed before used below. (YRG, 09/15/2011)
V2_prime = Va2_prime + Vc2_prime
V2 = Va2 + Vc2 + vsgd2 ! Modified by Eric Chiang ( July 2010 )
! 4.2 Bulk richardson number
! ----------------------
Pi = gravity * hs / (ths*V2)
rib_prime = (thvs_prime - thvg_prime &
- (thvs-thvg)/V2 * V2_prime &
- (thvs-thvg)/ths * ths_prime) * Pi
rib = (thvs - thvg) * Pi
! 5. CALCULATE PSI BASED UPON REGIME
! =======================================
iregime = int(regime)
select case (iregime)
! 5.1 Stable conditions (REGIME 1)
! ---------------------------
case (1);
psim_prime = 0.0
psimz_prime = 0.0
psim2_prime = 0.0
psim = -10.0*gzsoz0
psimz = -10.0*gz10oz0
psim2 = -10.0*gz2oz0
psim = max(psim,-10.0)
psimz = max(psimz,-10.0)
psim2 = max(psim2,-10.0)
psih_prime = psim_prime
psihz_prime = psimz_prime
psih2_prime = psim2_prime
psih = psim
psihz = psimz
psih2 = psim2
! 5.2 Mechanically driven turbulence (REGIME 2)
! ------------------------------------------
case (2);
Pi = - 1.0 / ((1.1 - 5.0*rib)*(1.1 - 5.0*rib))
psim_prime = 5.5 * gzsoz0 * rib_prime * Pi
psimz_prime = 5.5 * gz10oz0 * rib_prime * Pi
psim2_prime = 5.5 * gz2oz0 * rib_prime * Pi
Pi = (-5.0 * rib) / (1.1 - 5.0*rib)
psim = gzsoz0 * Pi
psimz = gz10oz0 * Pi
psim2 = gz2oz0 * Pi
if (psim >= -10.0) then
psim = psim
psim_prime = psim_prime
else
psim = -10.0
psim_prime = 0.0
end if
if (psimz >= -10.0) then
psimz = psimz
psimz_prime = psimz_prime
else
psimz = -10.0
psimz_prime = 0.0
end if
if (psim2 >= -10.0) then
psim2 = psim2
psim2_prime = psim2_prime
else
psim2 = -10.0
psim2_prime = 0.0
end if
psih_prime = psim_prime
psihz_prime = psimz_prime
psih2_prime = psim2_prime
psih = psim
psihz = psimz
psih2 = psim2
! 5.3 Unstable Forced convection (REGIME 3)
! -------------------------------------
case (3);
psim_prime = 0.0
psimz_prime = 0.0
psim2_prime = 0.0
psim = 0.0
psimz = 0.0
psim2 = 0.0
psih_prime = psim_prime
psihz_prime = psimz_prime
psih2_prime = psim2_prime
psih = psim
psihz = psimz
psih2 = psim2
! 5.4 Free convection (REGIME 4)
! --------------------------
case (4);
! Calculate psi m and pshi h using iteration method
psim_prime = 0.0
psih_prime = 0.0
psim = 0.0
psih = 0.0
cc = 2.0 * atan(1.0)
! do k = 1 , k_iteration
! 5.4.1 Calculate ust, m/L (mol), h/L (hol)
! --------------------------
! Friction speed
ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust
! Heat fux factor
mol = k_kar * (ths - thg)/(gzsoz0 - psih)
mol_prime = ((ths_prime - thg_prime) /(ths - thg) + &
psih_prime /(gzsoz0 - psih)) * mol
! Ratio of PBL height to Monin-Obukhov length
if (ust .LT. 0.01) then
hol_prime = rib_prime * gzsoz0
hol = rib * gzsoz0
else
hol = k_kar * gravity * hs * mol / (ths * ust * ust)
hol_prime = (mol_prime / mol - ths_prime / ths &
- 2.0* ust_prime / ust) * hol
end if
! 5.4.2 Calculate n, nz, R, Rz
! --------------------------
if (hol >= 0.0) then
hol_prime = 0.0
hol = 0.0
else
hol_prime = hol_prime
hol = hol
end if
if (hol >= -10.0) then
hol_prime = hol_prime
hol = hol
else
hol_prime = 0.0
hol = -10.0
end if
holz_prime = (h10 / hs) * hol_prime
holz = (h10 / hs) * hol
if (holz >= 0.0) then
holz_prime = 0.0
holz = 0.0
else
holz_prime = holz_prime
holz = holz
end if
if (holz >= -10.0) then
holz_prime = holz_prime
holz = holz
else
holz_prime = 0.0
holz = -10.0
end if
hol2_prime = (h2 / hs) * hol_prime
hol2 = (h2 / hs) * hol
if (hol2 >= 0.0) then
hol2_prime = 0.0
hol2 = 0.0
else
hol2_prime = hol2_prime
hol2 = hol2
end if
if (hol2 >= -10.0) then
hol2_prime = hol2_prime
hol2 = hol2
else
hol2_prime = 0.0
hol2 = -10.0
end if
! 5.4.3 Calculate Psim & psih
! --------------------------
! Using the continuous function:
xx_prime = -4.0* hol_prime /((1.0 - 16.0 * hol) ** 0.75)
xx = (1.0 - 16.0 * hol) ** 0.25
yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
yy = log((1.0+xx*xx)/2.0)
psim_prime = 2 * xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
psim = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
psih_prime = 2.0 * yy_prime
psih = 2.0 * yy
! Using the continuous function:
xx_prime = -4.0* holz_prime /((1.0 - 16.0 * holz) ** 0.75)
xx = (1.0 - 16.0 * holz) ** 0.25
yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
yy = log((1.0+xx*xx)/2.0)
psimz_prime = 2.0* xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
psimz = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
psihz_prime = 2.0 * yy_prime
psihz = 2.0 * yy
! Using the continuous function:
xx_prime = -4.0* hol2_prime /((1.0 - 16.0 * hol2) ** 0.75)
xx = (1.0 - 16.0 * hol2) ** 0.25
yy_prime = 2.0* xx * xx_prime /(1.0+xx*xx)
yy = log((1.0+xx*xx)/2.0)
psim2_prime = 2.0* xx_prime *(1.0/(1.0+xx)-1.0/(1+xx*xx)) + yy_prime
psim2 = 2.0 * log((1.0+xx)/2.0) + yy - 2.0 * atan(xx) + cc
psih2_prime = 2.0 * yy_prime
psih2 = 2.0 * yy
! end do
! 5.4.4 Define the limit value for psim & psih
! --------------------------
if (psim <= 0.9*gzsoz0) then
psim_prime = psim_prime
psim = psim
else
psim = 0.9*gzsoz0
psim_prime = 0.0
end if
if (psimz <= 0.9*gz10oz0) then
psimz_prime = psimz_prime
psimz = psimz
else
psimz_prime = 0.0
psimz = 0.9*gz10oz0
end if
if (psim2 <= 0.9*gz2oz0) then
psim2_prime = psim2_prime
psim2 = psim2
else
psim2_prime = 0.0
psim2 = 0.9*gz2oz0
end if
if (psih <= 0.9*gzsoz0) then
psih_prime = psih_prime
psih = psih
else
psih_prime = 0.0
psih = 0.9*gzsoz0
end if
if (psihz <= 0.9*gz10oz0) then
psihz_prime = psihz_prime
psihz = psihz
else
psihz_prime = 0.0
psihz = 0.9*gz10oz0
end if
if (psih2 <= 0.9*gz2oz0) then
psih2_prime = psih2_prime
psih2 = psih2
else
psih2_prime = 0.0
psih2 = 0.9*gz2oz0
end if
case default;
write(unit=message(1),fmt='(A,I2,A)') &
"Regime=",iregime," is invalid."
call da_error(__FILE__,__LINE__,message(1:1))
end select
! 6. CALCULATE PSI FOR WinD, TEMPERATURE AND MOISTURE
! =======================================
psiw_prime = - psim_prime
psiw = gzsoz0 - psim
psiz_prime = - psimz_prime
psiz = gz10oz0 - psimz
psit_prime = - psih_prime
psit = gzsoz0 - psih
psit2_prime = - psih2_prime
psit2 = gz2oz0 - psih2
ust = k_kar * sqrt(v2) /(gzsoz0 - psim)
ust_prime = (0.5/V2 * v2_prime + psim_prime /(gzsoz0 - psim)) * ust
psiq2_prime = k_kar*hs/(ka*(k_kar*ust*hs/ka + hs / zq0))*ust_prime
psiq_prime = psiq2_prime - psih_prime
psiq2_prime = psiq2_prime - psih2_prime
psiq = log(k_kar*ust*hs/ka + hs / zq0) - psih
psiq2 = log(k_kar*ust*h2/ka + h2 / zq0) - psih2
! 7. CALCULATE THE PERTURBATIONS for 10M WinD, 2M TEMPERATURE AND MOISTURE
! =======================================
Pi = psiz / psiw
u10_prime= (us_prime + us/psiz * psiz_prime - us/psiw * psiw_prime) * Pi
v10_prime= (vs_prime + vs/psiz * psiz_prime - vs/psiw * psiw_prime) * Pi
t2_prime = ((1.0-psit2/psit) * thg_prime + (ths_prime + &
(ths - thg)/psit2 * psit2_prime - &
(ths - thg)/psit * psit_prime) * psit2/psit &
+ rcp*(thg + (ths - thg)*psit2/psit)/psfc * psfc_prime) &
* (psfc/100000.0)**rcp
q2_prime = (1.0-psiq2/psiq) * qg_prime + psiq2/psiq * qs_prime + &
(qs -qg)*(psiq2/psiq) * (psiq2_prime/psiq2 - psiq_prime/psiq)
if (trace_use) call da_trace_exit("da_sfc_wtq_lin")
end subroutine da_sfc_wtq_lin