source: trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_asis.F90 @ 2613

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

      subroutine photolysis_asis(nlayer, ngrid,                    &
                                 lswitch, press, temp, sza,fractcol, tauref,   &
                                 zmmean, dist_sol, rmco2, rmo3, rmch4, v_phot, e_phot, nreact, three_prod)

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

      use comcstfi_mod
      use callkeys_mod
      use types_asis
      use chimiedata_h

      implicit none

!#include "chimiedata.h" 

!==========================================================================
!     input:
!==========================================================================
       
      integer, intent(in) :: nlayer ! number of atmospheric layers
      integer, intent(in) :: ngrid  ! number of atmospheric columns
      integer :: lswitch            ! interface level between chemistries
      real :: press(nlayer)         ! pressure (hPa)
      real :: temp(nlayer)          ! temperature (K)
      real :: sza                   ! solar zenith angle (deg)
      real :: fractcol              ! day fraction
      real :: tauref                ! optical depth at 7 hpa
      real :: zmmean(nlayer)        ! mean molecular mass (g/mol)
      real :: dist_sol              ! sun distance (AU) 
      real :: rmco2(nlayer)         ! co2 mixing ratio
      real :: rmo3(nlayer)          ! ozone mixing ratio
      real :: rmch4(nlayer)         ! ch4 mixing ratio
      integer, intent(in) :: nreact             ! number of reactions in reactions files
      logical, intent(in) :: three_prod(nreact) ! if the reaction have a three defferent products egal .true.

!==========================================================================
!     output: interpolated photodissociation rates (s-1)
!==========================================================================

      real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot
      real (kind = 8), dimension(nlayer,nb_phot_max) :: e_phot

!==========================================================================
!     local:
!==========================================================================

      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,     &
                 iozo, isza, itau, it, ich4, indch4, l, nb_phot 

      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
      real :: colch4(nlayer)               ! overhead CH4 column (molecule cm-2)
      real :: tauch4(nlayer)               ! estimation of optical depth by CH4
      real :: ch4_equ(nlayer)               ! equivalent constant mixing ratio for the same column of CH4
      real :: poids(2,2,2,2,2,2)          ! 6D interpolation weights 
      real :: tref                        ! temperature  at 1.9 hPa in the gcm (K)
      real :: table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
      real :: ch4ref                      ! ch4 mixing ratio at top of the atmosphere
      real :: cinf, csup, cicol, ciozo, cisza, citemp, citau, cich4
      real :: colo3min, dp, coef
      real :: ratio_o3(nlayer)
      real :: tau
      real :: j(nlayer,nd)
      real :: e(nlayer,nd)

!==========================================================================
!     day/night criterion
!==========================================================================

      if (sza <= 95.) then 

!==========================================================================
!     temperatures at 1.9 hPa in lookup table
!==========================================================================
      
!      table_temp(1) = 226.2
!      table_temp(2) = 206.2
!      table_temp(3) = 186.2
!      table_temp(4) = 169.8

!      table_temp(2) = 186.2
      table_temp(1) = 176.2

!==========================================================================
!     interpolation in solar zenith angle
!==========================================================================
 
      indsza = nsza - 1
      do isza = 1,nsza
         if (szatab(isza) >= sza) then
            indsza = isza - 1
            indsza = max(indsza, 1)
            exit
         end if
      end do

      if(nsza.eq.1) then
        cisza = 0.
        indsza=1
      else
        cisza = (sza - szatab(indsza))  &
             /(szatab(indsza + 1) - szatab(indsza))
      endif

!==========================================================================
!     interpolation in dust (tau)
!==========================================================================

      tau = min(tauref, tautab(ntau))
      tau = max(tau, tautab(1))

      indtau = ntau - 1
      do itau = 1,ntau
         if (tautab(itau) >= tau) then
            indtau = itau - 1
            indtau = max(indtau, 1)
            exit
         end if
      end do

      if(ntau.eq.1) then
        citau=0.
        indtau=1
      else
        citau = (tau - tautab(indtau))     &
             /(tautab(indtau + 1) - tautab(indtau))
      endif



