Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F	(revision 134)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F	(revision 334)
@@ -1,5 +1,6 @@
       subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
-     $                   zzlay,zday,pq,pdq,rice,
-     &                   dqchim,dqschim,dqcloud,dqscloud)
+     $                   zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,
+     $                   dqscloud,tauref,co2ice,
+     $                   pu,pdu,pv,pdv,surfdust,surfice)
 c
       implicit none
@@ -19,4 +20,5 @@
 c    update 03/05/2005 cosmetic changes (Franck Lefevre)
 c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
+c    update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre)
 c
 c   Arguments:
@@ -30,9 +32,16 @@
 c    pt(ngridmx,nlayermx)       Temperature (K)
 c    pdt(ngridmx,nlayermx)      Temperature tendency (K)
+c    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
+c    pdu(ngridmx,nlayermx)      u component tendency (K)
+c    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
+c    pdv(ngridmx,nlayermx)      v component tendency (K)
 c    dist_sol                   distance of the sun (AU)
 c    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
 c    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
 c    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
-c    rice(ngridmx,nlayermx)     Estimated ice crystal radius (m)
+c    tauref(ngridmx)            Optical depth at 7 hPa
+c    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
+c    surfdust(ngridmx,nlayermx) dust surface area (micron2/cm3)
+c    surfice(ngridmx,nlayermx)  ice surface area (micron2/cm3)
 c
 c  Output:
@@ -64,12 +73,20 @@
       real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
       real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
+      real    zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
       real    pt(ngridmx,nlayermx)       ! temperature
       real    pdt(ngridmx,nlayermx)      ! temperature tendency
+      real    pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
+      real    pdu(ngridmx,nlayermx)      ! u component tendency
+      real    pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
+      real    pdv(ngridmx,nlayermx)      ! v component tendency
       real    dist_sol                   ! distance of the sun (AU)
-      real    mu0(ngridmx)       ! cos of solar zenith angle (=1 when sun at zenith)
+      real    mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
       real    pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
       real    pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
       real    zday                       ! date (time since Ls=0, in martian days)
-      real    rice(ngridmx,nlayermx)     ! Estimated ice crystal radius (m)
+      real    tauref(ngridmx)            ! optical depth at 7 hPa
+      real    co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
+      real    surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3)
+      real    surfice(ngridmx,nlayermx)  !  ice surface area (micron2/cm3)
 
 c   outputs:
@@ -87,8 +104,15 @@
       integer  ig,l,i,iq
       integer  foundswitch, lswitch
+      integer  ig_vl1
+      real latvl1, lonvl1
       real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
                                      ! new mole fraction after
-      real colden(ngridmx,nqmx)      ! Column densities (cm-2)
       real zt(ngridmx,nlayermx)      ! temperature
+      real zu(ngridmx,nlayermx)      ! u component of the wind
+      real zv(ngridmx,nlayermx)      ! v component of the wind
+      real taucol                    ! optical depth at 7 hPa
+
+      logical depos                  ! switch for dry deposition
+      parameter (depos=.false.)
 c
 c  for each column of atmosphere:
@@ -100,11 +124,11 @@
       real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
       real szacol                 !  Solar zenith angle
+      real surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
+      real surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
       real jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
 c
 c for output:
 c
-      real zdens3d(ngridmx,nlayermx) !  Density  (cm-3)
       real jo3_3d(ngridmx,nlayermx)  !  Photodissociation rate O3->O1D (s-1)
-      real surfice(ngridmx,nlayermx) !  Surface of ice particules (um2/cm3)
 
       logical output              ! to issue calls to writediagfi and stats
@@ -112,8 +136,9 @@
 
       logical,save :: firstcall=.true.
-      integer,save :: nbq  ! number of tracers used in the chemistry
+      integer,save :: nbq       ! number of tracers used in the chemistry
       integer,save :: niq(nqmx) ! array storing the indexes of the tracers
 
 ! index of tracers:
+
       integer,save :: i_co2=0
       integer,save :: i_co=0
@@ -127,13 +152,9 @@
       integer,save :: i_ho2=0
       integer,save :: i_h2o2=0
+      integer,save :: i_ch4=0
       integer,save :: i_n2=0
       integer,save :: i_ar=0
       integer,save :: i_ice=0 ! water ice
       integer,save :: i_h2o=0 ! water vapour
