LCOV - code coverage report
Current view: | top level - physics/rrtmgp - radconstants.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 57 | 97 | 58.8 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 4 | 5 | 80.0 %

              Line data    Source code
           1             : module radconstants
           2             :
           3             : ! This module contains constants that are specific to the radiative transfer
           4             : ! code used in the RRTMGP model.
           5             :
           6             : use shr_kind_mod,         only: r8 => shr_kind_r8
           7             : use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp
           8             : use cam_abortutils,       only: endrun
           9             :
          10             : implicit none
          11             : private
          12             : save
          13             :
          14             : ! Number of bands in SW and LW.  These values must match data in the RRTMGP coefficients datasets.
          15             : ! But they are needed to allocate space in the physics buffer and need to be available before the
          16             : ! RRTMGP datasets are read.  So they are set as parameters here and checked in the
          17             : ! set_wavenumber_bands subroutine after the datasets are read.
          18             : integer, parameter, public :: nswbands = 14
          19             : integer, parameter, public :: nlwbands = 16
          20             :
          21             : ! Band limits (set from data in RRTMGP coefficient datasets)
          22             : real(r8), target :: wavenumber_low_shortwave(nswbands)
          23             : real(r8), target :: wavenumber_high_shortwave(nswbands)
          24             : real(r8), target :: wavenumber_low_longwave(nlwbands)
          25             : real(r8), target :: wavenumber_high_longwave(nlwbands)
          26             :
          27             : logical :: wavenumber_boundaries_set = .false.
          28             :
          29             : integer, public, protected :: nswgpts  ! number of SW g-points
          30             : integer, public, protected :: nlwgpts  ! number of LW g-points
          31             :
          32             : ! These are indices to specific bands for diagnostic output and COSP input.
          33             : integer, public, protected :: idx_sw_diag = -1     ! band contains 500-nm wave
          34             : integer, public, protected :: idx_nir_diag = -1    ! band contains 1000-nm wave
          35             : integer, public, protected :: idx_uv_diag = -1     ! band contains 400-nm wave
          36             : integer, public, protected :: idx_lw_diag = -1     ! band contains 1000 cm-1 wave (H20 window)
          37             : integer, public, protected :: idx_sw_cloudsim = -1 ! band contains 670-nm wave (for COSP)
          38             : integer, public, protected :: idx_lw_cloudsim = -1 ! band contains 10.5 micron wave (for COSP)
          39             :
          40             : ! GASES TREATED BY RADIATION (line spectra)
          41             : ! These names are recognized by RRTMGP.  They are in the coefficients files as
          42             : ! lower case strings.  These upper case names are used by CAM's namelist and
          43             : ! rad_constituents module.
          44             : integer, public, parameter :: gasnamelength = 5
          45             : integer, public, parameter :: nradgas = 8
          46             : character(len=gasnamelength), public, parameter :: gaslist(nradgas) &
          47             :    = (/'H2O  ','O3   ', 'O2   ', 'CO2  ', 'N2O  ', 'CH4  ', 'CFC11', 'CFC12'/)
          48             :
          49             : ! what is the minimum mass mixing ratio that can be supported by radiation implementation?
          50             : real(r8), public, parameter :: minmmr(nradgas) = epsilon(1._r8)
          51             :
          52             : public :: &
          53             :    set_wavenumber_bands,       &
          54             :    get_sw_spectral_boundaries, &
          55             :    get_lw_spectral_boundaries, &
          56             :    get_band_index_by_value,    &
          57             :    rad_gas_index
          58             :
          59             : !=========================================================================================
          60             : contains
          61             : !=========================================================================================
          62             :
          63        1536 : subroutine set_wavenumber_bands(kdist_sw, kdist_lw)
          64             :
          65             :    ! Set the low and high limits of the wavenumber grid for sw and lw.
          66             :    ! Values come from RRTMGP coefficients datasets, and are stored in the
          67             :    ! kdist objects.
          68             :    !
          69             :    ! Set band indices for bands containing specific wavelengths.
          70             :
          71             :    ! Arguments
          72             :    type(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw
          73             :    type(ty_gas_optics_rrtmgp), intent(in) :: kdist_lw
          74             :
          75             :    ! Local variables
          76             :    integer :: istat
          77        1536 :    real(r8), allocatable :: values(:,:)
          78             :
          79             :    character(len=128) :: errmsg
          80             :    character(len=*), parameter :: sub = 'set_wavenumber_bands'
          81             :    !----------------------------------------------------------------------------
          82             :
          83             :    ! Check that number of sw/lw bands in gas optics files matches the parameters.
          84        1536 :    if (kdist_sw%get_nband() /= nswbands) then
          85           0 :       write(errmsg,'(a,i4,a,i4)') 'number of sw bands in file, ', kdist_sw%get_nband(), &
          86           0 :          ", doesn't match parameter nswbands= ", nswbands
          87           0 :       call endrun(sub//': ERROR: '//trim(errmsg))
          88             :    end if
          89        1536 :    if (kdist_lw%get_nband() /= nlwbands) then
          90           0 :       write(errmsg,'(a,i4,a,i4)') 'number of lw bands in file, ', kdist_lw%get_nband(), &
          91           0 :          ", doesn't match parameter nlwbands= ", nlwbands
          92           0 :       call endrun(sub//': ERROR: '//trim(errmsg))
          93             :    end if
          94             :
          95        1536 :    nswgpts = kdist_sw%get_ngpt()
          96        1536 :    nlwgpts = kdist_lw%get_ngpt()
          97             :
          98             :    ! SW band bounds in cm^-1
          99        1536 :    allocate( values(2,nswbands), stat=istat )
         100        1536 :    if (istat/=0) then
         101           0 :       call endrun(sub//': ERROR allocating array: values(2,nswbands)')
         102             :    end if
         103        1536 :    values = kdist_sw%get_band_lims_wavenumber()
         104       23040 :    wavenumber_low_shortwave = values(1,:)
         105       23040 :    wavenumber_high_shortwave = values(2,:)
         106             :
         107             :    ! Indices into specific bands
         108        1536 :    idx_sw_diag     = get_band_index_by_value('sw', 500.0_r8, 'nm')
         109        1536 :    idx_nir_diag    = get_band_index_by_value('sw', 1000.0_r8, 'nm')
         110        1536 :    idx_uv_diag     = get_band_index_by_value('sw', 400._r8, 'nm')
         111        1536 :    idx_sw_cloudsim = get_band_index_by_value('sw', 0.67_r8, 'micron')
         112             :
         113        1536 :    deallocate(values)
         114             :
         115             :    ! LW band bounds in cm^-1
         116        1536 :    allocate( values(2,nlwbands), stat=istat )
         117        1536 :    if (istat/=0) then
         118           0 :       call endrun(sub//': ERROR allocating array: values(2,nlwbands)')
         119             :    end if
         120        1536 :    values = kdist_lw%get_band_lims_wavenumber()
         121       26112 :    wavenumber_low_longwave = values(1,:)
         122       26112 :    wavenumber_high_longwave = values(2,:)
         123             :
         124             :    ! Indices into specific bands
         125        1536 :    idx_lw_diag     = get_band_index_by_value('lw', 1000.0_r8, 'cm^-1')
         126        1536 :    idx_lw_cloudsim = get_band_index_by_value('lw', 10.5_r8, 'micron')
         127             :
         128        1536 :    wavenumber_boundaries_set = .true.
         129             :
         130        1536 : end subroutine set_wavenumber_bands
         131             :
         132             : !=========================================================================================
         133             :
         134        1536 : subroutine get_sw_spectral_boundaries(low_boundaries, high_boundaries, units)
         135             :
         136             :    ! provide spectral boundaries of each shortwave band
         137             :
         138             :    real(r8),    intent(out) :: low_boundaries(nswbands), high_boundaries(nswbands)
         139             :    character(*), intent(in) :: units ! requested units
         140             :
         141             :    character(len=*), parameter :: sub = 'get_sw_spectral_boundaries'
         142             :    !----------------------------------------------------------------------------
         143             :
         144        1536 :    if (.not. wavenumber_boundaries_set) then
         145           0 :       call endrun(sub//': ERROR, wavenumber boundaries not set. ')
         146             :    end if
         147             :
         148        1536 :    select case (units)
         149             :    case ('inv_cm','cm^-1','cm-1')
         150        1536 :       low_boundaries = wavenumber_low_shortwave
         151        1536 :       high_boundaries = wavenumber_high_shortwave
         152             :    case('m','meter','meters')
         153           0 :       low_boundaries = 1.e-2_r8/wavenumber_high_shortwave
         154           0 :       high_boundaries = 1.e-2_r8/wavenumber_low_shortwave
         155             :    case('nm','nanometer','nanometers')
         156           0 :       low_boundaries = 1.e7_r8/wavenumber_high_shortwave
         157           0 :       high_boundaries = 1.e7_r8/wavenumber_low_shortwave
         158             :    case('um','micrometer','micrometers','micron','microns')
         159           0 :       low_boundaries = 1.e4_r8/wavenumber_high_shortwave
         160           0 :       high_boundaries = 1.e4_r8/wavenumber_low_shortwave
         161             :    case('cm','centimeter','centimeters')
         162           0 :       low_boundaries  = 1._r8/wavenumber_high_shortwave
         163           0 :       high_boundaries = 1._r8/wavenumber_low_shortwave
         164             :    case default
         165        1536 :       call endrun(sub//': ERROR, requested spectral units not recognized: '//units)
         166             :    end select
         167             :
         168        1536 : end subroutine get_sw_spectral_boundaries
         169             :
         170             : !=========================================================================================
         171             :
         172           0 : subroutine get_lw_spectral_boundaries(low_boundaries, high_boundaries, units)
         173             :
         174             :    ! provide spectral boundaries of each longwave band
         175             :
         176             :    real(r8), intent(out) :: low_boundaries(nlwbands), high_boundaries(nlwbands)
         177             :    character(*), intent(in) :: units ! requested units
         178             :
         179             :    character(len=*), parameter :: sub = 'get_lw_spectral_boundaries'
         180             :    !----------------------------------------------------------------------------
         181             :
         182           0 :    if (.not. wavenumber_boundaries_set) then
         183           0 :       call endrun(sub//': ERROR, wavenumber boundaries not set. ')
         184             :    end if
         185             :
         186           0 :    select case (units)
         187             :    case ('inv_cm','cm^-1','cm-1')
         188           0 :       low_boundaries  = wavenumber_low_longwave
         189           0 :       high_boundaries = wavenumber_high_longwave
         190             :    case('m','meter','meters')
         191           0 :       low_boundaries  = 1.e-2_r8/wavenumber_high_longwave
         192           0 :       high_boundaries = 1.e-2_r8/wavenumber_low_longwave
         193             :    case('nm','nanometer','nanometers')
         194           0 :       low_boundaries  = 1.e7_r8/wavenumber_high_longwave
         195           0 :       high_boundaries = 1.e7_r8/wavenumber_low_longwave
         196             :    case('um','micrometer','micrometers','micron','microns')
         197           0 :       low_boundaries  = 1.e4_r8/wavenumber_high_longwave
         198           0 :       high_boundaries = 1.e4_r8/wavenumber_low_longwave
         199             :    case('cm','centimeter','centimeters')
         200           0 :       low_boundaries  = 1._r8/wavenumber_high_longwave
         201           0 :       high_boundaries = 1._r8/wavenumber_low_longwave
         202             :    case default
         203           0 :       call endrun(sub//': ERROR, requested spectral units not recognized: '//units)
         204             :    end select
         205             :
         206           0 : end subroutine get_lw_spectral_boundaries
         207             :
         208             : !=========================================================================================
         209             :
         210    10590848 : integer function rad_gas_index(gasname)
         211             :
         212             :    ! return the index in the gaslist array of the specified gasname
         213             :
         214             :    character(len=*),intent(in) :: gasname
         215             :    integer :: igas
         216             :
         217    10590848 :    rad_gas_index = -1
         218    43920396 :    do igas = 1, nradgas
         219    43920396 :       if (trim(gaslist(igas)).eq.trim(gasname)) then
         220    10590848 :          rad_gas_index = igas
         221    10590848 :          return
         222             :       endif
         223             :    enddo
         224           0 :    call endrun ("rad_gas_index: can not find gas with name "//gasname)
         225           0 : end function rad_gas_index
         226             :
         227             : !=========================================================================================
         228             :
         229        9216 : function get_band_index_by_value(swlw, targetvalue, units) result(ans)
         230             :
         231             :    ! Find band index for requested wavelength/wavenumber.
         232             :
         233             :    character(len=*), intent(in) :: swlw        ! sw or lw bands
         234             :    real(r8),         intent(in) :: targetvalue
         235             :    character(len=*), intent(in) :: units       ! units of targetvalue
         236             :    integer :: ans
         237             :
         238             :    ! local
         239        9216 :    real(r8), pointer, dimension(:) :: lowboundaries, highboundaries
         240             :    real(r8) :: tgt
         241             :    integer  :: nbnds, i
         242             :
         243             :    character(len=128) :: errmsg
         244             :    character(len=*), parameter :: sub = 'get_band_index_by_value'
         245             :    !----------------------------------------------------------------------------
         246             :
         247        6144 :    select case (swlw)
         248             :    case ('sw','SW','shortwave')
         249        6144 :       nbnds = nswbands
         250        6144 :       lowboundaries  => wavenumber_low_shortwave
         251        6144 :       highboundaries => wavenumber_high_shortwave
         252             :    case ('lw', 'LW', 'longwave')
         253        3072 :       nbnds = nlwbands
         254        3072 :       lowboundaries  => wavenumber_low_longwave
         255        3072 :       highboundaries => wavenumber_high_longwave
         256             :    case default
         257           0 :       call endrun('radconstants.F90: get_band_index_by_value: type of bands not recognized: '//swlw)
         258             :    end select
         259             :
         260             :    ! band info is in cm^-1 but target value may be other units,
         261             :    ! so convert targetvalue to cm^-1
         262        1536 :    select case (units)
         263             :    case ('inv_cm','cm^-1','cm-1')
         264        1536 :       tgt = targetvalue
         265             :    case('m','meter','meters')
         266           0 :       tgt = 1.0_r8 / (targetvalue * 1.e2_r8)
         267             :    case('nm','nanometer','nanometers')
         268        4608 :       tgt = 1.0_r8 / (targetvalue * 1.e-7_r8)
         269             :    case('um','micrometer','micrometers','micron','microns')
         270        3072 :       tgt = 1.0_r8 / (targetvalue * 1.e-4_r8)
         271             :    case('cm','centimeter','centimeters')
         272           0 :       tgt = 1._r8/targetvalue
         273             :    case default
         274        9216 :       call endrun('radconstants.F90: get_band_index_by_value: units not recognized: '//units)
         275             :    end select
         276             :
         277             :    ! now just loop through the array
         278        9216 :    ans = 0
         279       84480 :    do i = 1,nbnds
         280       84480 :       if ((tgt > lowboundaries(i)) .and. (tgt <= highboundaries(i))) then
         281             :          ans = i
         282             :          exit
         283             :       end if
         284             :    end do
         285             :
         286        9216 :    if (ans == 0) then
         287           0 :       write(errmsg,'(f10.3,a,a)') targetvalue, ' ', trim(units)
         288           0 :       call endrun(sub//': band not found containing wave: '//trim(errmsg))
         289             :    end if
         290             :
         291        9216 : end function get_band_index_by_value
         292             :
         293             : !=========================================================================================
         294             :
         295             : end module radconstants

Generated by: LCOV version 1.14