Index: /trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h	(revision 1429)
+++ /trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h	(revision 1430)
@@ -1,22 +1,25 @@
-! -----------------------------------------
-! ---  Data for photolysis lookup table ---
-! -----------------------------------------
+!--------------------------------------------
+!     data for photochemistry
+!--------------------------------------------
 
-      integer nd, nozo, nr, nsza, ntemp, ntau
+!--------------------------------------------
+!     dimensions of photolysis lookup table
+!--------------------------------------------
 
-      parameter (nd    = 29)
-      parameter (nozo  = 7)
-      parameter (nr    = nd + 28)
-      parameter (nsza  = 27)
-      parameter (ntemp = 4)
-      parameter (ntau  = 8)
+      integer, parameter :: nd    = 13  ! species
+      integer, parameter :: nz    = 143 ! altitude
+      integer, parameter :: nozo  = 7   ! ozone
+      integer, parameter :: nsza  = 27  ! solar zenith angle
+      integer, parameter :: ntemp = 4   ! temperature
+      integer, parameter :: ntau  = 8   ! dust
 
-      real kb
-      parameter (kb = 1.3806e-23)
+!--------------------------------------------
+
+      real, parameter :: kb = 1.3806e-23
 
       common/chimiedata/jphot,colairtab,table_ozo
 
-      real jphot(ntemp,nsza,0:200,nozo,ntau,nd)
-      real colairtab(0:200)
+      real jphot(ntemp,nsza,nz,nozo,ntau,nd)
+      real colairtab(nz)
       real szatab(nsza)
       real table_ozo(nozo)
@@ -30,2 +33,10 @@
 
       data tautab/0., 0.2, 0.4, 0.6, 0.8, 1., 2., 4./
+
+!--------------------------------------------
+!     number of reactions in ASIS solver
+!--------------------------------------------
+
+      integer, parameter :: nb_phot_max       = 18
+      integer, parameter :: nb_reaction_3_max = 6
+      integer, parameter :: nb_reaction_4_max = 30
Index: /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 1429)
+++ /trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F	(revision 1430)
@@ -133,6 +133,6 @@
       end do
 
-      call phot(nlayer, lswitch, press, temp, sza, tau, dist_sol, 
-     $          rmco2, rmo3, j)
+      call photolysis(nlayer, lswitch, press, temp, sza, tau, dist_sol,
+     $                rmco2, rmo3, j)
 
       j_o3_o1d = 5
@@ -245,5 +245,5 @@
       integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d,
      $           j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2,
-     $           j_ho2, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
+     $           j_ho2, j_no, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
      $           j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h,
      $           j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl,
@@ -295,26 +295,10 @@
       j_o3_o         =  6      ! o3 + hv     -> o2 + o
       j_h2o          =  7      ! h2o + hv    -> h + oh
-      j_hdo          =  8      ! hdo + hv    -> d + oh
-      j_h2o2         =  9      ! h2o2 + hv   -> oh + oh
-      j_ho2          =  10     ! ho2 + hv    -> oh + o
+      j_h2o2         =  8      ! h2o2 + hv   -> oh + oh
+      j_ho2          =  9      ! ho2 + hv    -> oh + o
+      j_no           =  10     ! no + hv     -> n + o
       j_no2          =  11     ! no2 + hv    -> no + o
