Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 1648)
+++ /trunk/LMDZ.MARS/README	(revision 1649)
@@ -2369,2 +2369,7 @@
 -  Moved CO2 clouds outuputs from physiq_mod.F to co2cloud.F
 -  Removed some CO2 clouds properties computation in physiq_mod.F (now done in post-treatment)
+
+== 24/01/2017 == JA
+-Modification of CO2 clouds routines : removal of some screen ouputs, modification for 3D of the writediagfi outputs in co2cloud.F, improvedCO2clouds.F and physiq_mod.F
+-callphys.def.co2clouds was added in deftank. It contains the necessary flags and constant for including co2clouds
+-traceur.def.co2clouds was added in deftank. It supplies the necessary tracers for running the co2 clouds scheme
Index: /trunk/LMDZ.MARS/deftank/callphys.def.co2clouds
===================================================================
--- /trunk/LMDZ.MARS/deftank/callphys.def.co2clouds	(revision 1649)
+++ /trunk/LMDZ.MARS/deftank/callphys.def.co2clouds	(revision 1649)
@@ -0,0 +1,168 @@
+##General options
+##~~~~~~~~~~~~~~~
+#Run with or without tracer transport ?
+tracer=.true.
+
+#Directory where external input files are:
+datadir=/home/audouard/GCM_MARS/datadir
+
+#Diurnal cycle ?  if diurnal=False, diurnal averaged solar heating
+diurnal=.true.
+
+#Seasonal cycle ? if season=False, Ls stays constant, to value set in "start"
+season = .true.   
+
+#write some more output on the screen ?
+lwrite = .false.    
+
+#Save statistics in file "stats.nc" ?
+callstats =.true.
+
+#Save EOF profiles in file "profiles" for Climate Database?
+calleofdump = .false.
+
+## Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+## ~~~~~~~~~~~~~
+#  =1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario,
+#  =4 Mars Year 24 from TES assimilation (old version of MY24; dust_tes.nc file)
+#  =6 "cold" (low dust) scenario ; =7 "warm" (high dust) scenario
+#  =24 Mars Year 24  ; =25 Mars Year 25 (year with a global dust storm) ; ...
+#  =30 Mars Year 30 
+iaervar = 24
+# Dust opacity at 610 Pa (when constant, i.e. for the iaervar=1 case)
+tauvis=0.1
+# Dust vertical distribution: 
+# (=0: old distrib. (Pollack90), =1: top set by "topdustref",
+#  =2: Viking scenario; =3 MGS scenario)
+iddist  = 1
+# Dust top altitude (km). (Matters only if iddist=1)
+topdustref = 130.
+
+## Physical Parameterizations :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call radiative transfer ?
+callrad   = .true.
+# call NLTE radiative schemes ?   matters only if callrad=T
+callnlte = .true.
+# NLTE 15um scheme to use.
+# 0-> Old scheme, static oxygen
+# 1-> Old scheme, dynamic oxygen
+# 2-> New scheme
+nltemodel = 0
+# call CO2 NIR absorption ?   matters only if callrad=T
+callnirco2 = .true.
+# NIR NLTE correction ? matters only if callnirco2=T
+nircorr=0
+# call turbulent vertical diffusion ?
+calldifv  = .true.
+# call convective adjustment ?
+calladj  = .true.
+# Thermals
+calltherm = .true.
+callrichsl = .true.
+# call CO2 condensation ?
+callcond =.true.
+# call thermal conduction in the soil ?
+callsoil = .true.
+# call Lott's gravity wave/subgrid topography scheme ?
+calllott  = .true.
+# Impose polar cap surface albedos as observed by TES?
+TESicealbedo = .true.
+## Coefficient for Northern cap albedoes
+TESice_Ncoef=1.6
+## Coefficient for Southern cap albedoes
+TESice_Scoef=1.6
+
+
+## Radiative transfer options :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# the rad.transfer is computed every "iradia" physical timestep
+iradia    = 1
+# Output of the exchange coefficient mattrix ? for diagnostic only
+callg2d  = .false.
+# Rayleigh scattering : (should be .false. for now)
+rayleigh  = .false.
+
+## Tracer (dust water, ice and/or chemical species) options (used if tracer=T):
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# DUST: Transported dust ? (if >0, use 'dustbin' dust bins)
+dustbin    = 2
+# DUST: Radiatively active dust ? (matters if dustbin>0)
+active  = .true.
+# DUST: use mass and number mixing ratios to predict dust size ?
+doubleq   = .true.
+# DUST: lifted by GCM surface winds ?
+lifting = .false.
+# DUST: lifted by dust devils ?
+callddevil = .false.
+# DUST: Scavenging by H2O/CO2 snowfall ?
+scavenging = .true.
+# DUST/WATERICE: Gravitationnal sedimentation ?
+sedimentation = .true.
+# WATERICE: Radiatively active transported atmospheric water ice ?
+activice   = .false.
+# WATER: Compute water cycle
+water = .true.
+# WATER: Microphysical scheme for water-ice clouds?
+microphys = .true.
+# WATER: parameter contact
+mteta = 0.95
+# WATER: Effective variance for sedimentation for the log-normal 
+#        distribution of ice particles ?
+nuice_sed=0.1
+# WATER: current permanent caps at both poles. True IS RECOMMENDED
+#        (with .true., North cap is a source of water and South pole
+#         is a cold trap)
+caps  = .true.
+# WATER: Water ice albedo ?
+albedo_h2o_ice = 0.35
+# WATER: Water ice thermal inertia
+inert_h2o_ice = 800
+# WATER: Frost thickness threshold for albedo
+frost_albedo_threshold = 0.005
+# PHOTOCHEMISTRY: include chemical species 
+photochem  = .false.
+
+# SCATTERERS: set number of scatterers. must be compliant with preceding options.
+naerkind = 1
+
+#include co2 clouds scheme?
+co2clouds=.true.
+# Microphysical scheme for co2 clouds?
+microphysco2 = .true.
+#Effective variance for sedimentation for the log-normal 
+#        distribution of CO2 clouds particles ?
+nuiceco2_sed=0.1
+
+## Thermospheric options (relevant if tracer=T) :
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call thermosphere ?
+callthermos = .false.
+# WATER: included without cycle (only if water=.false.)
+thermoswater = .false.
+# call thermal conduction ? (only if callthermos=.true.)
+callconduct = .false.
+# call molecular viscosity ? (only if callthermos=.true.)
+callmolvis = .false.
+# call molecular diffusion ? (only if callthermos=.true.)
+callmoldiff = .false.
+# call thermospheric photochemistry ? (only if callthermos=.true.)
+thermochem = .false.
+# call EUV heating ? (only if callthermos=.true.)
+calleuv=.false.
+#Method to include solar variability?
+#0-> Old method (fixed EUV input)  1-> Variability with E10.7 as observed
+solvarmod=1
+# date for solar flux calculation: (1985 < date < 2002)
+## (Solar min=1996.4 ave=1993.4 max=1990.6) ; Only used if solvarmod=0
+solarcondate = 1993.4
+#Solar variability as observed for MY? (must bebetween MY23 and MY30)
+# (only matters if solvarmod=1)
+solvaryear=24
+# value for the UV heating efficiency
+##(experimental values between 0.19 and 0.23, lower values may
+## be used to compensate for low 15 um cooling)
+euveff = 0.21       
+
+
+
Index: /trunk/LMDZ.MARS/deftank/traceur.def.co2clouds
===================================================================
--- /trunk/LMDZ.MARS/deftank/traceur.def.co2clouds	(revision 1649)
+++ /trunk/LMDZ.MARS/deftank/traceur.def.co2clouds	(revision 1649)
@@ -0,0 +1,11 @@
+10
+co2
+co2_ice
+dust_number
+dust_mass
+ccnco2_number
+ccnco2_mass
+h2o_ice
+h2o_vap
+ccn_mass
+ccn_number
Index: /trunk/LMDZ.MARS/libf/phymars/co2cloud.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1648)
+++ /trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1649)
@@ -256,5 +256,5 @@
        DO l=1, nlay
           DO ig=1,ngrid