-
-c
-c  scheme A: 1 ; scheme B: 2
-c
-      integer,parameter :: scheme=2
 c
 c=======================================================================
@@ -144,13 +165,6 @@
 c
          if (photochem) then
-            print*,'calchim: INIT CHEMISTRY'
-            if (scheme  .eq.  1) then
-               print*,'calchim: Scheme A : A METTRE A JOUR !!'
-               stop
-c              call init_chimie_A
-            else
-               print*,'calchim: Scheme B'
-               call init_chimie_B
-            end if
+            print*,'calchim: Read photolysis lookup table'
+            call read_phototable
          end if
 
@@ -245,4 +259,12 @@
            niq(nbq)=i_h2o2
          endif
+         i_ch4=igcm_ch4
+         if (i_ch4.eq.0) then
+           write(*,*) "calchim: Error; no CH4 tracer !!!"
+           stop
+         else
+           nbq=nbq+1
+           niq(nbq)=i_ch4
+         endif
          i_n2=igcm_n2
          if (i_n2.eq.0) then
@@ -278,22 +300,21 @@
          endif
 
-         write(*,*) 'calchim: found nbq=',nbq,' tracers'
-         write(*,*) '         i_co2=',i_co2
-         write(*,*) '         i_co=',i_co
-         write(*,*) '         i_o=',i_o
-         write(*,*) '         i_o1d=',i_o1d
-         write(*,*) '         i_o2=',i_o2
-         write(*,*) '         i_o3=',i_o3
-         write(*,*) '         i_h=',i_h
-         write(*,*) '         i_h2=',i_h2
-         write(*,*) '         i_oh=',i_oh
-         write(*,*) '         i_ho2=',i_ho2
-         write(*,*) '         i_h2o2=',i_h2o2
-         write(*,*) '         i_n2=',i_n2
-         write(*,*) '         i_ar=',i_ar
-         write(*,*) '         i_ice=',i_ice
-         write(*,*) '         i_h2o=',i_h2o
-!         write(*,*) '         niq(:)=',niq
-!         write(*,*) '     nqchem_min,nqmx=',nqchem_min,nqmx
+         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
+         write(*,*) '               i_co2  = ',i_co2
+         write(*,*) '               i_co   = ',i_co
+         write(*,*) '               i_o    = ',i_o
+         write(*,*) '               i_o1d  = ',i_o1d
+         write(*,*) '               i_o2   = ',i_o2
+         write(*,*) '               i_o3   = ',i_o3
+         write(*,*) '               i_h    = ',i_h
+         write(*,*) '               i_h2   = ',i_h2
+         write(*,*) '               i_oh   = ',i_oh
+         write(*,*) '               i_ho2  = ',i_ho2
+         write(*,*) '               i_h2o2 = ',i_h2o2
+         write(*,*) '               i_ch4  = ',i_ch4
+         write(*,*) '               i_n2   = ',i_n2
+         write(*,*) '               i_ar   = ',i_ar
+         write(*,*) '               i_ice  = ',i_ice
+         write(*,*) '               i_h2o  = ',i_h2o
          
          firstcall = .false.
@@ -302,8 +323,12 @@
 ! Initialize output tendencies to zero (to handle case of tracers which
 ! are not used in the chemistry (e.g. dust))
-      dqchim(:,:,:)=0
-      dqschim(:,:)=0
-
-
+
+      dqchim(:,:,:) = 0
+      dqschim(:,:)  = 0
+c
+!     latvl1= 22.27
+!     lonvl1= -47.94
+!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
+!    $        int(1.5+(lonvl1+180)*iim/360.)
 c
 c=======================================================================
@@ -317,10 +342,13 @@
          foundswitch = 0
          do l = 1,nlayermx
-            zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
-            do i=1,nbq
-              iq=niq(i) ! get tracer index
+            do i = 1,nbq
+              iq = niq(i) ! get tracer index
               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
-              zycol(l,iq) = zq(ig,l,iq) * mmean(ig,l)/mmol(iq)
-            enddo
+              zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
+            end do
+            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
+            zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
+            zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
+c
             zpress(l) = pplay(ig,l)/100.
             ztemp(l)  = zt(ig,l)
