var/da/da_tools/da_sfcprs.inc

   1 subroutine da_sfcprs (kts, kte, p, t, q, height, psfc, pslv, h_obs, t_obs)
   2
   3    ! Purpose: to compute the pressure at obs height (psfc) from the sea
   4    !          level pressure (pslv), obs height (h_obs), temperature
   5    !          (t_obs, optional), and the model profile (p,t,q).
   6
   7    implicit none
   8
   9    integer,        intent(in)  :: kts, kte
  10    real,           intent(in)  :: height(kts:kte)
  11    real,           intent(in)  :: p(kts:kte)
  12    real,           intent(in)  :: t(kts:kte)
  13    real,           intent(in)  :: q(kts:kte)
  14    real,           intent(in)  :: h_obs, pslv
  15    real,           intent(out) :: psfc
  16    real, optional, intent(in)  :: t_obs
  17
  18    real, parameter :: g = gravity
  19    real, parameter :: r =  gas_constant
  20    real, parameter :: rov2 = r / 2.0
  21    real, parameter :: gamma  = 6.5e-3   ! temperature lapse rate
  22    real, parameter :: gammarg= gamma*gas_constant/g
  23
  24    real    :: plmb(3)
  25
  26    integer :: k
  27    integer :: k500
  28    integer :: k700
  29    integer :: k850
  30
  31    logical :: l1
  32    logical :: l2
  33    logical :: l3
  34
  35    real    :: gamma78, gamma57
  36    real    :: ht
  37    real    :: p1
  38    real    :: t1
  39    real    :: t500
  40    real    :: t700
  41    real    :: t850
  42    real    :: tc
  43    real    :: tfixed
  44    real    :: tslv, tsfc
  45
  46    if (trace_use_dull) call da_trace_entry("da_sfcprs")
  47
  48    TC   = t_kelvin + 17.5
  49    PLMB = (/85000.0, 70000.0, 50000.0/)
  50
  51    !  Find the locations of the 850, 700 and 500 mb levels.
  52
  53 101 continue
  54    K850 = 0                              ! FinD K AT: P=850
  55    K700 = 0                              !            P=700
  56    K500 = 0                              !            P=500
  57
  58    do K = kte-1, kts, -1
  59       if      (p(k) > PLMB(1) .and. K850==0) then
  60          K850 = K + 1
  61       else if (p(k) > PLMB(2) .and. K700==0) then
  62          K700 = K + 1
  63       else if (p(k) > PLMB(3) .and. K500==0) then
  64          K500 = K + 1
  65       end if
  66
  67    end do
  68
  69    if ((K850 == 0) .OR. (K700 == 0) .OR. (K500 == 0)) then
  70       ! write(unit=message(1),fmt='(A,3f8.0,A)') &
  71       !    'Error in finding p level for ',PLMB(1), PLMB(2), PLMB(3),' Pa.'
  72       ! do K = 1, KX
  73       !    write(unit=message(k+1),fmt='(A,I3,A,F10.2)') 'K = ',K,'  PRESSURE = ',P(K)
  74       ! end do
  75       ! write(unit=message(kx+2),fmt='(A,2f8.0,A,f8.0,A)) ','Expected ',    &
  76       !    PLMB(1), PLMB(2),' and ',PLMB(3),' Pa. values, at least.'
  77       ! message(kx+3)="not enough levels"
  78       ! message(kx+4)="Change the pressure levels by -25mb"
  79       ! call da_error(__FILE__,__LINE__,message(1:kx+4))
  80       PLMB(1) = PLMB(1) - 2500.0
  81       PLMB(2) = PLMB(2) - 2500.0
  82       PLMB(3) = PLMB(3) - 2500.0
  83       goto 101
  84    end if
  85
  86    !  The 850 hPa level is called something special, and interpolated
  87    !  to cross points.  And then, we fill those last rows and columns.
  88
  89    ht = height(k850)
  90
  91    !  The variable ht is now -h_obs/ht(850 hPa).  The plot thickens.
  92
  93    ht = -h_obs / ht
  94
  95    !  Make an isothermal assumption to get a first guess at the surface
  96    !  pressure.  This is to tell us which levels to use for the lapse
  97    !  rates in a bit.
  98
  99    psfc = pslv * (pslv / p(k850)) ** ht
 100
 101    !  Get a pressure more than 100 hPa above the surface - P1.  The
 102    !  P1 is the top of the level that we will use for our lapse rate
 103    !  computations.
 104
 105    if      ((PSFC - p(k850) - 10000.0) >= 0.0) then
 106       P1 = p(k850)
 107    else if ((PSFC - p(k700)) >= 0.0) then
 108       P1 = PSFC - 10000.0
 109    else
 110       P1 = p(k500)
 111    end if
 112
 113    !  Compute virtual temperatures for k850, k700, and k500 layers.  Now
 114    !  you see WHY? we wanted Q on pressure levels, it all is beginning
 115    !  to make sense.
 116
 117    t850 = t(k850) * (1.0 + 0.608 * q(k850))
 118    t700 = t(k700) * (1.0 + 0.608 * q(k700))
 119    t500 = t(k500) * (1.0 + 0.608 * q(k500))
 120
 121    !  Compute two lapse rates between these three levels.  These are
 122    !  environmental values for each (i,j).
 123
 124    gamma78 = LOG(t850 / t700)  / LOG (p(k850) / p(k700))
 125    gamma57 = LOG(t700 / t500)  / LOG (p(k700) / p(k500))
 126
 127    if ((psfc - p(k850) - 10000.0) >= 0.0) then
 128       t1 = t850
 129    else if ((psfc - p(k850)) >= 0.0) then
 130       t1 = t700 * (p1 / p(k700)) ** gamma78
 131    else if ((psfc - p(k700)) >= 0.0) then
 132       t1 = t500 * (p1 / p(k500)) ** gamma57
 133    else
 134       t1 = t500
 135    end if
 136
 137    !  From our temperature way up in the air, we extrapolate down to
 138    !  the sea level to get a guess at the sea level temperature.
 139
 140    tslv = t1 * (pslv / p1) ** (gammarg)
 141
 142    !  The new surface temperature is computed from the with new sea level
 143    !  temperature, just using the elevation and a lapse rate.  This lapse
 144    !  rate is -6.5 K/km.
 145
 146    if (present (t_obs)) then
 147      tsfc = t_obs
 148    else
 149      tsfc = tslv - gamma * h_obs
 150    end if
 151
 152    !  A correction to the sea-level temperature, in case it is too warm.
 153
 154    TFIXED = TC - 0.005 * (TSFC - TC) ** 2
 155
 156    L1 = tslv .LT. tc
 157    L2 = tsfc .LE. tc
 158    L3 = .NOT. L1
 159    if (L2 .AND. L3) then
 160       tslv = tc
 161    else if ((.NOT. L2) .AND. L3) then
 162       tslv = tfixed
 163    end if
 164
 165    !  Finally, we can get to the surface pressure.
 166
 167    p1 = -h_obs * g / (rov2 * (tsfc + tslv))
 168    psfc = pslv * EXP(p1)
 169
 170    !  Surface pressure and sea-level pressure are the same at sea level.
 171
 172    if (ABS (h_obs) < 0.1) psfc = pslv
 173
 174    if (trace_use_dull) call da_trace_exit("da_sfcprs")
 175
 176 end subroutine da_sfcprs
 177
 178