Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 458)
+++ /trunk/LMDZ.MARS/README	(revision 459)
@@ -1302,2 +1302,17 @@
 
 >> tiny change : nuice_sed initialisation is now done in inifis. Also changed initracer and improvedclouds.
+
+
+== 08/12/2011 == EM
+>> Include changes and updates for photochemistry by FL:
+- aeronomars/calchim.F : change in units of surface density.
+- aeronomars/surfacearea.F : new routine to compute ice and dust surface area
+                             (m2/m3) available for heterogeneous reactions.
+- phymars/initracer.F : bug correction: initialize igcm_ch4 and change loop
+                        bounds (when guessing tracer names/properties with
+                        old input files).
+- phymars/watercloud.F : cleanup.
+- phymars/physiq.F : add call to surfacearea; photochemistry is now called
+                     after sedimentation to take into acount updated rdust
+                     and rice.
+ 
Index: /trunk/LMDZ.MARS/libf/aeronomars/calchim.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/calchim.F	(revision 458)
+++ /trunk/LMDZ.MARS/libf/aeronomars/calchim.F	(revision 459)
@@ -21,4 +21,5 @@
 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    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
 c
 c   Arguments:
@@ -42,6 +43,6 @@
 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    surfdust(ngridmx,nlayermx) dust surface area (m2/m3)
+c    surfice(ngridmx,nlayermx)  ice surface area (m2/m3)
 c
 c  Output:
@@ -87,6 +88,6 @@
       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)
+      real    surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3)
+      real    surfice(ngridmx,nlayermx)  !  ice surface area (m2/m3)
 
 c   outputs:
@@ -356,8 +357,8 @@
             zlocal(l) = zzlay(ig,l)/1000.
 c
-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           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
+c
+            surfdust1d(l) = surfdust(ig,l)*1.e-2
+            surfice1d(l)  = surfice(ig,l)*1.e-2
 c
 c           search for switch index between regions