@@ -328,5 +356,10 @@
             zlocal(l) = zzlay(ig,l)/1000.
 c
-c     search for switch index between regions
+c           surfdust1d and surfice1d in cm2/cm3
+c
+            surfdust1d(l) = surfdust(ig,l)*1.e-8
+            surfice1d(l)  = surfice(ig,l)*1.e-8
+c
+c           search for switch index between regions
 c
             if (photochem .and. thermochem) then
@@ -340,24 +373,11 @@
             end if
             if (.not.  thermochem) then
-               lswitch = min(33,nlayermx+1)
+               lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
             end if
 c
-c     ice surface area  in microns^2/cm^3
-c
-c     = 4 pi r^2 * [ zq * mugaz/NA / (rhoice*4/3 pi r^3) ] *zdens
-c     = 3/r * [ zq * mugaz/NA / rhoice ] *zdens
-c     with r in microns, rhoice = 0.92e-12 g microns^-3 and zdens in cm^-3
-c
-            if (water) then
-               zycol(l,i_ice) = (3.e-6/rice(ig,l))*zq(ig,l,i_ice)
-     $                           *(mugaz/6.022e23)*zdens(l)/0.92e-12
-c              write(*,*) "rice=",rice(ig,l)," m / zdens=",zdens(l),
-c    $        " cm-3 / icesurf=",zycol(l,nqmx-1)," microns^2/cm^3"
-               surfice(ig,l) = zycol(l,i_ice)
-            end if
-c
          end do ! of do l=1,nlayermx
 c
          szacol = acos(mu0(ig))*180./pi
+         taucol = tauref(ig)
 c
 c=======================================================================
@@ -365,19 +385,25 @@
 c=======================================================================
 c
-        if (photochem) then
-           if (scheme .eq. 1) then
-              print*,'Scheme A : A METTRE A JOUR !!'
-c             call photochemist_A(zycol,szacol,ptimestep,
-c    $                            zpress,ztemp,zdens,dist_sol)
-           else
-              call photochemist_B(lswitch,zycol,szacol,ptimestep,
-     $                            zpress,ztemp,zdens,dist_sol,jo3)
-           end if
-        end if
-
-        if (thermochem) then
-           call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
-     $                      zlocal,szacol,ptimestep,zday)
-        end if
+c        chemistry in lower atmosphere
+c
+         if (photochem) then
+            call photochemistry(lswitch,zycol,szacol,ptimestep,
+     $                          zpress,ztemp,zdens,dist_sol,
+     $                          surfdust1d,surfice1d,jo3,taucol)
+         end if
+c
+c        chemistry in upper atmosphere
+c
+         if (thermochem) then
+            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
+     $                       zlocal,szacol,ptimestep,zday)
+         end if
+c
+c        dry deposition
+c
+         if (depos) then
+            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev, 
+     $                      zu, zv, zt, zycol, ptimestep, co2ice)
+         end if
 c
 c=======================================================================
@@ -389,5 +415,4 @@
          if (water) then
             do l = 1,nlayermx
-!               dqchim(ig,l,nqmx-1) = 0.
                dqchim(ig,l,i_ice) = 0.
             end do
@@ -425,5 +450,4 @@
      $                                 - zq(ig,l,iq))/ptimestep
                    else if (iq.eq.i_o) then
-!                      i_o = iq
                       dqchim(ig,l,iq) = 0.
                    end if
@@ -440,15 +464,4 @@
           end do ! of do l = 1,nlayermx
 c
-c     dust: This is now taken care of as a first step at beginning of routine
-c
-!          if (nqchem_min .gt. 1) then
-!             do iq = 1,nqchem_min-1
-!                do l = 1,nlayermx
-!                   dqchim(ig,l,iq) = 0.
-!                end do
-!                dqschim(ig,iq) = 0.
-!             end do
-!          end if
-c
 c     condensation of h2o2
 c
@@ -456,44 +469,7 @@
      $                 ztemp,zycol,dqcloud,dqscloud)
 c
