!-------------------------------------------------------------------
! manages reading and interpolation of linoz data
! Created Dec 2020 by Juno Hsu (junoh@uci.edu) by modifying the original code by Francis Vitt (the new version can't read linoz v2 data because the named variables are hard-coded)
!             coded to read linoz v3 data: linv3_1845-2015_2010JPL_cmip_10deg_58km_c20201030.nc (monthly data every 5 year from 1845 up to 2015 generated by running the photochemial box, PRATMO (M. Prather, UCI)
!             format similar to its old v2 counterpart but for 4 species O3LNZ, N2OLNZ, NOYLNZ, CH4LNZ
!             in addition, the bottom surface level stores the surface concentration of the corresponding year for N2OLNZ and CH4LNZ, used in SFCSINK as a restoring value
!-------------------------------------------------------------------
module linoz_data
  use shr_kind_mod, only : r8 => shr_kind_r8
  use cam_abortutils,   only : endrun
  use spmd_utils,   only : masterproc
  use tracer_data,  only : trfld,trfile
  use cam_logfile,  only : iulog

  implicit none

  private  ! all unless made public
  save 

  public :: fields
  public :: linoz_data_init
  public :: linoz_data_adv
  public :: init_linoz_data_restart
  public :: write_linoz_data_restart
  public :: read_linoz_data_restart
  public :: has_linoz_data, has_linozv2_data, has_linozv3_data
  public :: linoz_data_defaultopts
  public :: linoz_data_setopts

  type(trfld), pointer :: fields(:) => null()
  type(trfile) :: file
  logical :: has_linozv2_data 
  logical :: has_linozv3_data 
  logical :: has_linoz_data 

  integer, parameter, public :: N_FLDS = 57
  integer :: number_flds

  character(len=256) :: filename = ''
  character(len=256) :: filelist = ''
  character(len=256) :: datapath = ''
  character(len=32)  :: datatype = 'CYCLICAL'
  logical            :: rmv_file = .false.
  integer            :: cycle_yr  = 0
  integer            :: fixed_ymd = 0
  integer            :: fixed_tod = 0

! The named variables in fields have to match those in the linoz nc file; if one of them is mis-spelled, the pointer will point to a wrong one without error messages.
  character(len=16), dimension(N_FLDS), parameter :: fld_names = & ! data field names
       (/'o3_clim         ','n2o_clim        ','noy_clim        ','ch4_clim        ',&!6 clim-terms
         'h2o_clim        ','t_clim          ','o3col_clim      ',&                                    
         'no3_PmL_clim    ','no3_dPmL_dO3    ','no3_dPmL_dN2O   ','no3_dPmL_dNOy   ',&!Taylar expan for o3
         'no3_dPmL_dCH4   ','no3_dPmL_dH2O   ','no3_dPmL_dT     ','no3_dPmL_dO3col ',&
         'pn2o_P_clim     ','pn2o_dP_dO3     ','pn2o_dP_dN2O    ','pn2o_dP_dNOy    ',&!for n2o pro
         'pn2o_dP_dCH4    ','pn2o_dP_dH2O    ','pn2o_dP_dT      ','pn2o_dP_dO3col  ',&
         'ln2o_L_clim     ','ln2o_dL_dO3     ','ln2o_dL_dN2O    ','ln2o_dL_dNOy    ',&!for n2o loss
         'ln2o_dL_dCH4    ','ln2o_dL_dH2O    ','ln2o_dL_dT      ','ln2o_dL_dO3col  ',&
         'pnoy_P_clim     ','pnoy_dP_dO3     ','pnoy_dP_dN2O    ','pnoy_dP_dNOy    ',&!for noy
         'pnoy_dP_dCH4    ','pnoy_dP_dH2O    ','pnoy_dP_dT      ','pnoy_dP_dO3col  ',&
         'lnoy_L_clim     ','lnoy_dL_dO3     ','lnoy_dL_dN2O    ','lnoy_dL_dNOy    ',&!for noy
         'lnoy_dL_dCH4    ','lnoy_dL_dH2O    ','lnoy_dL_dT      ','lnoy_dL_dO3col  ',&
         'nch4_PmL_clim   ','nch4_dPmL_dO3   ','nch4_dPmL_dN2O  ','nch4_dPmL_dNOy  ',&!for ch4
         'nch4_dPmL_dCH4  ','nch4_dPmL_dH2O  ','nch4_dPmL_dT    ','nch4_dPmL_dO3col',&
         'cariolle_pscs   ','o3lbs           '/)
  

!added o3lbs that stores prescribed o3 for CMIP surface layers evolving from 1850-2015 at L=end,and repeating at L=end-1 (237hPa) so the values won't be 
!significantly changed due to vertical interpolations from Linoz levels to model surface levels  

  character(len=16), dimension(N_FLDS), parameter :: fld_units = & ! data field names
       (/'vmr             ','vmr             ','vmr             ','vmr             ',&!clim
         'vmr             ','K               ','Dobson Units    ',                   &
         'vmr/s           ','vmr/vmr/s       ','vmr/vmr/s       ','vmr/vmr/s       ',&!net o3
         'vmr/vmr/s       ','vmr/vmr/s       ','vmr/K/s         ','vmr/DU/s        ',&
         'vmr/s           ','vmr/vmr/s       ','vmr/vmr/s       ','vmr/vmr/s       ',&!pn2o
         'vmr/vmr/s       ','vmr/vmr/s       ','vmr/K/s         ','vmr/DU/s        ',&
         '1/s             ','1/vmr/s         ','1/vmr/s         ','1/vmr/s         ',&!ln2o Taylor series expanding the frequency below 
         '1/vmr/s         ','1/vmr/s         ','1/K/s           ','1/DU/s          ',&
         '1/s             ','1/vmr/s         ','1/vmr/s         ','1/vmr/s         ',&!pnoy
         '1/vmr/s         ','1/vmr/s         ','1/K/s           ','1/DU/s          ',&
         '1/s             ','1/vmr/s         ','1/vmr/s         ','1/vmr/s         ',&!lnoy
         '1/vmr/s         ','1/vmr/s         ','1/K/s           ','1/DU/s          ',&
         '1/s             ','1/vmr/s         ','1/vmr/s         ','1/vmr/s         ',&!lnoy
         '1/vmr/s         ','1/vmr/s         ','1/K/s           ','1/DU/s          ',&
         '1/s             ','vmr             '/)

  integer :: index_map(N_FLDS)

  integer, public, parameter ::  o3_clim_ndx    =      1
  integer, public, parameter ::  n2o_clim_ndx   =      2
  integer, public, parameter ::  noy_clim_ndx   =      3
  integer, public, parameter ::  ch4_clim_ndx   =      4
  integer, public, parameter ::  h2o_clim_ndx   =      5
  integer, public, parameter ::  t_clim_ndx     =      6
  integer, public, parameter ::  o3col_clim_ndx =      7

  Integer, public, parameter :: nO3_PmL_clim_ndx  =    8
  integer, public, parameter :: nO3_dPmL_dO3_ndx  =    9
  integer, public, parameter :: nO3_dPmL_dN2O_ndx =   10
  integer, public, parameter :: nO3_dPmL_dNOy_ndx =   11
  integer, public, parameter :: nO3_dPmL_dCH4_ndx =   12
  integer, public, parameter :: nO3_dPmL_dH2O_ndx =   13
  integer, public, parameter :: nO3_dPmL_dT_ndx =     14
  integer, public, parameter :: nO3_dPmL_dO3col_ndx = 15

  integer, public, parameter :: pn2o_PmL_clim_ndx  =  16 
  integer, public, parameter :: pn2o_dPmL_dO3_ndx  =  17
  integer, public, parameter :: pn2o_dPmL_dN2O_ndx =  18  
  integer, public, parameter :: pn2o_dPmL_dNOy_ndx =  19
  integer, public, parameter :: pn2o_dPmL_dCH4_ndx =  20
  integer, public, parameter :: pn2o_dPmL_dH2O_ndx =  21
  integer, public, parameter :: pn2o_dPmL_dT_ndx   =  22
  integer, public, parameter :: pn2o_dPmL_dO3col_ndx = 23

  integer, public, parameter :: ln2o_PmL_clim_ndx  =   24
  integer, public, parameter :: ln2o_dPmL_dO3_ndx  =   25
  integer, public, parameter :: ln2o_dPmL_dN2O_ndx =   26
  integer, public, parameter :: ln2o_dPmL_dNOy_ndx =   27
  integer, public, parameter :: ln2o_dPmL_dCH4_ndx =   28
  integer, public, parameter :: ln2o_dPmL_dH2O_ndx =   29
  integer, public, parameter :: ln2o_dPmL_dT_ndx =     30
  integer, public, parameter :: ln2o_dPmL_dO3col_ndx = 31

  integer, public, parameter :: pnoy_PmL_clim_ndx  =   32
  integer, public, parameter :: pnoy_dPmL_dO3_ndx  =   33 
  integer, public, parameter :: pnoy_dPmL_dN2O_ndx =   34
  integer, public, parameter :: pnoy_dPmL_dNOy_ndx =   35
  integer, public, parameter :: pnoy_dPmL_dCH4_ndx =   36
  integer, public, parameter :: pnoy_dPmL_dH2O_ndx =   37
  integer, public, parameter :: pnoy_dPmL_dT_ndx =     38
  integer, public, parameter :: pnoy_dPmL_dO3col_ndx = 39

  integer, public, parameter :: lnoy_PmL_clim_ndx  =   40
  integer, public, parameter :: lnoy_dPmL_dO3_ndx  =   41
  integer, public, parameter :: lnoy_dPmL_dN2O_ndx =   42
  integer, public, parameter :: lnoy_dPmL_dNOy_ndx =   43
  integer, public, parameter :: lnoy_dPmL_dCH4_ndx =   44
  integer, public, parameter :: lnoy_dPmL_dH2O_ndx =   45
  integer, public, parameter :: lnoy_dPmL_dT_ndx =     46
  integer, public, parameter :: lnoy_dPmL_dO3col_ndx = 47

  integer, public, parameter :: nch4_PmL_clim_ndx  =   48
  integer, public, parameter :: nch4_dPmL_dO3_ndx  =   49
  integer, public, parameter :: nch4_dPmL_dN2O_ndx =   50
  integer, public, parameter :: nch4_dPmL_dNOy_ndx =   51
  integer, public, parameter :: nch4_dPmL_dCH4_ndx =   52
  integer, public, parameter :: nch4_dPmL_dH2O_ndx =   53
  integer, public, parameter :: nch4_dPmL_dT_ndx =     54
  integer, public, parameter :: nch4_dPmL_dO3col_ndx = 55
  integer, public, parameter :: cariolle_pscs_ndx =    56
  integer, public, parameter :: o3lbs_ndx         =    57

  
contains

!-------------------------------------------------------------------
!-------------------------------------------------------------------
  subroutine linoz_data_init()

    use tracer_data, only : trcdata_init
    use cam_history, only : addfld
    use ppgrid,      only : pver
    use error_messages, only: handle_err
    use ppgrid,         only: pcols, pver, begchunk, endchunk
    use physics_buffer, only : physics_buffer_desc

    implicit none

    integer :: ndx, istat, i
    
    if ( has_linoz_data ) then
       if ( masterproc ) then
          write(iulog,*) 'linoz_data_ini: linoz data :'//trim(filename)
       endif
    else
       return
    endif

    allocate(file%in_pbuf(size(fld_names)))
    file%in_pbuf(:) = .false.
    call trcdata_init( fld_names, filename, filelist, datapath, fields, file, &
                       rmv_file, cycle_yr, fixed_ymd, fixed_tod, datatype)
        
    number_flds = 0
    if (associated(fields)) number_flds = size( fields )

    if( number_flds < 1 ) then
       if ( masterproc ) then
          write(iulog,*) 'linoz_data_init: There are no linoz data'
          write(iulog,*) ' '
       endif
       return
    end if
!    write(iulog,*)'number_flds=',number_flds
    do i = 1,number_flds
       ndx = get_ndx( fields(i)%fldnam )
       index_map(i) = ndx
!       write(iulog,*)'i=',i,'ndx=',ndx,'index_map(i)=',index_map(i)
       if (ndx < 1) then
          call endrun('linoz_data_init: '//trim(fields(i)%fldnam)//' is not one of the named linoz data fields ')
       endif
!       write(iulog,*)'fld_names(i)=',fld_names(i)

       call addfld(fld_names(i), (/ 'lev' /), 'I', fld_units(i), 'linoz data' )
    enddo


  end subroutine linoz_data_init

!-------------------------------------------------------------------
!-------------------------------------------------------------------
  subroutine linoz_data_setopts(&
       linoz_data_file_in,      &
       linoz_data_filelist_in,  &
       linoz_data_path_in,      &
       linoz_data_type_in,      &
       linoz_data_rmfile_in,    &
       linoz_data_cycle_yr_in,  &
       linoz_data_fixed_ymd_in, &
       linoz_data_fixed_tod_in  &
       )

    implicit none

    character(len=*), intent(in), optional :: linoz_data_file_in
    character(len=*), intent(in), optional :: linoz_data_filelist_in
    character(len=*), intent(in), optional :: linoz_data_path_in
    character(len=*), intent(in), optional :: linoz_data_type_in
    logical,          intent(in), optional :: linoz_data_rmfile_in
    integer,          intent(in), optional :: linoz_data_cycle_yr_in
    integer,          intent(in), optional :: linoz_data_fixed_ymd_in
    integer,          intent(in), optional :: linoz_data_fixed_tod_in

    if ( present(linoz_data_file_in) ) then
       filename = linoz_data_file_in
     endif
    if ( present(linoz_data_filelist_in) ) then
       filelist = linoz_data_filelist_in
     endif
    if ( present(linoz_data_path_in) ) then
       datapath = linoz_data_path_in
     endif
    if ( present(linoz_data_type_in) ) then
       datatype = linoz_data_type_in
     endif
    if ( present(linoz_data_rmfile_in) ) then
       rmv_file = linoz_data_rmfile_in
     endif
    if ( present(linoz_data_cycle_yr_in) ) then
       cycle_yr = linoz_data_cycle_yr_in
     endif
    if ( present(linoz_data_fixed_ymd_in) ) then
       fixed_ymd = linoz_data_fixed_ymd_in
     endif
    if ( present(linoz_data_fixed_tod_in) ) then
       fixed_tod = linoz_data_fixed_tod_in
    endif
    
      has_linozv3_data= .false.
      has_linozv2_data= .false.
      has_linoz_data =  .false.
      if (len_trim(filename) > 0 .or. len_trim(filelist) >0 )then
         if ((filename(1:5) .eq. 'linv3') .or. (filelist(1:5) .eq. 'linv3'))then
            has_linozv3_data =.true.
            has_linozv2_data =.true.
            has_linoz_data=   .true.
         endif
      endif
!       write(iulog,*)'has_linoz_data=', has_linoz_data

  endsubroutine linoz_data_setopts

!-------------------------------------------------------------------
!-------------------------------------------------------------------
 

  subroutine linoz_data_defaultopts(   &
       linoz_data_file_out,     &
       linoz_data_filelist_out, &
       linoz_data_path_out,     &
       linoz_data_type_out,     &
       linoz_data_rmfile_out,   &
       linoz_data_cycle_yr_out, &
       linoz_data_fixed_ymd_out,&
       linoz_data_fixed_tod_out &
       ) 

    implicit none

    character(len=*), intent(out), optional :: linoz_data_file_out
    character(len=*), intent(out), optional :: linoz_data_filelist_out
    character(len=*), intent(out), optional :: linoz_data_path_out
    character(len=*), intent(out), optional :: linoz_data_type_out
    logical,          intent(out), optional :: linoz_data_rmfile_out
    integer,          intent(out), optional :: linoz_data_cycle_yr_out
    integer,          intent(out), optional :: linoz_data_fixed_ymd_out
    integer,          intent(out), optional :: linoz_data_fixed_tod_out

    if ( present(linoz_data_file_out) ) then
       linoz_data_file_out = filename
    endif
    if ( present(linoz_data_filelist_out) ) then
       linoz_data_filelist_out = filelist
    endif
    if ( present(linoz_data_path_out) ) then
       linoz_data_path_out = datapath
    endif
    if ( present(linoz_data_type_out) ) then
       linoz_data_type_out = datatype
    endif
    if ( present(linoz_data_rmfile_out) ) then
       linoz_data_rmfile_out = rmv_file
    endif
    if ( present(linoz_data_cycle_yr_out) ) then
       linoz_data_cycle_yr_out = cycle_yr
    endif
    if ( present(linoz_data_fixed_ymd_out) ) then
       linoz_data_fixed_ymd_out = fixed_ymd
    endif
    if ( present(linoz_data_fixed_tod_out) ) then
       linoz_data_fixed_tod_out = fixed_tod
    endif

  endsubroutine linoz_data_defaultopts

!-------------------------------------------------------------------
!-------------------------------------------------------------------
  subroutine linoz_data_adv( pbuf2d, state )

    use tracer_data,  only : advance_trcdata
    use physics_types,only : physics_state
    use ppgrid,       only : begchunk, endchunk
    use ppgrid,       only : pcols, pver
    use string_utils, only : to_lower, GLC
    use cam_history,  only : outfld
    use physconst,    only : boltz                ! J/K/molecule
    use physics_buffer, only : physics_buffer_desc

    implicit none

  ! args
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
    type(physics_state), intent(in):: state(begchunk:endchunk)                 

  ! local vars
    integer :: ind,c,ncol,i
    real(r8) :: to_mmr(pcols,pver)

    if( .not. has_linoz_data ) return

    call advance_trcdata( fields, file, state, pbuf2d  )
    
    ! set the tracer fields with the correct units
    do i = 1,number_flds
       ind = index_map(i)
       do c = begchunk,endchunk
          ncol = state(c)%ncol
          call outfld( fields(i)%fldnam, fields(i)%data(:ncol,:,c), ncol, state(c)%lchnk )
       enddo
    enddo

  end subroutine linoz_data_adv

!-------------------------------------------------------------------
!-------------------------------------------------------------------
  subroutine init_linoz_data_restart( piofile )
    use pio, only : file_desc_t
    use tracer_data, only : init_trc_restart
    implicit none
    type(file_desc_t),intent(inout) :: piofile     ! pio File pointer

    call init_trc_restart( 'linoz_data', piofile, file )

  end subroutine init_linoz_data_restart
!-------------------------------------------------------------------
  subroutine write_linoz_data_restart( PioFile )
    use tracer_data, only : write_trc_restart
    use pio, only : file_desc_t
    implicit none

    type(file_desc_T) :: piofile

    call write_trc_restart( piofile, file )

  end subroutine write_linoz_data_restart

!-------------------------------------------------------------------
!-------------------------------------------------------------------
  subroutine read_linoz_data_restart( PioFile )
    use tracer_data, only : read_trc_restart
    use pio, only : file_desc_t
    implicit none

    type(file_desc_T) :: piofile

    call read_trc_restart( 'linoz_data', piofile, file )

  end subroutine read_linoz_data_restart

!-------------------------------------------------------------------
!-------------------------------------------------------------------
  integer function get_ndx( name )

    implicit none
    character(len=*), intent(in) :: name

    integer :: i

    get_ndx = 0
    do i = 1,N_FLDS
      if ( trim(name) == trim(fld_names(i)) ) then
        get_ndx = i
        return
      endif
    enddo

  end function get_ndx

end module linoz_data
