<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_SFCPRS'><A href='../../html_code/tools/da_sfcprs.inc.html#DA_SFCPRS' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_sfcprs (kts, kte, p, t, q, height, psfc, pslv, h_obs, t_obs) 2,2
! Purpose: to compute the pressure at obs height (psfc) from the sea
! level pressure (pslv), obs height (h_obs), temperature
! (t_obs, optional), and the model profile (p,t,q).
implicit none
integer, intent(in) :: kts, kte
real, intent(in) :: height(kts:kte)
real, intent(in) :: p(kts:kte)
real, intent(in) :: t(kts:kte)
real, intent(in) :: q(kts:kte)
real, intent(in) :: h_obs, pslv
real, intent(out) :: psfc
real, optional, intent(in) :: t_obs
real, parameter :: g = gravity
real, parameter :: r = gas_constant
real, parameter :: rov2 = r / 2.0
real, parameter :: gamma = 6.5e-3 ! temperature lapse rate
real, parameter :: gammarg= gamma*gas_constant/g
real :: plmb(3)
integer :: k
integer :: k500
integer :: k700
integer :: k850
logical :: l1
logical :: l2
logical :: l3
real :: gamma78, gamma57
real :: ht
real :: p1
real :: t1
real :: t500
real :: t700
real :: t850
real :: tc
real :: tfixed
real :: tslv, tsfc
if (trace_use_dull) call da_trace_entry
("da_sfcprs")
TC = t_kelvin + 17.5
PLMB = (/85000.0, 70000.0, 50000.0/)
! Find the locations of the 850, 700 and 500 mb levels.
101 continue
K850 = 0 ! FinD K AT: P=850
K700 = 0 ! P=700
K500 = 0 ! P=500
do K = kte-1, kts, -1
if (p(k) > PLMB(1) .and. K850==0) then
K850 = K + 1
else if (p(k) > PLMB(2) .and. K700==0) then
K700 = K + 1
else if (p(k) > PLMB(3) .and. K500==0) then
K500 = K + 1
end if
end do
if ((K850 == 0) .OR. (K700 == 0) .OR. (K500 == 0)) then
! write(unit=message(1),fmt='(A,3f8.0,A)') &
! 'Error in finding p level for ',PLMB(1), PLMB(2), PLMB(3),' Pa.'
! do K = 1, KX
! write(unit=message(k+1),fmt='(A,I3,A,F10.2)') 'K = ',K,' PRESSURE = ',P(K)
! end do
! write(unit=message(kx+2),fmt='(A,2f8.0,A,f8.0,A)) ','Expected ', &
! PLMB(1), PLMB(2),' and ',PLMB(3),' Pa. values, at least.'
! message(kx+3)="not enough levels"
! message(kx+4)="Change the pressure levels by -25mb"
! call da_error(__FILE__,__LINE__,message(1:kx+4))
PLMB(1) = PLMB(1) - 2500.0
PLMB(2) = PLMB(2) - 2500.0
PLMB(3) = PLMB(3) - 2500.0
goto 101
end if
! The 850 hPa level is called something special, and interpolated
! to cross points. And then, we fill those last rows and columns.
ht = height(k850)
! The variable ht is now -h_obs/ht(850 hPa). The plot thickens.
ht = -h_obs / ht
! Make an isothermal assumption to get a first guess at the surface
! pressure. This is to tell us which levels to use for the lapse
! rates in a bit.
psfc = pslv * (pslv / p(k850)) ** ht
! Get a pressure more than 100 hPa above the surface - P1. The
! P1 is the top of the level that we will use for our lapse rate
! computations.
if ((PSFC - p(k850) - 10000.0) >= 0.0) then
P1 = p(k850)
else if ((PSFC - p(k700)) >= 0.0) then
P1 = PSFC - 10000.0
else
P1 = p(k500)
end if
! Compute virtual temperatures for k850, k700, and k500 layers. Now
! you see WHY? we wanted Q on pressure levels, it all is beginning
! to make sense.
t850 = t(k850) * (1.0 + 0.608 * q(k850))
t700 = t(k700) * (1.0 + 0.608 * q(k700))
t500 = t(k500) * (1.0 + 0.608 * q(k500))
! Compute two lapse rates between these three levels. These are
! environmental values for each (i,j).
gamma78 = LOG(t850 / t700) / LOG (p(k850) / p(k700))
gamma57 = LOG(t700 / t500) / LOG (p(k700) / p(k500))
if ((psfc - p(k850) - 10000.0) >= 0.0) then
t1 = t850
else if ((psfc - p(k850)) >= 0.0) then
t1 = t700 * (p1 / p(k700)) ** gamma78
else if ((psfc - p(k700)) >= 0.0) then
t1 = t500 * (p1 / p(k500)) ** gamma57
else
t1 = t500
end if
! From our temperature way up in the air, we extrapolate down to
! the sea level to get a guess at the sea level temperature.
tslv = t1 * (pslv / p1) ** (gammarg)
! The new surface temperature is computed from the with new sea level
! temperature, just using the elevation and a lapse rate. This lapse
! rate is -6.5 K/km.
if (present (t_obs)) then
tsfc = t_obs
else
tsfc = tslv - gamma * h_obs
end if
! A correction to the sea-level temperature, in case it is too warm.
TFIXED = TC - 0.005 * (TSFC - TC) ** 2
L1 = tslv .LT. tc
L2 = tsfc .LE. tc
L3 = .NOT. L1
if (L2 .AND. L3) then
tslv = tc
else if ((.NOT. L2) .AND. L3) then
tslv = tfixed
end if
! Finally, we can get to the surface pressure.
p1 = -h_obs * g / (rov2 * (tsfc + tslv))
psfc = pslv * EXP(p1)
! Surface pressure and sea-level pressure are the same at sea level.
if (ABS (h_obs) < 0.1) psfc = pslv
if (trace_use_dull) call da_trace_exit("da_sfcprs")
end subroutine da_sfcprs