subroutine sugas_corrk

!==================================================================
!     
!     Purpose
!     -------
!     Set up gaseous absorption parameters used by the radiation code.
!     This subroutine is a replacement for the old 'setrad', which contained
!     both absorption and scattering data.
!     
!     Authors
!     -------
!     Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009)
!     
!     Summary
!     -------
!     
!==================================================================

      use radinc_h
      use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
      use radcommon_h, only : tgasref,tgasmin,tgasmax
      use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
      use radcommon_h, only : wrefvar,gastype
      use datafile_mod, only: datadir

      implicit none

#include "callkeys.h"

!==================================================================

      logical file_ok

      integer n, nt, np, nh, ng, nw, m, i
      integer L_NGAUSScheck, L_NPREFcheck, L_NTREFcheck, L_REFVARcheck

      character(len=100) :: file_id
      character(len=100) :: file_path

      real*8 gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS)
      real*8 gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS)

      real*8 x, xi(4), yi(4), ans, p

      integer Nspecies

      real*8 Xtemp

!=======================================================================
!     Load variable species data, exit if we have wrong database
      file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/Q.dat'
      file_path=TRIM(datadir)//TRIM(file_id)
      ! check that the file exists
      inquire(FILE=file_path,EXIST=file_ok)
      if(.not.file_ok) then
         write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
         write(*,*)'was not found by sugas_corrk.F90, exiting.'
         call abort
      endif

      ! check that database matches varactive toggle
      open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
      read(111,*) Nspecies