-c     for outputs
-c
-          do i=1,nbq
-            iq=niq(i) ! get tracer index
-             colden(ig,iq) = 0.
-             do l = 1,nlayermx
-c
-c     column density converted in cm-2
-c     pplev en pa, mugaz en g.mol-1 et g en m.s-2
-c     not for ice
-c
-                if (iq.ne.i_h2o2) then
-                  colden(ig,iq) = colden(ig,iq) + zycol(l,iq)
-     $                         *6.022e22*(pplev(ig,l)-pplev(ig,l+1))
-     $                         /(mmean(ig,l)*g)
-                else   ! for H2O2, remove condensation from zycol
-                  colden(ig,iq) = colden(ig,iq) + (zycol(l,iq) +
-     $               dqcloud(ig,l,iq)*ptimestep*mmean(ig,l)/mmol(iq))
-     $                         *6.022e22*(pplev(ig,l)-pplev(ig,l+1))
-     $                         /(mmean(ig,l)*g)
-                end if
-c
-c     local densities, for outputs (put in zq)
-c     not for ice
-c
-                zq(ig,l,iq) = zycol(l,iq)*zdens(l)
-c                        for H2O2, remove condensation from zycol
-                if (iq.eq.i_h2o2) then
-                   zq(ig,l,iq) = zdens(l)*(zycol(l,iq) +
-     $               dqcloud(ig,l,iq)*ptimestep*mmean(ig,l)/mmol(iq))
-                end if
-             end do
-          end do
-c
 c     density and j(o3->o1d), for outputs
 c
-          zdens3d(ig,1) = zdens(1)
-          jo3_3d(ig,1) = jo3(1)
-          do l = 2,nlayermx
-             zdens3d(ig,l) = zdens(l)
+          do l = 1,nlayermx
              jo3_3d(ig,l) = jo3(l)
           end do
@@ -513,50 +489,10 @@
 
       if (output) then
-
          if (ngridmx .gt. 1) then
-c           call writediagfi(ngridmx,'dens','atm dens.','cm-3',3,zdens3d(1,1))
-c           call writediagfi(ngridmx,'jo3','j o3->o1d','s-1',3,jo3_3d(1,1))
-c           call writediagfi(ngridmx,'sice','ice surf.','um2/cm3',3,surfice(1,1))
-            do i=1,nbq
-              iq=niq(i) ! get tracer index
-               if (iq.ne.i_ice) then
-                 write(str20(1:20),'(a20)') noms(iq)
-                 call writediagfi(ngridmx,'n_'//trim(str20),'density',
-     $                             'cm-3',3,zq(1,1,iq))
-c                 call writediagfi(ngridmx,'dqch_'//str5,'density','cm-3',3,dqchim(1,1,iq))
-c                 if (noms(iq) .eq. "h2o2" .or. noms(iq) .eq. "h2o") then
-c                    call writediagfi(ngridmx,'cl_'//str5,'density','cm-3',3,dqcloud(1,1,iq))
-c                 end if
-                 call writediagfi(ngridmx,'c_'//trim(str20),
-     $                            'col. dens.','cm-2',2,colden(1,iq))
-               end if
-            end do
-c
-            if (callstats) then
-c
-c              convert to mole.cm-2 for the column densities
-c
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                  do ig = 1,ngridmx
-                     colden(ig,iq) = colden(ig,iq)/6.022e23
-                  end do   
-               end do 
-c
-c              call wstats(ngridmx,"jo3","jo3->o1d","s-1",3,jo3_3d)
-c
-               do i=1,nbq
-                 iq=niq(i) ! get tracer index
-                  if (iq.ne.i_ice) then
-                     write(str20(1:20),'(a20)') noms(iq)
-                     call wstats(ngridmx,"n_"//trim(str20),"density",
-     &                           "cm-3",3,zq(1,1,iq))
-                     call wstats(ngridmx,"c_"//trim(str20),"col. dens.",
-     &                           "mol cm-2",2,colden(1,iq))
-                  end if
-               end do ! of i=1,nbq
-            end if ! of if (callstats)
+            call writediagfi(ngridmx,'jo3','j o3->o1d',
+     $                       's-1',3,jo3_3d(1,1))
+            call wstats(ngridmx,'jo3','j o3->o1d',
+     $                       's-1',3,jo3_3d(1,1))
          end if ! of if (ngridmx.gt.1)
-c
       endif ! of if (output)
 c