!==========================================================================
!     air and ozone columns
!==========================================================================

!     co2 column at model top (molecule.cm-2)

!      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
!                       /(zmmean(lswitch-1)*g)
      col(lswitch-1) = 6.022e22*press(lswitch-1)*100./(zmmean(lswitch-1)*g)


!     ozone column at model top

      colo3(lswitch-1) = 0.
!     co2 and ozone columns for other levels (molecule.cm-2)

      do l = lswitch-2,1,-1
         dp = (press(l) - press(l+1))*100.
!         col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5   &
!                             *6.022e22*dp/(zmmean(l)*g)
         col(l) = col(l+1) +  6.022e22*dp/(zmmean(l)*g)
         col(l) = min(col(l), colairtab(1))
         colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 &
                                 *6.022e22*dp/(zmmean(l)*g)

      end do

!     ratio ozone column/minimal theoretical column (10 micron-atm) 

!     ro3 = 1.227e-10 is the o3 mixing ratio for a uniform
!     profile giving a column 10 micron-atmosphere at
!     a surface pressure of 1 bar.

      do l = 1,lswitch-1
!         colo3min    = col(l)*7.171e-10
         colo3min    = col(l)*1.227e-10*(g/9.81)*(mugaz/28)
         ratio_o3(l) = colo3(l)/colo3min
         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo))
         ratio_o3(l) = max(ratio_o3(l), 0.)
      end do

!      print*,'co3(1)=',colo3(1)
!      print*,'col(1)=',col(1)
!      print*,'ratio_o3(1)=',ratio_o3(1)
!      print*,'maxval(ratio_o3)=',maxval(ratio_o3(:))
!      print*,'maxval(ozo)=',table_ozo(nozo)/10.

!==========================================================================
!     ch4 dependence
!==========================================================================

!     1) search for temperature at 1.9 hPa (tref): vertical
!     interpolation

      ch4ref = rmch4(lswitch-2)
      colch4(lswitch-1) = 0.
      ch4_equ(lswitch-1) = 0.
      do l = lswitch-2,1,-1
         dp = (press(l) - press(l+1))*100.
         colch4(l) = colch4(l+1) + (rmch4(l+1) + rmch4(l))*0.5 &
                                 *6.022e22*dp/(zmmean(l)*g)
         ch4_equ(l)=colch4(l)/col(l)
!         tauch4(l)=1.8e-21*colch4(l)
!         if(tauch4(l).ge.1.0) exit
      end do
!      ch4ref = (rmch4(l+1)*(tauch4(l)-1)+rmch4(l)*(1-tauch4(l+1))) &
!                                 /(tauch4(l)-tauch4(l+1))

!     2) interpolation in CH4

!      ch4ref = min(ch4ref,table_ch4(nch4))
!      ch4ref = max(ch4ref,table_ch4(1))

!      indch4 = nch4 - 1
!      do ich4 = 1,nch4
!         if (table_ch4(ich4) >= ch4ref) then
!            indch4 = ich4 - 1
!            indch4 = max(indch4, 1)
!            exit
!         end if
!      end do
!      cich4 = (ch4ref - table_ch4(indch4))     &
!             /(table_ch4(indch4 + 1) - table_ch4(indch4))



!==========================================================================
!     temperature dependence
!==========================================================================

!     1) search for temperature at 1.9 hPa (tref): vertical interpolation

      tref = temp(1)
      do l = (lswitch-1)-1,1,-1
         if (press(l) > 1.9) then
            cinf = (press(l) - 1.9)/(press(l) - press(l+1))
            csup = 1. - cinf
            tref = cinf*temp(l+1) + csup*temp(l)
            exit
         end if
      end do 