-      j_ch4_ch3_h    =  12     ! ch4 + hv    -> ch3 + h
-      j_ch4_1ch2_h2  =  13     ! ch4 + hv    -> 1ch2 + h2
-      j_ch4_3ch2_h_h =  14     ! ch4 + hv    -> 3ch2 + h + h
-      j_ch4_ch_h2_h  =  15     ! ch4 + hv    -> ch + h2 + h
-      j_ch3o2h       =  16     ! ch3o2h + hv -> ch3o + oh
-      j_ch2o_hco     =  17     ! ch2o + hv   -> h + hco
-      j_ch2o_co      =  18     ! ch2o + hv   -> h2 + co
-      j_ch3oh        =  19     ! ch3oh + hv  -> ch3o + h
-      j_c2h6         =  20     ! c2h6 + hv   -> products
-      j_hcl          =  21     ! hcl + hv    -> h + cl
-      j_hocl         =  22     ! hocl + hv   -> oh + cl
-      j_clo          =  23     ! clo + hv    -> cl + o
-      j_so2          =  24     ! so2 + hv    -> so + o
-      j_so           =  25     ! so + hv     -> s + o
-      j_h2s          =  26     ! h2s + hv    -> hs + s
-      j_so3          =  27     ! so2 + hv    -> so2 + o
-      j_hno3         =  28     ! hno3 + hv   -> oh + no2
-      j_hno4         =  29     ! hno4 + hv   -> ho2 + no2
+      j_hno3         =  12     ! hno3 + hv   -> oh + no2
+      j_hno4         =  13     ! hno4 + hv   -> ho2 + no2
 
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -776,373 +760,4 @@
 c*****************************************************************
 
