var/da/da_physics/da_tpq_to_slp_adj.inc

   1 subroutine da_tpq_to_slp_adj (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
  17    implicit none
  18
  19    real,              intent(in)  :: TERR, PSFC
  20    real,              intent(in)  :: T(kms:kme), Q(kms:kme), P(kms:kme)
  21    real,              intent(out) :: T9(kms:kme), Q9(kms:kme), P9(kms:kme)
  22    real,              intent(in)  :: SLP9
  23    real,              intent(out) :: PSFC9
  24
  25    ! real :: SLP
  26
  27    integer              :: K, KLO, KHI
  28    real                 :: PL, T0, TS, XTERM,TLO, THI, TL
  29    real                 :: PL9,T09,TS9,TLO9,THI9,TL9,COEF1,COEF2
  30    real                 :: T08,TS8
  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_adj")
  38
  39    ! sea level pressure
  40
  41    xterm=gamma* gas_constant / gravity
  42
  43    ! compute pressure at pconst mb above surface (pl)
  44
  45    T9(:) = 0.0
  46    Q9(:) = 0.0
  47    P9(:) = 0.0
  48
  49    if (terr <= 0.0) then
  50       psfc9=slp9
  51       if (trace_use) call da_trace_exit("da_tpq_to_slp_adj")
  52       return
  53    end if
  54
  55    PL  = psfc - PCONST
  56    klo = 0
  57
  58    ! find 2 levels on sigma surfaces surrounding pl at each i,j
  59
  60    do k=kts, kte-1
  61       if ((p(k) >= pl) .and. (p(k+1) < pl)) then
  62          khi = k+1
  63          klo = k
  64          exit
  65       end if
  66    end do
  67
  68    if (klo < 1) then
  69       write(unit=message(1),fmt='(A,F11.3,A)') &
  70          'Error finding pressure level ',PCONST,' Mb above the surface'
  71       write(unit=message(2),fmt='(A,F11.3,2X,A,F11.3)') &
  72          'PL=',PL,'  PSFC=',psfc
  73       call da_error(__FILE__,__LINE__,message(1:2))
  74    end if
  75
  76    ! get temperature at pl (tl), extrapolate t at surface (ts)
  77    ! and t at sea level (t0) with 6.5 k/km lapse rate
  78
  79    TLO   = t(KLO) * (EPS+q(KLO))/(EPS*(1.0+q(KLO)))
  80    THI   = t(KHI) * (EPS+q(KHI))/(EPS*(1.0+q(KHI)))
  81    COEF1 = ALOG(PL/p(KHI))
  82    COEF2 = ALOG(p(KLO)/p(KHI))
  83    TL    = THI-(THI-TLO)*COEF1/COEF2
  84    TS    = TL*(psfc/PL)**XTERM
  85    TS8   = TS
  86    T0    = TS +GAMMA*terr
  87    T08   = T0
  88
  89    ! Correct sea level temperature if too hot
  90
  91    if (t0 >= tc) then
  92       if (ts <= tc) then
  93          t0 = tc
  94       else
  95          t0 = tc-0.005*(ts-tc)**2
  96       end if
  97    end if
  98
  99    psfc9 = slp9*EXP(2.0*gravity*TERR/(gas_constant*(TS+T0)))
 100    TS9 = -slp9*psfc*EXP(2.0*gravity*TERR/(gas_constant*(TS+T0)))*  &
 101        2.0*gravity*TERR/(gas_constant*(TS+T0)**2)
 102    T09 = -slp9*psfc*EXP(2.0*gravity*TERR/(gas_constant*(TS+T0)))*  &
 103        2.0*gravity*TERR/(gas_constant*(TS+T0)**2)
 104
 105    if ( t08 >= tc ) then
 106       if ( ts8 > tc ) then
 107          ts9=ts9-0.01*(ts-tc)*t09
 108       end if
 109       t09=0.0
 110    end if
 111
 112    TS9     = TS9+T09
 113    TL9     = TS9*(psfc/PL)**XTERM
 114    psfc9   = psfc9+TS9*XTERM*(TL/PL)*(psfc/PL)**(XTERM-1)
 115    PL9     = -TS9*XTERM*(TL*psfc/(PL*PL))*(psfc/PL)**(XTERM-1)
 116    THI9    = (1.0-COEF1/COEF2)*TL9
 117    TLO9    = COEF1/COEF2*TL9
 118    PL9     = PL9-(THI-TLO)/(COEF2*PL)*TL9
 119    p9(KHI) = p9(KHI)+((THI-TLO)/(COEF2*p(KHI))*(1.0-COEF1/COEF2))*TL9
 120    p9(KLO) = p9(KLO)+(THI-TLO)*COEF1/(COEF2*COEF2*p(KLO))*TL9
 121    t9(KHI) = t9(KHI)+THI9* (EPS+q(KHI))/(EPS*(1.0+q(KHI)))
 122    q9(KHI) = q9(KHI)+THI9*t(KHI)*(1.0-EPS)/(EPS*(1.0+q(KHI))**2)
 123    t9(KLO) = t9(KLO)+TLO9* (EPS+q(KLO))/(EPS*(1.0+q(KLO)))
 124    q9(KLO) = q9(KLO)+TLO9*t(KLO)*(1.0-EPS)/(EPS*(1.0+q(KLO))**2)
 125
 126    psfc9=psfc9+PL9
 127
 128    if (trace_use) call da_trace_exit("da_tpq_to_slp_adj")
 129
 130 end subroutine da_tpq_to_slp_adj
 131
 132