-             Niceco2=tempo_traceurs(ig,l,igcm_co2_ice)
+             Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30)
              Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number),
      &            1.e-30)
@@ -271,19 +271,19 @@
                 rdust(ig,l)=max(rdust(ig,l),1.e-9)    
              end if
-             rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 
-     &            + Qccnco2*rho_dust)
-     &            / (Niceco2 + Qccnco2)
+c             rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 
+c     &            + Qccnco2*rho_dust)
+c     &            / (Niceco2 + Qccnco2)
              riceco2(ig,l)= Niceco2*3.0/
      &            (4.0*rho_ice_co2*pi*Nccnco2)
      &            +rdust(ig,l)*rdust(ig,l)*rdust(ig,l)
              riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
-             write(*,*) "in co2clouds, rice = ",riceco2(ig,l)
-             write(*,*) "in co2clouds, rho = ",rhocloudco2t(ig,l)
-
-              call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
+c             write(*,*) "in co2clouds, rice = ",riceco2(ig,l)
+c             write(*,*) "in co2clouds, rho = ",rhocloudco2t(ig,l)
+
+             call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
      &            tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l)) 
-              write(*,*) "in co2clouds, rice update = ",riceco2(ig,l)
-              write(*,*) "in co2clouds, rho update = "
-     &             ,rhocloudco2t(ig,l)
+c              write(*,*) "in co2clouds, rice update = ",riceco2(ig,l)
+c              write(*,*) "in co2clouds, rho update = "
+c     &             ,rhocloudco2t(ig,l)
 
               rsedcloudco2(ig,l)=max(riceco2(ig,l)*
@@ -291,5 +291,5 @@
      &            rdust(ig,l))
              rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),5.e-4)