!     2) interpolation in temperature

      tref = min(tref,table_temp(1))
      tref = max(tref,table_temp(ntemp))

      if(ntemp.eq.1) then
        citemp = 1.
        indtemp = 1
      else
        do it = 2, ntemp
           if (table_temp(it) <= tref) then
              citemp = (log(tref) - log(table_temp(it)))              &
                      /(log(table_temp(it-1)) - log(table_temp(it)))
              indtemp = it - 1
              exit
           end if
        end do
      endif



!==========================================================================
!     loop over vertical levels
!==========================================================================

      do l = 1,lswitch-1

!     interpolation in air column

         indcol = nz - 1
         do icol = 1,nz
            if (colairtab(icol) < col(l)) then
               indcol = icol - 1
               exit
            end if
         end do
         cicol = (log(col(l)) - log(colairtab(indcol + 1)))              &
                /(log(colairtab(indcol)) - log(colairtab(indcol + 1)))

!     interpolation in ozone column

         indozo = nozo - 1
         do iozo = 1,nozo
            if (table_ozo(iozo) >= ratio_o3(l)) then
               indozo = iozo - 1
               indozo = max(indozo, 1)
               exit
            end if
         end do

      if(nozo.eq.1) then
         ciozo = 0.
      else
         ciozo = (ratio_o3(l) - table_ozo(indozo))             &
                /(table_ozo(indozo + 1) - table_ozo(indozo))
      endif

!     2) interpolation in CH4

      ch4ref = min(ch4_equ(l),table_ch4(nch4))
      ch4ref = max(ch4ref,table_ch4(1))

      indch4 = nch4 - 1
      do ich4 = 1,nch4
         if (table_ch4(ich4) >= ch4ref) then
            indch4 = ich4 - 1
            indch4 = max(indch4, 1)
            exit
         end if
      end do
      if(nch4.eq.1) then
        cich4=0.
        indch4=1
      else
        cich4 = (ch4ref - table_ch4(indch4))     &
               /(table_ch4(indch4 + 1) - table_ch4(indch4))
      endif

!     5-dimensional interpolation weights

!     poids(temp,sza,co2,o3,tau,ch4)

         poids(1,1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)
         poids(1,1,1,2,1,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)*(1.-cich4)       
         poids(1,1,2,1,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)
         poids(1,1,2,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4) 
         poids(1,2,1,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)               
         poids(1,2,1,2,1,1) = citemp*cisza*cicol*ciozo*(1.-citau)*(1.-cich4)                    
         poids(1,2,2,1,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)          
         poids(1,2,2,2,1,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4)               
