source: trunk/LMDZ.GENERIC/libf/aeronostd/read_phototable.F90 @ 3746

!***********************************************************************

      subroutine read_phototable

!***********************************************************************
!
!   subject:
!   --------
!
!   read photolysis lookup table
!
!   VERSION: 8/10/2014
!
!   Author:   Franck Lefevre
!   update 06/03/2021 dimension set in table + CH4 dimension + photoheat (Yassin Jaziri)
!
!   Arguments:
!   ----------
!
!   The output variable is jphot/ephot and is put in common chimiedata.
!
!***********************************************************************

      use ioipsl_getincom 
      use ioipsl_getin_p_mod, only: getin_p
      use datafile_mod
      use callkeys_mod
      use chimiedata_h
      implicit none

!#include "chimiedata.h"

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     local:
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      integer :: fic, ij, iozo, isza, itemp, iz, itau, ich4, ierr
      real    :: xsza, dummy

      character(len = 128) :: phototable ! photolysis table file name
      character(len = 128) :: ephototable ! energie photolysis table file name

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! set photolysis table input file name
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



! look for a " phototable= ..." option in def files
     write(*,*) "Directory where external input files are:"
      phototable = "jmars.20140930" ! default
     call getin_p("phototable",phototable) ! default path
     write(*,*) " phototable = ",trim(phototable)


      fic = 81

      open(fic, form = 'formatted', status = 'old',                &
           file =trim(datadir)//"/"//trim(phototable),iostat=ierr)

      if (ierr /= 0) THEN
        write(*,*)'Error : cannot open photolysis lookup table ', trim(phototable)
        write(*,*)'It should be in :',trim(datadir),'/'
        write(*,*)'1) You can change this directory in callphys.def'
        write(*,*)'   with:'
        write(*,*)'   datadir=/path/to/the/directory'
        write(*,*)'2) You can change the input phototable file name in'
        write(*,*)'   callphys.def with:'
        write(*,*)'   phototable=filename'
        stop
      end if

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! read photolys table
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      print*, 'read photolysis lookup table ',trim(phototable)

      read(fic,*) nz, nsza, nozo, dummy, ntemp, ntau, nch4, nd

      allocate(jphot(ntemp,nsza,nz,nozo,ntau,nch4,nd))
      allocate(ephot(ntemp,nsza,nz,nozo,ntau,nch4,nd))
      allocate(colairtab(nz))
      allocate(szatab(nsza))
      allocate(table_ozo(nozo))
      allocate(tautab(ntau))
      allocate(table_ch4(nch4))

      do itau = 1,ntau
         do itemp = 1,ntemp
            do iozo = 1,nozo
               do ich4 =1,nch4
                  do isza = 1,nsza
                     do iz = nz,1,-1
                        read(fic,*) colairtab(iz), szatab(isza),table_ozo(iozo), dummy, dummy, dummy, table_ch4(ich4)
                        read(fic,'(7e11.4)') (jphot(itemp,isza,iz,iozo,itau,ich4,ij), ij= 1,nd)
                        do ij = 1,nd
                           if (jphot(itemp,isza,iz,iozo,itau,ich4,ij) == 1.e-30) then
                              jphot(itemp,isza,iz,iozo,itau,ich4,ij) = 0.
                           end if
                        end do
                     end do
                  end do
               end do
            end do
         end do
      end do

      print*, 'lookup table...ok'
      close(fic)

      if (photoheat) then 
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! set energie photolysis table input file name
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



! look for a " ephototable= ..." option in def files
     write(*,*) "Directory where external input files are:"
      ephototable = "zdtearth_N+_CO2_0.0004_O2_0.21_G11" ! default
     call getin_p("ephototable",ephototable) ! default path
     write(*,*) " ephototable = ",trim(ephototable)


      fic = 81

      open(fic, form = 'formatted', status = 'old',                &
           file =trim(datadir)//"/"//trim(ephototable),iostat=ierr)

      if (ierr /= 0) THEN
        write(*,*)'Error : cannot open energie photolysis lookup table ', trim(ephototable)
        write(*,*)'It should be in :',trim(datadir),'/'
        write(*,*)'1) You can change this directory in callphys.def'
        write(*,*)'   with:'
        write(*,*)'   datadir=/path/to/the/directory'
        write(*,*)'2) You can change the input ephototable file name in'
        write(*,*)'   callphys.def with:'
        write(*,*)'   ephototable=filename'
        stop
      end if

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! read energie photolys table
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      print*, 'read energie photolysis lookup table ',trim(ephototable)

      do itau = 1,ntau
         do itemp = 1,ntemp
            do iozo = 1,nozo
               do ich4 =1,nch4
                  do isza = 1,nsza
                     do iz = nz,1,-1
                        read(fic,*) colairtab(iz), xsza,table_ozo(iozo), dummy, dummy, dummy, table_ch4(ich4)
                        read(fic,'(7e11.4)') (ephot(itemp,isza,iz,iozo,itau,ich4,ij), ij= 1,nd)
                        do ij = 1,nd
                           if (ephot(itemp,isza,iz,iozo,itau,ich4,ij) == 1.e-30) then
                              ephot(itemp,isza,iz,iozo,itau,ich4,ij) = 0.
                           end if
                        end do
                     end do
                  end do
               end do
            end do
         end do
      end do

      print*, 'lookup etable...ok'
      close(fic)
      end if

      return
      end