!$$$ module documentation block
! . . . .
! module: sfc_model
! prgmmr: lee/parrish org: np23 date: 2005-12-22
!
! abstract: This module containes routines (including forward and tangent
! linear) to apply a boundary layer model for use when
! assimilating near surface conventional observations
!
! program history log:
! 2004 s.j.lee
! 2005-12-22 parrish - reformat and introduce into then current GSI release
! 2006-04-07 park - change the sensible temperature of the lowest 1st and 2nd sigma level
! and ground virtual temperature on SFC_WTQ_FWD
! 2006-04-07 park - change the sensible temperature perturbation of the lowest sigma level
! and use the T2 on sfc_wtq_Lin
! 2006-05-05 park - use the constants module
! 2006-07-14 parrish - modify so input/output pressure is in units of cb, and remove ln(psfc) references
! 2006-07-28 derber - use r1000 from constants module
! 2006-09-20 derber - add t sensible or virtual dependency on iqtflg
! 2006-09-28 treadon - add f10 (10m wind factor) to subroutine argument list
! 2007-04-03 treadon - replace 1_1 on "pi =( -five * rib ) / (1_1 - five * rib )" line with r1_1
!
! subroutines included:
!
! variable definitions:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
SUBROUTINE SFC_WTQ_FWD (psfc_in,tg,ps_in,tvs,qs,us,vs, &
ps2_in,tvs2,qs2, hs, roughness, iland, &
f10, u10, v10, t2, q2, regime, iqtflg)
!$$$ subprogram documentation block
! . . . .
! subprogram: SFC_WTQ_FWD
!
! prgrmmr:
!
! abstract: Calculate the 10m wind, 2m (virtual) temperature and moisture based on the
! similarity theory/
!
! The unit for input pressure variables psfc_in, ps_in, ps2_in is cb.
! The unit for pressure : psfc, ps, ps2 is Pa.
! The unit for temperature: tg, tvs, tvs2, tv2,t2 is K.
! The unit for moisture : qs, qs2, 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.
!
! Reference:
! ---------
! Detail MM5/physics/pbl_sfc/mrfpbl/MRFPBL.F
!
! program history log:
! NOTE: changed code so input layer temps are virtual, and are internally
! converted to sensible as required. this corrects error where
! ths is potential temperature for lowest half sigma layer, but
! before this change is in fact potential virtual temperature.
! also, output t2 should be potential temperature, so changed t2 to tv2
! 2008-04-14 safford - completed standard documentation block, rm unused vars
!
! Input argument list:
! psfc, tg : surface pressure and ground temperature
! ps, tvs, qs, us, vs, hs : model variable at lowlest half sigma leve1
! ps2, tvs2, qs2 : model variables at the second lowest half
! sigma level
! iqtflg : flag true if output is virtual temperature
! otherwise sensible temperature
!
! Output argument list:
! regime : PBL regime
! u10, v10 : 10-m high observed wind components
! t2 , q2 : 2-m high observed virtual (or sensible) temperature and
! mixing ratio
! Constant argument list:
!
! hs : height at the lowest half sigma level (above surface)
! roughness : roughness
! xland : land-water-mask (<=0.5 water, >0.5 land)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: grav,fv,rd_over_cp,zero,quarter,one,two,four,five,r1000,r10,r100,r0_01
IMPLICIT NONE
real(r_kind) ,intent(in ) :: ps_in,tvs,qs,us,vs
real(r_kind) ,intent(in ) :: ps2_in,tvs2,qs2,psfc_in,tg
real(r_kind) ,intent(in ) :: hs,roughness
integer(i_kind),intent(in ) :: iland
logical ,intent(in ) :: iqtflg
integer(i_kind),intent( out) :: regime
real(r_kind) ,intent( out) :: f10,u10,v10,q2,t2
! Maximum number of iterations in computing psim, psih
! INTEGER(i_kind), PARAMETER :: k_iteration = 10
! INTEGER(i_kind), 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(r_kind), PARAMETER :: h10 = 10.0_r_kind
real(r_kind), parameter :: h2 = two
!
! Default roughness over the land
REAL(r_kind), PARAMETER :: zint0 = 0.01_r_kind
!
! Von Karman constant
REAL(r_kind), PARAMETER :: k_kar = 0.4_r_kind
!
! Working variables
real(r_kind) :: psfc,ps,ps2
REAL(r_kind) :: Vc2, Va2, V2
REAL(r_kind) :: rib, xx, yy, cc
REAL(r_kind) :: psiw, psiz, mol, ust, hol, holz, hol2
REAL(r_kind) :: psim, psimz, psim2, psih, psihz, psih2
REAL(r_kind) :: psit, psit2, psiq, psiq2
REAL(r_kind) :: gzsoz0, gz10oz0, gz2oz0
REAL(r_kind) :: eg, qg, tvg, ts, ts2 ! change tvs, tvs2 to ts, ts2, add t2
REAL(r_kind) :: ths, thg, thvs, thvg, thvs2
REAL(r_kind) :: zq0, z0
REAL(r_kind), PARAMETER :: ka = 2.4E-5_r_kind
! local
real(r_kind),parameter :: r16 = 16.0_r_kind
real(r_kind),parameter :: r1_1 = 1.1_r_kind
real(r_kind),parameter :: r0_9 = 0.9_r_kind
real(r_kind),parameter :: r0_2 = 0.2_r_kind
!-----------------------------------------------------------------------------!
! convert input pressure variables to mb from cb.
psfc = r10*psfc_in
ps = r10*ps_in
ps2 = r10*ps2_in
! 1 Compute the roughness length based upon season and land use
! =====================================
! 1.1 Define the rouhness length
! -----------------
z0 = roughness/r100 ! unit change : originally cm --> m
if( z0 < 0.0001_r_kind) z0 = 0.0001_r_kind
! 1.2 Define the rouhgness length for moisture
! -----------------
if ( iland == 0 ) then
zq0 = z0
else
zq0 = zint0
endif
! 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
! ---------------------------------------------------------
ts = tvs / (one + fv * qs)
! 2.2 Compute Virtual temperature on the second lowest half sigma level
! -----------------------------------------------------------------
ts2 = tvs2 / (one + fv * qs2)
! 2.3 Convert ground virtual temperature assuming it's saturated
! ----------------------------------------------------------
call DA_TP_To_Qs( tg, psfc, eg, qg )
tvg = tg * (one + fv * qg)
! 3. Compute the potential temperature
! ======================================
! 3.1 Potential temperature on the lowest half sigma level
! ----------------------------------------------------
ths = ts * (r1000 / ps) ** rd_over_cp
! 3.2 Potential temperature at the ground
! -----------------------------------
thg = tg * (r1000 / psfc) ** rd_over_cp
! 4. Virtual potential temperature
! ==================================
! 4.1 Virtual potential temperature on the lowest half sigma level
! ------------------------------------------------------------
thvs = tvs * (r1000 / ps) ** rd_over_cp
! 4.2 Virtual potential temperature on the second lowest half sigma level
! ------------------------------------------------------------------
thvs2 = tvs2 * (r1000 / ps2) ** rd_over_cp
! 4.3 Virtual potential temperature at ground
! ---------------------------------------
thvg = tvg * (r1000 / psfc) ** rd_over_cp
! 5. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
! =================================================
! 5.1 Velocity
! --------
!
! Wind speed:
Va2 = us*us + vs*vs
!
! Convective velocity:
if ( thvg >= thvs ) then
Vc2 = four * (thvg - thvs)
else
Vc2 = zero
endif
!
V2 = .000001_r_kind + Va2 + Vc2 ! can't let V2 = 0, so add small value on, corresponding to 1mm/sec wind
! 5.2 Bulk richardson number
! ----------------------
rib = (grav * hs / ths) * (thvs - thvg) / V2
! 6. CALCULATE PSI BASED UPON REGIME
! =======================================
! 6.1 Stable conditions (REGIME 1)
! ---------------------------
IF (rib >= r0_2) THEN
regime = 1
psim = -r10*gzsoz0
psimz = -r10*gz10oz0
psim2 = -r10*gz2oz0
psim = max(psim,-r10)
psimz = max(psimz,-r10)
psim2 = max(psim2,-r10)
psih = psim
psihz = psimz
psih2 = psim2
! 6.2 Mechanically driven turbulence (REGIME 2)
! ------------------------------------------
ELSE IF ((rib < r0_2) .AND. (rib > zero)) THEN
regime = 2
psim = ( -five * rib ) * gzsoz0 / (r1_1 - five*rib)
psimz = ( -five * rib ) * gz10oz0 / (r1_1 - five*rib)
psim2 = ( -five * rib ) * gz2oz0 / (r1_1 - five*rib)
psim = max(psim,-r10)
psimz = max(psimz,-r10)
psim2 = max(psim2,-r10)
psih = psim
psihz = psimz
psih2 = psim2
! 6.3 Unstable Forced convection (REGIME 3)
! -------------------------------------
ELSE IF ((rib == zero) .or. (rib<zero .and. thvs2>thvs)) THEN
regime = 3
psim = zero
psimz = zero
psim2 = zero
psih = psim
psihz = psimz
psih2 = psim2
! 6.4 Free convection (REGIME 4)
! --------------------------
ELSE
regime = 4
! Calculate psi m and pshi h using iteration method
psim = zero
psih = zero
cc = two * atan(one)
! do k = 1 , k_iteration
! 6.4.1 Calculate ust, m/L (mol), h/L (hol)
! --------------------------
! Friction speed
ust = k_kar * sqrt(v2) /( gzsoz0 - psim)
! Heat flux factor
mol = k_kar * (ths - thg )/( gzsoz0 - psih)
! Ratio of PBL height to Monin-Obukhov length
if ( ust < r0_01 ) then
hol = rib * gzsoz0
else
hol = k_kar * grav * hs * mol / ( ths * ust * ust )
endif
! 6.4.2 Calculate n, nz, R, Rz
! --------------------------
hol = min(hol,zero)
hol = max(hol,-r10)
holz = (h10 / hs) * hol
holz = min(holz,zero)
holz = max(holz,-r10)
hol2 = (h2 / hs) * hol
hol2 = min(hol2,zero)
hol2 = max(hol2,-r10)
! 6.4.3 Calculate Psim & psih
! --------------------------
! Using the look-up table:
! nn = int( -r100 * hol )
! rr = ( -r100 * hol ) - nn
! psim = psimtb(nn) + rr * ( psimtb(nn+1) - psimtb(nn))
! psih = psihtb(nn) + rr * ( psihtb(nn+1) - psihtb(nn))
! Using the continuous function:
xx = (one - r16 * hol) ** quarter
yy = log((one+xx*xx)/two)
psim = two * log((one+xx)/two) + yy - two * atan(xx) + cc
psih = two * yy
! Using the look-up table:
! nz = int( -r100 * holz )
! rz = ( -r100 * holz ) - nz
! psimz = psimtb(nz) + rz * ( psimtb(nz+1) - psimtb(nz))
! psihz = psihtb(nz) + rz * ( psihtb(nz+1) - psihtb(nz))
! Using the continuous function:
xx = (one - r16 * holz) ** quarter
yy = log((one+xx*xx)/two)
psimz = two * log((one+xx)/two) + yy - two * atan(xx) + cc
psihz = two * yy
! Using the look-up table:
! n2 = int( -r100 * hol2 )
! r2 = ( -r100 * hol2 ) - n2
! psim2 = psimtb(n2) + r2 * ( psimtb(n2+1) - psimtb(n2))
! psih2 = psihtb(n2) + r2 * ( psihtb(n2+1) - psihtb(n2))
! Using the continuous function:
xx = (one - r16 * hol2) ** quarter
yy = log((one+xx*xx)/two)
psim2 = two * log((one+xx)/two) + yy - two * atan(xx) + cc
psih2 = two * yy
! enddo
! 6.4.4 Define the limit value for psim & psih
! --------------------------
psim = min(psim,r0_9*gzsoz0)
psimz = min(psimz,r0_9*gz10oz0)
psim2 = min(psim2,r0_9*gz2oz0)
psih = min(psih,r0_9*gzsoz0)
psihz = min(psihz,r0_9*gz10oz0)
psih2 = min(psih2,r0_9*gz2oz0)
ENDIF ! Regime
! 7. CALCULATE PSI FOR WIND, TEMPERATURE AND MOISTURE
! =======================================
psiw = gzsoz0 - psim
psiz = gz10oz0 - psimz
psit = gzsoz0 - psih
psit2 = gz2oz0 - psih2
! Friction speed
ust = k_kar * sqrt(v2) /( gzsoz0 - psim)
psiq = log(k_kar*ust*hs/ka + hs / zq0 ) - psih
psiq2 = log(k_kar*ust*h2/ka + h2 / zq0 ) - psih2
! 8. CALCULATE 10M WIND, 2M TEMPERATURE AND MOISTURE
! =======================================
u10 = us * psiz / psiw
v10 = vs * psiz / psiw
f10 = psiz / psiw
t2 = ( thg + ( ths - thg )*psit2/psit)*(psfc/r1000)**rd_over_cp
q2 = qg + (qs - qg)*psiq2/psiq
if(iqtflg)then
t2 = t2 * (one + fv * q2)
end if
return
END SUBROUTINE SFC_WTQ_FWD
SUBROUTINE DA_TP_To_Qs( t, p, es, qs )
!$$$ subprogram documentation block
! . . . .
! subprogram: DA_TP_To_Qs
!
! prgrmmr:
!
! abstract: Convert T/p to saturation specific humidity.
!
! METHOD: qs = es_alpha * es / ( p - ( 1 - rd_over_rv ) * es ).
! Use Rogers & Yau (1989) formula: es = a exp( bTc / (T_c + c) ).
!
! program history log:
! 2000-10-03 Barker - Creation of F90 version.
! 2008-04-14 safford - added standard documentation block
!
! input argument list:
! t - Temperature.
! p - Pressure.
!
! output argument list:
! es - Sat. vapour pressure.
! qs - Sat. specific humidity.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind
use constants, only: eps,omeps,t0c,r0_01
IMPLICIT NONE
REAL(r_kind), INTENT(IN ) :: t ! Temperature.
REAL(r_kind), INTENT(IN ) :: p ! Pressure.
REAL(r_kind), INTENT( OUT) :: es ! Sat. vapour pressure.
REAL(r_kind), INTENT( OUT) :: qs ! Sat. specific humidity.
! Saturation Vapour Pressure Constants(Rogers & Yau, 1989)
REAL(r_kind), PARAMETER :: es_alpha = 611.2_r_kind
REAL(r_kind), PARAMETER :: es_beta = 17.67_r_kind
REAL(r_kind), PARAMETER :: es_gamma = 243.5_r_kind
REAL(r_kind) :: t_c ! T in degreesC.
!------------------------------------------------------------------------------
! [1.0] Initialise:
!------------------------------------------------------------------------------
t_c = t - t0c
!------------------------------------------------------------------------------
! [2.0] Calculate saturation vapour pressure:
!------------------------------------------------------------------------------
es = r0_01 * es_alpha * exp( es_beta * t_c / ( t_c + es_gamma ) ) ! change to mb
!------------------------------------------------------------------------------
! [3.0] Calculate saturation specific humidity:
!------------------------------------------------------------------------------
qs = eps * es / ( p - omeps * es )
return
END SUBROUTINE DA_TP_To_Qs
SUBROUTINE sfc_wtq_Lin(psfc_in, tg, ps_in, tvs, qs, us, vs, regime, &
psfc_prime_in,tg_prime,sig1,ts_in,qs_prime, &
us_prime,vs_prime, hs,roughness,iland, &
u10_prime,v10_prime,t2_prime,q2_prime,iqtflg)
!$$$ subprogram documentation block
! . . . .
! subprogram: sfc_wtq_Lin
!
! prgrmmr:
!
! abstract: Calculate the 10m wind, 2m (virtual) temperature and moisture based on the
! similarity theory/
!
! Reference:
! ---------
! Detail MM5/physics/pbl_sfc/mrfpbl/MRFPBL.F
!
!
! program history log:
! NOTE: changed code so input layer temps are virtual, and are internally
! converted to sensible as required. this corrects error where
! ths is potential temperature for lowest half sigma layer, but
! before this change is in fact potential virtual temperature.
! also, output t2 should be potential temperature, so changed t2 to tv2
! 2008-04-14 safford - added standard documentation block, rm unused vars
!
! Input argument list(basic state):
!
! psfc_in, tg : surface pressure and ground temperature
! ps_in, tvs, qs, us, vs, hs : model variable at lowlest half sigma leve1
! regime : PBL regime
! iqtflg : if true return virtual temperature
! : otherwise return sensible temperature
!
! Input argument list(pertubation):
!
! psfc_prime_in, tg_prime : Surface pressure and ground temperature
! ps_prime, tvs_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 argument list(pertubation):
!
! u10_prime, v10_prime : 10-m high observed wind components
! tv2_prime , q2_prime : 2-m high observed temperature and mixing ratio
!
! attributes:
! language: f90
! machine ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: grav,fv,rd_over_cp,zero,quarter,half,one,two,four,five,r1000,r10,r100,r0_01
IMPLICIT NONE
INTEGER(i_kind), INTENT (in ) :: regime,iland
REAL(r_kind) , INTENT (in ) :: ps_in , tvs , qs , us, vs, psfc_in, tg ! change ts to tvs
REAL(r_kind) , INTENT (in ) :: sig1, ts_in, qs_prime , & ! change ts_prime to tvs_prime
us_prime, vs_prime, psfc_prime_in, tg_prime
REAL(r_kind) , INTENT (in ) :: hs, roughness
REAL(r_kind) , INTENT ( out) :: u10_prime, v10_prime, t2_prime, q2_prime
logical , intent (in ) :: iqtflg
! Maximum number of iterations in computing psim, psih
! INTEGER(i_kind), 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(r_kind), PARAMETER :: h10 = 10.0_r_kind
real(r_kind), parameter :: h2 = two
!
! Default roughness over the land
REAL(r_kind), PARAMETER :: zint0 = 0.01_r_kind
!
! Von Karman constant
REAL(r_kind), PARAMETER :: k_kar = 0.4_r_kind
!
! Working variables
real(r_kind) :: psfc,ps
REAL(r_kind) :: Vc2, Va2, V2
REAL(r_kind) :: rib, xx, yy, cc, Pi
REAL(r_kind) :: psiw, psiz, mol, ust, hol, holz, hol2
REAL(r_kind) :: psim, psimz, psim2, psih, psihz, psih2
REAL(r_kind) :: psit, psit2, psiq, psiq2
REAL(r_kind) :: gzsoz0, gz10oz0, gz2oz0
REAL(r_kind) :: eg, qg, tvg, ts ! change tvs to ts
REAL(r_kind) :: ths, thg, thvs, thvg
REAL(r_kind) :: zq0, z0
real(r_kind) ps_prime,psfc_prime
REAL(r_kind) :: Vc2_prime, Va2_prime, V2_prime
REAL(r_kind) :: rib_prime, xx_prime, yy_prime
REAL(r_kind) :: psiw_prime, psiz_prime, mol_prime, ust_prime, &
hol_prime, holz_prime, hol2_prime
REAL(r_kind) :: psim_prime, psimz_prime, psim2_prime, &
psih_prime, psihz_prime, psih2_prime
REAL(r_kind) :: psit_prime, psit2_prime, &
psiq_prime, psiq2_prime
REAL(r_kind) :: qg_prime, tvg_prime, ts_prime, tvs_prime
REAL(r_kind) :: ths_prime, thg_prime, thvs_prime, thvg_prime ! add t2_prime
real(r_kind) t2,q2
REAL(r_kind), PARAMETER :: ka = 2.4E-5_r_kind
! local
real(r_kind),parameter :: r16 = 16.0_r_kind
real(r_kind),parameter :: r5_5 = 5.5_r_kind
real(r_kind),parameter :: r1_1 = 1.1_r_kind
real(r_kind),parameter :: r0_9 = 0.9_r_kind
real(r_kind),parameter :: r0_75 = 0.75_r_kind
!-----------------------------------------------------------------------------!
! convert input pressure variables from cb to mb
psfc = r10*psfc_in
ps = r10*ps_in
psfc_prime = r10*psfc_prime_in
! 1 Compute the roughness length based upon season and land use
! =====================================
! 1.1 Define the rouhness length
! -----------------
z0 = roughness/r100 ! unit change : originally cm --> m
if( z0 < 0.0001_r_kind) z0 = 0.0001_r_kind
! 1.2 Define the rouhgness length for moisture
! -----------------
if ( iland == 0 ) then
zq0 = z0
else
zq0 = zint0
endif
! 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 * (one + fv * qs) + fv * ts * qs_prime
! tvs = ts * (one + fv * qs)
ts = tvs / (one + fv * qs) ! input now tvs, and need also ts
if(iqtflg)then
ts_prime = (ts_in + fv*(qs*ts_in-tvs*qs_prime))/(one+fv*qs)**2
tvs_prime=ts_in
else
ts_prime = ts_in
tvs_prime = (ts_in*(one+fv*qs)**2+fv*tvs*qs_prime)/(one+fv*qs)
end if
! 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_Lin( tg, psfc, eg, tg_prime, psfc_prime, qg_prime )
qg_prime = qg_prime * qg
tvg_prime = tg_prime * (one + fv * qg) + fv * tg * qg_prime
tvg = tg * (one + fv * qg)
! 3. Compute the potential temperature and virtual potential temperature
! =======================================================================
! 3.1 Potential temperature on the lowest half sigma level
! ----------------------------------------------------
ps_prime=sig1*psfc_prime
Pi = (r1000 / ps) ** rd_over_cp
ths_prime = (ts_prime - ps_prime * rd_over_cp * ts/ps ) * Pi
ths = ts * Pi
! 3.2 Virtual potential temperature on the lowest half sigma level
! ------------------------------------------------------------
thvs_prime = (tvs_prime - ps_prime * rd_over_cp * tvs/ps) * Pi
thvs = tvs * Pi
! 3.3 Potential temperature at the ground
! -----------------------------------
Pi = (r1000 / psfc) ** rd_over_cp
thg_prime = (tg_prime - psfc_prime * rd_over_cp * tg/psfc) * Pi
thg = tg * Pi
! 3.4 Virtual potential temperature at ground
! ---------------------------------------
thvg_prime = (tvg_prime - psfc_prime * rd_over_cp * tvg/psfc) * Pi
thvg = tvg * Pi
! 4. BULK RICHARDSON NUMBER AND MONI-OBUKOV LENGTH
! =================================================
! 4.1 Velocity
! --------
!
! Wind speed:
Va2_prime = two*us*us_prime + two*vs*vs_prime
Va2 = us*us + vs*vs
!
! Convective velocity:
if ( thvg >= thvs ) then
Vc2_prime = four * (thvg_prime - thvs_prime)
Vc2 = four * (thvg - thvs)
else
Vc2_prime = zero
Vc2 = zero
endif
!
V2_prime = Va2_prime+ Vc2_prime
V2 = .000001_r_kind + Va2 + Vc2
! 4.2 Bulk richardson number
! ----------------------
Pi = grav * 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
! =======================================
SELECT CASE ( regime )
! 5.1 Stable conditions (REGIME 1)
! ---------------------------
CASE ( 1 );
psim_prime = zero
psimz_prime = zero
psim2_prime = zero
psim = -r10*gzsoz0
psimz = -r10*gz10oz0
psim2 = -r10*gz2oz0
psim = max(psim,-r10)
psimz = max(psimz,-r10)
psim2 = max(psim2,-r10)
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 = - one / ((r1_1 - five*rib)*(r1_1 - five*rib))
psim_prime = r5_5 * gzsoz0 * rib_prime * Pi
psimz_prime = r5_5 * gz10oz0 * rib_prime * Pi
psim2_prime = r5_5 * gz2oz0 * rib_prime * Pi
Pi = ( -five * rib ) / (r1_1 - five*rib)
psim = gzsoz0 * Pi
psimz = gz10oz0 * Pi
psim2 = gz2oz0 * Pi
if ( psim >= -r10 ) then
psim = psim
psim_prime = psim_prime
else
psim = -r10
psim_prime = zero
endif
if ( psimz >= -r10 ) then
psimz = psimz
psimz_prime = psimz_prime
else
psimz = -r10
psimz_prime = zero
endif
if ( psim2 >= -r10 ) then
psim2 = psim2
psim2_prime = psim2_prime
else
psim2 = -r10
psim2_prime = zero
endif
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 = zero
psimz_prime = zero
psim2_prime = zero
psim = zero
psimz = zero
psim2 = zero
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 = zero
psih_prime = zero
psim = zero
psih = zero
cc = two * atan(one)
! 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 = (half/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 < r0_01 ) then
hol_prime = rib_prime * gzsoz0
hol = rib * gzsoz0
else
hol = k_kar * grav * hs * mol / ( ths * ust * ust )
hol_prime = ( mol_prime / mol - ths_prime / ths &
- two* ust_prime / ust ) * hol
endif
! 5.4.2 Calculate n, nz, R, Rz
! --------------------------
if ( hol >= zero ) then
hol_prime = zero
hol = zero
else
hol_prime = hol_prime
hol = hol
endif
if ( hol >= -r10 ) then
hol_prime = hol_prime
hol = hol
else
hol_prime = zero
hol = -r10
endif
holz_prime = (h10 / hs) * hol_prime
holz = (h10 / hs) * hol
if ( holz >= zero ) then
holz_prime = zero
holz = zero
else
holz_prime = holz_prime
holz = holz
endif
if ( holz >= -r10 ) then
holz_prime = holz_prime
holz = holz
else
holz_prime = zero
holz = -r10
endif
hol2_prime = (h2 / hs) * hol_prime
hol2 = (h2 / hs) * hol
if ( hol2 >= zero ) then
hol2_prime = zero
hol2 = zero
else
hol2_prime = hol2_prime
hol2 = hol2
endif
if ( hol2 >= -r10 ) then
hol2_prime = hol2_prime
hol2 = hol2
else
hol2_prime = zero
hol2 = -r10
endif
! 5.4.3 Calculate Psim & psih
! --------------------------
! Using the continuous function:
xx_prime = -four* hol_prime /((one - r16 * hol) ** r0_75)
xx = (one - r16 * hol) ** quarter
yy_prime = two* xx * xx_prime /(one+xx*xx)
yy = log((one+xx*xx)/two)
psim_prime = 2 * xx_prime *(one/(one+xx)-one/(1+xx*xx)) + yy_prime
psim = two * log((one+xx)/two) + yy - two * atan(xx) + cc
psih_prime = two * yy_prime
psih = two * yy
! Using the continuous function:
xx_prime = -four* holz_prime /((one - r16 * holz) ** r0_75)
xx = (one - r16 * holz) ** quarter
yy_prime = two* xx * xx_prime /(one+xx*xx)
yy = log((one+xx*xx)/two)
psimz_prime = two* xx_prime *(one/(one+xx)-one/(1+xx*xx)) + yy_prime
psimz = two * log((one+xx)/two) + yy - two * atan(xx) + cc
psihz_prime = two * yy_prime
psihz = two * yy
! Using the continuous function:
xx_prime = -four* hol2_prime /((one - r16 * hol2) ** r0_75)
xx = (one - r16 * hol2) ** quarter
yy_prime = two* xx * xx_prime /(one+xx*xx)
yy = log((one+xx*xx)/two)
psim2_prime = two* xx_prime *(one/(one+xx)-one/(1+xx*xx)) + yy_prime
psim2 = two * log((one+xx)/two) + yy - two * atan(xx) + cc
psih2_prime = two * yy_prime
psih2 = two * yy
! enddo
! 5.4.4 Define the limit value for psim & psih
! --------------------------
if ( psim <= r0_9*gzsoz0 ) then
psim_prime = psim_prime
psim = psim
else
psim = r0_9*gzsoz0
psim_prime = zero
endif
if ( psimz <= r0_9*gz10oz0 ) then
psimz_prime = psimz_prime
psimz = psimz
else
psimz_prime = zero
psimz = r0_9*gz10oz0
endif
if ( psim2 <= r0_9*gz2oz0 ) then
psim2_prime = psim2_prime
psim2 = psim2
else
psim2_prime = zero
psim2 = r0_9*gz2oz0
endif
if ( psih <= r0_9*gzsoz0 ) then
psih_prime = psih_prime
psih = psih
else
psih_prime = zero
psih = r0_9*gzsoz0
endif
if ( psihz <= r0_9*gz10oz0 ) then
psihz_prime = psihz_prime
psihz = psihz
else
psihz_prime = zero
psihz = r0_9*gz10oz0
endif
if ( psih2 <= r0_9*gz2oz0 ) then
psih2_prime = psih2_prime
psih2 = psih2
else
psih2_prime = zero
psih2 = r0_9*gz2oz0
endif
CASE DEFAULT;
write(unit=*, fmt='(/a,i2,a/)') "Regime=",regime," is invalid."
stop "sfc_wtq_Lin"
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 = (half/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 = ( (one-psit2/psit) * thg_prime + ( ths_prime + &
(ths - thg)/psit2 * psit2_prime - &
(ths - thg)/psit * psit_prime ) * psit2/psit &
+ rd_over_cp*( thg + ( ths - thg )*psit2/psit)/psfc * psfc_prime ) &
* (psfc/r1000)**rd_over_cp
q2_prime = (one-psiq2/psiq) * qg_prime + psiq2/psiq * qs_prime + &
(qs -qg)*(psiq2/psiq) * (psiq2_prime/psiq2 - psiq_prime/psiq)
t2 = ( thg + ( ths - thg )*psit2/psit)*(psfc/r1000)**rd_over_cp
q2 = qg + (qs - qg)*psiq2/psiq
if(iqtflg)then
t2_prime=t2_prime*(one+fv*q2)+fv*t2*q2_prime
end if
END SUBROUTINE sfc_wtq_Lin
SUBROUTINE DA_TP_To_Qs_Lin( t, p, es, t_prime, p_prime, &
qs_prime_over_qs )
!$$$ subprogram documentation block
! . . . .
! subprogram: DA_TP_To_Qs_Lin
!
! prgrmmr:
!
! abstract: Convert es/p/es_prime to saturation specific humidity increment.
! METHOD: qs~ = qs * ( p es'/es - p' ) / ( p - (1-rd_over_rv) es ).
! Use Rogers & Yau (1989) formula: es = a exp( bTc / (T_c + c) ).
!
! program history log:
! 2000-10-03 Barker - Creation of F90 version.
! 2008-04-14 safford - add standard documentation block, rm unused uses
!
! input argument list:
! t - temperature
! p - pressure
! es - Sat. vapour pressure
! t_prime - temperature increment
! p_prime - pressure increment
!
! output argument list:
! qs_prime_over_qs - qs~/qs
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds,only: r_kind
use constants, only: omeps,t0c
IMPLICIT NONE
REAL(r_kind), INTENT(IN ) :: t ! Temperature.
REAL(r_kind), INTENT(IN ) :: p ! Pressure.
REAL(r_kind), INTENT(IN ) :: es ! Sat. vapour pressure.
REAL(r_kind), INTENT(IN ) :: t_prime ! Temperature increment.
REAL(r_kind), INTENT(IN ) :: p_prime ! Pressure increment.
REAL(r_kind), INTENT( OUT) :: qs_prime_over_qs ! qs~/qs.
! Saturation Vapour Pressure Constants(Rogers & Yau, 1989)
! REAL(r_kind), PARAMETER :: es_alpha = 611.2_r_kind
REAL(r_kind), PARAMETER :: es_beta = 17.67_r_kind
REAL(r_kind), PARAMETER :: es_gamma = 243.5_r_kind
REAL(r_kind), PARAMETER :: es_gammakelvin = es_gamma-t0c
REAL(r_kind), PARAMETER :: es_gammabeta = es_gamma*es_beta
REAL(r_kind) :: temp ! Temporary value.
REAL(r_kind) :: es_prime_over_es ! es~/es
!------------------------------------------------------------------------------
! [1.0] Initialise:
!------------------------------------------------------------------------------
temp = t + es_gammakelvin
!------------------------------------------------------------------------------
! [2.0] Calculate saturation vapour pressure increment:
!------------------------------------------------------------------------------
es_prime_over_es = es_gammabeta * t_prime / ( temp * temp )
!------------------------------------------------------------------------------
! [3.0] Calculate saturation specific humidity increment:
!------------------------------------------------------------------------------
qs_prime_over_qs = ( p * es_prime_over_es - p_prime ) / &
( p - omeps * es )
END SUBROUTINE DA_TP_To_Qs_Lin
subroutine get_tlm_tsfc(tlm_tsfc,psges2,tgges,prsltmp2, &
tvtmp,qtmp,utmp,vtmp,hsges,roges,msges,regime,iqtflg)
!$$$ subprogram documentation block
! . . . .
! subprogram: get_tlm_tsfc
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-14 safford - add documentation block
!
! input argument list:
! psges2 -
! tgges -
! prsltmp2 -
! tvtmp -
! qtmp -
! utmp -
! vtmp -
! hsges -
! roges -
! msges -
! regime -
! iqtflg -
!
! output argument list:
! tlm_tsfc -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: zero,one
implicit none
real(r_kind) ,intent( out) :: tlm_tsfc(6)
real(r_kind) ,intent(in ) :: psges2,tgges,prsltmp2,tvtmp,qtmp,utmp,vtmp,hsges,roges
integer(i_kind),intent(in ) :: regime,msges
logical ,intent(in ) :: iqtflg
real(r_kind) perturb(6),u10_prime,v10_prime,q2_prime
real(r_kind) sig1
integer(i_kind) i
sig1=prsltmp2/psges2
do i=1,6
! 1 -- psfc_prime (perturbation of psfc, psfc in cb)
! 2 -- tg_prime
! 3 -- ts_prime (sensible or virtual depending on iqtflg)
! 4 -- qs_prime
! 5 -- us_prime
! 6 -- vs_prime
perturb=zero
perturb(i)=one
call sfc_wtq_lin(psges2,tgges,prsltmp2,tvtmp,qtmp,utmp,vtmp,regime, &
perturb(1),perturb(2),sig1,perturb(3),perturb(4),perturb(5),perturb(6), &
hsges,roges,msges,u10_prime,v10_prime,tlm_tsfc(i),q2_prime,iqtflg)
end do
end subroutine get_tlm_tsfc