!         poids(2,1,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)     
!         poids(2,1,1,2,1,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)*(1.-cich4)          
!         poids(2,1,2,1,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)
!         poids(2,1,2,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4)     
!         poids(2,2,1,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)*(1.-cich4)          
!         poids(2,2,1,2,1,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)*(1.-cich4)               
!         poids(2,2,2,1,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*(1.-cich4)     
!         poids(2,2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)*(1.-cich4)          
!!
!         poids(1,1,1,1,2,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau*(1.-cich4)               
!         poids(1,1,1,2,2,1) = citemp*(1.-cisza)*cicol*ciozo*citau*(1.-cich4)                    
!         poids(1,1,2,1,2,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)          
!         poids(1,1,2,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau*(1.-cich4)               
!         poids(1,2,1,1,2,1) = citemp*cisza*cicol*(1.-ciozo)*citau*(1.-cich4)                    
!         poids(1,2,1,2,2,1) = citemp*cisza*cicol*ciozo*citau*(1.-cich4)                         
!         poids(1,2,2,1,2,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)               
!         poids(1,2,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*citau*(1.-cich4)                    
!         poids(2,1,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau*(1.-cich4)          
!         poids(2,1,1,2,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau*(1.-cich4)               
!         poids(2,1,2,1,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)     
!         poids(2,1,2,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau*(1.-cich4)          
!         poids(2,2,1,1,2,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau*(1.-cich4)               
!         poids(2,2,1,2,2,1) = (1.-citemp)*cisza*cicol*ciozo*citau*(1.-cich4)                    
!         poids(2,2,2,1,2,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau*(1.-cich4)          
!         poids(2,2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau*(1.-cich4)               
!      
         poids(1,1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*cich4               
         poids(1,1,1,2,1,2) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)*cich4                    
         poids(1,1,2,1,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4          
         poids(1,1,2,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*cich4               
         poids(1,2,1,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)*cich4                    
         poids(1,2,1,2,1,2) = citemp*cisza*cicol*ciozo*(1.-citau)*cich4                         
         poids(1,2,2,1,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4               
         poids(1,2,2,2,1,2) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)*cich4                    
!         poids(2,1,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)*cich4          
!         poids(2,1,1,2,1,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)*cich4               
!         poids(2,1,2,1,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4     
!         poids(2,1,2,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)*cich4          
!         poids(2,2,1,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)*cich4               
!         poids(2,2,1,2,1,2) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)*cich4                    
!         poids(2,2,2,1,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)*cich4          
!         poids(2,2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)*cich4               
!!
!         poids(1,1,1,1,2,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau*cich4                    
!         poids(1,1,1,2,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau*cich4                         
!         poids(1,1,2,1,2,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*cich4               
!         poids(1,1,2,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau*cich4                    
!         poids(1,2,1,1,2,2) = citemp*cisza*cicol*(1.-ciozo)*citau*cich4                         
!         poids(1,2,1,2,2,2) = citemp*cisza*cicol*ciozo*citau*cich4                              
!         poids(1,2,2,1,2,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau*cich4                    
!         poids(1,2,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau*cich4                         
!         poids(2,1,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau*cich4               
!         poids(2,1,1,2,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau*cich4                    
!         poids(2,1,2,1,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau*cich4          
!         poids(2,1,2,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau*cich4               
!         poids(2,2,1,1,2,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau*cich4                    
!         poids(2,2,1,2,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau*cich4                         
!         poids(2,2,2,1,2,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau*cich4               
!         poids(2,2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau*cich4                    








!     4-dimensional interpolation in the lookup table

         do ij = 1,nd
            j(l,ij) =                                                                         & 
            poids(1,1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,indch4,ij)           &
          + poids(1,1,1,2,1,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,indch4,ij)         &
          + poids(1,1,2,1,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,indch4,ij)         &
          + poids(1,1,2,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4,ij)       &
          + poids(1,2,1,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,indch4,ij)         &
          + poids(1,2,1,2,1,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4,ij)       &
          + poids(1,2,2,1,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4,ij)       &
          + poids(1,2,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)     &
!          + poids(2,1,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,indch4,ij)         &
!          + poids(2,1,1,2,1,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4,ij)       &
!          + poids(2,1,2,1,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4,ij)       &
!          + poids(2,1,2,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4,ij)     &
!          + poids(2,2,1,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4,ij)       &
!          + poids(2,2,1,2,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4,ij)     &
!          + poids(2,2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4,ij)     &
!          + poids(2,2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)   &
!!
!          + poids(1,1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo,indtau+1,indch4,ij)         &
!          + poids(1,1,1,2,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4,ij)       &
!          + poids(1,1,2,1,2,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4,ij)       &
!          + poids(1,1,2,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)     &
!          + poids(1,2,1,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4,ij)       &
!          + poids(1,2,1,2,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)     &
!          + poids(1,2,2,1,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)     &
!          + poids(1,2,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)   &
!          + poids(2,1,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4,ij)       &
!          + poids(2,1,1,2,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4,ij)     &
!          + poids(2,1,2,1,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4,ij)     &
!          + poids(2,1,2,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)   &
!          + poids(2,2,1,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4,ij)     &
!          + poids(2,2,1,2,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)   &
!          + poids(2,2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)   &
!          + poids(2,2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)
! CH4
          + poids(1,1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau,indch4+1,ij)         &
          + poids(1,1,1,2,1,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau,indch4+1,ij)       &
          + poids(1,1,2,1,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau,indch4+1,ij)       &
          + poids(1,1,2,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4+1,ij)     &
          + poids(1,2,1,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau,indch4+1,ij)       &
          + poids(1,2,1,2,1,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4+1,ij)     &
          + poids(1,2,2,1,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4+1,ij)     &
          + poids(1,2,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) 
!          + poids(2,1,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau,indch4+1,ij) &
!          + poids(2,1,1,2,1,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4+1,ij) &
!          + poids(2,1,2,1,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4+1,ij) &
!          + poids(2,1,2,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4+1,ij) &
!          + poids(2,2,1,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4+1,ij) &
!          + poids(2,2,1,2,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4+1,ij) &
!          + poids(2,2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4+1,ij) &
!          + poids(2,2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) &
!!
!          + poids(1,1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,1,1,2,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(1,1,2,1,2,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,1,2,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(1,2,1,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,2,1,2,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(1,2,2,1,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,2,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,1,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,1,1,2,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,1,2,1,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,1,2,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,2,1,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,2,1,2,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij)

         end do

      if (photoheat) then
         do ij = 1,nd
            e(l,ij) =                                                                         & 
            poids(1,1,1,1,1,1)*ephot(indtemp,indsza,indcol,indozo,indtau,indch4,ij)           &
          + poids(1,1,1,2,1,1)*ephot(indtemp,indsza,indcol,indozo+1,indtau,indch4,ij)         &
          + poids(1,1,2,1,1,1)*ephot(indtemp,indsza,indcol+1,indozo,indtau,indch4,ij)         &
          + poids(1,1,2,2,1,1)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4,ij)       &
          + poids(1,2,1,1,1,1)*ephot(indtemp,indsza+1,indcol,indozo,indtau,indch4,ij)         &
          + poids(1,2,1,2,1,1)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4,ij)       &
          + poids(1,2,2,1,1,1)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4,ij)       &
          + poids(1,2,2,2,1,1)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)     &
!          + poids(2,1,1,1,1,1)*ephot(indtemp+1,indsza,indcol,indozo,indtau,indch4,ij)         &
!          + poids(2,1,1,2,1,1)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4,ij)       &
!          + poids(2,1,2,1,1,1)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4,ij)       &
!          + poids(2,1,2,2,1,1)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4,ij)     &
!          + poids(2,2,1,1,1,1)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4,ij)       &
!          + poids(2,2,1,2,1,1)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4,ij)     &
!          + poids(2,2,2,1,1,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4,ij)     &
!          + poids(2,2,2,2,1,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4,ij)   &
!!
!          + poids(1,1,1,1,2,1)*ephot(indtemp,indsza,indcol,indozo,indtau+1,indch4,ij)         &
!          + poids(1,1,1,2,2,1)*ephot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4,ij)       &
!          + poids(1,1,2,1,2,1)*ephot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4,ij)       &
!          + poids(1,1,2,2,2,1)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)     &
!          + poids(1,2,1,1,2,1)*ephot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4,ij)       &
!          + poids(1,2,1,2,2,1)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)     &
!          + poids(1,2,2,1,2,1)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)     &
!          + poids(1,2,2,2,2,1)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)   &
!          + poids(2,1,1,1,2,1)*ephot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4,ij)       &
!          + poids(2,1,1,2,2,1)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4,ij)     &
!          + poids(2,1,2,1,2,1)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4,ij)     &
!          + poids(2,1,2,2,2,1)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4,ij)   &
!          + poids(2,2,1,1,2,1)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4,ij)     &
!          + poids(2,2,1,2,2,1)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4,ij)   &
!          + poids(2,2,2,1,2,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4,ij)   &
!          + poids(2,2,2,2,2,1)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4,ij)
! CH4
          + poids(1,1,1,1,1,2)*ephot(indtemp,indsza,indcol,indozo,indtau,indch4+1,ij)         &
          + poids(1,1,1,2,1,2)*ephot(indtemp,indsza,indcol,indozo+1,indtau,indch4+1,ij)       &
          + poids(1,1,2,1,1,2)*ephot(indtemp,indsza,indcol+1,indozo,indtau,indch4+1,ij)       &
          + poids(1,1,2,2,1,2)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau,indch4+1,ij)     &
          + poids(1,2,1,1,1,2)*ephot(indtemp,indsza+1,indcol,indozo,indtau,indch4+1,ij)       &
          + poids(1,2,1,2,1,2)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau,indch4+1,ij)     &
          + poids(1,2,2,1,1,2)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau,indch4+1,ij)     &
          + poids(1,2,2,2,1,2)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) 
!          + poids(2,1,1,1,1,2)*ephot(indtemp+1,indsza,indcol,indozo,indtau,indch4+1,ij) &
!          + poids(2,1,1,2,1,2)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau,indch4+1,ij) &
!          + poids(2,1,2,1,1,2)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau,indch4+1,ij) &
!          + poids(2,1,2,2,1,2)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau,indch4+1,ij) &
!          + poids(2,2,1,1,1,2)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau,indch4+1,ij) &
!          + poids(2,2,1,2,1,2)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau,indch4+1,ij) &
!          + poids(2,2,2,1,1,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau,indch4+1,ij) &
!          + poids(2,2,2,2,1,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,indch4+1,ij) &
!!
!          + poids(1,1,1,1,2,2)*ephot(indtemp,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,1,1,2,2,2)*ephot(indtemp,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(1,1,2,1,2,2)*ephot(indtemp,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,1,2,2,2,2)*ephot(indtemp,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(1,2,1,1,2,2)*ephot(indtemp,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,2,1,2,2,2)*ephot(indtemp,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(1,2,2,1,2,2)*ephot(indtemp,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(1,2,2,2,2,2)*ephot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,1,1,1,2,2)*ephot(indtemp+1,indsza,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,1,1,2,2,2)*ephot(indtemp+1,indsza,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,1,2,1,2,2)*ephot(indtemp+1,indsza,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,1,2,2,2,2)*ephot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,2,1,1,2,2)*ephot(indtemp+1,indsza+1,indcol,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,2,1,2,2,2)*ephot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,indch4+1,ij) &
!          + poids(2,2,2,1,2,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,indch4+1,ij) &
!          + poids(2,2,2,2,2,2)*ephot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,indch4+1,ij)

         end do
       end if

