
      module mo_seto2

      use shr_kind_mod, only : r8 => shr_kind_r8
      use cam_abortutils,   only : endrun
      use spmd_utils,   only : masterproc
      use cam_logfile,  only : iulog

      implicit none

      private
      public :: o2_xsect_inti
      public :: set_o2_xsect

      save

      integer :: nsrc
      integer :: ngast
      integer :: nla
      integer :: nwint
      real(r8), allocatable    :: wlint(:)
      real(r8), allocatable    :: xso2int(:)
      real(r8), allocatable    :: wlla(:)
      real(r8), allocatable    :: wlgast(:)

      contains

      subroutine o2_xsect_inti( o2_xsect_file )
!-----------------------------------------------------------------------------
!   purpose:
!   compute equivalent optical depths for o2 absorption, parameterized in
!   the sr bands and the lyman-alpha line.
!-----------------------------------------------------------------------------
!   parameters:
!   nz      - integer, number of specified altitude levels in the working (i)
!             grid
!   z       - real(r8), specified altitude working grid (km)                  (i)
!   nw      - integer, number of specified intervals + 1 in working       (i)
!             wavelength grid
!   wl      - real(r8), vector of lower limits of wavelength intervals in     (i)
!             working wavelength grid
!   cz      - real(r8), number of air molecules per cm^2 at each specified    (i)
!             altitude layer
!   zen     - real(r8), solar zenith angle                                    (i)
!   dto2    - real(r8), optical depth due to o2 absorption at each specified  (o)
!             vertical layer at each specified wavelength
!   xso2    - real(r8), molecular absorption cross section in sr bands at     (o)
!             each specified altitude and wavelength.  includes herzberg
!             continuum.
!-----------------------------------------------------------------------------

      use mo_params,     only : deltax
      use mo_inter,      only : inter2
      use mo_inter,      only : inter_inti
      use mo_wavelen,    only : nw, wl
      use ioFileMod,     only : getfil
      use cam_pio_utils, only : cam_pio_openfile
      use pio, only : file_desc_t, pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, &
           pio_get_var, pio_closefile, pio_nowrite

      implicit none

!-----------------------------------------------------------------------------
!	... dummy arguments
!-----------------------------------------------------------------------------
      character(len=*), intent(in) :: o2_xsect_file

!-----------------------------------------------------------------------------
!	... local variables
!-----------------------------------------------------------------------------
      type(file_desc_t) :: ncid
      integer :: dimid
      integer :: vid
      integer :: astat
      integer :: iret
      integer :: i, wn, n
      integer :: wrk_ind(4)
      real(r8), allocatable :: x1(:)
      real(r8), allocatable :: y1(:)
      character(len=256) :: filespec
      character(len=256) :: locfn
      integer :: ierr
!-----------------------------------------------------------------------------
! 	... cross section data for use outside the sr-bands (combined from
!           brasseur and solomon and the jpl 1994 recommendation)
!-----------------------------------------------------------------------------

!-----------------------------------------------------------------------------
! 	... read o2 cross section data outside sr-bands
!-----------------------------------------------------------------------------
!	... o2 absorption cross sections:
!           from 116 nm to 245 nm, including schumann-runge continumm
!           from brasseur and solomon 1986.
!-----------------------------------------------------------------------------
      filespec = trim( o2_xsect_file )
      call getfil( filespec, locfn, 0 )
      call cam_pio_openfile( ncid, trim( locfn ), PIO_NOWRITE )
!---------------------------------------------------------------------------
! 	... get the dimensions
!---------------------------------------------------------------------------
      ierr = pio_inq_dimid( ncid, 'nosr', dimid )
      ierr = pio_inq_dimlen( ncid, dimid, nsrc )
      ierr = pio_inq_dimid( ncid, 'ngast', dimid )
      ierr = pio_inq_dimlen( ncid, dimid, ngast )
      ierr = pio_inq_dimid( ncid, 'nla', dimid )
      ierr = pio_inq_dimlen( ncid, dimid, nla )
