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)
      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             :     
      61             :    integer :: i, ii, k, nlev
      62             :    character(len=*), parameter :: routine = 'std_atm_pres'
      63             :    !----------------------------------------------------------------------------
      64             :     
      65           0 :    nlev = size(height)
      66           0 :    do k = 1, nlev
      67           0 :       if (height(k) < 0.0_r8) then
      68             :          ! Extrapolate below mean sea level using troposphere lapse rate.
      69             :          ii = 1
      70             :       else
      71             :          ! find region containing height
      72           0 :          find_region: do i = nreg, 1, -1
      73           0 :             if (height(k) >= hb(i)) then
      74             :                ii = i
      75             :                exit find_region
      76             :             end if
      77             :          end do find_region
      78             :       end if
      79             :       
      80           0 :       if (lb(ii) /= 0._r8) then
      81           0 :          pstd(k) = pb(ii) * ( tb(ii) / (tb(ii) + lb(ii)*(height(k) - hb(ii)) ) )**(c1/lb(ii))
      82             :       else
      83           0 :          pstd(k) = pb(ii) * exp( -c1*(height(k) - hb(ii))/tb(ii) )
      84             :       end if
      85             :       
      86             :    end do
      87             : 
      88           0 : end subroutine std_atm_pres
      89             : 
      90             : !=========================================================================================
      91             : 
      92          12 : subroutine std_atm_height(pstd, height)
      93             :     
      94             :    ! arguments
      95             :    real(r8), intent(in)   :: pstd(:)   ! std pressure in Pa
      96             :    real(r8), intent(out)  :: height(:) ! height above sea level in meters
      97             :     
      98             :    integer :: i, ii, k, nlev
      99             :    logical :: found_region
     100             :    character(len=*), parameter :: routine = 'std_atm_height'
     101             :    !----------------------------------------------------------------------------
     102             :     
     103          12 :    nlev = size(height)
     104         208 :    do k = 1, nlev
     105             :       
     106         196 :       if (pstd(k) <= pb(nreg)) then
     107             :          ii = nreg
     108         196 :       else if (pstd(k) > pb(1)) then
     109             :          ii = 1
     110             :       else
     111             :          ! find region containing pressure
     112         514 :          find_region: do i = 2, nreg
     113         514 :             if (pstd(k) > pb(i)) then
     114         196 :                ii = i - 1
     115         196 :                exit find_region
     116             :             end if
     117             :          end do find_region
     118             :       end if
     119             : 
     120         208 :       if (lb(ii) /= 0._r8) then
     121         156 :          height(k) = hb(ii) + (tb(ii)/lb(ii)) * ( (pb(ii)/pstd(k))**(lb(ii)/c1) - 1._r8 )
     122             :       else
     123          40 :          height(k) = hb(ii) + (tb(ii)/c1)*log(pb(ii)/pstd(k))
     124             :       end if
     125             :    end do
     126             : 
     127          12 : end subroutine std_atm_height
     128             : 
     129             : !=========================================================================================
     130             : 
     131           0 : subroutine std_atm_temp(height, temp)
     132             :     
     133             :    ! arguments
     134             :    real(r8), intent(in)   :: height(:) ! std pressure in Pa
     135             :    real(r8), intent(out)  :: temp(:)   ! temperature
     136             :     
     137             :    ! local vars
     138             :    integer :: i, ii, k, nlev
     139             :    character(len=*), parameter :: routine = 'std_atm_temp'
     140             :    !----------------------------------------------------------------------------
     141             :     
     142           0 :    nlev = size(height)
     143           0 :    do k = 1, nlev
     144           0 :       if (height(k) < 0.0_r8) then
     145             :          ii = 1
     146             :       else
     147             :          ! find region containing height
     148           0 :          find_region: do i = nreg, 1, -1
     149           0 :             if (height(k) >= hb(i)) then
     150             :                ii = i
     151             :                exit find_region
     152             :             end if
     153             :          end do find_region
     154             :       end if
     155             : 
     156           0 :       if (lb(ii) /= 0._r8) then
     157           0 :          temp(k) = tb(ii) + lb(ii)*(height(k) - hb(ii))
     158             :       else
     159           0 :          temp(k) = tb(ii)
     160             :       end if
     161             :       
     162             :    end do
     163             : 
     164           0 : end subroutine std_atm_temp
     165             : 
     166             : end module std_atm_profile