-      subroutine phot(nlayer,
-     $                lswitch, press, temp, sza, tauref, dist_sol, 
-     $                rmco2, rmo3, j)
-
-c*****************************************************************
-
-      USE comcstfi_h
-      implicit none
-
-#include "chimiedata.h" 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     inputs:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-       
-      integer, intent(in) :: nlayer ! number of atmospheric layers
-      integer :: lswitch            ! interface level between chemistries
-      real :: press(nlayer)         ! pressure (hPa)
-      real :: temp(nlayer)          ! temperature (K)
-      real :: sza                   ! solar zenith angle (deg)
-      real :: tauref                ! optical depth at 7 hpa
-      real :: dist_sol              ! sun distance (AU) 
-      real :: rmco2(nlayer)         ! co2 mixing ratio
-      real :: rmo3(nlayer)          ! ozone mixing ratio
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     output:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1) 
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp, 
-     $           iozo, isza, itau, it, l
-
-      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
-      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
-      real :: poids(2,2,2,2,2)            ! 5D 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 :: cinf, csup, cicol, 
-     $        ciozo, cisza, citemp, citau
-      real :: colo3min, dp, coef
-      real :: ratio_o3(nlayer)
-      real :: tau
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     day/night criterion
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      if (sza .le. 95.) then 
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     temperatures at 1.9 hPa in lookup table
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      
-      table_temp(1) = 226.2
-      table_temp(2) = 206.2
-      table_temp(3) = 186.2
-      table_temp(4) = 169.8
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     interpolation in solar zenith angle
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      indsza = nsza - 1
-      do isza = 1,nsza
-         if (szatab(isza) .ge. sza) then
-            indsza = min(indsza,isza - 1)
-            indsza = max(indsza, 1)
-         end if
-      end do
-      cisza = (sza - szatab(indsza))
-     $       /(szatab(indsza + 1) - szatab(indsza))
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     interpolation in dust (tau)
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      tau = min(tauref, tautab(ntau))
-      tau = max(tau, tautab(1))
-
-      indtau = ntau - 1
-      do itau = 1,ntau
-         if (tautab(itau) .ge. tau) then
-            indtau = min(indtau,itau - 1)
-            indtau = max(indtau, 1)
-         end if
-      end do
-      citau = (tau - tautab(indtau))
-     $       /(tautab(indtau + 1) - tautab(indtau))
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     co2 and ozone columns
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-c     co2 column at model top (molecule.cm-2)
-
-      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.
-     $                 /(mugaz*g)
-
-c     ozone column at model top
-
-      colo3(lswitch-1) = 0.
-
-c     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/(mugaz*g)
-         col(l) = min(col(l), colairtab(0))
-         colo3(l) = colo3(l+1)
-     $            + (rmo3(l+1) + rmo3(l))*0.5
-     $              *6.022e22*dp/(mugaz*g)
-      end do
-
-c     ratio ozone column/minimal theoretical column (0.1 micron-atm) 
-c
-c     ro3 = 7.171e-10 is the o3 mixing ratio for a uniform
-c     profile giving a column 0.1 micron-atmosphere at
-c     a surface pressure of 10 hpa.
-
-      do l = 1,lswitch-1
-         colo3min    = col(l)*7.171e-10
-         ratio_o3(l) = colo3(l)/colo3min
-         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.)
-         ratio_o3(l) = max(ratio_o3(l), 1.)
-      end do
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     temperature dependence
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-c     1) search for temperature at 1.9 hPa (tref): vertical interpolation
-
-      tref = temp(1)
-      do l = (lswitch-1)-1,1,-1
-         if (press(l) .gt. 1.9) then
-            cinf = (press(l) - 1.9)
-     $            /(press(l) - press(l+1))
-            csup = 1. - cinf
-            tref = cinf*temp(l+1) + csup*temp(l)
-            go to 10
-         end if
-      end do 
- 10   continue
-
-c     2) interpolation in temperature
-
-      tref = min(tref,table_temp(1))
-      tref = max(tref,table_temp(ntemp))
-
-      do it = 2, ntemp
-         if (table_temp(it) .le. tref) then
-            citemp = (log(tref) - log(table_temp(it)))
-     $              /(log(table_temp(it-1)) - log(table_temp(it)))
-            indtemp = it - 1
-            goto 20
-         end if
-      end do
-  20  continue
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     loop over vertical levels
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do l = 1,lswitch-1
-
-c     interpolation in air column
-
-         do icol = 0,200
-            if (colairtab(icol) .lt. col(l)) then
-               cicol = (log(col(l)) - log(colairtab(icol)))
-     $                /(log(colairtab(icol-1)) - log(colairtab(icol)))
-               indcol = icol - 1
-               goto 30
-            end if
-         end do
- 30      continue
-
-cc    interpolation in ozone column
-
-         indozo = nozo - 1
-         do iozo = 1,nozo
-            if (table_ozo(iozo)*10. .ge. ratio_o3(l)) then
-               indozo = min(indozo, iozo - 1)
-               indozo = max(indozo, 1)
-            end if
-         end do
-         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)
-     $          /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)
-
-cc    4-dimensional interpolation weights
-
-c     poids(temp,sza,co2,o3,tau)
-
-         poids(1,1,1,1,1) = citemp
-     $                     *(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
-         poids(1,1,1,2,1) = citemp
-     $                     *(1.-cisza)*cicol*ciozo*(1.-citau)
-         poids(1,1,2,1,1) = citemp
-     $                     *(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(1,1,2,2,1) = citemp
-     $                     *(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
-         poids(1,2,1,1,1) = citemp
-     $                     *cisza*cicol*(1.-ciozo)*(1.-citau)
-         poids(1,2,1,2,1) = citemp
-     $                     *cisza*cicol*ciozo*(1.-citau)
-         poids(1,2,2,1,1) = citemp
-     $                     *cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(1,2,2,2,1) = citemp
-     $                     *cisza*(1.-cicol)*ciozo*(1.-citau)
-         poids(2,1,1,1,1) = (1.-citemp)
-     $                     *(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
-         poids(2,1,1,2,1) = (1.-citemp)
-     $                     *(1.-cisza)*cicol*ciozo*(1.-citau)
-         poids(2,1,2,1,1) = (1.-citemp)
-     $                     *(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(2,1,2,2,1) = (1.-citemp)
-     $                     *(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
-         poids(2,2,1,1,1) = (1.-citemp)
-     $                     *cisza*cicol*(1.-ciozo)*(1.-citau)
-         poids(2,2,1,2,1) = (1.-citemp)
-     $                     *cisza*cicol*ciozo*(1.-citau)
-         poids(2,2,2,1,1) = (1.-citemp)
-     $                     *cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
-         poids(2,2,2,2,1) = (1.-citemp)
-     $                     *cisza*(1.-cicol)*ciozo*(1.-citau)
-
-         poids(1,1,1,1,2) = citemp
-     $                     *(1.-cisza)*cicol*(1.-ciozo)*citau
-         poids(1,1,1,2,2) = citemp
-     $                     *(1.-cisza)*cicol*ciozo*citau
-         poids(1,1,2,1,2) = citemp
-     $                     *(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
-         poids(1,1,2,2,2) = citemp
-     $                     *(1.-cisza)*(1.-cicol)*ciozo*citau
-         poids(1,2,1,1,2) = citemp
-     $                     *cisza*cicol*(1.-ciozo)*citau
-         poids(1,2,1,2,2) = citemp
-     $                     *cisza*cicol*ciozo*citau
-         poids(1,2,2,1,2) = citemp
-     $                     *cisza*(1.-cicol)*(1.-ciozo)*citau
-         poids(1,2,2,2,2) = citemp
-     $                     *cisza*(1.-cicol)*ciozo*citau
-         poids(2,1,1,1,2) = (1.-citemp)
-     $                     *(1.-cisza)*cicol*(1.-ciozo)*citau
-         poids(2,1,1,2,2) = (1.-citemp)
-     $                     *(1.-cisza)*cicol*ciozo*citau
-         poids(2,1,2,1,2) = (1.-citemp)
-     $                     *(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
-         poids(2,1,2,2,2) = (1.-citemp)
-     $                     *(1.-cisza)*(1.-cicol)*ciozo*citau
-         poids(2,2,1,1,2) = (1.-citemp)
-     $                     *cisza*cicol*(1.-ciozo)*citau
-         poids(2,2,1,2,2) = (1.-citemp)
-     $                     *cisza*cicol*ciozo*citau
-         poids(2,2,2,1,2) = (1.-citemp)
-     $                     *cisza*(1.-cicol)*(1.-ciozo)*citau
-         poids(2,2,2,2,2) = (1.-citemp)
-     $                     *cisza*(1.-cicol)*ciozo*citau
-
-cc    4-dimensional interpolation in the lookup table
-
-         do ij = 1,nd
-           j(l,ij) =
-     $     poids(1,1,1,1,1)
-     $     *jphot(indtemp,indsza,indcol,indozo,indtau,ij)
-     $   + poids(1,1,1,2,1)
-     $     *jphot(indtemp,indsza,indcol,indozo+1,indtau,ij)
-     $   + poids(1,1,2,1,1)
-     $     *jphot(indtemp,indsza,indcol+1,indozo,indtau,ij)
-     $   + poids(1,1,2,2,1)
-     $     *jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij)
-     $   + poids(1,2,1,1,1)
-     $     *jphot(indtemp,indsza+1,indcol,indozo,indtau,ij)
-     $   + poids(1,2,1,2,1)
-     $     *jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij)
-     $   + poids(1,2,2,1,1)
-     $     *jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij)
-     $   + poids(1,2,2,2,1)
-     $     *jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij)
-     $   + poids(2,1,1,1,1)
-     $     *jphot(indtemp+1,indsza,indcol,indozo,indtau,ij)
-     $   + poids(2,1,1,2,1)
-     $     *jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij)
-     $   + poids(2,1,2,1,1)
-     $     *jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij)
-     $   + poids(2,1,2,2,1)
-     $     *jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij)
-     $   + poids(2,2,1,1,1)
-     $     *jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij)
-     $   + poids(2,2,1,2,1)
-     $     *jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij)
-     $   + poids(2,2,2,1,1)
-     $     *jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij)
-     $   + poids(2,2,2,2,1)
-     $     *jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)
-
-     $   + poids(1,1,1,1,2)
-     $     *jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)
-     $   + poids(1,1,1,2,2)
-     $     *jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij)
-     $   + poids(1,1,2,1,2)
-     $     *jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij)
-     $   + poids(1,1,2,2,2)
-     $     *jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij)
-     $   + poids(1,2,1,1,2)
-     $     *jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij)
-     $   + poids(1,2,1,2,2)
-     $     *jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij)
-     $   + poids(1,2,2,1,2)
-     $     *jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij)
-     $   + poids(1,2,2,2,2)
-     $     *jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij)
-     $   + poids(2,1,1,1,2)
-     $     *jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij)
-     $   + poids(2,1,1,2,2)
-     $     *jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij)
-     $   + poids(2,1,2,1,2)
-     $     *jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij)
-     $   + poids(2,1,2,2,2)
-     $     *jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij)
-     $   + poids(2,2,1,1,2)
-     $     *jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij)
-     $   + poids(2,2,1,2,2)
-     $     *jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij)
-     $   + poids(2,2,2,1,2)
-     $     *jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij)
-     $   + poids(2,2,2,2,2)
-     $     *jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)
-         end do
-
-cc    correction for sun distance
-
-         do ij = 1,nd
-            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
-         end do 
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     end of loop over vertical levels
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      end do
-
-      else
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     night
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      do ij = 1,nd
-         do l = 1,lswitch-1
-            j(l,ij) = 0.
-         end do
-      end do
-
-      end if
-
-      return
-      end
-
-c*****************************************************************
-
       subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
 