!---------------------------------------------------------------------------
! 	... allocate arrays
!---------------------------------------------------------------------------
      allocate( wlint(nsrc), xso2int(nsrc), x1(nsrc), y1(nsrc), stat=astat )
      if( astat /= 0 ) then
         write(iulog,*) 'o2_xsect_inti: failed to allocate wlint ... y1; error = ',astat
         call endrun
      end if
      allocate( wlgast(ngast), wlla(nla), stat=astat )
      if( astat /= 0 ) then
         write(iulog,*) 'o2_xsect_inti: failed to allocate wlgast, wlla; error = ',astat
         call endrun
      end if
!---------------------------------------------------------------------------
! 	... read the wave bin coordinates
!---------------------------------------------------------------------------
      ierr = pio_inq_varid( ncid, 'wl_src', vid )
      ierr = pio_get_var( ncid, vid, x1 )
      ierr = pio_inq_varid( ncid, 'xs_src', vid )
      ierr = pio_get_var( ncid, vid, y1 )
      ierr = pio_inq_varid( ncid, 'wl_gast', vid )
      ierr = pio_get_var( ncid, vid, wlgast )
      ierr = pio_inq_varid( ncid, 'wl_lym', vid )
      ierr = pio_get_var( ncid, vid, wlla )
      call pio_closefile(  ncid )
!-----------------------------------------------------------------------------
! 	... put together the internal grid by "pasting" the lyman-alpha grid and 
!           kockarts grid into the combination of brasseur/solomon and jpl grid
!-----------------------------------------------------------------------------
      wlint(1:9) = x1(1:9)
      nwint = 9
      wlint(nwint+1:nwint+2) = wlla(1:2)
      nwint = 11
      wlint(nwint+1:nwint+36) = x1(12:47)
      nwint = 47
      wlint(nwint+1:nwint+ngast) = wlgast(1:ngast)
      nwint = nwint + ngast
      wlint(nwint+1:nwint+41) = x1(65:105)
      nwint = nwint + 41
      wrk_ind(1:4) = (/ nsrc, ngast, nla, nwint /)
!-----------------------------------------------------------------------------
! 	... initialize interpolation module
!-----------------------------------------------------------------------------
      call inter_inti( nw+1, wl, nsrc, wlint )
!-----------------------------------------------------------------------------
! 	... interpolate brasseur/solomon and jpl data onto internal grid
!-----------------------------------------------------------------------------
      call inter2( nsrc, wlint, xso2int, nsrc, x1, y1, iret )
      deallocate( x1, y1 )


      end subroutine o2_xsect_inti

      subroutine set_o2_xsect( z, nw, wl, cz, &
                               vcol, scol, dto2, xso2 )
!-----------------------------------------------------------------------------
!   purpose:
!   compute equivalent optical depths for o2 absorption, parameterized in
!   the sr bands and the lyman-alpha line.
!-----------------------------------------------------------------------------
!   parameters:
!   nz      - integer, number of specified altitude levels in the working (i)
!             grid
!   z       - real(r8), specified altitude working grid (km)
!   nw      - integer, number of specified intervals + 1 in working
!             wavelength grid
!   wl      - real(r8), vector of lower limits of wavelength intervals in
!             working wavelength grid
!   cz      - real(r8), number of air molecules per cm^2 at each specified
!             altitude layer
!   zen     - real(r8), solar zenith angle
!   dto2    - real(r8), optical depth due to o2 absorption at each specified  (o)
!             vertical layer at each specified wavelength
!   xso2    - real(r8), molecular absorption cross section in sr bands at     (o)
!             each specified altitude and wavelength.  includes herzberg
!             continuum.
!-----------------------------------------------------------------------------
!   edit history:
!   02/98  included lyman-alpha parameterization
!   03/97  fix dto2 problem at top level (nz)
!   02/97  changed offset for grid-end interpolation to relative number
!          (x * (1 +- deltax))
!   08/96  modified for early exit, no redundant read of data and smaller
!          internal grid if possible;  internal grid uses user grid points
!          whenever possible
!   07/96  modified to work on internal grid and interpolate final values
!          onto the user-defined grid
!-----------------------------------------------------------------------------

      use mo_params,  only : kw
      use mo_wavelen, only : delw_bin
      use mo_inter,   only : inter3
      use mo_schu,    only : schu
      use mo_lymana,  only : lymana
      use ppgrid,     only : pver, pverp

      implicit none