Index: /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 459)
+++ /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 459)
@@ -0,0 +1,98 @@
+      subroutine surfacearea(ngrid, nlay, pplay, pzlay, pt, pq, nq,
+     $                       rdust, rice, tau, tauscaling,
+     $                       surfdust, surfice)
+
+      implicit none
+
+!==========================================================================
+!     calculation of the ice and dust surface area (m2/m3)
+!     available for heterogeneous reactions
+!
+!     Franck Lefevre
+!     version 1.0 november 2011
+!==========================================================================
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "callkeys.h"
+#include "tracer.h"
+#include "dimradmars.h"
+#include "chimiedata.h"
+#include "conc.h"
+
+! input
+
+      integer ngrid, nlay
+      integer nq                   ! number of tracers
+      real pplay(ngrid,nlay)       ! pressure at mid-layers (Pa)
+      real pzlay(ngrid,nlay)       ! altitude at mid-layers (m)
+      real pt(ngrid,nlay)          ! temperature at mid-layers (K)
+      real pq(ngrid,nlay,nq)       ! tracers (kg/kg)
+      real rdust(ngrid,nlay)       ! dust geometric mean radius (m)
+      real rice(ngrid,nlay)        ! ice mass mean radius (m)
+      real tau(ngrid,naerkind)     ! column dust optical depth at each point
+      real tauscaling(ngrid)       ! conversion factor for dust amount
+
+! output
+
+      real surfdust(ngrid,nlay)    ! dust surface area (m2/m3)
+      real surfice(ngrid,nlay)     ! water-ice surface area (m2/m3)
+
+! local
+
+      integer l, ig
+      real rho                     ! density (kg/m3)
+      real ccntyp                  ! typical dust number density (#/kg)
+                                   ! (microphys = false)
+      real rdusttyp                ! typical dust radius (m)
+                                   ! (microphys = false)
+
+!==========================================================================
+
+      if (microphys) then ! improvedclouds
+         do l = 1,nlay
+            do ig = 1,ngrid
+               rho = pplay(ig,l)/(rnew(ig,l)*pt(ig,l))
+               surfdust(ig,l) = pq(ig,l,igcm_dust_number)*rho
+     $                          *tauscaling(ig)
+     $                          *4.*pi*rdust(ig,l)**2
+               surfice(ig,l)  = pq(ig,l,igcm_ccn_number)*rho
+     $                          *tauscaling(ig)
+     $                          *4.*pi*rice(ig,l)**2
+            end do
+         end do
+      else            ! simpleclouds
+         do l = 1,nlay
+            do ig = 1,ngrid
+               rho = pplay(ig,l)/(rnew(ig,l)*pt(ig,l))
+               rdusttyp = max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
+               ccntyp = 1.3e+8*max(tau(ig,1),0.001)/0.1
+     $                  *exp(-pzlay(ig,l)/10000.)
+               ccntyp = ccntyp/ccn_factor
+               if (rice(ig,l) .gt. rdust(ig,l)) then
+                  surfdust(ig,l) = ccntyp*(ccn_factor - 1.)*rho
+     $                             *4.*pi*rdusttyp**2
+                  surfice(ig,l)  = ccntyp*4.*pi*rice(ig,l)**2
+               else
+                  surfdust(ig,l) = ccntyp*ccn_factor*rho
+     $                             *4.*pi*rdusttyp**2
+                  surfice(ig,l)  = 0.
+               end if
+            end do
+         end do
+      end if         ! of microphys
+
+! write diagnostics in micron2/cm3
+
+      call wstats(ngrid,"surfdust", "Dust surface area",
+     $            "micron2 cm-3",3,surfdust*1.e6)
+      call writediagfi(ngrid,"surfdust", "Dust cloud surface area",
+     $            "micron2 cm-3",3,surfdust*1.e6)
+      call wstats(ngrid,"surfice", "Ice cloud surface area",
+     $            "micron2 cm-3",3,surfice*1.e6)
+      call writediagfi(ngrid,"surfice", "Ice cloud surface area",
+     $            "micron2 cm-3",3,surfice*1.e6)
+
+      return
+      end
Index: /trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 458)
+++ /trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 459)
@@ -139,5 +139,5 @@
             nqchem_start=dustbin+1
           end if
-          do iq=nqchem_start,nqchem_start+ncomp-1
+          do iq=nqchem_start,nqchem_start+ncomp-2
             noms(iq)=nomchem(iq-nqchem_start+1)
             mmol(iq)=mmolchem(iq-nqchem_start+1)
@@ -229,4 +229,5 @@
       igcm_ho2=0
       igcm_h2o2=0
+      igcm_ch4=0
       igcm_n2=0
       igcm_ar=0
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 458)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 459)
@@ -201,6 +201,6 @@
       real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
       real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
-      REAL surfdust(ngridmx,nlayermx) ! dust surface area (micron2/cm3, if photochemistry)
-      REAL surfice(ngridmx,nlayermx)  !  ice surface area (micron2/cm3, if photochemistry)
+      REAL surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3, if photochemistry)
+      REAL surfice(ngridmx,nlayermx)  ! ice surface area  (m2/m3, if photochemistry)
 
 c     Variables used by the slope model
@@ -1022,58 +1022,5 @@
          END IF  ! of IF (water)
 
-
-c   7b. Chemical species
-c     ------------------
-
-#ifndef MESOSCALE
-c        --------------
-c        photochemistry :
-c        --------------
-         IF (photochem .or. thermochem) then
-!NB: Photochemistry includes condensation of H2O2
-            PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.'
-            call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
-     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
-     $                   zdqcloud,zdqscloud,tauref,co2ice,
-     $                   pu,pdu,pv,pdv,surfdust,surfice)
-
-           ! increment values of tracers:
-           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
-                      ! tracers is zero anyways
-             DO l=1,nlayer
-               DO ig=1,ngrid
-                 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
-               ENDDO
-             ENDDO
-           ENDDO ! of DO iq=1,nq
-           ! add condensation tendency for H2O2
-           if (igcm_h2o2.ne.0) then
-             DO l=1,nlayer
-               DO ig=1,ngrid
-                 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
-     &                                +zdqcloud(ig,l,igcm_h2o2)
-               ENDDO
-             ENDDO
-           endif
-
-           ! increment surface values of tracers:
-           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
-                      ! tracers is zero anyways
-             DO ig=1,ngrid
-               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
-             ENDDO
-           ENDDO ! of DO iq=1,nq
-           ! add condensation tendency for H2O2
-           if (igcm_h2o2.ne.0) then
-             DO ig=1,ngrid
-               dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
-     &                                +zdqscloud(ig,igcm_h2o2)
-             ENDDO
-           endif
-
-         END IF  ! of IF (photochem.or.thermochem)
-#endif
-
-c   7c. Aerosol particles
+c   7b. Aerosol particles
 c     -------------------
 
