reorganize temperature based cross section

src/cross_sections/CMakeLists.txt

@@ -31,4 +31,6 @@ target_sources(tuvx_object
     rayliegh.F90
 )
 
+add_subdirectory(util)
+
 ################################################################################

src/cross_sections/util/CMakeLists.txt

@@ -0,0 +1,10 @@
+################################################################################
+# utilities for cross section parameterizations
+
+target_sources(tuvx_object
+  PRIVATE
+    temperature_parameterization.F90
+    temperature_range.F90  
+)
+
+################################################################################

src/cross_sections/util/temperature_range.F90

@@ -0,0 +1,159 @@
+! Copyright (C) 2020-4 National Center for Atmospheric Research
+! SPDX-License-Identifier: Apache-2.0
+
+module tuvx_temperature_range
+! Defines a temperature range for use in temperature-based cross
+! section parameterizations
+
+  ! Including musica_config at the module level to avoid an ICE
+  ! with Intel 2022.1 compiler
+  use musica_config,                   only : config_t
+  use musica_constants,                only : dk => musica_dk
+
+  implicit none
+
+  private
+  public :: temperature_range_t
+
+
+  !> Range for temperature-based calculations
+  type :: temperature_range_t
+    !> Minimum temperature [K] for inclusion in range
+    real(kind=dk) :: min_temperature_ = 0.0_dk
+    !> Maximum temperature [K] for include in range
+    real(kind=dk) :: max_temperature_ = huge(1.0_dk)
+    !> Indicates whether to use a fixed temperature for the
+    !! parameterization calculation. If FALSE, the actual
+    !! temperature is used.
+    logical :: is_fixed_ = .false.
+    !> Fixed temperature [K] to use in paramterization calculation
+    !!
+    !! Is only used if is_fixed == TRUE
+    real(kind=dk) :: fixed_temperature_ = 0.0_dk
+  contains
+    !> Returns the number of bytes required to pack the range onto a
+    !! character buffer
+    procedure :: pack_size => temperature_range_pack_size
+    !> Packs the range onto a character buffer
+    procedure :: mpi_pack => temperature_range_mpi_pack
+    !> Unpacks a range from a character buffer
+    procedure :: mpi_unpack => temperature_range_mpi_unpack
+  end type temperature_range_t
+
+  !> Constructor for temperature_range_t
+  interface temperature_range_t
+    module procedure :: temperature_range_constructor
+  end interface temperature_range_t
+
+contains
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  function temperature_range_constructor( config ) result( this )
+    ! Constructs temperature range objects
+
+    use musica_assert,                 only : assert_msg
+    use musica_config,                 only : config_t
+    use musica_string,                 only : string_t
+
+    type(temperature_range_t)               :: this
+    type(config_t),           intent(inout) :: config
+
+    character(len=*), parameter :: my_name = "temperature range constructor"
+    type(string_t) :: required_keys(0), optional_keys(3)
+    logical :: found
+
+    optional_keys(1) = "minimum"
+    optional_keys(2) = "maximum"
+    optional_keys(3) = "fixed value"
+    call assert_msg( 355912601,                                               &
+                     config%validate( required_keys, optional_keys ),         &
+                     "Bad configuration for temperature range" )
+
+    call config%get( "minimum", this%min_temperature_, my_name,               &
+                     default = 0.0_dk )
+    call config%get( "maximum", this%max_temperature_, my_name,               &
+                     default = huge(1.0_dk) )
+    call config%get( "fixed value", this%fixed_temperature_, my_name,         &
+                     found = found )
+    this%is_fixed_ = found
+
+  end function temperature_range_constructor
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  integer function temperature_range_pack_size( this, comm )                  &
+      result( pack_size )
+    ! Returns the size of a character buffer required to pack the range
+
+    use musica_mpi,                    only : musica_mpi_pack_size
+
+    class(temperature_range_t), intent(in) :: this ! temperature range to be packed
+    integer,                    intent(in) :: comm ! MPI communicator
+
+#ifdef MUSICA_USE_MPI
+    pack_size = musica_mpi_pack_size( this%min_temperature_,   comm ) +       &
+                musica_mpi_pack_size( this%max_temperature_,   comm ) +       &
+                musica_mpi_pack_size( this%is_fixed_,          comm ) +       &
+                musica_mpi_pack_size( this%fixed_temperature_, comm )
+#else
+    pack_size = 0
+#endif
+
+  end function temperature_range_pack_size
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  subroutine temperature_range_mpi_pack( this, buffer, position, comm )
+    ! Packs the temperature range onto a character buffer
+
+    use musica_assert,                 only : assert
+    use musica_mpi,                    only : musica_mpi_pack
+
+    class(temperature_range_t), intent(in)    :: this      ! temperature range to be packed
+    character,                  intent(inout) :: buffer(:) ! memory buffer
+    integer,                    intent(inout) :: position  ! current buffer position
+    integer,                    intent(in)    :: comm      ! MPI communicator
+
+#ifdef MUSICA_USE_MPI
+    integer :: prev_pos
+
+    prev_pos = position
+    call musica_mpi_pack( buffer, position, this%min_temperature_,   comm )
+    call musica_mpi_pack( buffer, position, this%max_temperature_,   comm )
+    call musica_mpi_pack( buffer, position, this%is_fixed_,          comm )
+    call musica_mpi_pack( buffer, position, this%fixed_temperature_, comm )
+    call assert( 409699380, position - prev_pos <= this%pack_size( comm ) )
+#endif
+
+  end subroutine temperature_range_mpi_pack
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  subroutine temperature_range_mpi_unpack( this, buffer, position, comm )
+    ! Unpacks a temperature range from a character buffer
+
+    use musica_assert,                 only : assert
+    use musica_mpi,                    only : musica_mpi_unpack
+
+    class(temperature_range_t), intent(out)   :: this      ! temperature range to be unpacked
+    character,                  intent(inout) :: buffer(:) ! memory buffer
+    integer,                    intent(inout) :: position  ! current buffer position
+    integer,                    intent(in)    :: comm      ! MPI communicator
+
+#ifdef MUSICA_USE_MPI
+    integer :: prev_pos
+
+    prev_pos = position
+    call musica_mpi_unpack( buffer, position, this%min_temperature_,   comm )
+    call musica_mpi_unpack( buffer, position, this%max_temperature_,   comm )
+    call musica_mpi_unpack( buffer, position, this%is_fixed_,          comm )
+    call musica_mpi_unpack( buffer, position, this%fixed_temperature_, comm )
+    call assert( 164457757, position - prev_pos <= this%pack_size( comm ) )
+#endif
+
+  end subroutine temperature_range_mpi_unpack
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+end module tuvx_temperature_range