-             write(*,*) 'Rsedcloud = ',rsedcloudco2(ig,l)
+c             write(*,*) 'Rsedcloud = ',rsedcloudco2(ig,l)
              !write(*,*) 'Rhocloudco2 = ',rhocloudco2t(ig,l)
 
@@ -657,6 +657,6 @@
      &       pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)*
      &       tauscaling(ig),1.e-30)
-        rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 + Qccnco2*rho_dust)
-     &       / (Niceco2 + Qccnco2)
+c        rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 + Qccnco2*rho_dust)
+c     &       / (Niceco2 + Qccnco2)
 c        rhocloudco2(ig,l) = min(max(rhocloudco2t,rho_ice_co2),rho_dust)
 
@@ -669,15 +669,15 @@
 
         riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
-        write(*,*) "In co2cloud, after loop, riceco2 =",riceco2(ig,l)
-        write(*,*) "In co2cloud, after loop, rhoco2 ="
-     &       ,rhocloudco2t(ig,l)
+c        write(*,*) "In co2cloud, after loop, riceco2 =",riceco2(ig,l)
+c        write(*,*) "In co2cloud, after loop, rhoco2 ="
+c     &       ,rhocloudco2t(ig,l)
 
         call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
      &             tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
 
-        write(*,*) "In co2cloud, after loop and update, riceco2 ="
-     &       ,riceco2(ig,l)
-        write(*,*) "In co2cloud, after loop and update, rhoco2 ="
-     &       ,rhocloudco2t(ig,l)
+c        write(*,*) "In co2cloud, after loop and update, riceco2 ="
+c     &       ,riceco2(ig,l)
+c        write(*,*) "In co2cloud, after loop and update, rhoco2 ="
+c     &       ,rhocloudco2t(ig,l)
 
         if ( Niceco2 
@@ -752,11 +752,11 @@
      &              (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
                   
-               write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
+c               write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
             enddo
          enddo
-       call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",1,
+       call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
      &        satuco2)
        call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
-     &        ,1,riceco2)
+     &        ,3,riceco2)
 ! or output in diagfi.nc (for testphys1d)
 c         call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
@@ -764,5 +764,5 @@
 c     &                       'K JA',1,pt)
          
-      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",1,
+      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
      &   rsedcloudco2)
       
Index: /trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F	(revision 1648)
+++ /trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F	(revision 1649)
@@ -302,6 +302,6 @@
 
            IF ( satu .ge. 1d0 ) THEN ! if there is condensation
-              write(*,*)
-              write(*,*) "l, pco2, satu= ",l,pco2,satu
+c              write(*,*)
+c              write(*,*) "l, pco2, satu= ",l,pco2,satu
 c              Masse_atm=mmean(ig,l)*1.e-3*pplay(ig,l)/rgp/zt(ig,l) !Kg par couche
 
@@ -309,5 +309,5 @@
               call updaterccn(zq(ig,l,igcm_dust_mass),
      &             zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
-              write(*,*) "Improved, l,Rdust = ",l,rdust(ig,l)
+c              write(*,*) "Improved, l,Rdust = ",l,rdust(ig,l)
 
               rdust(ig,l)= zq(ig,l,igcm_dust_mass)
@@ -315,15 +315,15 @@
      &             / zq(ig,l,igcm_dust_number)
               rdust(ig,l)= rdust(ig,l)**(1./3.)
-              write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l)
+c              write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l)
             rdust(ig,l)=max(1.e-9,rdust(ig,l))
             rdust(ig,l)=min(2e-6,rdust(ig,l))
-              write(*,*) "Improved3, l,Rdust = ",l,rdust(ig,l)
+c              write(*,*) "Improved3, l,Rdust = ",l,rdust(ig,l)
 
 c       Expand the dust moments into a binned distribution
               Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
               No = zq(ig,l,igcm_dust_number)* tauscaling(ig)