@@ -1130,4 +1077,62 @@
            ENDDO
          END IF   ! of IF (sedimentation)
+
+c   7c. Chemical species
+c     ------------------
+
+#ifndef MESOSCALE
+c        --------------
+c        photochemistry :
+c        --------------
+         IF (photochem .or. thermochem) then
+
+!           dust and ice surface area
+            call surfacearea(ngrid, nlayer, pplay, zzlay, pt, pq, nq, 
+     $                       rdust, rice, tau, tauscaling, 
+     $                       surfdust, surfice)
+!           call photochemistry
+            call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
+     $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
+     $                   zdqcloud,zdqscloud,tauref,co2ice,
+     $                   pu,pdu,pv,pdv,surfdust,surfice)
+
+           ! increment values of tracers:
+           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
+                      ! tracers is zero anyways
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq)
+               ENDDO
+             ENDDO
+           ENDDO ! of DO iq=1,nq
+           
+           ! add condensation tendency for H2O2
+           if (igcm_h2o2.ne.0) then
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2)
+     &                                +zdqcloud(ig,l,igcm_h2o2)
+               ENDDO
+             ENDDO
+           endif
+
+           ! increment surface values of tracers:
+           DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry
+                      ! tracers is zero anyways
+             DO ig=1,ngrid
+               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
+             ENDDO
+           ENDDO ! of DO iq=1,nq
+
+           ! add condensation tendency for H2O2
+           if (igcm_h2o2.ne.0) then
+             DO ig=1,ngrid
+               dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
+     &                                +zdqscloud(ig,igcm_h2o2)
+             ENDDO
+           endif
+           
+         END IF  ! of IF (photochem.or.thermochem)
+#endif
 
 
@@ -1743,4 +1748,6 @@
 c             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
 c     &                        'm',3,rdust*ref_r0)
+c             call WRITEDIAGFI(ngridmx,'rdust','rdust',
+c    &                        'm',3,rdust)
 c             call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
 c     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
Index: /trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 458)
+++ /trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 459)
@@ -122,34 +122,4 @@
       end do
 
-c=======================================================================
-
-!!!!!!!!!! FOR PHOTOCHEMISTRY, REIMPLEMENT output of surfdust/surfice
-!!      if (photochem) then
-!!c        computation of dust and ice surface area (micron2/cm3)
-!!c        for heterogeneous chemistry
-!!
-!!         do l = 1,nlay
-!!            do ig = 1,ngrid
-!!c                                                                      
-!!c              npart: number density of ccn in #/cm3    
-!!c                                                     
-!!               npart(ig,l) = 1.e-6*ccn(ig,l)  
-!!     $                       *masse(ig,l)/epaisseur(ig,l) 
-!!c
-!!c              dust and ice surface area
-!!c                                                                  
-!!               surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2
-!!c                                                                  
-!!               if (rice(ig,l) .ge. rdust(ig,l)) then                   
-!!                  surfice(ig,l)  = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2
-!!                  surfdust(ig,l) = 0.
-!!               else                                                    
-!!                  surfice(ig,l) = 0.                                 
-!!               end if                                             
-!!            end do      ! of do ig=1,ngrid
-!!         end do         ! of do l=1,nlay
-!!      end if            ! of photochem
-
-
       RETURN
       END
