Line data    Source code

       1             : module std_atm_profile
       2             : 
       3             : !-------------------------------------------------------------------------------
       4             : !
       5             : ! The barometric formula for U.S. Standard Atmosphere is valid up to 86 km.
       6             : ! see https://en.wikipedia.org/wiki/Barometric_formula.
       7             : !
       8             : ! N.B.  The extension above 86 km is using data from Hanli.  It is not complete
       9             : !       since the hardcoded parameter (c1) needs adjustment above 86 km.
      10             : !
      11             : !-------------------------------------------------------------------------------
      12             : 
      13             : use shr_kind_mod,        only: r8 => shr_kind_r8
      14             : use cam_logfile,         only: iulog
      15             : use cam_abortutils,      only: endrun
      16             : 
      17             : implicit none
      18             : private
      19             : save
      20             : 
      21             : public :: &
      22             :    std_atm_pres,   & ! compute pressure given height
      23             :    std_atm_height, & ! compute height given pressure
      24             :    std_atm_temp      ! compute temperature given height
      25             : 
      26             : ! Parameters for barometric formula for U.S. Standard Atmosphere.
      27             : 
      28             : integer, parameter  :: nreg = 15  ! number of regions
      29             : 
      30             : real(r8), parameter :: hb(nreg) = & ! height at bottom of layer (m)
      31             :      (/0.0_r8, 1.1e4_r8, 2.0e4_r8, 3.2e4_r8, 4.7e4_r8, 5.1e4_r8, 7.1e4_r8, 8.6e4_r8, &
      32             :      9.1e4_r8, 1.1e5_r8, 1.2e5_r8, 1.5e5_r8, 2.0e5_r8, 3.0e5_r8, 7.e5_r8/)
      33             : 
      34             : real(r8), parameter :: pb(nreg) = & ! standard pressure (Pa)
      35             :      (/101325._r8, 22632.1_r8, 5474.89_r8, 868.02_r8, 110.91_r8, 66.94_r8, 3.96_r8, 3.7e-1_r8,  &
      36             :      1.5e-1_r8, 7.1e-3_r8, 2.5e-3_r8, 4.5e-4_r8, 8.47e-5_r8, 8.77e-6_r8, 3.19e-8_r8/)
      37             : 
      38             : real(r8), parameter :: tb(nreg) = & ! standard temperature (K)
      39             :      (/288.15_r8, 216.65_r8, 216.65_r8, 228.65_r8, 270.65_r8, 270.65_r8, 214.65_r8, 186.87_r8,  &
      40             :      186.87_r8, 240._r8, 360._r8, 634.39_r8, 854.56_r8, 976.01_r8, 1.e3_r8/)
      41             : 
      42             : real(r8), parameter :: lb(nreg) = & ! temperature lapse rate (K/m)
      43             :      (/-0.0065_r8, 0.0_r8, 0.001_r8, 0.0028_r8, 0.0_r8, -0.0028_r8, -0.001852_r8, 0.0_r8,       &
      44             :      2.796e-3_r8, 0.012_r8, 9.15e-3_r8, 4.4e-3_r8, 1.21e-3_r8, 6.e-5_r8, 0.0_r8/)
      45             : 
      46             : real(r8), parameter :: rg = 8.3144598_r8 ! universal gas constant (J/mol/K)
      47             : real(r8), parameter :: g0 = 9.80665_r8   ! gravitational acceleration (m/s^2)
      48             : real(r8), parameter :: mw = 0.0289644_r8 ! molar mass of dry air (kg/mol)
      49             : real(r8), parameter :: c1 = g0*mw/rg
      50             : 
      51             : !=========================================================================================
      52             : CONTAINS
      53             : !=========================================================================================
      54             : 
      55           0 : subroutine std_atm_pres(height, pstd, user_specified_ps)
      56             : 
      57             :    ! arguments
      58             :    real(r8),           intent(in)  :: height(:) ! height above sea level in meters
      59             :    real(r8),           intent(out) :: pstd(:)   ! std pressure in Pa
      60             :    real(r8), optional, intent(in)  :: user_specified_ps
      61             : 
      62             :    integer  :: i, ii, k, nlev
      63             :    integer  :: ierr
      64             :    real(r8) :: pb_local(nreg)
      65             : 
      66             :    character(len=*), parameter :: routine = 'std_atm_pres'
      67             :    !----------------------------------------------------------------------------
      68             : 
      69             :    ! Initialize local standard pressure values array
      70           0 :    pb_local = pb
      71             : 
      72             :    ! Set new surface pressure value if provided by the caller
      73           0 :    if (present(user_specified_ps)) then
      74           0 :       pb_local(1) = user_specified_ps
      75             :    end if
      76             : 
      77           0 :    nlev = size(height)
      78           0 :    do k = 1, nlev
      79           0 :       if (height(k) < 0.0_r8) then
      80             :          ! Extrapolate below mean sea level using troposphere lapse rate.
      81             :          ii = 1
      82             :       else
      83             :          ! find region containing height
      84           0 :          find_region: do i = nreg, 1, -1
      85           0 :             if (height(k) >= hb(i)) then
      86             :                ii = i
      87             :                exit find_region
      88             :             end if
      89             :          end do find_region
      90             :       end if
      91             : 
      92           0 :       if (lb(ii) /= 0._r8) then
      93           0 :          pstd(k) = pb_local(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
      94             :       else
      95           0 :          pstd(k) = pb_local(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
      96             :       end if
      97             : 
      98             :    end do
      99           0 : end subroutine std_atm_pres
     100             : 
     101             : !=========================================================================================
     102             : 
     103          12 : subroutine std_atm_height(pstd, height)
     104             : 
     105             :    ! arguments
     106             :    real(r8), intent(in)   :: pstd(:)   ! std pressure in Pa
     107             :    real(r8), intent(out)  :: height(:) ! height above sea level in meters
     108             : 
     109             :    integer :: i, ii, k, nlev
     110             :    logical :: found_region
     111             :    character(len=*), parameter :: routine = 'std_atm_height'
     112             :    !----------------------------------------------------------------------------
     113             : 
     114          12 :    nlev = size(height)
     115         208 :    do k = 1, nlev
     116             : 
     117         196 :       if (pstd(k) <= pb(nreg)) then
     118             :          ii = nreg
     119         196 :       else if (pstd(k) > pb(1)) then
     120             :          ii = 1
     121             :       else
     122             :          ! find region containing pressure
     123         514 :          find_region: do i = 2, nreg
     124         514 :             if (pstd(k) > pb(i)) then
     125         196 :                ii = i - 1
     126         196 :                exit find_region
     127             :             end if
     128             :          end do find_region
     129             :       end if
     130             : 
     131         208 :       if (lb(ii) /= 0._r8) then
     132         156 :          height(k) = hb(ii) + (tb(ii)/lb(ii)) * ( (pb(ii)/pstd(k))**(lb(ii)/c1) - 1._r8 )
     133             :       else
     134          40 :          height(k) = hb(ii) + (tb(ii)/c1)*log(pb(ii)/pstd(k))
     135             :       end if
     136             :    end do
     137             : 
     138          12 : end subroutine std_atm_height
     139             : 
     140             : !=========================================================================================
     141             : 
     142           0 : subroutine std_atm_temp(height, temp)
     143             : 
     144             :    ! arguments
     145             :    real(r8), intent(in)   :: height(:) ! std pressure in Pa
     146             :    real(r8), intent(out)  :: temp(:)   ! temperature
     147             : 
     148             :    ! local vars
     149             :    integer :: i, ii, k, nlev
     150             :    character(len=*), parameter :: routine = 'std_atm_temp'
     151             :    !----------------------------------------------------------------------------
     152             : 
     153           0 :    nlev = size(height)
     154           0 :    do k = 1, nlev
     155           0 :       if (height(k) < 0.0_r8) then
     156             :          ii = 1
     157             :       else
     158             :          ! find region containing height
     159           0 :          find_region: do i = nreg, 1, -1
     160           0 :             if (height(k) >= hb(i)) then
     161             :                ii = i
     162             :                exit find_region
     163             :             end if
     164             :          end do find_region
     165             :       end if
     166             : 
     167           0 :       if (lb(ii) /= 0._r8) then
     168           0 :          temp(k) = tb(ii) + lb(ii)*(height(k) - hb(ii))
     169             :       else
     170           0 :          temp(k) = tb(ii)
     171             :       end if
     172             : 
     173             :    end do
     174             : 
     175           0 : end subroutine std_atm_temp
     176             : 
     177             : end module std_atm_profile