Index: /trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90	(revision 1430)
+++ /trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90	(revision 1430)
@@ -0,0 +1,298 @@
+!==========================================================================
+
+      subroutine photolysis(nlayer,                              &
+                            lswitch, press, temp, sza, tauref,   &
+                            dist_sol, rmco2, rmo3, j)
+
+!==========================================================================
+
+      use comcstfi_h
+
+      implicit none
+
+#include "chimiedata.h" 
+
+!==========================================================================
+!     input:
+!==========================================================================
+       
+      integer, intent(in) :: nlayer ! number of atmospheric layers
+      integer :: lswitch            ! interface level between chemistries
+      real :: press(nlayer)         ! pressure (hPa)
+      real :: temp(nlayer)          ! temperature (K)
+      real :: sza                   ! solar zenith angle (deg)
+      real :: tauref                ! optical depth at 7 hpa
+      real :: dist_sol              ! sun distance (AU) 
+      real :: rmco2(nlayer)         ! co2 mixing ratio
+      real :: rmo3(nlayer)          ! ozone mixing ratio
+
+!==========================================================================
+!     output: interpolated photodissociation rates (s-1)
+!==========================================================================
+
+      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1) 
+
+!==========================================================================
+!     local:
+!==========================================================================
+
+      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,     &
+                 iozo, isza, itau, it, l
+
+      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
+      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
+      real :: poids(2,2,2,2,2)            ! 5D 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 :: cinf, csup, cicol, ciozo, cisza, citemp, citau
+      real :: colo3min, dp, coef
+      real :: ratio_o3(nlayer)
+      real :: tau
+
+!==========================================================================
+!     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
+
+!==========================================================================
+!     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
+      cisza = (sza - szatab(indsza))  &
+             /(szatab(indsza + 1) - szatab(indsza))
+
+!==========================================================================
+!     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
+      citau = (tau - tautab(indtau))     &
+             /(tautab(indtau + 1) - tautab(indtau))
+
+!==========================================================================
+!     co2 and ozone columns
+!==========================================================================
+
+!     co2 column at model top (molecule.cm-2)
+
+      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
+                       /(mugaz*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/(mugaz*g)
+         col(l) = min(col(l), colairtab(1))
+         colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 &
+                                 *6.022e22*dp/(mugaz*g)
+      end do
+
+!     ratio ozone column/minimal theoretical column (0.1 micron-atm) 
+
+!     ro3 = 7.171e-10 is the o3 mixing ratio for a uniform
+!     profile giving a column 0.1 micron-atmosphere at
+!     a surface pressure of 10 hpa.
+
+      do l = 1,lswitch-1
+         colo3min    = col(l)*7.171e-10
+         ratio_o3(l) = colo3(l)/colo3min
+         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.)
+         ratio_o3(l) = max(ratio_o3(l), 1.)
+      end do
+
+!==========================================================================
+!     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))
+
+      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
+
+!==========================================================================
+!     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)*10. >= ratio_o3(l)) then
+               indozo = iozo - 1
+               indozo = max(indozo, 1)
+               exit
+            end if
+         end do
+         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)             &
+                /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)
+
+!     4-dimensional interpolation weights
+
+!     poids(temp,sza,co2,o3,tau)
+
+         poids(1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
+         poids(1,1,1,2,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)
+         poids(1,1,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
+         poids(1,1,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
+         poids(1,2,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)
+         poids(1,2,1,2,1) = citemp*cisza*cicol*ciozo*(1.-citau)
+         poids(1,2,2,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
+         poids(1,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)
+         poids(2,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
+         poids(2,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)
+         poids(2,1,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
+         poids(2,1,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
+         poids(2,2,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)
+         poids(2,2,1,2,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)
+         poids(2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
+         poids(2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)
+!
+         poids(1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau
+         poids(1,1,1,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau
+         poids(1,1,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
+         poids(1,1,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau
+         poids(1,2,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*citau
+         poids(1,2,1,2,2) = citemp*cisza*cicol*ciozo*citau
+         poids(1,2,2,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau
+         poids(1,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau
+         poids(2,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau
+         poids(2,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau
+         poids(2,1,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
+         poids(2,1,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau
+         poids(2,2,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau
+         poids(2,2,1,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau
+         poids(2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau
+         poids(2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau
+
+!     4-dimensional interpolation in the lookup table
+
+         do ij = 1,nd
+            j(l,ij) =                                                                & 
+            poids(1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,ij)           &
+          + poids(1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,ij)         &
+          + poids(1,1,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,ij)         &
+          + poids(1,1,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij)       &
+          + poids(1,2,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,ij)         &
+          + poids(1,2,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij)       &
+          + poids(1,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij)       &
+          + poids(1,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij)     &
+          + poids(2,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,ij)         &
+          + poids(2,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij)       &
+          + poids(2,1,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij)       &
+          + poids(2,1,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij)     &
+          + poids(2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij)       &
+          + poids(2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij)     &
+          + poids(2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij)     &
+          + poids(2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)   &
+!
+          + poids(1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)         &
+          + poids(1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij)       &
+          + poids(1,1,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij)       &
+          + poids(1,1,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij)     &
+          + poids(1,2,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij)       &
+          + poids(1,2,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij)     &
+          + poids(1,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij)     &
+          + poids(1,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij)   &
+          + poids(2,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij)       &
+          + poids(2,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij)     &
+          + poids(2,1,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij)     &
+          + poids(2,1,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij)   &
+          + poids(2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij)     &
+          + poids(2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij)   &
+          + poids(2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij)   &
+          + poids(2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)
+         end do
+
+!     correction for sun distance
+
+         do ij = 1,nd
+            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
+         end do 
+
+!==========================================================================
+!     end of loop over vertical levels
+!==========================================================================
+
+      end do
+
+      else
+
+!==========================================================================
+!     night
+!==========================================================================
+
+         j(:,:) = 0.
+
+      end if
+
+      return
+      end
Index: unk/LMDZ.MARS/libf/aeronomars/read_phototable.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/read_phototable.F	(revision 1429)
+++ 	(revision )
@@ -1,93 +1,0 @@
-      subroutine read_phototable
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c   subject:
-c   --------
-c
-c   read photolysis lookup table
-c
-c   VERSION: 17/03/2011
-c
-c   Author:   Franck Lefevre
-c
-c   Arguments:
-c   ----------
-c
-c   The output variable is jphot and is put in common chimiedata.
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-! to use  'getin'
-      USE ioipsl_getincom 
-
-      implicit none
-c
-#include "dimensions.h"
-#include "chimiedata.h"
-#include "datafile.h"
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     local:
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-      integer fic, ij, iozo, isza, itemp, iz, itau, ierr
-c
-      real xsza
-      character(len=128) :: phototable ! photolysis table file name
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-! set photolysis table input file name
-      phototable="jmars.20111014" ! default
-      ! look for a " phototable= ..." option in def files
-      call getin("phototable",phototable)
-
-      fic = 81
-
-      open(fic, form = 'formatted', status = 'old',
-     &file =trim(datafile)//"/"//trim(phototable),iostat=ierr)
-
-      IF (ierr.NE.0) THEN
-        write(*,*)'Error : cannot open photolysis lookup table ',
-     &             trim(phototable)
-        write(*,*)'(in aeronomars/init_chemistry.F)'
-        write(*,*)'It should be in :',trim(datafile),'/'
-        write(*,*)'1) You can change this directory path in physiq.def'
-        write(*,*)'   with:'
-        write(*,*)'   datadir=/path/to/the/directory'
-        write(*,*)'2) You can change the input phototable file name in'
-        write(*,*)'   physiq.def with:'
-        write(*,*)'   phototable=filename'
-        STOP
-      ENDIF
-
-c
-      print*, 'read photolysis lookup table ',trim(phototable)
-c
-      do itau = 1,ntau
-         do itemp = 1,ntemp
-            do iozo = 1,nozo
-               do isza = 1,nsza
-                  do iz = 200,0,-1
-                     read(fic,*) colairtab(iz), xsza, table_ozo(iozo)
-c                    print*, colairtab(iz), xsza, table_ozo(iozo)
-                     read(fic,'(7e11.4)')
-     $               (jphot(itemp,isza,iz,iozo,itau,ij), ij= 1,nd)
-                     do ij = 1,nd
-                        if (jphot(itemp,isza,iz,iozo,itau,ij)
-     $                     .eq. 1.e-30) then
-                           jphot(itemp,isza,iz,iozo,itau,ij) = 0.
-                        end if
-                     end do
-                  end do
-               end do
-            end do
-         end do
-      end do
-c
-      print*, 'lookup table...ok'
-      close(fic)
-c
-      return
-      end
Index: /trunk/LMDZ.MARS/libf/aeronomars/read_phototable.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/read_phototable.F90	(revision 1430)
+++ /trunk/LMDZ.MARS/libf/aeronomars/read_phototable.F90	(revision 1430)
@@ -0,0 +1,94 @@
+!***********************************************************************
+
+      subroutine read_phototable
+
+!***********************************************************************
+!
+!   subject:
+!   --------
+!
+!   read photolysis lookup table
+!
+!   VERSION: 8/10/2014
+!
+!   Author:   Franck Lefevre
+!
+!   Arguments:
+!   ----------
+!
+!   The output variable is jphot and is put in common chimiedata.
+!
+!***********************************************************************
+
+      use ioipsl_getincom 
+
+      implicit none
+
+#include "chimiedata.h"
+#include "datafile.h"
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!     local:
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      integer :: fic, ij, iozo, isza, itemp, iz, itau, ierr
+      real    :: xsza
+
+      character(len = 128) :: phototable ! photolysis table file name
+
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! set photolysis table input file name
+!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      phototable = "jmars.20140930" ! default
+
+! look for a " phototable= ..." option in def files
+
+      call getin("phototable",phototable)
+
+      fic = 81
+
+      open(fic, form = 'formatted', status = 'old',                &
+           file =trim(datafile)//"/"//trim(phototable),iostat=ierr)
+
+      if (ierr /= 0) THEN
+        write(*,*)'Error : cannot open photolysis lookup table ', trim(phototable)
+        write(*,*)'It should be in :',trim(datafile),'/'
+        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)
+
+      do itau = 1,ntau
+         do itemp = 1,ntemp
+            do iozo = 1,nozo
+               do isza = 1,nsza
+                  do iz = nz,1,-1
+                     read(fic,*) colairtab(iz), xsza, table_ozo(iozo)
+                     read(fic,'(7e11.4)') (jphot(itemp,isza,iz,iozo,itau,ij), ij= 1,nd)
+                     do ij = 1,nd
+                        if (jphot(itemp,isza,iz,iozo,itau,ij) == 1.e-30) then
+                           jphot(itemp,isza,iz,iozo,itau,ij) = 0.
+                        end if
+                     end do
+                  end do
+               end do
+            end do
+         end do
+      end do
+
+      print*, 'lookup table...ok'
+      close(fic)
+
+      return
+      end