!-----------------------------------------------------------------------------
!	... dummy arguments
!-----------------------------------------------------------------------------
      integer, intent(in)    :: nw
      real(r8), intent(in)   :: wl(kw)
      real(r8), intent(in)   :: cz(pverp)
      real(r8), intent(in)   :: z(pverp)
      real(r8), intent(in)   :: vcol(pverp)
      real(r8), intent(in)   :: scol(pverp)
      real(r8), intent(out)  :: dto2(pver,nw)
      real(r8), intent(out)  :: xso2(nw,pverp)

!-----------------------------------------------------------------------------
!	... local variables
!-----------------------------------------------------------------------------
      integer  :: wn, k, igast
      integer  :: astat
      real(r8) :: secchi(pverp)

!-----------------------------------------------------------------------------
! 	... o2 optical depth and equivalent cross section on kockarts grid
!-----------------------------------------------------------------------------
      real(r8), allocatable :: dto2k(:,:)
      real(r8), allocatable :: xso2k(:,:)
!-----------------------------------------------------------------------------
! 	... o2 optical depth and equivalent cross section in the lyman-alpha region
!-----------------------------------------------------------------------------
      real(r8), allocatable :: dto2la(:,:)
      real(r8), allocatable :: xso2la(:,:)
!-----------------------------------------------------------------------------
! 	... temporary one-dimensional storage for optical depth and cross section values
!           xxtmp  - on internal grid
!           xxuser - on user defined grid
!-----------------------------------------------------------------------------
      real(r8), dimension(2*kw) :: dttmp, xstmp
      real(r8) :: dtuser(kw)
      real(r8) :: xsuser(kw)
      real(r8) :: o2col(pverp)

      real(r8) :: x, y
      real(r8) :: delo2

!-----------------------------------------------------------------------------
!	... allocate local variables
!-----------------------------------------------------------------------------
      allocate( dto2k(pver,ngast-1), xso2k(pverp,ngast-1), stat=astat )
      if( astat /= 0 ) then
         write(iulog,*) 'set_o2_xsect: failed to allocate dto2k,xso2k; error = ',astat
         call endrun
      end if
      allocate( dto2la(pver,nla-1), xso2la(pverp,nla-1), stat=astat )
      if( astat /= 0 ) then
         write(iulog,*) 'set_o2_xsect: failed to allocate dto2k,xso2k; error = ',astat
         call endrun
      end if
!-----------------------------------------------------------------------------
! 	... check, whether user grid is in the o2 absorption band at all...
!           if not, set cross section and optical depth values to zero and return
!-----------------------------------------------------------------------------
      dto2(:pver,:nw)  = 0._r8
      xso2(:nw,:pverp) = 0._r8
      if( wl(1) > 243._r8 ) then
         return
      end if

!-----------------------------------------------------------------------------
! 	... sec xhi or chapman calculation
!           for zen > 95 degrees, use zen = 95.  (this is only to compute effective o2
!           cross sections. still, better than setting dto2 = 0. as was done up to 
!           version 4.0) sm 1/2000
!           in future could replace with mu2(iz) (but mu2 is also wavelength-depenedent)
!           or imporved chapman function 
!-----------------------------------------------------------------------------

!-----------------------------------------------------------------------------
! 	... slant o2 column 
!-----------------------------------------------------------------------------
      o2col(1:pverp) = 0.2095_r8 * scol(1:pverp)

!-----------------------------------------------------------------------------
! 	... effective secant of solar zenith angle.  use 2.0 if no direct sun. 
!           for nz, use value at nz-1
!-----------------------------------------------------------------------------
      secchi(1:pver) = scol(1:pver)/vcol(1:pver)
      where( secchi(1:pver) == 0._r8 )
         secchi(1:pver) = 2._r8
      endwhere
      secchi(pverp) = secchi(pver)

!-----------------------------------------------------------------------------
! 	... if necessary:
!           kockarts parameterization of the sr bands, output values of o2
!           optical depth and o2 equivalent cross section are on his grid
!-----------------------------------------------------------------------------
      if( wl(1) < wlgast(ngast) .and. wl(nw+1) > wlgast(1) ) then
           call schu( o2col, secchi, dto2k, xso2k )
      else
         dto2k(:,:) = 0._r8
         xso2k(:,:) = 0._r8
      end if