-              write(*,*) "dust number, mass = ",
-     &             zq(ig,l,igcm_dust_number)* tauscaling(ig),
-     &             zq(ig,l,igcm_dust_mass)* tauscaling(ig)
+c              write(*,*) "dust number, mass = ",
+c     &             zq(ig,l,igcm_dust_number)* tauscaling(ig),
+c     &             zq(ig,l,igcm_dust_mass)* tauscaling(ig)
 c              write(*,*) "No, Mo = ",No, Mo 
               Rn = rdust(ig,l)
@@ -461,10 +461,10 @@
 c + enlever les CCN a la distri de dust
 
-              write(*,*) "new dust_mass, number =",
-     &             zq(ig,l,igcm_dust_mass)* tauscaling(ig), 
-     &             zq(ig,l,igcm_dust_number)*tauscaling(ig)
-              write(*,*) "new ccn mass, number =",
-     &             zq(ig,l,igcm_ccnco2_mass)* tauscaling(ig)
-     &             ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)
+c              write(*,*) "new dust_mass, number =",
+c     &             zq(ig,l,igcm_dust_mass)* tauscaling(ig), 
+c     &             zq(ig,l,igcm_dust_number)*tauscaling(ig)
+c              write(*,*) "new ccn mass, number =",
+c     &             zq(ig,l,igcm_ccnco2_mass)* tauscaling(ig)
+c     &             ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)
            
            ENDIF                ! of is satu >1
@@ -497,5 +497,5 @@
      &             zq(ig,l,igcm_dust_number) .lt. 0. ) THEN
                
-                 write(*,*) "before growth CCN N,M = "          
+                 write(*,*) "CO2 clouds before growth CCN N,M = "          
      &                ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)
      &                ,zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig)
@@ -522,6 +522,6 @@
           ! WATCH OUT: CO2 nuclei is supposed to be dust
           ! only when deriving rhocloud (otherwise would need to keep info on  water embedded in co2) - listo
-              write(*,*) "Rdust before growth = ",rdust(ig,l)
-              write(*,*) "Riceco2 before growth = ",riceco2(ig,l)
+c              write(*,*) "Rdust before growth = ",rdust(ig,l)
+c              write(*,*) "Riceco2 before growth = ",riceco2(ig,l)
 
               !! Niceco2,Qccnco2,Nccnco2
@@ -531,5 +531,5 @@
               call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
      &             tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 
-              write(*,*) "Riceco2 update before growth = ",riceco2(ig,l)
+c              write(*,*) "Riceco2 update before growth = ",riceco2(ig,l)
 
               No   = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)
@@ -550,5 +550,5 @@
      &             riceco2(ig,l)*rho_ice_co2)*Ic_rice
               dMice =  No * Ic_rice * ptimestep ! Kg par kg d'air, <0 si croissance !           
-              write(*,*) "dMicev0 in improved = " , dMice
+c              write(*,*) "dMicev0 in improved = " , dMice
 
              if (dMice .gt. 0) dMice = min(dMice,zq0(ig,l,igcm_co2_ice))
@@ -561,5 +561,5 @@
 c              write(*,*) "riceco2+dr/dt = ", riceco2(ig,l)
 
-              write(*,*) "dMice in improved = " , dMice
+c              write(*,*) "dMice in improved = " , dMice
               
               
@@ -574,5 +574,5 @@
      &             *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
      &             *rdust(ig,l) )**(1.0/3.0)
-              write(*,*) "new riceco2 = ",riceco2(ig,l) 
+c              write(*,*) "new riceco2 = ",riceco2(ig,l) 
                
               !! Niceco2,Qccnco2,Nccnco2
@@ -582,5 +582,5 @@
               call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
      &             tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 
-              write(*,*) "new riceco2 updaterad= ",riceco2(ig,l) 
+c              write(*,*) "new riceco2 updaterad= ",riceco2(ig,l) 
 
 ! latent heat release        
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1648)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1649)
@@ -1165,9 +1165,5 @@
      &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
      &           rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2)
-            
-            
-            call WRITEdiagfi(ngrid,"rhocloudco2","rho cloud co2","kk",1,
-     &           rhocloudco2)
-            
+                
 
 c Temperature variation due to latent heat release
@@ -2282,4 +2278,11 @@
 c        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
 c    &                   'w.m-2',3,zdtlw)
+  
+         if (co2clouds) then
+            do iq=1, nq
+               call WRITEDIAGFI(ngrid,trim(noms(iq)),
+     &              trim(noms(iq)),'kg/kg',3,zq(:,:,iq))
+            end do
+         endif
             if (.not.activice) then
                CALL WRITEDIAGFI(ngrid,'tauTESap',