!     correction for sun distance

         do ij = 1,nd
!            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
            j(l,ij) = j(l,ij)*(1.0/dist_sol)**2.
            if (photoheat) e(l,ij) = e(l,ij)*(1.0/dist_sol)**2.
  
            ! Only during daylight.
            if((ngrid.eq.1))then
                  j(l,ij)= j(l,ij)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
                  if (photoheat) e(l,ij)= e(l,ij)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
            elseif(diurnal .eqv. .false.) then
                  j(l,ij)= j(l,ij) * fractcol
                  if (photoheat) e(l,ij)= e(l,ij) * fractcol
            endif


         end do 



!==========================================================================
!     end of loop over vertical levels
!==========================================================================

      end do

      else

!==========================================================================
!     night
!==========================================================================

         j(:,:) = 0.
         if (photoheat) e(:,:) = 0.

      end if

! fill v_phot array

      v_phot(:,:) = 0.
      e_phot(:,:) = 0.

! Order of photolysis reaction has to be the same in monoreact file and the phototable file
      ij      = 0
      nb_phot = 0
      do while(nb_phot<nb_phot_hv_max)
        ij      = ij + 1
        nb_phot = nb_phot + 1
        do l = 1,lswitch-1
          v_phot(l,nb_phot) = j(l,ij)
          if (photoheat) e_phot(l,nb_phot) = e(l,ij)
        end do
        if (three_prod(ij)) then
          nb_phot = nb_phot + 1
          do l = 1,lswitch-1
            v_phot(l,nb_phot) = j(l,ij)
            if (photoheat) e_phot(l,nb_phot) = e(l,ij)
          end do
        end if
      end do

      return
      end

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