LCOV - code coverage report
Current view: | top level - utils - std_atm_profile.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 13 | 36 | 36.1 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 3 | 33.3 %

              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

Generated by: LCOV version 1.14