!-----------------------------------------------------------------------------
! 	... lyman-alpha parameterization, output values of o2 opticaldepth
!           and o2 effective (equivalent) cross section
!-----------------------------------------------------------------------------
      if( wl(1) <= wlla(nla) .and. wl(nw+1) >= wlla(1) ) then
         call lymana( o2col, secchi, dto2la, xso2la )
      else
         dto2la(:,:) = 0._r8
         xso2la(:,:) = 0._r8
      end if

!-----------------------------------------------------------------------------
! 	... loop through the altitude levels
!-----------------------------------------------------------------------------
level_loop : &
      do k = 1,pverp
         igast = 0
!-----------------------------------------------------------------------------
! 	... loop through the internal wavelength grid
!-----------------------------------------------------------------------------
         do wn = 1,nwint-1
!-----------------------------------------------------------------------------
! 	... if outside kockarts grid and outside lyman-alpha, use the 
!           jpl/brasseur+solomon data, if inside
!           kockarts grid, use the parameterized values from the call to schu,
!           if inside lyman-alpha, use the paraemterized values from call to lymana
!-----------------------------------------------------------------------------
            if( wlint(wn+1) <= wlgast(1) .or. wlint(wn) >= wlgast(ngast) ) then
              if( wlint(wn+1) <= wlla(1) .or. wlint(wn) >= wlla(nla) ) then
                 xstmp(wn) = xso2int(wn)
              else
                 xstmp(wn) = xso2la(k,1)
              end if
            else
               igast = igast + 1
               xstmp(wn) = xso2k(k,igast)
            end if
!-----------------------------------------------------------------------------
! 	... compute the area in each bin (for correct interpolation purposes only!)
!-----------------------------------------------------------------------------
            xstmp(wn) = xstmp(wn) * (wlint(wn+1) - wlint(wn))
         end do
!-----------------------------------------------------------------------------
! 	... interpolate o2 cross section from the internal grid onto the user grid
!-----------------------------------------------------------------------------
         call inter3( nw+1, wl, xsuser, nwint, wlint, xstmp )
         xso2(:nw,k) = xsuser(:nw) * delw_bin(:nw)
      end do level_loop

      do k = 1,pver
         igast = 0
         delo2 = .2095_r8 * cz(k)    ! vertical o2 column
!-----------------------------------------------------------------------------
! 	... loop through the internal wavelength grid
!-----------------------------------------------------------------------------
         do wn = 1,nwint-1
!-----------------------------------------------------------------------------
! 	... if outside kockarts grid and outside lyman-alpha, use the 
!           jpl/brasseur+solomon data, if inside
!           kockarts grid, use the parameterized values from the call to schu,
!           if inside lyman-alpha, use the paraemterized values from call to lymana
!-----------------------------------------------------------------------------
            if( wlint(wn+1) <= wlgast(1) .or. wlint(wn) >= wlgast(ngast) ) then
              if( wlint(wn+1) <= wlla(1) .or. wlint(wn) >= wlla(nla) ) then
                 dttmp(wn) = xso2int(wn) * delo2
              else
                 dttmp(wn) = dto2la(k,1)
              end if
            else
               igast = igast + 1
               dttmp(wn) = dto2k(k,igast)
            end if
!-----------------------------------------------------------------------------
! 	... compute the area in each bin (for correct interpolation purposes only!)
!-----------------------------------------------------------------------------
            dttmp(wn) = dttmp(wn) * (wlint(wn+1) - wlint(wn))
         end do
!-----------------------------------------------------------------------------
! 	... interpolate o2 optical depth from the internal grid onto the user grid
!-----------------------------------------------------------------------------
         call inter3( nw+1, wl, dtuser, nwint, wlint, dttmp )
         dto2(k,:nw) = dtuser(:nw) * delw_bin(:nw)
      end do

      deallocate( dto2k, xso2k, dto2la, xso2la )

      end subroutine set_o2_xsect

      end module mo_seto2