!     if(Nspecies.eq.1 .and. (varactive.or.varfixed))then
!        print*,'You have varactive/fixed=.true. but the database [',  &
!                    corrkdir(1:LEN_TRIM(corrkdir)),           &
!               '] has no variable species, exiting.'
!        call abort
!     elseif(Nspecies.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
!        print*,'You have varactive and varfixed=.false. and the database [', &
!                    corrkdir(1:LEN_TRIM(corrkdir)),           &
!               '] has a variable species.'
!        call abort
!     elseif(Nspecies.gt.3 .or. Nspecies.lt.1)then
!        print*,Nspecies,' species in database [',               &
!                    corrkdir(1:LEN_TRIM(corrkdir)),           &
!               '], radiative code cannot handle this.'
!        call abort
!     endif

      do n=1,Nspecies
         read(111,*) gastype(n)
         print*,'Gas ',n,' is ',gastype(n)
      enddo

      ! check the array size is correct, load the coefficients
      open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
      read(111,*) L_REFVARcheck
      if(.not.(L_REFVARcheck.eq.L_REFVAR)) then      
         print*,'The size of your radiative transfer mixing ratio array does '
         print*,'not match the value given in Q.dat, exiting.'
         call abort
      endif
      read(111,*) wrefvar
      close(111)

      ! display the values
      print*,'Variable gas mixing ratios:'
      do n=1,L_REFVAR
         print*,n,'.',wrefvar(n),' kg/kg'
      end do
      print*,''


!=======================================================================
!     Set the weighting in g-space

      file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/g.dat'
      file_path=TRIM(datadir)//TRIM(file_id)

      ! check that the file exists
      inquire(FILE=file_path,EXIST=file_ok)
      if(.not.file_ok) then
         write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
         write(*,*)'was not found by sugas_corrk.F90, exiting.'
         call abort
      endif
      
      ! check the array size is correct, load the coefficients
      open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
      read(111,*) L_NGAUSScheck
      if(.not.(L_NGAUSScheck.eq.L_NGAUSS)) then
         print*,'The size of your radiative transfer g-space array does '
         print*,'not match the value given in g.dat, exiting.'
         call abort
      endif
      read(111,*) gweight
      close(111)

      ! display the values
      print*,'Correlated-k g-space grid:'
      do n=1,L_NGAUSS
         print*,n,'.',gweight(n)
      end do
      print*,''

!=======================================================================
!     Set the reference pressure and temperature arrays.  These are
!     the pressures and temperatures at which we have k-coefficients.

!-----------------------------------------------------------------------
! pressure

      file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/p.dat'
      file_path=TRIM(datadir)//TRIM(file_id)

      ! check that the file exists
      inquire(FILE=file_path,EXIST=file_ok)
      if(.not.file_ok) then
         write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
         write(*,*)'was not found by sugas_corrk.F90, exiting.'
         call abort
      endif
      
      ! check the array size is correct, load the coefficients
      open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
      read(111,*) L_NPREFcheck
      if(.not.(L_NPREFcheck.eq.L_NPREF)) then
         print*,'file_path = ',file_path, L_NPREFcheck,L_NPREF
         print*,'check in radinc_h.F90'
         print*,'The size of your radiative transfer pressure array does '
         print*,'not match the value given in p.dat, exiting.'
         call abort
      endif
      read(111,*) pgasref
      close(111)

      ! display the values
      print*,'Correlated-k pressure grid (mBar):'
      do n=1,L_NPREF
         print*,n,'. 1 x 10^',pgasref(n),' mBar'
      end do
      print*,''

      ! save the min / max matrix values
      pgasmin = 10.0**pgasref(1)
      pgasmax = 10.0**pgasref(L_NPREF)

      ! interpolate to finer grid
      do n=1,L_NPREF-1
         Xtemp=pgasref(n+1)-pgasref(n)
         do m=1,5
!            pfgasref((n-1)*5+m) = pgasref(n)+(m-1)*0.2
       
            pfgasref((n-1)*5+m) = pgasref(n)+(m-1)*(Xtemp/5)

         end do
      end do
      pfgasref(L_PINT) = pgasref(L_NPREF)
      ! Warning! this may need to be generalised if we want to use uneven grids!
!      print*,'pgasref',pgasref
!      print*,'pfgasref',pfgasref


!-----------------------------------------------------------------------
! temperature

      file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/T.dat'
      file_path=TRIM(datadir)//TRIM(file_id)

      ! check that the file exists
      inquire(FILE=file_path,EXIST=file_ok)
      if(.not.file_ok) then
         write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
         write(*,*)'was not found by sugas_corrk.F90, exiting.'
         call abort
      endif

      ! check the array size is correct, load the coefficients
      open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
      read(111,*) L_NTREFcheck
      if(.not.(L_NTREFcheck.eq.L_NTREF)) then
         print*,'The size of your radiative transfer temperature array does '
         print*,'not match the value given in T.dat, exiting.'
         call abort
      endif
      read(111,*) tgasref
      close(111)

      ! display the values
      print*,'Correlated-k temperature grid:'
      do n=1,L_NTREF
         print*,n,'.',tgasref(n),' K'
      end do

      ! save the min / max matrix values
      tgasmin = tgasref(1)
      tgasmax = tgasref(L_NTREF)

!=======================================================================
!     Get gaseous k-coefficients and interpolate onto finer pressure grid

         ! VISIBLE
      if (callgasvis) then
         file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 
         file_path=TRIM(datadir)//TRIM(file_id)

         ! check that the file exists
         inquire(FILE=file_path,EXIST=file_ok)
         if(.not.file_ok) then
            write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
            write(*,*)'was not found by sugas_corrk.F90.'
            write(*,*)'Are you sure you have absorption data for these bands?'
            call abort
         endif
         
         open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
         read(111,*) gasv8
         close(111)
         
      else
         print*,'Visible gaseous absorption is set to zero.'
      endif

      ! INFRA-RED
      file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 
      file_path=TRIM(datadir)//TRIM(file_id)
      ! check that the file exists
      inquire(FILE=file_path,EXIST=file_ok)
      if(.not.file_ok) then
         write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
         write(*,*)'was not found by sugas_corrk.F90.'
         write(*,*)'Are you sure you have absorption data for these bands?'
         call abort
      endif
      
      open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
      read(111,*) gasi8
      close(111)
      
      do nw=1,L_NSPECTI
         fzeroi(nw) = 0.
      end do
      do nw=1,L_NSPECTV
         fzerov(nw) = 0.
      end do

!     Take log10 of the values - this is what we will interpolate.
!     Smallest value is 1.0E-200.

      do nt=1,L_NTREF
         do np=1,L_NPREF
            do nh=1,L_REFVAR
               do ng = 1,L_NGAUSS

                  do nw=1,L_NSPECTV
                     if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
                        gasv8(nt,np,nh,nw,ng) = &
                            log10(gasv8(nt,np,nh,nw,ng))
                     else
                        gasv8(nt,np,nh,nw,ng) = -200.0
                     end if
                  end do

                  do nw=1,L_NSPECTI
                     if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
                        gasi8(nt,np,nh,nw,ng) = &
                            log10(gasi8(nt,np,nh,nw,ng))
                     else
                        gasi8(nt,np,nh,nw,ng) = -200.0
                     end if
                  end do
                  
               end do
            end do
         end do
      end do

!     Interpolate the values:  first the longwave

      do nt=1,L_NTREF
         do nh=1,L_REFVAR
            do nw=1,L_NSPECTI
               do ng=1,L_NGAUSS

!     First, the initial interval

                  n = 1 
                  do m=1,5
                     x     = pfgasref(m)
                     xi(1) = pgasref(n)
                     xi(2) = pgasref(n+1)
                     xi(3) = pgasref(n+2)
                     xi(4) = pgasref(n+3)
                     yi(1) = gasi8(nt,n,nh,nw,ng)
                     yi(2) = gasi8(nt,n+1,nh,nw,ng)
                     yi(3) = gasi8(nt,n+2,nh,nw,ng)
                     yi(4) = gasi8(nt,n+3,nh,nw,ng)
                     call lagrange(x,xi,yi,ans)
                     gasi(nt,m,nh,nw,ng) = 10.0**ans
                  end do 
                  
                  do n=2,L_NPREF-2
                     do m=1,5
                        i     = (n-1)*5+m
                        x     = pfgasref(i)
                        xi(1) = pgasref(n-1)
                        xi(2) = pgasref(n)
                        xi(3) = pgasref(n+1)
                        xi(4) = pgasref(n+2)
                        yi(1) = gasi8(nt,n-1,nh,nw,ng)
                        yi(2) = gasi8(nt,n,nh,nw,ng)
                        yi(3) = gasi8(nt,n+1,nh,nw,ng)
                        yi(4) = gasi8(nt,n+2,nh,nw,ng)
                        call lagrange(x,xi,yi,ans)
                        gasi(nt,i,nh,nw,ng) = 10.0**ans
                     end do 
                  end do

!     Now, get the last interval

                  n = L_NPREF-1                 
                  do m=1,5
                     i     = (n-1)*5+m
                     x     = pfgasref(i)
                     xi(1) = pgasref(n-2)
                     xi(2) = pgasref(n-1)
                     xi(3) = pgasref(n)
                     xi(4) = pgasref(n+1)
                     yi(1) = gasi8(nt,n-2,nh,nw,ng)
                     yi(2) = gasi8(nt,n-1,nh,nw,ng)
                     yi(3) = gasi8(nt,n,nh,nw,ng)
                     yi(4) = gasi8(nt,n+1,nh,nw,ng)
                     call lagrange(x,xi,yi,ans)
                     gasi(nt,i,nh,nw,ng) = 10.0**ans
                  end do  

!     Fill the last pressure point

                  gasi(nt,L_PINT,nh,nw,ng) = &
                       10.0**gasi8(nt,L_NPREF,nh,nw,ng)

               end do
            end do
         end do
      end do

!     Interpolate the values:  now the shortwave

      do nt=1,L_NTREF
         do nh=1,L_REFVAR
            do nw=1,L_NSPECTV
               do ng=1,L_NGAUSS

!     First, the initial interval

                  n = 1 
                  do m=1,5
                     x     = pfgasref(m)
                     xi(1) = pgasref(n)
                     xi(2) = pgasref(n+1)
                     xi(3) = pgasref(n+2)
                     xi(4) = pgasref(n+3)
                     yi(1) = gasv8(nt,n,nh,nw,ng)
                     yi(2) = gasv8(nt,n+1,nh,nw,ng)
                     yi(3) = gasv8(nt,n+2,nh,nw,ng)
                     yi(4) = gasv8(nt,n+3,nh,nw,ng)
                     call lagrange(x,xi,yi,ans)
                     gasv(nt,m,nh,nw,ng) = 10.0**ans
                  end do 
                  
                  do n=2,L_NPREF-2
                     do m=1,5
                        i     = (n-1)*5+m
                        x     = pfgasref(i)
                        xi(1) = pgasref(n-1)
                        xi(2) = pgasref(n)
                        xi(3) = pgasref(n+1)
                        xi(4) = pgasref(n+2)
                        yi(1) = gasv8(nt,n-1,nh,nw,ng)
                        yi(2) = gasv8(nt,n,nh,nw,ng)
                        yi(3) = gasv8(nt,n+1,nh,nw,ng)
                        yi(4) = gasv8(nt,n+2,nh,nw,ng)
                        call lagrange(x,xi,yi,ans)
                        gasv(nt,i,nh,nw,ng) = 10.0**ans
                     end do 
                  end do

!     Now, get the last interval

                  n = L_NPREF-1
                  do m=1,5
                     i     = (n-1)*5+m
                     x     = pfgasref(i)
                     xi(1) = pgasref(n-2)
                     xi(2) = pgasref(n-1)
                     xi(3) = pgasref(n)
                     xi(4) = pgasref(n+1)
                     yi(1) = gasv8(nt,n-2,nh,nw,ng)
                     yi(2) = gasv8(nt,n-1,nh,nw,ng)
                     yi(3) = gasv8(nt,n,nh,nw,ng)
                     yi(4) = gasv8(nt,n+1,nh,nw,ng)
                     call lagrange(x,xi,yi,ans)
                     gasv(nt,i,nh,nw,ng) = 10.0**ans
                  end do  

!     Fill the last pressure point

                  gasv(nt,L_PINT,nh,nw,ng) = &
                      10.0**gasv8(nt,L_NPREF,nh,nw,ng)
                  
               end do
            end do
         end do
      end do

      return
    end subroutine sugas_corrk