test/unit/tuv_doug/JCALC/CMakeLists.txt

@@ -4,6 +4,7 @@
 target_sources(tuv_doug
   PRIVATE
     XSQY_BRO.f
+    XSQY_BRONO2.f
     XSQY_CF2CL2.f
     XSQY_CFC113.f
     XSQY_CFC114.f
@@ -21,6 +22,7 @@ target_sources(tuv_doug
     XSQY_HCFC141b.f
     XSQY_HCFC142b.f
     XSQY_HNO3.f
+    XSQY_HO2NO2.f
     XSQY_N2O5.f
 )
 

test/unit/tuv_doug/JCALC/XSQY_BRONO2.f

@@ -0,0 +1,176 @@
+      subroutine XSQY_BRONO2(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn)
+!-----------------------------------------------------------------------------!
+!   purpose:                                                                  !
+!   provide product (cross section) x (quantum yield) for brono2 photolysis:  !
+!         BrONO2 + hv -> products                                             !
+!                                                                             !
+!   cross section: jpl 06 recommendation                                      !
+!   quantum yield: jpl 06 recommendation                                      !
+!-----------------------------------------------------------------------------!
+!   parameters:                                                               !
+!   nw     - integer, number of specified intervals + 1 in working        (i) !
+!            wavelength grid                                                  !
+!   wl     - real, vector of lower limits of wavelength intervals in      (i) !
+!            working wavelength grid                                          !
+!   wc     - real, vector of center points of wavelength intervals in     (i) !
+!            working wavelength grid                                          !
+!   nz     - integer, number of altitude levels in working altitude grid  (i) !
+!   tlev   - real, temperature (k) at each specified altitude level       (i) !
+!   airlev - real, air density (molec/cc) at each altitude level          (i) !
+!   j      - integer, counter for number of weighting functions defined  (io) !
+!   sq     - real, cross section x quantum yield (cm^2) for each          (o) !
+!            photolysis reaction defined, at each defined wavelength and      !
+!            at each defined altitude level                                   !
+!   jlabel - character*60, string identifier for each photolysis reaction (o) !
+!            defined                                                          !
+!-----------------------------------------------------------------------------!
+!   edit history:                                                             !
+!   07/27/07: Doug Kinnison                                                   !
+!-----------------------------------------------------------------------------!
+      implicit none
+      include 'params'
+
+!-----------------------------------------------------------------------------!
+!     ... input                                                               !
+!-----------------------------------------------------------------------------!
+      real, intent(in) :: wl(kw)
+      real, intent(in) :: wc(kw)
+      real, intent(in) :: tlev(kz)
+      real, intent(in) :: airlev(kz)
+
+      integer, intent(in) :: nz
+      integer, intent(in) :: nw
+
+      character*80, intent(in)  ::  pn
+      character*60, intent(out) :: jlabel(kj)
+
+      real, intent(out) :: sq(kj,kz,kw)
+
+!-----------------------------------------------------------------------------!
+!     ... input/output                                                        !
+!-----------------------------------------------------------------------------!
+      integer, intent(inout) :: j
+
+!-----------------------------------------------------------------------------!
+!     ... local                                                               !
+!-----------------------------------------------------------------------------!
+      integer kdata
+      parameter(kdata=200)
+      integer i, iw, n, n1, idum, ierr, iz
+      real x1  (kdata), y1(kdata)
+      real xin (kdata)
+      real a1  (kdata), a2(kdata)
+      real ytmp(nz,kdata)
+      real ytd (nz,kw)
+      real yg1 (kw)
+      real tin (nz)
+      real qy1
+      real qy2
+!-----------------------------------------------
+!     ... tin set to tlev 
+!-----------------------------------------------
+      tin(:) = tlev(:)
+
+!-----------------------------------------------
+!     ... jlabel(j) = 'BrONO2 + hv -> Br + NO3'
+!     ... jlabel(j) = 'BrONO2 + hv -> BrO + NO2'
+!-----------------------------------------------
+      j = j+1
+      jlabel(j) = 'BrONO2 + hv -> Br + NO3'
+
+!-----------------------------------------------
+!     ... cross sections from JPL06 
+!-----------------------------------------------
+      open(kin,
+     &     file=TRIM(pn)//'XS_BRONO2_JPL06.txt',status='old')
+
+      read(kin,*) idum, n
+      do i = 1, idum-2
+        read(kin,*)
+      enddo
+
+      do iw = 1, n
+        read(kin,*) xin(iw), y1(iw)
+      enddo
+      close(kin)
+
+!-----------------------------------------------
+!     ... Read in temperature dep coeff
+!-----------------------------------------------
+      open(kin,
+     &     file=TRIM(pn)//'XS_BRONO2_td_JPL06.txt',status='old')
+
+      read(kin,*) idum, n
+      do i = 1, idum-2
+        read(kin,*)
+      enddo
+
+      do iw = 1, n
+        read(kin,*) xin(iw), a1(iw), a2(iw)
+      enddo
+      close(kin)
+
+!-----------------------------------------------
+!     ... derive T-dep (200-296K)
+!-----------------------------------------------
+      do  iz = 1, nz
+        do iw = 1 , n
+          if ((tin(iz) .GE. 200.) .AND. (tin(iz) .LE. 296.)) Then
+            ytmp(iz,iw) = y1(iw) * 
+     &                   ( 1. + 
+     &                     A1(iw)* (tin(iz)-296.) +  
+     &                     A2(iw)*((tin(iz)-296.)**2) 
+     &                    )      
+          endif
+          if (tin(iz) .LT. 200.) then
+             ytmp(iz,iw) = y1(iw) * 
+     &                   ( 1. + 
+     &                     A1(iw)* (200.-296.) +  
+     &                     A2(iw)*((200.-296.)**2) 
+     &                    )      
+          endif
+          if (tin(iz) .GT. 296.) then
+             ytmp(iz,iw) = y1(iw) 
+          endif
+        enddo
+      enddo
+
+!-----------------------------------------------
+!     ... Interpolate
+!-----------------------------------------------
+      do iz = 1, nz
+        n1 = n
+        y1 = ytmp(iz,:)
+        x1 = xin
+
+        call addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
+        call addpnt(x1,y1,kdata,n1,          0.,0.)
+        call addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
+        call addpnt(x1,y1,kdata,n1,        1e38,0.)
+        call inter2(nw,wl,yg1,n1,x1,y1,ierr)
+        ytd(iz,:) = yg1(:)
+
+        if (ierr .ne. 0) then
+          write(*,*) ierr, jlabel(j)
+          stop
+        endif
+      enddo
+
+!----------------------------------------------
+!     ...Quantum yields JPL06
+!     ...This recommendation is only for >300nm
+!        However, it is used at all wavelengths
+!----------------------------------------------
+      qy1 = 0.85
+      qy2 = 0.15
+      do iw = 1, nw-1
+         do iz = 1, nz
+            sq(j,iz,iw)   = qy1 * ytd(iz,iw)
+            sq(j+1,iz,iw) = qy2 * ytd(iz,iw)
+         enddo
+      enddo
+
+      j = j+1
+      jlabel(j) = 'BrONO2 + hv -> BrO + NO2'
+
+      end subroutine XSQY_BRONO2