var/da/da_physics/da_tpq_to_slp_lin.inc

   1 subroutine da_tpq_to_slp_lin ( T, Q, P, TERR, PSFC, T9, Q9, P9, PSFC9, SLP9)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: computes sea level pressure from the rule
   5    !              t1/t2=(p1/p2)**(gamma*r/g).
   6    !
   7    !     input       t        temperature
   8    !                 q        mixing ratio
   9    !                 p        pressure
  10    !                 terr     terrain
  11    !                 psfc     surface pressure
  12    !
  13    !     output      slp      sea level pressure
  14    !-----------------------------------------------------------------------
  15
  16    implicit none
  17
  18    real, intent(in)    :: terr, psfc, psfc9
  19    real, intent(in)    :: t(kms:kme)
  20    real, intent(in)    :: q(kms:kme)
  21    real, intent(in)    :: p(kms:kme)
  22    real, intent(in)    :: t9(kms:kme)
  23    real, intent(in)    :: q9(kms:kme)
  24    real, intent(in)    :: p9(kms:kme)
  25    ! real                :: slp
  26    real, intent(out)   :: slp9
  27
  28    integer :: k, klo, khi
  29    real    :: pl, t0, ts, xterm,tlo, thi, tl
  30    real    :: pl9,t09,ts9,tlo9,thi9,tl9,coef1,coef2
  31
  32    real, parameter      :: gamma  = 6.5e-3
  33    real, parameter      :: tc     = t_kelvin+17.5
  34    real, parameter      :: pconst = 10000.0
  35    real, parameter      :: eps    = 0.622
  36
  37    if (trace_use) call da_trace_entry("da_tpq_to_slp_lin")
  38
  39    !     ... sea level pressure
  40
  41    xterm=gamma* gas_constant / gravity
  42
  43    ! compute pressure at pconst mb above surface (pl)
  44
  45    if (terr <= 0.0) then
  46       slp9 = psfc9
  47       ! slp = psfc
  48       if (trace_use) call da_trace_exit("da_tpq_to_slp_lin")
  49       return
  50    end if
  51
  52
  53    pl9  = psfc9
  54    pl  = psfc - pconst
  55    klo = 0
  56
  57    ! find 2 levels on sigma surfaces surrounding pl at each i,j
  58
  59    k_loop: do k=kts, kte-1
  60       if ((p(k) >= pl) .and. (p(k+1) < pl)) then
  61          khi = k+1
  62          klo = k
  63          exit k_loop
  64       end if
  65    end do k_loop
  66
  67    if (klo < 1) then
  68       write(unit=message(1),fmt='(A,F11.3,A)') &
  69            'error finding pressure level ',pconst,' mb above the surface'
  70       write(unit=message(2),fmt='(A,F11.3,2X,A,F11.3)') 'PL=',PL,'  PSFC=',psfc
  71       call da_error(__FILE__,__LINE__,message(1:2))
  72    end if
  73
  74    ! get temperature at pl (tl), extrapolate t at surface (ts)
  75    ! and t at sea level (t0) with 6.5 k/km lapse rate
  76
  77    tlo9=t9(klo) * (eps+q(klo))/(eps*(1.0+q(klo))) + &
  78         q9(klo)*t(klo)*(1.0-eps)/(eps*(1.0+q(klo))**2)
  79    tlo=t(klo) * (eps+q(klo))/(eps*(1.0+q(klo)))
  80    thi9=t9(khi) * (eps+q(khi))/(eps*(1.0+q(khi)))+   &
  81         q9(khi)*t(khi)*(1.0-eps)/(eps*(1.0+q(khi))**2)
  82    thi=t(khi) * (eps+q(khi))/(eps*(1.0+q(khi)))
  83    coef1=alog(pl/p(khi))
  84    coef2=alog(p(klo)/p(khi))
  85    tl9=(1.0-coef1/coef2)*thi9+coef1/coef2*tlo9       &
  86        -(thi-tlo)/(coef2*pl)*pl9                 &
  87        +((thi-tlo)/(coef2*p(khi))*(1.0-coef1/coef2))*p9(khi)   &
  88        +(thi-tlo)*coef1/(coef2*coef2*p(klo))*p9(klo)
  89    tl=thi-(thi-tlo)*coef1/coef2
  90    ts9=tl9*(psfc/pl)**xterm+psfc9*xterm*(tl/pl)*(psfc/pl)**  &
  91        (xterm-1)-pl9*xterm*(tl*psfc/(pl*pl))*(psfc/pl)**(xterm-1)
  92    ts=tl*(psfc/pl)**xterm
  93    t09=ts9
  94    t0=ts +gamma*terr
  95
  96    ! correct sea level temperature if too hot
  97
  98    if ( t0 >= tc ) then
  99       if ( ts <= tc ) then
 100          t09 = 0.0
 101          t0 = tc
 102       else
 103          t09 = -0.01*(ts-tc)*ts9
 104          t0 = tc-0.005*(ts-tc)**2
 105       end if
 106    end if
 107
 108    ! compute sea level pressure
 109
 110    slp9=psfc9*exp(2.0*gravity*terr/(gas_constant*(ts+t0)))  &
 111           -psfc*exp(2.0*gravity*terr/(gas_constant*(ts+t0)))*  &
 112           2.0*gravity*terr/(gas_constant*(ts+t0)**2)*(ts9+t09)
 113    ! slp=psfc*exp(2.0*gravity*terr/(gas_constant*(ts+t0)))
 114
 115    if (trace_use) call da_trace_exit("da_tpq_to_slp_lin")
 116
 117 end subroutine da_tpq_to_slp_lin
 118
 119