Line data    Source code

       1             : module cam_diagnostic_utils
       2             : 
       3             : ! Collection of routines used for diagnostic calculations.
       4             : 
       5             : use shr_kind_mod, only: r8 => shr_kind_r8
       6             : use ppgrid,       only: pcols, pver
       7             : 
       8             : 
       9             : implicit none
      10             : private
      11             : save
      12             : 
      13             : public :: &
      14             :    cpslec      ! compute sea level pressure
      15             : 
      16             : !===============================================================================
      17             : contains
      18             : !===============================================================================
      19             : 
      20     1495368 : subroutine cpslec(ncol, pmid, phis, ps, t, psl, gravit, rair)
      21             : 
      22             : !----------------------------------------------------------------------- 
      23             : ! 
      24             : ! Compute sea level pressure.
      25             : ! 
      26             : ! Uses ECMWF formulation Algorithm: See section 3.1.b in NCAR NT-396 "Vertical 
      27             : ! Interpolation and Truncation of Model-Coordinate Data
      28             : !
      29             : !-----------------------------------------------------------------------
      30             : 
      31             :   !-----------------------------Arguments---------------------------------
      32             :   integer , intent(in) :: ncol             ! longitude dimension
      33             : 
      34             :   real(r8), intent(in) :: pmid(pcols,pver) ! Atmospheric pressure (pascals)
      35             :   real(r8), intent(in) :: phis(pcols)      ! Surface geopotential (m**2/sec**2)
      36             :   real(r8), intent(in) :: ps(pcols)        ! Surface pressure (pascals)
      37             :   real(r8), intent(in) :: T(pcols,pver)    ! Vertical slice of temperature (top to bot)
      38             :   real(r8), intent(in) :: gravit           ! Gravitational acceleration
      39             :   real(r8), intent(in) :: rair             ! gas constant for dry air
      40             : 
      41             :   real(r8), intent(out):: psl(pcols)       ! Sea level pressures (pascals)
      42             : 
      43             :   !-----------------------------Parameters--------------------------------
      44             :   real(r8), parameter :: xlapse = 6.5e-3_r8   ! Temperature lapse rate (K/m)
      45             : 
      46             :   !-----------------------------Local Variables---------------------------
      47             :   integer  :: i             ! Loop index
      48             :   real(r8) :: alpha         ! Temperature lapse rate in terms of pressure ratio (unitless)
      49             :   real(r8) :: Tstar         ! Computed surface temperature
      50             :   real(r8) :: TT0           ! Computed temperature at sea-level
      51             :   real(r8) :: alph          ! Power to raise P/Ps to get rate of increase of T with pressure
      52             :   real(r8) :: beta          ! alpha*phis/(R*T) term used in approximation of PSL
      53             :   !-----------------------------------------------------------------------
      54             : 
      55     1495368 :   alpha = rair*xlapse/gravit
      56    24969168 :   do i=1,ncol
      57    24969168 :      if ( abs(phis(i)/gravit) < 1.e-4_r8 )then
      58    23473800 :         psl(i)=ps(i)
      59             :      else
      60           0 :         Tstar=T(i,pver)*(1._r8+alpha*(ps(i)/pmid(i,pver)-1._r8)) ! pg 7 eq 5
      61             : 
      62           0 :         TT0=Tstar + xlapse*phis(i)/gravit                  ! pg 8 eq 13
      63             : 
      64           0 :         if ( Tstar<=290.5_r8 .and. TT0>290.5_r8 ) then           ! pg 8 eq 14.1
      65           0 :            alph=rair/phis(i)*(290.5_r8-Tstar)  
      66           0 :         else if (Tstar>290.5_r8  .and. TT0>290.5_r8) then        ! pg 8 eq 14.2
      67           0 :            alph=0._r8
      68           0 :            Tstar= 0.5_r8 * (290.5_r8 + Tstar)  
      69             :         else  
      70           0 :            alph=alpha  
      71           0 :            if (Tstar<255._r8) then  
      72           0 :               Tstar= 0.5_r8 * (255._r8 + Tstar)                  ! pg 8 eq 14.3
      73             :            endif
      74             :         endif
      75             : 
      76           0 :         beta = phis(i)/(rair*Tstar)
      77           0 :         psl(i)=ps(i)*exp( beta*(1._r8-alph*beta/2._r8+((alph*beta)**2)/3._r8))
      78             :      end if
      79             :   enddo
      80             : 
      81     1495368 : end subroutine cpslec
      82             : 
      83             : !===============================================================================
      84             : 
      85             : end module cam_diagnostic_utils