subroutine compute_fact10(u,v,t,q,ps,prsi1,prsi2,skint,z0rl,islimsk,f10m)
!$$$ subprogram documentation block
! . . . .
! subprogram: compute_fact10 compute 10m wind factor
! prgmmr: treadon org: np23 date: 2006-09-28
!
! abstract: Use GFS surface physics routines to compute 10m wind factor
!
! program history log:
! 2006-09-28 treadon - initial routine
! 2008-06-05 safford - rm unused vars and uses, comment out unused params
!
! input argument list:
! u - u wind component (2d field, 1st model layer)
! v - v wind component
! t - sensible temperature
! q - specific humidity
! ps - surface pressure
! prsi1 - surface pressure pressure
! prsi2 - interface pressure at top of first layer
! skint - skin temperature
! z0rl - surface roughness
! islimsk - land/sea/ice mask
!
! output argument list:
! f10m - 10m wind factor
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: grav,zero,quarter,half,one,two,four,five,&
fv,rd,rd_over_cp,r1000
implicit none
! Passed Variables
real(r_kind) ,intent(in ) :: u,v,t,q,ps,skint,z0rl
integer(i_kind),intent(in ) :: islimsk
real(r_kind) ,intent( out) :: f10m
real(r_kind) ,intent(in ) :: prsi1,prsi2
! Local Variables
real(r_kind):: prsl,prkl
real(r_kind):: prki1,prki2
real(r_kind):: q0,tem,del,rkap,rkapi,rkapp1
real(r_kind)::ustar,wind,z0,z0max,ztmax,tv1,z1, &
rat,thv1,theta1,tvs,dtv,rb,fm,fh,hlinf, &
hl1,pm,ph,pm10,hl12,ph2,fm10
real(r_kind):: psurf,ps1
real(r_kind) restar,aa0,bb0,fhs,fms,hl0,hlt,adtv,bb,aa,hl0inf,hltinf,&
hl110,olinf
! Local Parameters
real(r_kind),parameter:: charnok=0.014_r_kind
! real(r_kind),parameter:: ca=0.4_r_kind
real(r_kind),parameter:: alpha=five
real(r_kind),parameter:: a0=-3.975_r_kind
real(r_kind),parameter:: a1=12.32_r_kind
real(r_kind),parameter:: b1=-7.755_r_kind
real(r_kind),parameter:: b2=6.041_r_kind
real(r_kind),parameter:: a0p=-7.941_r_kind
real(r_kind),parameter:: a1p=24.75_r_kind
real(r_kind),parameter:: b1p=-8.705_r_kind
real(r_kind),parameter:: b2p=7.899_r_kind
real(r_kind),parameter:: vis=1.4e-5_r_kind
! real(r_kind),parameter:: aa1=-1.076_r_kind
! real(r_kind),parameter:: bb1=0.7045_r_kind
! real(r_kind),parameter:: cc1=-0.05808_r_kind
! real(r_kind),parameter:: bb2=-0.1954_r_kind
! real(r_kind),parameter:: cc2=0.009999_r_kind
! real(r_kind),parameter:: rnu=1.51e-5_r_kind
! real(r_kind),parameter:: arnu=0.135_r_kind*rnu
real(r_kind),parameter:: ten=10.0_r_kind
! INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIED
! PSURF IS IN PASCALS
! WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1
! SURFACE ROUGHNESS LENGTH IS CONVERTED TO M FROM CM
rkap = rd_over_cp
RKAPI = one / RKAP
RKAPP1 = one + RKAP
rat=zero ; restar=zero ; ustar=zero
del = prsi1-prsi2
tem = rkapp1*del
prki1 = (prsi1*0.01_r_kind)**rkap
prki2 = (prsi2*0.01_r_kind)**rkap
prkl = (prki1*prsi1-prki2*prsi2)/tem
prsl = 100.0_r_kind*prkl**rkapi
psurf=ps*r1000
ps1=prsl*r1000
wind= sqrt( u*u + v*v )
wind=max(wind,one)
q0=max(q,1.e-8_r_kind)
theta1=t*(prki1/prkl)
tv1 =t*(one+fv*q0)
thv1=theta1 *(one+fv*q0)
tvs =skint*(one+fv*q0)
z0=0.01_r_kind*z0rl
z1=-rd*tv1*log(ps1/psurf)/grav
! COMPUTE STABILITY DEPENDENT EXCHANGE COEFFICIENTS
if (islimsk == 0) then
ustar = sqrt(grav * z0 / charnok)
end if
! COMPUTE STABILITY INDICES (RB AND HLINF)
z0max = min(z0,one*z1)
ztmax = z0max
if (islimsk == 0) then
restar=ustar*z0max/vis
restar=max(restar,1.e-6_r_kind)
! Rat taken from Zeng, Zhao and Dickinson 1997
rat = 2.67_r_kind * restar**quarter - 2.57_r_kind
rat = min(rat,7.0_r_kind)
ztmax = z0max * exp(-rat)
end if
dtv = thv1 - tvs
adtv = abs(dtv)
adtv = max(adtv,1.e-3_r_kind)
dtv = sign(one,dtv)*adtv
rb = grav*dtv*z1 / (half*(thv1+tvs)*wind*wind)
rb=max(rb,-5.e3_r_kind)
fm=log((z0max+z1) / z0max)
fh=log((ztmax+z1) / ztmax)
hlinf=rb*fm*fm / fh
fm10=log((z0max+ten) / z0max)
! STABLE CASE
if (dtv >= zero) then
hl1=hlinf
end if
if ((dtv >= zero) .AND. (hlinf > quarter)) then
hl0inf=z0max*hlinf/z1
hltinf=ztmax*hlinf/z1
aa=sqrt(one + four*alpha*hlinf)
aa0=sqrt(one + four*alpha*hl0inf)
bb=aa
bb0=sqrt(one + four*alpha*hltinf)
pm=aa0 - aa + log((aa+one)/(aa0+one))
ph=bb0 - bb + log((bb+one)/(bb0+one))
fms=fm-pm
fhs=fh-ph
hl1=fms*fms*rb/fhs
end if
! SECOND ITERATION
if (dtv >= zero) then
hl0=z0max*hl1/z1
hlt=ztmax*hl1/z1
aa=sqrt(one + four*alpha*hl1)
aa0=sqrt(one + four*alpha*hl0)
bb=aa
bb0=sqrt(one + four*alpha*hlt)
pm=aa0 - aa + log((aa+one)/(aa0+one))
ph=bb0 - bb + log((bb+one)/(bb0+one))
hl110=hl1*ten/z1
aa=sqrt(one + four*alpha*hl110)
pm10=aa0 - aa + log((aa+one)/(aa0+one))
hl12=hl1*two/z1
bb=sqrt(one + four*alpha*hl12)
ph2=bb0 - bb + log((bb+one)/(bb0+one))
end if
! UNSTABLE CASE
! CHECK FOR UNPHYSICAL OBUKHOV LENGTH
if (dtv < zero) then
olinf = z1/hlinf
if ( abs(olinf) <= z0max*50.0_r_kind ) then
hlinf = -z1/(50.0_r_kind*z0max)
end if
end if
! GET PM AND PH
if (dtv < zero .AND. hlinf >= (-half)) then
hl1=hlinf
pm=(a0 + a1*hl1)*hl1/(one + b1*hl1 + b2*hl1*hl1)
ph=(a0p + a1p*hl1)*hl1/(one + b1p*hl1 + b2*hl1*hl1)
hl110=hl1*ten/z1
pm10=(a0 + a1*hl110)*hl110/(one + b1*hl110 + b2*hl110*hl110)
hl12=hl1*two/z1
ph2=(a0p + a1p*hl12)*hl12/(one + b1p*hl12 + b2p*hl12*hl12)
end if
if (dtv < zero .AND. hlinf < (-half)) then
hl1=-hlinf
pm=log(hl1) + two*hl1**(-quarter) - 0.8776_r_kind
ph=log(hl1) + half*hl1**(-half) + 1.386_r_kind
hl110=hl1*ten/z1
pm10=log(hl110) + two*hl110**(-quarter) - 0.8776_r_kind
hl12=hl1*two/z1
ph2=log(hl12) + half*hl12**(-half) + 1.386_r_kind
end if
! FINISH THE EXCHANGE COEFFICIENT COMPUTATION TO PROVIDE FM AND FH
fm=fm-pm
fh=fh-ph
fm10=fm10-pm10
f10m=fm10/fm
return
end subroutine compute_fact10