Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/README
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/README	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/README	(revision 214)
@@ -0,0 +1,29 @@
+> tar xvf LMDZ.MARS.source.tar
+ajouter perso
+change resol in compile
+check NETCDF NCDFLIB NCDFINC
+> compile
+> mkdir myGCM
+> cd myGCM
+lier z2sig.def (voir what_I_need)
+lier un start_archive avec bon nombre traceurs
+> ln -sf ../gcm.e
+> ln -sf ../newstart.e
+> ln -sf ../deftank/run0 
+lier un callphys.def
+lier un run.def
+cp num_run
+> newstart.e
+> mv restart.nc start0.nc
+> mv restartfi.nc startfi0.nc
+LANCER AVEC
+> nohup run_mcd_3days > log_run_mcd_3days &
+ATTENDRE LONGTEMPS ...
+> mv tmpsave ../STARTBASE_64_48_32_t2
+> rm -rf start.nc start0.nc startfi.nc startfi0.nc
+
+lier le run.def simu
+modifier physiq.F pour sortir des donnees 
+recompiler
+link diagfi.nc vers un espace temp
+ex: ln -sf $MODELS/TMPDIR/GCMINI/diagfi.nc
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile	(revision 214)
@@ -0,0 +1,68 @@
+#! /bin/bash
+
+## changer eventuellement le callphys.def egalement
+
+##################
+tracers=2
+nx=64
+ny=48
+nz=32
+##################
+tracers=1
+nx=64
+ny=48
+nz=25
+##################
+tracers=1
+nx=64
+ny=48
+nz=18
+##################
+#tracers=1
+#nx=64
+#ny=48
+#nz=32
+##################
+tracers=5
+nx=64
+ny=48
+nz=25
+##################
+#tracers=2
+#nx=64
+#ny=48
+#nz=25
+##################
+
+
+
+
+# GCM environment variables
+#--------------------------
+export LMDGCM=$PWD
+export LIBOGCM=$PWD/libo
+
+echo $LMDGCM
+echo $LIBOGCM
+echo $NCDFLIB
+echo $NCDFINC 
+
+#rm -rf libo
+mkdir libo
+touch libo/log_compile_newstart
+touch libo/log_compile_gcm
+
+echo bricolo
+\rm libf/grid/dimensions.h 2> /dev/null
+
+echo compile newstart
+echo check in $PWD/libo/log_compile_newstart
+makegcm -t ${tracers} -p mars -d ${nx}x${ny}x${nz} newstart | tee libo/log_compile_newstart | grep 'warnings' >> libo/log_compile_newstart
+
+echo compile gcm
+echo check in $PWD/libo/log_compile_gcm
+makegcm -t ${tracers} -p mars -d ${nx}x${ny}x${nz} gcm | tee libo/log_compile_gcm | grep 'warnings' >> libo/log_compile_gcm
+
+echo well, OK.
+
+
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/diff_avec_premiere_version_new
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/diff_avec_premiere_version_new	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/diff_avec_premiere_version_new	(revision 214)
@@ -0,0 +1,63 @@
+21,25d20
+< #tracers=1
+< #nx=64
+< #ny=48
+< #nz=32
+< ##################
+Only in /donnees/aslmd/SVN/trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_new/myGCM: .svn
+diff myGCM/callphys.def /donnees/aslmd/SVN/trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_new/myGCM/callphys.def
+5d4
+< ##true???
+17c16
+< stats =.false.
+---
+> stats =.true.
+116,122d114
+< ## albedoes of CO2 ice caps from TES ?
+< TESicealbedo=.true.
+< ## Coefficient for Northern cap albedoes
+< TESice_Ncoef=1.4
+< ## Coefficient for Southern cap albedoes
+< TESice_Scoef=1.4
+< 
+diff myGCM/launch_gcm /donnees/aslmd/SVN/trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_new/myGCM/launch_gcm
+15,19c15
+< #startbase='../../STARTS/'
+< #startbase='/d2/emlmd/run_LMDZ.MARS.BETA.TASI/tau1/'
+< #startbase='/tmp7/aslmd/flush/'
+< #startbase='/tmp7/aslmd/flush2/'
+< startbase='/d6/vblmd/MERIDIANI_EXOMARS/starts_en_cours/tau05/'
+---
+> startbase='../../STARTS/'
+26,29d21
+< if [[ "$daygcm" == 0 ]]
+< then
+< daygcm=669 
+< fi
+5d4
+< ##true???
+17c16
+< stats =.false.
+---
+> stats =.true.
+116,122d114
+< ## albedoes of CO2 ice caps from TES ?
+< TESicealbedo=.true.
+< ## Coefficient for Northern cap albedoes
+< TESice_Ncoef=1.4
+< ## Coefficient for Southern cap albedoes
+< TESice_Scoef=1.4
+< 
+15,19c15
+< #startbase='../../STARTS/'
+< #startbase='/d2/emlmd/run_LMDZ.MARS.BETA.TASI/tau1/'
+< #startbase='/tmp7/aslmd/flush/'
+< #startbase='/tmp7/aslmd/flush2/'
+< startbase='/d6/vblmd/MERIDIANI_EXOMARS/starts_en_cours/tau05/'
+---
+> startbase='../../STARTS/'
+26,29d21
+< if [[ "$daygcm" == 0 ]]
+< then
+< daygcm=669 
+< fi
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def	(revision 214)
@@ -0,0 +1,122 @@
+##General options
+##~~~~~~~~~~~~~~~
+#Run with or without tracer transport ?
+tracer=.false.
+
+#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" ?
+stats =.false.
+
+#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 (same as =24 for now)
+#  =24 Mars Year 24 from TES assimilation (ie: MCD reference case)
+#  =25 Mars Year 25 from TES assimilation (ie: a year with a global dust storm)
+#  =26 Mars Year 26 from TES assimilation
+iaervar = 24 
+# Dust opacity at 700 Pa (when constant, i.e. for the iaervar=1 case)
+tauvis=0.2
+# Dust vertical distribution: 
+# (=0: old distrib. (Pollack90), =1: top set by "topdustref",
+#  =2: Viking scenario; =3 MGS scenario)
+iddist  = 3
+# Dust top altitude (km). (Matters only if iddist=1)
+topdustref = 55.
+
+## Physical Parameterizations :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call radiative transfer ?
+callrad   = .true.
+# call NLTE radiative schemes ?   matters only if callrad=T
+callnlte = .true.
+# call CO2 NIR absorption ?   matters only if callrad=T
+callnirco2 = .true.
+# call turbulent vertical diffusion ?
+calldifv  = .true.
+# call convective adjustment ?
+calladj  = .true.
+# call CO2 condensation ?
+callcond =.true.
+# call thermal conduction in the soil ?
+callsoil = .true.
+# call Lott's gravity wave/subgrid topography scheme ?
+calllott  = .true.
+
+## 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    = 0
+# DUST: Radiatively active dust ? (matters if dustbin>0)
+active  = .false.
+# DUST: use mass and number mixing ratios to predict dust size ?
+# (must also have dustbin=1)
+doubleq   = .false.
+# DUST: lifted by GCM surface winds ?
+lifting = .false.
+# DUST: lifted by dust devils ?
+callddevil = .false.
+# DUST: Scavenging by CO2 snowfall ?
+scavenging = .false.
+# DUST/WATERICE: Gravitationnal sedimentation ?
+sedimentation = .false.
+# WATERICE: Radiatively active transported atmospheric water ice ?
+activice   = .false.
+# WATER: Compute water cycle
+water = .false.
+# 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.
+# PHOTOCHEMISTRY: include chemical species 
+photochem  = .false.
+
+## 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 EUV heating ? (only if callthermos=.true.)
+calleuv=.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.
+# date for solar flux calculation: (1985 < date < 2002)
+## (Solar min=1996.4 ave=1993.4 max=1990.6)
+solarcondate = 1993.4
+       
+
+## albedoes of CO2 ice caps from TES ?
+TESicealbedo=.true.
+## Coefficient for Northern cap albedoes
+TESice_Ncoef=1.4
+## Coefficient for Southern cap albedoes
+TESice_Scoef=1.4
+
+
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def.csttau
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def.csttau	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def.csttau	(revision 214)
@@ -0,0 +1,122 @@
+##General options
+##~~~~~~~~~~~~~~~
+#Run with or without tracer transport ?
+tracer=.true.
+
+#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" ?
+stats =.false.
+
+#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 (same as =24 for now)
+#  =24 Mars Year 24 from TES assimilation (ie: MCD reference case)
+#  =25 Mars Year 25 from TES assimilation (ie: a year with a global dust storm)
+#  =26 Mars Year 26 from TES assimilation
+iaervar = 1 
+# Dust opacity at 700 Pa (when constant, i.e. for the iaervar=1 case)
+tauvis=0.5
+# Dust vertical distribution: 
+# (=0: old distrib. (Pollack90), =1: top set by "topdustref",
+#  =2: Viking scenario; =3 MGS scenario)
+iddist  = 3
+# Dust top altitude (km). (Matters only if iddist=1)
+topdustref = 55.
+
+## Physical Parameterizations :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call radiative transfer ?
+callrad   = .true.
+# call NLTE radiative schemes ?   matters only if callrad=T
+callnlte = .true.
+# call CO2 NIR absorption ?   matters only if callrad=T
+callnirco2 = .true.
+# call turbulent vertical diffusion ?
+calldifv  = .true.
+# call convective adjustment ?
+calladj  = .true.
+# call CO2 condensation ?
+callcond =.true.
+# call thermal conduction in the soil ?
+callsoil = .true.
+# call Lott's gravity wave/subgrid topography scheme ?
+calllott  = .true.
+
+## 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    = 0
+# DUST: Radiatively active dust ? (matters if dustbin>0)
+active  = .false.
+# DUST: use mass and number mixing ratios to predict dust size ?
+# (must also have dustbin=1)
+doubleq   = .false.
+# DUST: lifted by GCM surface winds ?
+lifting = .false.
+# DUST: lifted by dust devils ?
+callddevil = .false.
+# DUST: Scavenging by CO2 snowfall ?
+scavenging = .false.
+# DUST/WATERICE: Gravitationnal sedimentation ?
+sedimentation = .false.
+# WATERICE: Radiatively active transported atmospheric water ice ?
+activice   = .false.
+# WATER: Compute water cycle
+water = .false.
+# 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.
+# PHOTOCHEMISTRY: include chemical species 
+photochem  = .false.
+
+## 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 EUV heating ? (only if callthermos=.true.)
+calleuv=.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.
+# date for solar flux calculation: (1985 < date < 2002)
+## (Solar min=1996.4 ave=1993.4 max=1990.6)
+solarcondate = 1993.4
+       
+
+## albedoes of CO2 ice caps from TES ?
+TESicealbedo=.true.
+## Coefficient for Northern cap albedoes
+TESice_Ncoef=1.4
+## Coefficient for Southern cap albedoes
+TESice_Scoef=1.4
+
+
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def.water
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def.water	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/callphys.def.water	(revision 214)
@@ -0,0 +1,94 @@
+General options
+~~~~~~~~~~~~~~~
+tracer    (Run with or without tracer transport ?)
+T
+diurnal   (Diurnal cycle ?  if diurnal=F, diurnal averaged solar heating)
+T
+season    (Seasonal cycle ? if season=F, Ls stays constant like in "start")
+T
+lwrite    (want some more output on the screen ?) 
+F
+stats     (Saving statistics in file "stats" ?)
+F
+calleofdump (Saving EOF profiles in file "profiles" for Climate Database ?)
+F
+Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+~~~~~~~~~~~~~
+iaervar  (=1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario)
+4        (=4 Mars Year 24 from TES assimilation) 
+iddist  (Dust vertical distribution: =0: old distrib. (Pollack90) 
+3       (=1: top set by "topdustref"; =2: Viking scenario; =3 MGS scenario )
+topdustref (Dust top altitude (km). Matter only if iddist=1)
+55.
+Physical Parameterizations :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+callrad   (call radiative transfer ?)
+T
+callnlte (call NLTE radiative schemes ?   matter only if callrad=T)
+T
+callnirco2 (call CO2 NIR absorption ?   matter only if callrad=T)
+T
+calldifv  (call turbulent vertical diffusion ?)
+T
+calladj   (call convective adjustment ?)
+T
+callcond  (call CO2 condensation ?)
+T
+callsoil  (call thermal conduction in the soil ?)
+T
+calllott  (call Lott's gravity wave/subgrid topography scheme ?)
+T
+Radiative transfer options :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+iradia    (the rad.transfer is computed every "iradia" physical timestep)
+1
+callg2d   (Output of the exchange coefficient mattrix ? for diagnostic only)
+F
+rayleigh  (Rayleigh scattering : should be =F for now)
+F
+Tracer (dust water, ice and/or chemical species) options (use if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dustbin    (DUST: Transported dust ? (if >0, uses q(1) to q(dustbin))
+0
+active     (DUST: Radiatively active dust ? (uses q(1) to q(dustbin))
+F
+doubleq    (DUST: needs dustbin=1, use mass (q(1)) and number (q(2)) mixing ratio to predict dust size ?)
+F
+lifting    (DUST: lifted by GCM surface winds ?)
+F
+dustdevil  (DUST: lifted by dust devils ?)
+F
+scavenging (DUST: Scavenging by CO2 snowfall ?)
+F
+sedimentation (DUST/WATERICE: Gravitationnal sedimentation ?)
+T              (set to T to simulate water cycle with clouds)                
+iceparty   (WATERICE: includes water ice: q(nqmx-1)). "water" must be T !
+T
+activice   (WATERICE: Radiatively active transported atmospheric water ice ?)
+F
+water      (WATER: Compute water cycle using q(nqmx) )
+T
+caps       (WATER: current permanent caps at both poles. T IS RECOMMENDED
+T          (caps=T means Ncap is a source of water, S pole is a cold trap)
+photochem  (PHOTOCHEMISTRY: chemical species included)
+F
+Thermospheric options (relevant if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+callthermos  (call thermosphere ?)
+F
+thermoswater  (WATER: included without cycle	only if water=F)
+F
+callconduct  (call thermal conduction ?     matter only if callthermos=T)
+F
+calleuv  (call EUV heating ?                matter only if callthermos=T)
+F
+callmolvis  (call molecular viscosity ?     matter only if callthermos=T)
+F
+callmoldiff  (call molecular diffusion ?    matter only if callthermos=T)
+F
+thermochem  (call thermospheric photochemistry ?  matter only if callthermos=T)
+F
+solarcondate (date for solar flux calculation: 1985 < date < 2002))
+1993.4       (Solar min=1996.4 ave=1993.4 max=1990.6)
+
+
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/launch_gcm
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/launch_gcm	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/launch_gcm	(revision 214)
@@ -0,0 +1,74 @@
+#! /bin/bash
+
+##############################
+# use: - launch_gcm          #
+#      - echo 2 | launch_gcm #
+##############################
+
+# Author : A. Spiga - Spring/Winter 2008
+
+################## penser a recompiler si on change ceci
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_32_t2'
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/old_STARTBASE_64_48_25'
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_25'
+#startbase='/tmp7/emlmd/run_LMDZ.MARS.BETA.fit.widor/RUN/'
+#startbase='../../STARTS/'
+#startbase='/d2/emlmd/run_LMDZ.MARS.BETA.TASI/tau1/'
+#startbase='/tmp7/aslmd/flush/'
+#startbase='/tmp7/aslmd/flush2/'
+#startbase='/d6/vblmd/MERIDIANI_EXOMARS/starts_en_cours/tau05/'
+startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_25_t2dust_JBM/'
+##################
+
+#
+# get simulation start day and compute start/startfi required date (simulation day - 1)
+#
+echo '>>>' choose start sol ? ; read daygcm
+if [[ "$daygcm" == 0 ]]
+then
+daygcm=669 
+fi
+daygcm=$(expr ${daygcm} - 1)
+
+#
+# find the correct start/startfi files among 3-day startbase
+#
+ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null
+if [[ "$?" == 0 ]]
+then
+    echo found ! ${daygcm}
+else
+    #echo not found ! ${daygcm}	
+    daygcm=$(expr ${daygcm} - 1)
+    ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null
+    if [[ "$?" == 0 ]]     # probably a 'while' loop would be better if more than 3
+    then
+	echo found ! ${daygcm}
+    else
+	#echo not found ! ${daygcm}    
+	daygcm=$(expr ${daygcm} - 1)
+        ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null	
+        echo found ! ${daygcm}
+    fi	
+fi	
+
+#
+# link start/startfi files
+#
+ln -sf ${startbase}/start${daygcm}.nc start.nc
+ln -sf ${startbase}/startfi${daygcm}.nc startfi.nc
+
+#
+# clean the place
+#
+rm -rf restart.nc 2> /dev/null 
+rm -rf restartfi.nc 2> /dev/null
+rm -rf log_gcm 2> /dev/null
+
+#
+# run !
+#
+echo running the GCM - check in ${PWD}/log_gcm 
+nohup gcm.e > log_gcm
+
+exit
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/num_run.STARTBASE
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/num_run.STARTBASE	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/num_run.STARTBASE	(revision 214)
@@ -0,0 +1,1 @@
+0
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run.def
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run.def	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run.def	(revision 214)
@@ -0,0 +1,112 @@
+#
+#----------------------------------------------------------------------- 
+#Parametres de controle du run:                                          
+#------------------------------                                          
+
+# Nombre de jours d'integration                                         
+     nday=20
+
+# nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step = 960
+
+# periode pour le pas Matsuno (en pas)                                  
+  iperiod=5
+
+# periode de sortie des variables de controle (en pas)                  
+  iconser=120
+
+# periode d'ecriture du fichier histoire (en jour)                      
+    iecri=100
+
+# periode de stockage fichier histmoy (en jour)                         
+ periodav=60.
+
+# periode de la dissipation (en pas)                                    
+  idissip=5
+
+# choix de l'operateur de dissipation (star ou  non star )              
+ lstardis=.true.
+
+# avec ou sans coordonnee hybrides
+ hybrid=.true.
+
+# nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv=1
+
+# nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot=2
+
+# nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh=2
+
+# temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv= 3500.
+
+# temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot=4500.
+
+# temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp=4500.
+
+# coefficient pour gamdissip                                            
+  coefdis=0.
+
+# choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats=.false.
+
+# avec ou sans physique                                                 
+   physic=.true.
+
+# periode de la physique (en pas)                                       
+  iphysiq=20
+
+# choix d'une grille reguliere                                          
+  grireg=.true.
+
+# frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy=40
+
+# longitude en degres du centre du zoom                                 
+   clon=63.
+
+# latitude en degres du centre du zoom                                  
+   clat=0.
+
+# facteur de grossissement du zoom,selon longitude                      
+  grossismx=1.
+
+# facteur de grossissement du zoom ,selon latitude                      
+ grossismy=1.
+
+#  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb=.false.
+
+# extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx= 0.
+
+# extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy=0.
+
+#  raideur du zoom en  X
+    taux=2.
+
+#  raideur du zoom en  Y
+    tauy=2.
+
+#  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus= .false.
+
+# Avec sponge layer
+  callsponge  = .true.
+ 
+# Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge= 2
+
+# Sponge:  hauteur de sponge (km)
+  hsponge= 130
+
+# Sponge:  tetasponge (secondes)
+  tetasponge = 50000
+
+# some definitions for the physics, in file 'callphys.def'
+INCLUDEDEF=callphys.def
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run.def.64x48x25
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run.def.64x48x25	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run.def.64x48x25	(revision 214)
@@ -0,0 +1,112 @@
+#
+#----------------------------------------------------------------------- 
+#Parametres de controle du run:                                          
+#------------------------------                                          
+
+# Nombre de jours d'integration                                         
+     nday=10
+
+# nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step = 960
+
+# periode pour le pas Matsuno (en pas)                                  
+  iperiod=5
+
+# periode de sortie des variables de controle (en pas)                  
+  iconser=120
+
+# periode d'ecriture du fichier histoire (en jour)                      
+    iecri=100
+
+# periode de stockage fichier histmoy (en jour)                         
+ periodav=60.
+
+# periode de la dissipation (en pas)                                    
+  idissip=5
+
+# choix de l'operateur de dissipation (star ou  non star )              
+ lstardis=.true.
+
+# avec ou sans coordonnee hybrides
+ hybrid=.true.
+
+# nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv=1
+
+# nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot=2
+
+# nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh=2
+
+# temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv= 4000.
+
+# temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot=5000.
+
+# temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp=5000.
+
+# coefficient pour gamdissip                                            
+  coefdis=0.
+
+# choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats=.false.
+
+# avec ou sans physique                                                 
+   physic=.true.
+
+# periode de la physique (en pas)                                       
+  iphysiq=20
+
+# choix d'une grille reguliere                                          
+  grireg=.true.
+
+# frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy=40
+
+# longitude en degres du centre du zoom                                 
+   clon=63.
+
+# latitude en degres du centre du zoom                                  
+   clat=0.
+
+# facteur de grossissement du zoom,selon longitude                      
+  grossismx=1.
+
+# facteur de grossissement du zoom ,selon latitude                      
+ grossismy=1.
+
+#  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb=.false.
+
+# extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx= 0.
+
+# extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy=0.
+
+#  raideur du zoom en  X
+    taux=2.
+
+#  raideur du zoom en  Y
+    tauy=2.
+
+#  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus= .false.
+
+# Avec sponge layer
+  callsponge  = .true.
+ 
+# Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge= 2
+
+# Sponge:  hauteur de sponge (km)
+  hsponge= 90
+
+# Sponge:  tetasponge (secondes)
+  tetasponge = 50000
+
+# some definitions for the physics, in file 'callphys.def'
+INCLUDEDEF=callphys.def
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run_mcd_3days
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run_mcd_3days	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/run_mcd_3days	(revision 214)
@@ -0,0 +1,60 @@
+#! /bin/bash
+
+#
+# - please link or copy start0.nc and startfi.nc (or start_archive.nc)
+# - please set 3 in the first parameter of run.def
+# - please set 0 in the first line of num_run 
+# - please set daylim below (669 for the whole year)
+#
+
+##########
+daylim=669
+##########
+
+rm -rf tmpsave
+mkdir tmpsave
+
+runnum=0 ## 0 whatever, this is just the way one counts successive runs
+runnum_obsolete=0
+daynum=0 #61 #0 
+
+while [ ${daynum} -lt ${daylim} ]  
+do
+	# run GCM during three days
+	echo run starting at day ${daynum}
+	echo '****'
+	time run0 
+
+	# save generated starting files 
+	echo saving files
+	echo '****'
+	runnum=$(expr ${runnum} + 1)
+        save=$(expr ${daynum} + 2)
+	cp start${runnum}.nc tmpsave/start${save}.nc
+        cp startfi${runnum}.nc tmpsave/startfi${save}.nc
+	mv lrun${runnum} tmpsave/lrun_for_start${save}
+
+	# increment day to day+3 for while loop assessement
+	daynum=$(expr ${daynum} + 3)
+
+        # get rid of what was generated two steps before
+	echo cleaning
+	echo '****'
+	if [ ${runnum} -gt 1 ]	#because we want to keep start0.nc and startfi0.nc
+	then
+	runnum_obsolete=$(expr ${runnum} - 1)
+	rm -rf start${runnum_obsolete}.nc
+	rm -rf startfi${runnum_obsolete}.nc
+	rm -rf diagfi${runnum_obsolete}.nc
+	fi
+done
+
+# clean the last files
+rm -rf start${runnum}.nc
+rm -rf startfi${runnum}.nc
+rm -rf diagfi${runnum}.nc
+
+# end
+echo -- end --
+echo please check results in file tmpsave
+echo if you want to keep the results, rename this folder
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/traceur.def
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/traceur.def	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/traceur.def	(revision 214)
@@ -0,0 +1,2 @@
+1
+co2
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/what_I_need.info
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/what_I_need.info	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/myGCM/what_I_need.info	(revision 214)
@@ -0,0 +1,9 @@
+total 24K
+-rw-r--r--  1 aslmd lmdjus 3.4K Mar 11 11:50 callphys.def
+-rw-r--r--  1 aslmd lmdjus 3.4K Mar 11 11:45 callphys.def~
+lrwxrwxrwx  1 aslmd lmdjus   18 Mar 11 11:43 gcm.e -> ../LMDZ.MARS/gcm.e
+lrwxrwxrwx  1 aslmd lmdjus   23 Mar 11 11:43 newstart.e -> ../LMDZ.MARS/newstart.e
+-rw-r--r--  1 aslmd lmdjus 3.2K Mar 11 11:49 run.def
+lrwxrwxrwx  1 aslmd lmdjus   57 Mar 11 11:44 start_archive.nc -> /u/forget/WWW/datagcm/Starts/start_archive_64x48x25_t2.nc
+-rw-r--r--  1 aslmd lmdjus    0 Mar 11 11:53 what_I_need
+lrwxrwxrwx  1 aslmd lmdjus   30 Mar 11 11:52 z2sig.def -> ../LMDZ.MARS/deftank/z2sig.def
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F	(revision 214)
@@ -0,0 +1,1485 @@
+      SUBROUTINE physiq(ngrid,nlayer,nq,
+     $            firstcall,lastcall,
+     $            pday,ptime,ptimestep,
+     $            pplev,pplay,pphi,
+     $            pu,pv,pt,pq,
+     $            pw,
+     $            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
+
+
+      IMPLICIT NONE
+c=======================================================================
+c
+c   subject:
+c   --------
+c
+c   Organisation of the physical parametrisations of the LMD 
+c   martian atmospheric general circulation model.
+c
+c   The GCM can be run without or with tracer transport
+c   depending on the value of Logical "tracer" in file  "callphys.def"
+c   Tracers may be water vapor, ice OR chemical species OR dust particles
+c
+c   SEE comments in initracer.F about numbering of tracer species...
+c
+c   It includes:
+c
+c      1. Initialization:
+c      1.1 First call initializations
+c      1.2 Initialization for every call to physiq
+c      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
+c      2. Compute radiative transfer tendencies
+c         (longwave and shortwave) for CO2 and aerosols.
+c      3. Gravity wave and subgrid scale topography drag :
+c      4. Vertical diffusion (turbulent mixing):
+c      5. Convective adjustment
+c      6. Condensation and sublimation of carbon dioxide.
+c      7.  TRACERS :
+c       7a. water and water ice
+c       7b. call for photochemistry when tracers are chemical species
+c       7c. other scheme for tracer (dust) transport (lifting, sedimentation)
+c       7d. updates (CO2 pressure variations, surface budget)
+c      8. Contribution to tendencies due to thermosphere
+c      9. Surface and sub-surface temperature calculations
+c     10. Write outputs :
+c           - "startfi", "histfi" (if it's time)
+c           - Saving statistics (if "callstats = .true.")
+c           - Dumping eof (if "calleofdump = .true.")
+c           - Output any needed variables in "diagfi" 
+c     11. Diagnostic: mass conservation of tracers
+c 
+c   author: 
+c   ------- 
+c           Frederic Hourdin	15/10/93
+c           Francois Forget		1994
+c           Christophe Hourdin	02/1997 
+c           Subroutine completly rewritten by F.Forget (01/2000)
+c           Introduction of the photochemical module: S. Lebonnois (11/2002)
+c           Introduction of the thermosphere module: M. Angelats i Coll (2002)
+c           Water ice clouds: Franck Montmessin (update 06/2003)
+c           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
+c             Nb: See callradite.F for more information.
+c           
+c   arguments:
+c   ----------
+c
+c   input:
+c   ------
+c    ecri                  period (in dynamical timestep) to write output
+c    ngrid                 Size of the horizontal grid.
+c                          All internal loops are performed on that grid.
+c    nlayer                Number of vertical layers.
+c    nq                    Number of advected fields
+c    firstcall             True at the first call
+c    lastcall              True at the last call
+c    pday                  Number of days counted from the North. Spring
+c                          equinoxe.
+c    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
+c    ptimestep             timestep (s)
+c    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
+c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
+c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
+c    pu(ngrid,nlayer)      u component of the wind (ms-1)
+c    pv(ngrid,nlayer)      v component of the wind (ms-1)
+c    pt(ngrid,nlayer)      Temperature (K)
+c    pq(ngrid,nlayer,nq)   Advected fields
+c    pudyn(ngrid,nlayer)    \ 
+c    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
+c    ptdyn(ngrid,nlayer)     / corresponding variables
+c    pqdyn(ngrid,nlayer,nq) /
+c    pw(ngrid,?)           vertical velocity
+c
+c   output:
+c   -------
+c
+c    pdu(ngrid,nlayermx)        \
+c    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
+c    pdt(ngrid,nlayermx)         /  variables due to physical processes.
+c    pdq(ngrid,nlayermx,nqmx)   /
+c    pdpsrf(ngrid)             /
+c    tracerdyn                 call tracer in dynamical part of GCM ?
+
+c
+c=======================================================================
+c
+c    0.  Declarations :
+c    ------------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comgeomfi.h"
+#include "surfdat.h"
+#include "comsoil.h"
+#include "comdiurn.h"
+#include "callkeys.h"
+#include "comcstfi.h"
+#include "planete.h"
+#include "comsaison.h"
+#include "control.h"
+#include "dimradmars.h"
+#include "comg1d.h"
+#include "tracer.h"
+#include "nlteparams.h"
+
+#include "chimiedata.h"
+#include "watercap.h"
+#include "param.h"
+#include "param_v3.h"
+#include "conc.h"
+
+#include "netcdf.inc"
+
+
+
+c Arguments :
+c -----------
+
+c   inputs:
+c   -------
+      INTEGER ngrid,nlayer,nq
+      REAL ptimestep
+      REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
+      REAL pphi(ngridmx,nlayer)
+      REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
+      REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
+      REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
+      REAL zh(ngridmx,nlayermx)      ! potential temperature (K)
+      LOGICAL firstcall,lastcall
+
+      REAL pday
+      REAL ptime 
+      logical tracerdyn
+
+c   outputs:
+c   --------
+c     physical tendencies
+      REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
+      REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
+      REAL pdpsrf(ngridmx) ! surface pressure tendency
+
+
+c Local saved variables:
+c ----------------------
+c     aerosol (dust or ice) extinction optical depth  at reference wavelength 
+c     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
+c      aerosol optical properties  :
+      REAL aerosol(ngridmx,nlayermx,naerkind)
+
+      INTEGER day_ini  ! Initial date of the run (sol since Ls=0) 
+      INTEGER icount     ! counter of calls to physiq during the run.
+      REAL tsurf(ngridmx)            ! Surface temperature (K)
+      REAL tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
+      REAL co2ice(ngridmx)           ! co2 ice surface layer (kg.m-2)  
+      REAL albedo(ngridmx,2)         ! Surface albedo in each solar band
+      REAL emis(ngridmx)             ! Thermal IR surface emissivity
+      REAL dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
+      REAL fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
+      REAL fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
+      REAL capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
+      REAL fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
+      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
+      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
+      INTEGER ig_vl1                 ! Grid Point near VL1   (for diagnostic) 
+
+c     Variables used by the water ice microphysical scheme:
+      REAL rice(ngridmx,nlayermx)    ! Water ice geometric mean radius (m)
+      REAL nuice(ngridmx,nlayermx)   ! Estimated effective variance
+                                     !   of the size distribution
+c     Albedo of deposited surface ice
+      REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
+
+      SAVE day_ini, icount
+      SAVE aerosol, tsurf,tsoil
+      SAVE co2ice,albedo,emis, q2
+      SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
+      SAVE ig_vl1
+
+      REAL stephan   
+      DATA stephan/5.67e-08/  ! Stephan Boltzman constant
+      SAVE stephan
+
+c Local variables :
+c -----------------
+
+      REAL CBRT
+      EXTERNAL CBRT
+
+      CHARACTER*80 fichier 
+      INTEGER l,ig,ierr,igout,iq,i, tapphys
+
+      REAL fluxsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2)
+      REAL fluxsurf_sw(ngridmx,2)    !incident SW (solar) surface flux (W.m-2)
+      REAL fluxtop_lw(ngridmx)       !Outgoing LW (IR) flux to space (W.m-2)
+      REAL fluxtop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
+      REAL tauref(ngridmx)           ! Reference column optical depth at 700 Pa
+                                     ! (used if active=F) 
+      REAL tau(ngridmx,naerkind)     ! Column dust optical depth at each point
+      REAL zls                       !  solar longitude (rad)
+      REAL zday                      ! date (time since Ls=0, in martian days)
+      REAL zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
+      REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
+      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
+
+c     Tendancies due to various processes:
+      REAL dqsurf(ngridmx,nqmx)
+      REAL zdtlw(ngridmx,nlayermx)     ! (K/s)
+      REAL zdtsw(ngridmx,nlayermx)     ! (K/s)
+      REAL cldtlw(ngridmx,nlayermx)     ! (K/s) LW heating rate for clear area
+      REAL cldtsw(ngridmx,nlayermx)     ! (K/s) SW heating rate for clear area
+      REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
+      REAL zdtnlte(ngridmx,nlayermx)   ! (K/s)
+      REAL zdtsurf(ngridmx)            ! (K/s)
+      REAL zdtcloud(ngridmx,nlayermx)
+      REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
+      REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhadj(ngridmx,nlayermx)                           ! (K/s)
+      REAL zdtgw(ngridmx,nlayermx)                            ! (K/s)
+      REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
+      REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
+      REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
+
+      REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
+      REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
+      REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
+      REAL zdqadj(ngridmx,nlayermx,nqmx)
+      REAL zdqc(ngridmx,nlayermx,nqmx)
+      REAL zdqcloud(ngridmx,nlayermx,nqmx)
+      REAL zdqscloud(ngridmx,nqmx)
+      REAL zdqchim(ngridmx,nlayermx,nqmx)
+      REAL zdqschim(ngridmx,nqmx)
+
+      REAL zdteuv(ngridmx,nlayermx)    ! (K/s)
+      REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
+      REAL zdumolvis(ngridmx,nlayermx)
+      REAL zdvmolvis(ngridmx,nlayermx)
+      real zdqmoldiff(ngridmx,nlayermx,nqmx)
+
+c     Local variable for local intermediate calcul:
+      REAL zflubid(ngridmx)
+      REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
+      REAL zdum1(ngridmx,nlayermx)
+      REAL zdum2(ngridmx,nlayermx)
+      REAL ztim1,ztim2,ztim3, z1,z2
+      REAL ztime_fin
+      REAL zdh(ngridmx,nlayermx)
+      INTEGER length
+      PARAMETER (length=100)
+
+c local variables only used for diagnostic (output in file "diagfi" or "stats")
+c -----------------------------------------------------------------------------
+      REAL ps(ngridmx), zt(ngridmx,nlayermx)
+      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
+      REAL zq(ngridmx,nlayermx,nqmx)
+      REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
+      character*2 str2
+      character*5 str5
+      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
+      REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
+                                   !   (particules kg-1)
+      real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
+      real qtot1,qtot2 ! total aerosol mass
+      integer igmin, lmin
+      logical tdiag
+
+      real co2col(ngridmx)        ! CO2 column
+      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
+      REAL zstress(ngrid), cd
+      real hco2(nqmx),tmean, zlocal(nlayermx)
+      real rho(ngridmx,nlayermx)  ! density
+      real vmr(ngridmx,nlayermx)  ! volume mixing ratio
+      REAL mtot(ngridmx)          ! Total mass of water vapor (kg/m2)
+      REAL icetot(ngridmx)        ! Total mass of water ice (kg/m2)
+      REAL rave(ngridmx)          ! Mean water ice effective radius (m)
+      REAL opTES(ngridmx,nlayermx)! abs optical depth at 825 cm-1
+      REAL tauTES(ngridmx)        ! column optical depth at 825 cm-1
+      REAL Qabsice                ! Water ice absorption coefficient
+
+
+      REAL time_phys
+
+c=======================================================================
+
+c 1. Initialisation:
+c -----------------
+
+c  1.1   Initialisation only at first call
+c  ---------------------------------------
+      IF (firstcall) THEN
+
+c        variables set to 0
+c        ~~~~~~~~~~~~~~~~~~
+         call zerophys(ngrid*nlayer*naerkind,aerosol)
+         call zerophys(ngrid*nlayer,dtrad)
+         call zerophys(ngrid,fluxrad)
+
+c        read startfi 
+c        ~~~~~~~~~~~~
+
+! Read netcdf initial physical parameters.
+         CALL phyetat0 ("startfi.nc",0,0,
+     &         nsoilmx,nq,
+     &         day_ini,time_phys,
+     &         tsurf,tsoil,emis,q2,qsurf,co2ice)
+
+         if (pday.ne.day_ini) then
+           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
+     &                "physics and dynamics"
+           write(*,*) "dynamics day: ",pday
+           write(*,*) "physics day:  ",day_ini
+           stop
+         endif
+
+         write (*,*) 'In physiq day_ini =', day_ini
+
+c        Initialize albedo and orbital calculation
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         CALL surfini(ngrid,co2ice,qsurf,albedo)
+         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+
+c        initialize soil 
+c        ~~~~~~~~~~~~~~~
+         IF (callsoil) THEN
+            CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
+     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+         ELSE
+            PRINT*,
+     &     'PHYSIQ WARNING! Thermal conduction in the soil turned off'
+            DO ig=1,ngrid
+               capcal(ig)=1.e5
+               fluxgrd(ig)=0.
+            ENDDO
+         ENDIF
+         icount=1
+
+
+c        initialize tracers
+c        ~~~~~~~~~~~~~~~~~~
+         tracerdyn=tracer
+         IF (tracer) THEN
+            CALL initracer(qsurf,co2ice)
+         ENDIF  ! end tracer
+
+c        Determining gridpoint near Viking Lander 1 (used for diagnostic only)
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if(ngrid.ne.1) then
+           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.)
+           write(*,*) 'Viking Lander 1 GCM point: lat,lon',
+     &              lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi
+         end if 
+
+c        Initialize thermospheric parameters
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if (callthermos) call param_read
+
+c        Initialize R and Cp as constant
+
+         if (.not.callthermos .and. .not.photochem) then
+                 do l=1,nlayermx
+                  do ig=1,ngridmx
+                   rnew(ig,l)=r
+                   cpnew(ig,l)=cpp
+                   mmean(ig,l)=mugaz
+                   enddo
+                  enddo  
+         endif         
+
+        IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
+          write(*,*)"physiq: water_param Surface ice alb:",alb_surfice
+        ENDIF
+                   
+      ENDIF        !  (end of "if firstcall")
+
+c ---------------------------------------------------
+c 1.2   Initializations done at every physical timestep:
+c ---------------------------------------------------
+c
+      IF (ngrid.NE.ngridmx) THEN
+         PRINT*,'STOP in PHYSIQ'
+         PRINT*,'Probleme de dimensions :'
+         PRINT*,'ngrid     = ',ngrid
+         PRINT*,'ngridmx   = ',ngridmx
+         STOP
+      ENDIF
+
+c     Initialize various variables
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      call zerophys(ngrid*nlayer, pdv)
+      call zerophys(ngrid*nlayer, pdu)
+      call zerophys(ngrid*nlayer, pdt)
+      call zerophys(ngrid*nlayer*nq, pdq)
+      call zerophys(ngrid, pdpsrf)
+      call zerophys(ngrid, zflubid)
+      call zerophys(ngrid, zdtsurf)
+      call zerophys(ngrid*nq, dqsurf)
+      igout=ngrid/2+1 
+
+
+      zday=pday+ptime ! compute time, in sols (and fraction thereof)
+
+c     Compute Solar Longitude (Ls) :
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (season) then
+         call solarlong(zday,zls)
+      else
+         call solarlong(float(day_ini),zls)
+      end if
+
+c     Compute geopotential at interlayers
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c     ponderation des altitudes au niveau des couches en dp/p
+
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zzlay(ig,l)=pphi(ig,l)/g
+         ENDDO
+      ENDDO
+      DO ig=1,ngrid
+         zzlev(ig,1)=0.
+         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+      ENDDO
+      DO l=2,nlayer
+         DO ig=1,ngrid
+            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
+            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
+            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+         ENDDO
+      ENDDO
+
+
+!     Potential temperature calculation not the same in physiq and dynamic
+
+c     Compute potential temperature
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO l=1,nlayer
+         DO ig=1,ngrid 
+            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
+            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c    1.2.5 Compute mean mass, cp, and R
+c    --------------------------------
+
+      if(photochem.or.callthermos) then
+         call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
+      endif
+c-----------------------------------------------------------------------
+c    2. Compute radiative tendencies :
+c------------------------------------
+
+
+      IF (callrad) THEN
+         IF( MOD(icount-1,iradia).EQ.0) THEN
+
+c          Local Solar zenith angle
+c          ~~~~~~~~~~~~~~~~~~~~~~~~
+           CALL orbite(zls,dist_sol,declin)
+
+           IF(diurnal) THEN
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
+     s         ztim1,ztim2,ztim3, mu0,fract)
+
+           ELSE
+               CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
+           ENDIF
+
+c          NLTE cooling from CO2 emission
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+           IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
+
+c          Find number of layers for LTE radiation calculations
+           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
+     &          CALL nlthermeq(ngrid,nlayer,pplev,pplay)
+
+c          Note: Dustopacity.F has been transferred to callradite.F
+         
+c          Call main radiative transfer scheme
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c          Transfer through CO2 (except NIR CO2 absorption)
+c            and aerosols (dust and water ice)
+
+c          Radiative transfer
+c          ------------------
+           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
+     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
+     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
+     &     tauref,tau,aerosol,ccn,rdust,rice,nuice)
+
+c          CO2 near infrared absorption
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           call zerophys(ngrid*nlayer,zdtnirco2)
+           if (callnirco2) then
+              call nirco2abs (ngrid,nlayer,pplay,dist_sol,
+     .                       mu0,fract,declin, zdtnirco2)
+           endif
+
+c          Radiative flux from the sky absorbed by the surface (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           DO ig=1,ngrid
+               fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
+     $         +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
+     $         +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
+           ENDDO
+
+
+c          Net atmospheric radiative heating rate (K.s-1)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           IF(callnlte) THEN
+              CALL blendrad(ngrid, nlayer, pplay,
+     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
+           ELSE
+              DO l=1,nlayer
+                 DO ig=1,ngrid
+                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
+     &                          +zdtnirco2(ig,l)
+                  ENDDO
+              ENDDO
+           ENDIF
+
+
+
+        ENDIF ! of if(mod(icount-1,iradia).eq.0)
+
+c    Transformation of the radiative tendencies:
+c    -------------------------------------------
+
+c          Net radiative surface flux (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+           DO ig=1,ngrid
+               zplanck(ig)=tsurf(ig)*tsurf(ig)
+               zplanck(ig)=emis(ig)*
+     $         stephan*zplanck(ig)*zplanck(ig)
+               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
+           ENDDO
+
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
+            ENDDO
+         ENDDO
+
+      ENDIF ! of IF (callrad)
+
+c-----------------------------------------------------------------------
+c    3. Gravity wave and subgrid scale topography drag :
+c    -------------------------------------------------
+
+
+      IF(calllott)THEN
+
+        CALL calldrag_noro(ngrid,nlayer,ptimestep,
+     &                 pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
+ 
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
+          ENDDO
+        ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c    4. Vertical diffusion (turbulent mixing):
+c    -----------------------------------------
+c
+      IF (calldifv) THEN
+
+
+         DO ig=1,ngrid
+            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
+         ENDDO
+
+         CALL zerophys(ngrid*nlayer,zdum1)
+         CALL zerophys(ngrid*nlayer,zdum2)
+         do l=1,nlayer
+            do ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            enddo
+         enddo
+         
+c        Calling vdif (Martian version WITH CO2 condensation)
+         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
+     $        ptimestep,capcal,lwrite,
+     $        pplay,pplev,zzlay,zzlev,z0,
+     $        pu,pv,zh,pq,tsurf,emis,qsurf,
+     $        zdum1,zdum2,zdh,pdq,zflubid,
+     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
+     &        zdqdif,zdqsdif)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
+               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
+
+               zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
+
+            ENDDO
+         ENDDO
+
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
+         ENDDO
+
+         if (tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+           DO iq=1, nq
+              DO ig=1,ngrid
+                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+              ENDDO
+           ENDDO
+         end if ! of if (tracer)
+
+      ELSE    
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+
+     s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+         ENDDO
+      ENDIF ! of IF (calldifv)
+
+
+c-----------------------------------------------------------------------
+c   5. Dry convective adjustment:
+c   -----------------------------
+
+      IF(calladj) THEN
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            ENDDO
+         ENDDO
+         CALL zerophys(ngrid*nlayer,zduadj)
+         CALL zerophys(ngrid*nlayer,zdvadj)
+         CALL zerophys(ngrid*nlayer,zdhadj)
+
+         CALL convadj(ngrid,nlayer,nq,ptimestep,
+     $                pplay,pplev,zpopsk,
+     $                pu,pv,zh,pq,
+     $                pdu,pdv,zdh,pdq,
+     $                zduadj,zdvadj,zdhadj,
+     $                zdqadj)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
+               pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
+
+               zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
+            ENDDO
+         ENDDO
+
+         if(tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         end if
+      ENDIF ! of IF(calladj)
+
+c-----------------------------------------------------------------------
+c   6. Carbon dioxide condensation-sublimation:
+c   -------------------------------------------
+
+      IF (callcond) THEN
+         CALL newcondens(ngrid,nlayer,nq,ptimestep,
+     $              capcal,pplay,pplev,tsurf,pt,
+     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
+     $              co2ice,albedo,emis,
+     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
+     $	            fluxsurf_sw,zls) 
+
+         DO l=1,nlayer
+           DO ig=1,ngrid
+             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
+             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
+             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
+           ENDDO
+         ENDDO
+         DO ig=1,ngrid
+            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
+         ENDDO
+
+         IF (tracer) THEN
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         ENDIF ! of IF (tracer)
+
+      ENDIF  ! of IF (callcond)
+
+c-----------------------------------------------------------------------
+c   7. Specific parameterizations for tracers 
+c:   -----------------------------------------
+
+      if (tracer) then 
+
+c   7a. Water and ice
+c     ---------------
+
+c        ---------------------------------------
+c        Water ice condensation in the atmosphere
+c        ----------------------------------------
+         IF (water) THEN
+
+           call watercloud(ngrid,nlayer,ptimestep,
+     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
+     &                pq,pdq,zdqcloud,zdqscloud,zdtcloud,
+     &                nq,naerkind,tau,
+     &                ccn,rdust,rice,nuice)
+           if (activice) then
+c Temperature variation due to latent heat release
+           DO l=1,nlayer
+             DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
+             ENDDO
+           ENDDO
+           endif
+
+! increment water vapour and ice atmospheric tracers tendencies
+           IF (water) THEN
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                 pdq(ig,l,igcm_h2o_vap)=pdq(ig,l,igcm_h2o_vap)+
+     &                                   zdqcloud(ig,l,igcm_h2o_vap)
+                 pdq(ig,l,igcm_h2o_ice)=pdq(ig,l,igcm_h2o_ice)+
+     &                                   zdqcloud(ig,l,igcm_h2o_ice)
+               ENDDO
+             ENDDO
+           ENDIF ! of IF (water) THEN
+! Increment water ice surface tracer tendency
+         DO ig=1,ngrid
+           dqsurf(ig,igcm_h2o_ice)=dqsurf(ig,igcm_h2o_ice)+
+     &                               zdqscloud(ig,igcm_h2o_ice)
+         ENDDO
+         
+         END IF  ! of IF (water)
+
+c   7b. Chemical species
+c     ------------------
+
+c        --------------
+c        photochemistry :
+c        --------------
+         IF (photochem .or. thermochem) then
+          call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
+     &      zzlay,zday,pq,pdq,rice,
+     &      zdqchim,zdqschim,zdqcloud,zdqscloud)
+!NB: Photochemistry includes condensation of H2O2
+
+           ! 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)
+
+c   7c. Aerosol particles
+c     -------------------
+
+c        ----------
+c        Dust devil :
+c        ----------
+         IF(callddevil) then 
+           call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
+     &                zdqdev,zdqsdev)
+ 
+           if (dustbin.ge.1) then
+              do iq=1,nq
+                 DO l=1,nlayer
+                    DO ig=1,ngrid
+                       pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
+                    ENDDO
+                 ENDDO
+              enddo
+              do iq=1,nq
+                 DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
+                 ENDDO
+              enddo
+           endif  ! of if (dustbin.ge.1)
+
+         END IF ! of IF (callddevil)
+
+c        ------------- 
+c        Sedimentation :   acts also on water ice
+c        ------------- 
+         IF (sedimentation) THEN 
+           !call zerophys(ngrid*nlayer*nq, zdqsed)
+           zdqsed(1:ngrid,1:nlayer,1:nq)=0
+           !call zerophys(ngrid*nq, zdqssed)
+           zdqssed(1:ngrid,1:nq)=0
+
+           call callsedim(ngrid,nlayer, ptimestep,
+     &                pplev,zzlev, pt, rdust, rice,
+     &                pq, pdq, zdqsed, zdqssed,nq)
+
+           DO iq=1, nq
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
+               ENDDO
+             ENDDO
+           ENDDO
+           DO iq=1, nq
+             DO ig=1,ngrid
+                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
+             ENDDO
+           ENDDO
+         END IF   ! of IF (sedimentation)
+
+c   7d. Updates
+c     ---------
+
+        DO iq=1, nq
+          DO ig=1,ngrid
+
+c       ---------------------------------
+c       Updating tracer budget on surface
+c       ---------------------------------
+            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
+
+          ENDDO  ! (ig)
+        ENDDO    ! (iq)
+
+      endif !  of if (tracer) 
+
+
+c-----------------------------------------------------------------------
+c   8. THERMOSPHERE CALCULATION
+c-----------------------------------------------------------------------
+
+      if (callthermos) then
+        call thermosphere(pplev,pplay,dist_sol,
+     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
+     &     pt,pq,pu,pv,pdt,pdq,
+     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
+
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
+     &                         +zdteuv(ig,l)
+            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
+            DO iq=1, nq
+              pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
+            ENDDO
+          ENDDO
+        ENDDO
+
+      endif ! of if (callthermos)
+
+c-----------------------------------------------------------------------
+c   9. Surface  and sub-surface soil temperature
+c-----------------------------------------------------------------------
+c
+c
+c   9.1 Increment Surface temperature:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      DO ig=1,ngrid
+         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
+      ENDDO
+
+c  Prescribe a cold trap at south pole (except at high obliquity !!)
+c  Temperature at the surface is set there to be the temperature
+c  corresponding to equilibrium temperature between phases of CO2
+
+      IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
+         if (caps.and.(obliquit.lt.27.)) then
+           ! NB: Updated surface pressure, at grid point 'ngrid', is
+           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
+           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
+     &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
+         endif
+c       -------------------------------------------------------------
+c       Change of surface albedo (set to 0.4) in case of ground frost
+c       everywhere except on the north permanent cap and in regions
+c       covered by dry ice. 
+c              ALWAYS PLACE these lines after newcondens !!!
+c       -------------------------------------------------------------
+         do ig=1,ngrid
+           if ((co2ice(ig).eq.0).and.
+     &        (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
+              albedo(ig,1) = alb_surfice
+              albedo(ig,2) = alb_surfice
+           endif
+         enddo  ! of do ig=1,ngrid
+      ENDIF  ! of IF (tracer.AND.water.AND.(ngridmx.NE.1))
+
+c
+c   9.2 Compute soil temperatures and subsurface heat flux:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (callsoil) THEN
+         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
+     &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+      ENDIF
+
+c-----------------------------------------------------------------------
+c  10. Write output files
+c  ----------------------
+
+c    -------------------------------
+c    Dynamical fields incrementation
+c    -------------------------------
+c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
+      ! temperature, zonal and meridional wind
+      DO l=1,nlayer
+        DO ig=1,ngrid
+          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
+        ENDDO
+      ENDDO
+
+      ! tracers
+      DO iq=1, nq
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
+          ENDDO
+        ENDDO
+      ENDDO
+
+      ! surface pressure
+      DO ig=1,ngrid
+          ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
+      ENDDO
+
+      ! pressure
+      DO l=1,nlayer
+        DO ig=1,ngrid
+             zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
+             zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
+        ENDDO
+      ENDDO
+
+      ! Density 
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
+         ENDDO
+      ENDDO
+
+c    Compute surface stress : (NB: z0 is a common in planete.h)
+c     DO ig=1,ngrid
+c        cd = (0.4/log(zzlay(ig,1)/z0))**2
+c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
+c     ENDDO
+
+c     Sum of fluxes in solar spectral bands (for output only)
+      DO ig=1,ngrid
+	     fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
+	     fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
+      ENDDO
+c ******* TEST ******************************************************
+      ztim1 = 999
+      DO l=1,nlayer
+        DO ig=1,ngrid
+           if (pt(ig,l).lt.ztim1) then
+               ztim1 = pt(ig,l)
+               igmin = ig
+               lmin = l 
+           end if
+        ENDDO
+      ENDDO
+      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then	
+        write(*,*) 'PHYSIQ: stability WARNING :'
+        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
+     &              'ig l =', igmin, lmin
+      end if
+c *******************************************************************
+
+c     ---------------------
+c     Outputs to the screen 
+c     ---------------------
+
+      IF (lwrite) THEN
+         PRINT*,'Global diagnostics for the physics'
+         PRINT*,'Variables and their increments x and dx/dt * dt'
+         WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
+         WRITE(*,'(2f10.5,2f15.5)')
+     s   tsurf(igout),zdtsurf(igout)*ptimestep,
+     s   pplev(igout,1),pdpsrf(igout)*ptimestep
+         WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
+         WRITE(*,'(i4,6f10.5)') (l,
+     s   pu(igout,l),pdu(igout,l)*ptimestep,
+     s   pv(igout,l),pdv(igout,l)*ptimestep,
+     s   pt(igout,l),pdt(igout,l)*ptimestep,
+     s   l=1,nlayer)
+      ENDIF ! of IF (lwrite)
+
+      IF (ngrid.NE.1) THEN
+         print*,'Ls =',zls*180./pi,
+     &   ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
+     &   ' tau(Viking1) =',tau(ig_vl1,1)
+
+
+c        -------------------------------------------------------------------
+c        Writing NetCDF file  "RESTARTFI" at the end of the run
+c        -------------------------------------------------------------------
+c        Note: 'restartfi' is stored just before dynamics are stored
+c              in 'restart'. Between now and the writting of 'restart',
+c              there will have been the itau=itau+1 instruction and
+c              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
+c              thus we store for time=time+dtvr
+
+         IF(lastcall) THEN
+            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) 
+            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
+     .              ptimestep,pday,
+     .              ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
+     .              area,albedodat,inertiedat,zmea,zstd,zsig,
+     .              zgam,zthe)
+         ENDIF
+
+c        -------------------------------------------------------------------
+c        Calculation of diagnostic variables written in both stats and
+c          diagfi files
+c        -------------------------------------------------------------------
+
+         if (tracer) then
+           if (water) then
+
+             call zerophys(ngrid,mtot)
+             call zerophys(ngrid,icetot)
+             call zerophys(ngrid,rave)
+             call zerophys(ngrid,tauTES)
+             do ig=1,ngrid 
+               do l=1,nlayermx
+                 mtot(ig) = mtot(ig) + 
+     &                      zq(ig,l,igcm_h2o_vap) * 
+     &                      (pplev(ig,l) - pplev(ig,l+1)) / g
+                 icetot(ig) = icetot(ig) + 
+     &                        zq(ig,l,igcm_h2o_ice) * 
+     &                        (pplev(ig,l) - pplev(ig,l+1)) / g
+                 rave(ig) = rave(ig) + 
+     &                      zq(ig,l,igcm_h2o_ice) *
+     &                      (pplev(ig,l) - pplev(ig,l+1)) / g * 
+     &                      rice(ig,l) * (1.+nuice_ref)
+c                Computing abs optical depth at 825 cm-1 in each
+c                  layer to simulate NEW TES retrieval
+                 Qabsice = min(
+     &             max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2
+     &                        )
+                 opTES(ig,l)= 0.75 * Qabsice * 
+     &             zq(ig,l,igcm_h2o_ice) *
+     &             (pplev(ig,l) - pplev(ig,l+1)) / g
+     &             / (rho_ice * rice(ig,l) * (1.+nuice_ref))
+                 tauTES(ig)=tauTES(ig)+ opTES(ig,l) 
+               enddo
+               rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
+               if (icetot(ig)*1e3.lt.0.01) rave(ig)=0.
+             enddo
+
+           endif ! of if (water)
+         endif ! of if (tracer)
+
+c        -----------------------------------------------------------------
+c        WSTATS: Saving statistics
+c        -----------------------------------------------------------------
+c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+c        which can later be used to make the statistic files of the run:
+c        "stats")          only possible in 3D runs !
+
+         
+         IF (callstats) THEN
+
+           call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
+           call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+           call wstats(ngrid,"co2ice","CO2 ice cover",
+     &                "kg.m-2",2,co2ice)
+           call wstats(ngrid,"fluxsurf_lw",
+     &                "Thermal IR radiative flux to surface","W.m-2",2,
+     &                fluxsurf_lw)
+           call wstats(ngrid,"fluxsurf_sw",
+     &                "Solar radiative flux to surface","W.m-2",2,
+     &                fluxsurf_sw_tot)
+           call wstats(ngrid,"fluxtop_lw",
+     &                "Thermal IR radiative flux to space","W.m-2",2,
+     &                fluxtop_lw)
+           call wstats(ngrid,"fluxtop_sw",
+     &                "Solar radiative flux to space","W.m-2",2,
+     &                fluxtop_sw_tot)
+           call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
+           call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
+           call wstats(ngrid,"v","Meridional (North-South) wind",
+     &                "m.s-1",3,zv)
+           call wstats(ngrid,"w","Vertical (down-up) wind",
+     &                "m.s-1",3,pw)
+           call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
+           call wstats(ngrid,"pressure","Pressure","Pa",3,pplay)
+c          call wstats(ngrid,"q2",
+c    &                "Boundary layer eddy kinetic energy",
+c    &                "m2.s-2",3,q2)
+c          call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
+c    &                emis)
+c          call wstats(ngrid,"ssurf","Surface stress","N.m-2",
+c    &                2,zstress)
+c          call wstats(ngrid,"sw_htrt","sw heat.rate",
+c    &                 "W.m-2",3,zdtsw)
+c          call wstats(ngrid,"lw_htrt","lw heat.rate",
+c    &                 "W.m-2",3,zdtlw)
+
+           if (tracer) then
+             if (water) then
+               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
+     &                  *mugaz/mmol(igcm_h2o_vap)
+               call wstats(ngrid,"vmr_h2ovapor",
+     &                    "H2O vapor volume mixing ratio","mol/mol",
+     &                    3,vmr)
+               vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+     &                  *mugaz/mmol(igcm_h2o_ice)
+               call wstats(ngrid,"vmr_h2oice",
+     &                    "H2O ice volume mixing ratio","mol/mol",
+     &                    3,vmr)
+               call wstats(ngrid,"h2o_ice_s",
+     &                    "surface h2o_ice","kg/m2",
+     &                    2,qsurf(1,igcm_h2o_ice))
+
+               call wstats(ngrid,"mtot",
+     &                    "total mass of water vapor","kg/m2",
+     &                    2,mtot)
+               call wstats(ngrid,"icetot",
+     &                    "total mass of water ice","kg/m2",
+     &                    2,icetot)
+               call wstats(ngrid,"reffice",
+     &                    "Mean reff","m",
+     &                    2,rave)
+c              call wstats(ngrid,"rice",
+c    &                    "Ice particle size","m",
+c    &                    3,rice)
+c              If activice is true, tauTES is computed in aeropacity.F;
+               if (.not.activice) then
+                 call wstats(ngrid,"tauTESap",
+     &                      "tau abs 825 cm-1","",
+     &                      2,tauTES)
+               endif
+
+             endif ! of if (water)
+
+             if (thermochem.or.photochem) then
+                do iq=1,nq
+                   if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
+     .                (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
+     .                (noms(iq).eq."h2").or.
+     .                (noms(iq).eq."o3")) then
+                        do l=1,nlayer
+                          do ig=1,ngrid
+                            vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
+                          end do
+                        end do
+                        call wstats(ngrid,"vmr_"//trim(noms(iq)),
+     .                     "Volume mixing ratio","mol/mol",3,vmr)
+                   endif
+                enddo
+             endif ! of if (thermochem.or.photochem)
+
+           endif ! of if (tracer)
+
+           IF(lastcall) THEN
+             write (*,*) "Writing stats..."
+             call mkstats(ierr)
+           ENDIF
+
+         ENDIF !if callstats
+
+c        (Store EOF for Mars Climate database software)
+         IF (calleofdump) THEN
+            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
+         ENDIF
+
+c        ==========================================================
+c        WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
+c          any variable for diagnostic (output with period
+c          "ecritphy", set in "run.def")
+c        ==========================================================
+c        WRITEDIAGFI can ALSO be called from any other subroutines
+c        for any variables !!
+        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     &                  emis)
+         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     &                  tsurf)
+         call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     &                  co2ice)
+c         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
+c     &                  fluxsurf_lw)
+c         call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
+c     &                  fluxsurf_sw_tot)
+c         call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
+c     &                  fluxtop_lw)
+c         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
+c     &                  fluxtop_sw_tot)
+         call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+c        call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
+c         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
+c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
+c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
+c        call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
+c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
+c    &                  zstress)
+c        call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
+c    &                   'w.m-2',3,zdtsw)
+c        call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
+c    &                   'w.m-2',3,zdtlw)
+
+!!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL
+        call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
+     &                       "K",3,tsoil)
+        call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
+     &                       "K",3,inertiedat)
+!!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL
+
+c        ----------------------------------------------------------
+c        Outputs of the CO2 cycle
+c        ----------------------------------------------------------
+
+         if (tracer.and.(igcm_co2.ne.0)) then
+!          call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
+!    &                     "kg/kg",2,zq(1,1,igcm_co2))
+!          call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
+!    &                     "kg/kg",3,zq(1,1,igcm_co2))
+        
+         ! Compute co2 column
+         call zerophys(ngrid,co2col)
+         do l=1,nlayermx
+           do ig=1,ngrid
+             co2col(ig)=co2col(ig)+
+     &                  zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g
+           enddo
+         enddo
+c         call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
+c     &                  co2col)
+         endif ! of if (tracer.and.(igcm_co2.ne.0))
+
+c        ----------------------------------------------------------
+c        Outputs of the water cycle
+c        ----------------------------------------------------------
+         if (tracer) then
+           if (water) then
+
+            !!!! waterice = q01, voir readmeteo.F90
+            call WRITEDIAGFI(ngridmx,'q01',noms(iq),
+     &                      'kg/kg',3,
+     &                       zq(1:ngridmx,1:nlayermx,igcm_h2o_ice))
+            !!!! watervapor = q02, voir readmeteo.F90
+            call WRITEDIAGFI(ngridmx,'q02',noms(iq),
+     &                      'kg/kg',3,
+     &                       zq(1:ngridmx,1:nlayermx,igcm_h2o_vap))
+
+c            call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
+c     &                     'kg.m-2',2,qsurf(1,iq))
+
+
+c             CALL WRITEDIAGFI(ngridmx,'mtot',
+c     &                       'total mass of water vapor',
+c     &                       'kg/m2',2,mtot)
+c             CALL WRITEDIAGFI(ngridmx,'icetot',
+c     &                       'total mass of water ice',
+c     &                       'kg/m2',2,icetot)
+cc            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+cc    &                *mugaz/mmol(igcm_h2o_ice)
+cc            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
+cc    &                       'mol/mol',3,vmr)
+c             CALL WRITEDIAGFI(ngridmx,'reffice',
+c     &                       'Mean reff',
+c     &                       'm',2,rave)
+cc            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
+cc    &                       'm',3,rice)
+cc            If activice is true, tauTES is computed in aeropacity.F;
+c             if (.not.activice) then
+c               CALL WRITEDIAGFI(ngridmx,'tauTESap',
+c     &                         'tau abs 825 cm-1',
+c     &                         '',2,tauTES)
+c             endif
+c             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
+c     &                       'surface h2o_ice',
+c     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+           endif !(water)
+
+
+           if (water.and..not.photochem) then
+             iq=nq
+c            write(str2(1:2),'(i2.2)') iq
+c            call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
+c    &                       'kg.m-2',2,zdqscloud(1,iq))
+c            call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
+c    &                       'kg/kg',3,zdqchim(1,1,iq))
+c            call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
+c    &                       'kg/kg',3,zdqdif(1,1,iq))
+c            call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
+c    &                       'kg/kg',3,zdqadj(1,1,iq))
+c            call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
+c    &                       'kg/kg',3,zdqc(1,1,iq))
+           endif  !(water.and..not.photochem)
+         endif
+
+c        ----------------------------------------------------------
+c        Outputs of the dust cycle
+c        ----------------------------------------------------------
+
+c        call WRITEDIAGFI(ngridmx,'tauref',
+c    &                    'Dust ref opt depth','NU',2,tauref)
+
+         if (tracer.and.(dustbin.ne.0)) then
+c          call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
+           if (doubleq) then
+cc            call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
+cc    &                       'kg.m-2',2,qsurf(1,1))
+cc            call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
+cc    &                       'N.m-2',2,qsurf(1,2))
+cc            call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
+cc    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
+cc            call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
+cc    &                       'kg.m-2.s-1',2,zdqssed(1,1))
+c             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
+c     &                        'm',3,rdust*ref_r0)
+             call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
+     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
+            call WRITEDIAGFI(ngridmx,'dustN','Dust number',
+     &                        'part/kg',3,pq(1,1,igcm_dust_number))
+           else
+             do iq=1,dustbin
+               write(str2(1:2),'(i2.2)') iq
+               call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
+     &                         'kg/kg',3,zq(1,1,iq))
+               call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
+     &                         'kg.m-2',2,qsurf(1,iq))
+             end do
+           endif ! (doubleq)
+c          if (submicron) then
+c            call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr',
+c    &                        'kg/kg',3,pq(1,1,igcm_dust_submicron))
+c          endif ! (submicron)
+         end if  ! (tracer.and.(dustbin.ne.0))
+
+c        ----------------------------------------------------------
+c        Output in netcdf file "diagsoil.nc" for subterranean
+c          variables (output every "ecritphy", as for writediagfi)
+c        ----------------------------------------------------------
+
+         ! Write soil temperature
+!        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
+!    &                     3,tsoil)
+         ! Write surface temperature
+!        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
+!    &                     2,tsurf)
+
+c        ==========================================================
+c        END OF WRITEDIAGFI
+c        ==========================================================
+
+      ELSE     ! if(ngrid.eq.1)
+
+         print*,'Ls =',zls*180./pi,
+     &  '  tauref(700 Pa) =',tauref
+c      ----------------------------------------------------------------------
+c      Output in grads file "g1d" (ONLY when using testphys1d)
+c      (output at every X physical timestep)
+c      ----------------------------------------------------------------------
+c
+c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
+c         CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
+c         CALL writeg1d(ngrid,1,ps,'ps','Pa')
+         
+c         CALL writeg1d(ngrid,nlayer,zt,'T','K')
+c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
+
+! or output in diagfi.nc (for testphys1d)
+         call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
+         call WRITEDIAGFI(ngridmx,'temp','Temperature',
+     &                       'K',1,zt)
+
+         if(tracer) then
+c           CALL writeg1d(ngrid,1,tau,'tau','SI')
+            do iq=1,nq
+c              CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 
+               call WRITEDIAGFI(ngridmx,trim(noms(iq)),
+     &              trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
+            end do
+         end if
+
+         zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
+
+         do l=2,nlayer-1
+            tmean=zt(1,l)
+            if(zt(1,l).ne.zt(1,l-1))
+     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
+              zlocal(l)= zlocal(l-1)
+     &        -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g
+         enddo
+         zlocal(nlayer)= zlocal(nlayer-1)-
+     &                   log(pplay(1,nlayer)/pplay(1,nlayer-1))*
+     &                   rnew(1,nlayer)*tmean/g
+
+      END IF       ! if(ngrid.ne.1)
+
+      icount=icount+1
+      RETURN
+      END
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F.newold
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F.newold	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F.newold	(revision 214)
@@ -0,0 +1,1351 @@
+      SUBROUTINE physiq(ngrid,nlayer,nq,
+     $            firstcall,lastcall,
+     $            pday,ptime,ptimestep,
+     $            pplev,pplay,pphi,
+     $            pu,pv,pt,pq,
+     $            pw,
+     $            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
+
+
+      IMPLICIT NONE
+c=======================================================================
+c
+c   subject:
+c   --------
+c
+c   Organisation of the physical parametrisations of the LMD 
+c   martian atmospheric general circulation model.
+c
+c   The GCM can be run without or with tracer transport
+c   depending on the value of Logical "tracer" in file  "callphys.def"
+c   Tracers may be water vapor, ice OR chemical species OR dust particles
+c
+c   SEE comments in initracer.F about numbering of tracer species...
+c
+c   It includes:
+c
+c      1. Initialization:
+c      1.1 First call initializations
+c      1.2 Initialization for every call to physiq
+c      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
+c      2. Compute radiative transfer tendencies
+c         (longwave and shortwave) for CO2 and dust.
+c      3. Gravity wave and subgrid scale topography drag :
+c      4. Vertical diffusion (turbulent mixing):
+c      5. Convective adjustment
+c      6. Condensation and sublimation of carbon dioxide.
+c      7.  TRACERS :
+c       7a. water and water ice
+c       7b. call for photochemistry when tracers are chemical species
+c       7c. other scheme for tracer (dust) transport (lifting, sedimentation)
+c       7d. updates (CO2 pressure variations, surface budget)
+c      8. Contribution to tendencies due to thermosphere
+c      9. Surface and sub-surface temperature calculations
+c     10. Write outputs :
+c           - "startfi", "histfi" (if it's time)
+c           - Saving statistics (if "callstats = .true.")
+c           - Dumping eof (if "calleofdump = .true.")
+c           - Output any needed variables in "diagfi" 
+c      10. Diagnostic: mass conservation of tracers
+c 
+c   author: 
+c   ------- 
+c           Frederic Hourdin	15/10/93
+c           Francois Forget		1994
+c           Christophe Hourdin	02/1997 
+c           Subroutine completly rewritten by F.Forget (01/2000)
+c           Introduction of the photochemical module: S. Lebonnois (11/2002)
+c           Introduction of the thermosphere module: M. Angelats i Coll (2002)
+c           Water ice clouds: Franck Montmessin (update 06/2003)
+c           
+c
+c
+c   arguments:
+c   ----------
+c
+c   input:
+c   ------
+c    ecri                  period (in dynamical timestep) to write output
+c    ngrid                 Size of the horizontal grid.
+c                          All internal loops are performed on that grid.
+c    nlayer                Number of vertical layers.
+c    nq                    Number of advected fields
+c    firstcall             True at the first call
+c    lastcall              True at the last call
+c    pday                  Number of days counted from the North. Spring
+c                          equinoxe.
+c    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
+c    ptimestep             timestep (s)
+c    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
+c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
+c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
+c    pu(ngrid,nlayer)      u component of the wind (ms-1)
+c    pv(ngrid,nlayer)      v component of the wind (ms-1)
+c    pt(ngrid,nlayer)      Temperature (K)
+c    pq(ngrid,nlayer,nq)   Advected fields
+c    pudyn(ngrid,nlayer)    \ 
+c    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
+c    ptdyn(ngrid,nlayer)     / corresponding variables
+c    pqdyn(ngrid,nlayer,nq) /
+c    pw(ngrid,?)           vertical velocity
+c
+c   output:
+c   -------
+c
+c    pdu(ngrid,nlayermx)        \
+c    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
+c    pdt(ngrid,nlayermx)         /  variables due to physical processes.
+c    pdq(ngrid,nlayermx)        /
+c    pdpsrf(ngrid)             /
+c    tracerdyn                 call tracer in dynamical part of GCM ?
+
+c
+c=======================================================================
+c
+c    0.  Declarations :
+c    ------------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comgeomfi.h"
+#include "surfdat.h"
+#include "comdiurn.h"
+#include "callkeys.h"
+#include "comcstfi.h"
+#include "planete.h"
+#include "comsaison.h"
+#include "control.h"
+#include "dimradmars.h"
+#include "comg1d.h"
+#include "tracer.h"
+#include "nlteparams.h"
+
+#include "chimiedata.h"
+#include "watercap.h"
+#include "fisice.h"
+#include "param.h"
+#include "param_v3.h"
+#include "conc.h"
+
+#include "netcdf.inc"
+
+
+
+c Arguments :
+c -----------
+
+c   inputs:
+c   -------
+      INTEGER ngrid,nlayer,nq
+      REAL ptimestep
+      REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
+      REAL pphi(ngridmx,nlayer)
+      REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
+      REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
+      REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
+      REAL zh(ngridmx,nlayermx)      ! potential temperature (K)
+      LOGICAL firstcall,lastcall
+
+      REAL pday
+      REAL ptime 
+      logical tracerdyn
+
+c   outputs:
+c   --------
+c     physical tendencies
+      REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
+      REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
+      REAL pdpsrf(ngridmx) ! surface pressure tendency
+
+
+c Local saved variables:
+c ----------------------
+c     aerosol (dust or ice) extinction optical depth  at reference wavelength 
+c     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
+c      aerosol optical properties  :
+      REAL aerosol(ngridmx,nlayermx,naerkind) 
+
+      INTEGER day_ini  ! Initial date of the run (sol since Ls=0) 
+      INTEGER icount     ! counter of calls to physiq during the run.
+      REAL tsurf(ngridmx)            ! Surface temperature (K)
+      REAL tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
+      REAL co2ice(ngridmx)           ! co2 ice surface layer (kg.m-2)  
+      REAL albedo(ngridmx,2)         ! Surface albedo in each solar band
+      REAL emis(ngridmx)             ! Thermal IR surface emissivity
+      REAL dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
+      REAL fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
+      REAL fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
+      REAL capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
+      REAL fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
+      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
+      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
+      INTEGER ig_vl1                 ! Grid Point near VL1   (for diagnostic) 
+
+      SAVE day_ini, icount
+      SAVE aerosol, tsurf,tsoil
+      SAVE co2ice,albedo,emis, q2
+      SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
+      SAVE ig_vl1
+
+      REAL stephan   
+      DATA stephan/5.67e-08/  ! Stephan Boltzman constant
+      SAVE stephan
+
+c Local variables :
+c -----------------
+
+
+      CHARACTER*80 fichier 
+      INTEGER l,ig,ierr,igout,iq, tapphys
+      INTEGER iqmin                     ! Used if iceparty engaged
+
+      REAL fluxsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2)
+      REAL fluxsurf_sw(ngridmx,2)    !incident SW (solar) surface flux (W.m-2)
+      REAL fluxtop_lw(ngridmx)       !Outgoing LW (IR) flux to space (W.m-2)
+      REAL fluxtop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
+c     for clear area (uncovered by clouds) :
+      REAL clsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2) 
+      REAL clsurf_sw(ngridmx,2)    !incident SW (solar) surface flux (W.m-2)
+      REAL cltop_lw(ngridmx)       !Outgoing LW (IR) flux to space (W.m-2)
+      REAL cltop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
+      REAL tauref(ngridmx)           ! Reference column optical depth at 700 Pa
+                                     ! (used if active=F) 
+      REAL tau(ngridmx,naerkind)     ! Column dust optical depth at each point
+      REAL zls                       !  solar longitude (rad)
+      REAL zday                      ! date (time since Ls=0, in martian days)
+      REAL zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
+      REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
+      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
+
+
+c     Tendancies due to various processes:
+      REAL dqsurf(ngridmx,nqmx)
+      REAL zdtlw(ngridmx,nlayermx)     ! (K/s)
+      REAL zdtsw(ngridmx,nlayermx)     ! (K/s)
+      REAL cldtlw(ngridmx,nlayermx)     ! (K/s) LW heating rate for clear area
+      REAL cldtsw(ngridmx,nlayermx)     ! (K/s) SW heating rate for clear area
+      REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
+      REAL zdtnlte(ngridmx,nlayermx)   ! (K/s)
+      REAL zdtsurf(ngridmx)            ! (K/s)
+      REAL zdtcloud(ngridmx,nlayermx)
+      REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
+      REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhadj(ngridmx,nlayermx)                           ! (K/s)
+      REAL zdtgw(ngridmx,nlayermx)                            ! (K/s)
+      REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
+      REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
+      REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
+
+      REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
+      REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
+      REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
+      REAL zdqadj(ngridmx,nlayermx,nqmx)
+      REAL zdqc(ngridmx,nlayermx,nqmx)
+      REAL zdqcloud(ngridmx,nlayermx,nqmx)
+      REAL zdqscloud(ngridmx,nqmx)
+      REAL zdqchim(ngridmx,nlayermx,nqmx)
+      REAL zdqschim(ngridmx,nqmx)
+
+      REAL zdteuv(ngridmx,nlayermx)    ! (K/s)
+      REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
+      REAL zdumolvis(ngridmx,nlayermx)
+      REAL zdvmolvis(ngridmx,nlayermx)
+      real zdqmoldiff(ngridmx,nlayermx,nqmx)
+
+c     Local variable for local intermediate calcul:
+      REAL zflubid(ngridmx)
+      REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
+      REAL zdum1(ngridmx,nlayermx)
+      REAL zdum2(ngridmx,nlayermx)
+      REAL ztim1,ztim2,ztim3, z1,z2
+      REAL ztime_fin
+      REAL zdh(ngridmx,nlayermx)
+      REAL pclc_min ! minimum of the cloud fraction over the whole domain
+      INTEGER length
+      PARAMETER (length=100)
+
+c local variables only used for diagnostic (output in file "diagfi" or "stats")
+c -----------------------------------------------------------------------------
+      REAL ps(ngridmx), zt(ngridmx,nlayermx)
+      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
+      REAL zq(ngridmx,nlayermx,nqmx)
+      REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
+      character*2 str2
+      character*5 str5
+      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
+      real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
+      real qtot1,qtot2 ! total aerosol mass
+      integer igmin, lmin
+      logical tdiag
+
+
+      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
+      REAL zstress(ngrid), cd
+      real hco2(nqmx),tmean, zlocal(nlayermx)
+      real rho(ngridmx,nlayermx) ! density
+      real vmr(ngridmx,nlayermx) ! volume mixing ratio
+
+      REAL time_phys
+     
+c=======================================================================
+
+c 1. Initialisation:
+c -----------------
+
+c  1.1   Initialisation only at first call
+c  ---------------------------------------
+      IF (firstcall) THEN
+
+c        variables set to 0
+c        ~~~~~~~~~~~~~~~~~~
+         call zerophys(ngrid*nlayer*naerkind,aerosol)
+         call zerophys(ngrid*nlayer,dtrad)
+         call zerophys(ngrid,fluxrad)
+
+c        read startfi 
+c        ~~~~~~~~~~~~
+
+! Read netcdf initial physical parameters.
+         CALL phyetat0 ("startfi.nc",0,0,
+     &         nsoilmx,nq,
+     &         day_ini,time_phys,
+     &         tsurf,tsoil,emis,q2,qsurf,co2ice)
+
+         if (pday.ne.day_ini) then
+           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
+     &                "physics and dynamics"
+           write(*,*) "dynamics day: ",pday
+           write(*,*) "physics day:  ",day_ini
+           stop
+         endif
+
+         write (*,*) 'In physiq day_ini =', day_ini
+
+c        Initialize albedo and orbital calculation
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         CALL surfini(ngrid,co2ice,qsurf,albedo)
+         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+
+c        initialize soil 
+c        ~~~~~~~~~~~~~~~
+         IF (callsoil) THEN
+            CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
+     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+         ELSE
+            PRINT*,'WARNING! Thermal conduction in the soil turned off'
+            DO ig=1,ngrid
+               capcal(ig)=1.e5
+               fluxgrd(ig)=0.
+            ENDDO
+         ENDIF
+         icount=1
+
+
+c        initialize tracers
+c        ~~~~~~~~~~~~~~~~~~
+         tracerdyn=tracer
+         IF (tracer) THEN
+            CALL initracer(qsurf,co2ice)
+         ENDIF  ! end tracer
+
+c        Determining gridpoint near Viking Lander 1 (used for diagnostic only)
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if(ngrid.ne.1) then
+           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.)
+           write(*,*) 'Viking Lander 1 GCM point: lat,lon',
+     &              lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi
+         end if 
+
+c        Initialize thermospheric parameters
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if (callthermos) call param_read
+
+c        Initialize R and Cp as constant
+
+         if (.not.callthermos .and. .not.photochem) then
+                 do l=1,nlayermx
+                  do ig=1,ngridmx
+                   rnew(ig,l)=r
+                   cpnew(ig,l)=cpp
+                   mmean(ig,l)=mugaz
+                   enddo
+                  enddo  
+         endif         
+                   
+      ENDIF        !  (end of "if firstcall")
+
+c ---------------------------------------------------
+c 1.2   Initializations done at every physical timestep:
+c ---------------------------------------------------
+c
+      IF (ngrid.NE.ngridmx) THEN
+         PRINT*,'STOP in PHYSIQ'
+         PRINT*,'Probleme de dimensions :'
+         PRINT*,'ngrid     = ',ngrid
+         PRINT*,'ngridmx   = ',ngridmx
+         STOP
+      ENDIF
+
+c     Initialize various variables
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      call zerophys(ngrid*nlayer, pdv)
+      call zerophys(ngrid*nlayer, pdu)
+      call zerophys(ngrid*nlayer, pdt)
+      call zerophys(ngrid*nlayer*nq, pdq)
+      call zerophys(ngrid, pdpsrf)
+      call zerophys(ngrid, zflubid)
+      call zerophys(ngrid, zdtsurf)
+      call zerophys(ngrid*nq, dqsurf)
+      igout=ngrid/2+1 
+
+
+      zday=pday+ptime ! compute time, in sols (and fraction thereof)
+
+c     Compute Solar Longitude (Ls) :
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (season) then
+         call solarlong(zday,zls)
+      else
+         call solarlong(float(day_ini),zls)
+      end if
+
+c     Compute geopotential at interlayers
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c     ponderation des altitudes au niveau des couches en dp/p
+
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zzlay(ig,l)=pphi(ig,l)/g
+         ENDDO
+      ENDDO
+      DO ig=1,ngrid
+         zzlev(ig,1)=0.
+         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+      ENDDO
+      DO l=2,nlayer
+         DO ig=1,ngrid
+            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
+            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
+            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+         ENDDO
+      ENDDO
+
+
+!     Potential temperature calculation not the same in physiq and dynamic
+
+c     Compute potential temperature
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO l=1,nlayer
+         DO ig=1,ngrid 
+            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
+            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c    1.2.5 Compute mean mass, cp, and R
+c    --------------------------------
+
+      if(photochem.or.callthermos) then
+         call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
+      endif
+c-----------------------------------------------------------------------
+c    2. Compute radiative tendencies :
+c------------------------------------
+
+
+      IF (callrad) THEN
+         IF( MOD(icount-1,iradia).EQ.0) THEN
+
+c          Local Solar zenith angle
+c          ~~~~~~~~~~~~~~~~~~~~~~~~
+           CALL orbite(zls,dist_sol,declin)
+
+           IF(diurnal) THEN
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
+     s         ztim1,ztim2,ztim3, mu0,fract)
+
+           ELSE
+               CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
+           ENDIF
+
+c          NLTE cooling from CO2 emission
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+           IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
+
+c          Find number of layers for LTE radiation calculations
+           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
+     &          CALL nlthermeq(ngrid,nlayer,pplev,pplay)
+
+
+c          Atmospheric dust opacity and aerosol distribution:
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           CALL dustopacity(ngrid,nlayer,nq,zday,pplay,pplev,zls,pq,
+     $      tauref,tau,aerosol)
+         
+c          Call main radiative transfer scheme
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c          Transfer through dust and CO2, except NIR CO2 absorption
+
+           CALL callradite(icount,ngrid,nlayer,aerosol,albedo,
+     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
+     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw)
+
+c          CO2 near infrared absorption
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           call zerophys(ngrid*nlayer,zdtnirco2)
+           if (callnirco2) then
+              call nirco2abs (ngrid,nlayer,pplay,dist_sol,
+     .                       mu0,fract,declin, zdtnirco2)
+           endif
+
+c          Radiative flux from the sky absorbed by the surface (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           DO ig=1,ngrid
+               fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
+     $         +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
+     $         +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
+           ENDDO
+
+
+c          Net atmospheric radiative heating rate (K.s-1)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           IF(callnlte) THEN
+              CALL blendrad(ngrid, nlayer, pplay,
+     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
+           ELSE
+              DO l=1,nlayer
+                 DO ig=1,ngrid
+                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
+     &                          +zdtnirco2(ig,l)
+                  ENDDO
+              ENDDO
+           ENDIF
+
+
+
+        ENDIF ! of if(mod(icount-1,iradia).eq.0)
+
+
+
+c    Transformation of the radiative tendencies:
+c    -------------------------------------------
+
+c          Net radiative surface flux (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+           DO ig=1,ngrid
+               zplanck(ig)=tsurf(ig)*tsurf(ig)
+               zplanck(ig)=emis(ig)*
+     $         stephan*zplanck(ig)*zplanck(ig)
+               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
+           ENDDO
+
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
+            ENDDO
+         ENDDO
+
+      ENDIF ! of IF (callrad)
+
+c-----------------------------------------------------------------------
+c    3. Gravity wave and subgrid scale topography drag :
+c    -------------------------------------------------
+
+
+      IF(calllott)THEN
+
+        CALL calldrag_noro(ngrid,nlayer,ptimestep,
+     &                 pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
+ 
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
+          ENDDO
+        ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c    4. Vertical diffusion (turbulent mixing):
+c    -----------------------------------------
+c
+      IF (calldifv) THEN
+
+
+         DO ig=1,ngrid
+            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
+         ENDDO
+
+         CALL zerophys(ngrid*nlayer,zdum1)
+         CALL zerophys(ngrid*nlayer,zdum2)
+         do l=1,nlayer
+            do ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            enddo
+         enddo
+
+c        Calling vdif (Martian version WITH CO2 condensation)
+         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
+     $        ptimestep,capcal,lwrite,
+     $        pplay,pplev,zzlay,zzlev,z0,
+     $        pu,pv,zh,pq,tsurf,emis,qsurf,
+     $        zdum1,zdum2,zdh,pdq,zflubid,
+     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
+     &        zdqdif,zdqsdif)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
+               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
+
+               zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
+
+            ENDDO
+         ENDDO
+
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
+         ENDDO
+
+         if (tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+           DO iq=1, nq
+              DO ig=1,ngrid
+                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+              ENDDO
+           ENDDO
+         end if ! of if (tracer)
+
+      ELSE    
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+
+     s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+         ENDDO
+      ENDIF ! of IF (calldifv)
+
+
+c-----------------------------------------------------------------------
+c   5. Dry convective adjustment:
+c   -----------------------------
+
+      IF(calladj) THEN
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            ENDDO
+         ENDDO
+         CALL zerophys(ngrid*nlayer,zduadj)
+         CALL zerophys(ngrid*nlayer,zdvadj)
+         CALL zerophys(ngrid*nlayer,zdhadj)
+
+         CALL convadj(ngrid,nlayer,nq,ptimestep,
+     $                pplay,pplev,zpopsk,
+     $                pu,pv,zh,pq,
+     $                pdu,pdv,zdh,pdq,
+     $                zduadj,zdvadj,zdhadj,
+     $                zdqadj)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
+               pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
+
+               zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
+            ENDDO
+         ENDDO
+
+         if(tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         end if
+      ENDIF ! of IF(calladj)
+
+c-----------------------------------------------------------------------
+c   6. Carbon dioxide condensation-sublimation:
+c   -------------------------------------------
+
+      IF (callcond) THEN
+         CALL newcondens(ngrid,nlayer,nq,ptimestep,
+     $              capcal,pplay,pplev,tsurf,pt,
+     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
+     $              co2ice,albedo,emis,
+     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
+     $	            fluxsurf_sw) 
+
+         DO l=1,nlayer
+           DO ig=1,ngrid
+             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
+             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
+             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
+           ENDDO
+         ENDDO
+         DO ig=1,ngrid
+            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
+         ENDDO
+
+         IF (tracer) THEN
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         ENDIF ! of IF (tracer)
+
+      ENDIF  ! of IF (callcond)
+
+c-----------------------------------------------------------------------
+c   7. Specific parameterizations for tracers 
+c:   -----------------------------------------
+
+      if (tracer) then 
+
+c   7a. Water and ice
+c     ---------------
+
+c        ---------------------------------------
+c        Water ice condensation in the atmosphere
+c        ----------------------------------------
+         IF (water) THEN
+           call watercloud(ngrid,nlayer, ptimestep,
+     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
+     &                pq,pdq,zdqcloud,qsurf,zdqscloud,zdtcloud,
+     &                nq,naerkind,tau,icount,zls)
+
+           if (activice) then
+c Temperature variation due to latent heat release
+           DO l=1,nlayer
+             DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
+             ENDDO
+           ENDDO
+           endif
+
+           IF (iceparty) then
+             iqmin=nq-1
+           ELSE
+             iqmin=nq
+           ENDIF
+
+           DO iq=iqmin,nq
+             DO l=1,nlayer
+                DO ig=1,ngrid
+                   pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
+                ENDDO
+             ENDDO
+             DO ig=1,ngrid
+                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
+             ENDDO
+           ENDDO
+
+         END IF  ! of IF (water)
+
+c   7b. Chemical species
+c     ------------------
+
+c        --------------
+c        photochemistry :
+c        --------------
+         IF (photochem .or. thermochem) then
+          call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
+     .      zzlay,zday,pq,pdq,zdqchim,zdqschim,zdqcloud,zdqscloud)
+           
+c Photochemistry includes condensation of H2O2
+
+           do iq=nqchem_min,nq
+            if (noms(iq).eq."h2o2") then
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
+               ENDDO
+             ENDDO
+            else
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
+               ENDDO
+             ENDDO
+            endif
+           ENDDO
+           do iq=nqchem_min,nq
+            if (noms(iq).eq."h2o2") then
+               DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
+               ENDDO
+            else
+               DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
+               ENDDO
+            endif
+           ENDDO
+
+         END IF  ! of IF (photochem.or.thermochem)
+
+c   7c. Aerosol particles
+c     -------------------
+
+c        ----------
+c        Dust devil :
+c        ----------
+         IF(callddevil) then 
+           call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
+     &                zdqdev,zdqsdev)
+ 
+           if (dustbin.ge.1) then
+              do iq=1,nq
+                 DO l=1,nlayer
+                    DO ig=1,ngrid
+                       pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
+                    ENDDO
+                 ENDDO
+              enddo
+              do iq=1,nq
+                 DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
+                 ENDDO
+              enddo
+           endif  ! of if (dustbin.ge.1)
+
+         END IF ! of IF (callddevil)
+
+c        ------------- 
+c        Sedimentation :   acts also on water ice
+c        ------------- 
+         IF (sedimentation) THEN 
+           call zerophys(ngrid*nlayer*nq, zdqsed)
+           call zerophys(ngrid*nq, zdqssed)
+
+           if(doubleq) then
+            call callsedim2q(ngrid,nlayer, ptimestep,
+     &                pplev,zzlev, pt,
+     &                pq, pdq, zdqsed, zdqssed,nq)
+           else
+            call callsedim(ngrid,nlayer, ptimestep,
+     &                pplev,zzlev, pt,
+     &                pq, pdq, zdqsed, zdqssed,nq)
+           end if
+
+
+           DO iq=1, nq
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
+               ENDDO
+             ENDDO
+           ENDDO
+           DO iq=1, nq
+             DO ig=1,ngrid
+                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
+             ENDDO
+           ENDDO
+         END IF   ! of IF (sedimentation)
+
+c   7d. Updates
+c     ---------
+
+        DO iq=1, nq
+          DO ig=1,ngrid
+
+c       ---------------------------------
+c       Updating tracer budget on surface
+c       ---------------------------------
+            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
+
+          ENDDO  ! (ig)
+        ENDDO    ! (iq)
+
+      endif !  of if (tracer) 
+
+
+c-----------------------------------------------------------------------
+c   8. THERMOSPHERE CALCULATION
+c-----------------------------------------------------------------------
+
+      if (callthermos) then
+        call thermosphere(pplev,pplay,dist_sol,
+     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
+     &     pt,pq,pu,pv,pdt,pdq,
+     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
+c           do iq=nqchem_min,nq
+c           write(*,*) 'thermo iq,pq',iq,pq(690,1,iq)
+c           enddo
+
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
+     &                         +zdteuv(ig,l)
+            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
+            DO iq=1, nq
+              pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
+            ENDDO
+          ENDDO
+        ENDDO
+
+
+      endif
+c-----------------------------------------------------------------------
+c   9. Surface  and sub-surface soil temperature
+c-----------------------------------------------------------------------
+c
+c
+c   9.1 Increment Surface temperature:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      DO ig=1,ngrid
+         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
+      ENDDO
+
+c  Prescribe a cold trap at south pole (except at high obliquity !!)
+c  Temperature at the surface is set there to be the temperature
+c  corresponding to equilibrium temperature between phases of CO2
+
+      IF (tracer.AND.water.AND.ngridmx.NE.1) THEN
+         if (caps.and.(obliquit.lt.27.)) then
+           ! NB: Updated surface pressure, at grid point 'ngrid', is
+           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
+           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
+     &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
+         endif
+c       -------------------------------------------------------------
+c       Change of surface albedo (set to 0.4) in case of ground frost
+c       everywhere except on the north permanent cap and in regions
+c       covered by dry ice. 
+c              ALWAYS PLACE these lines after newcondens !!!
+c       -------------------------------------------------------------
+         do ig=1,ngrid
+
+c       -------------- Version temporaire fit TES 2008 ------------
+         if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
+              albedo(ig,1)=0.45                                     
+              albedo(ig,2)=0.45                                     
+          endif
+
+c         if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
+c           if (.not.watercaptag(ig)) then
+c             albedo(ig,1)=0.4
+c             albedo(ig,2)=0.4
+c           endif
+c         endif
+c       -------------- version Francois ---------------
+c          if (co2ice(ig).eq.0.and.
+c    &        ((qsurf(ig,nqmx).gt.0.005).or.(watercaptag(ig)))) then
+c              albedo(ig,1)=max(0.4,albedodat(ig))
+c              albedo(ig,2)=albedo(ig,1)
+c          endif
+c       ---------------------------------------------
+         enddo  ! of do ig=1,ngrid
+      ENDIF  ! of IF (tracer.AND.water.AND.ngridmx.NE.1)
+
+c
+c   9.2 Compute soil temperatures and subsurface heat flux:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (callsoil) THEN
+         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
+     &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+      ENDIF
+
+c-----------------------------------------------------------------------
+c  10. Write output files
+c  ----------------------
+
+c    -------------------------------
+c    Dynamical fields incrementation
+c    -------------------------------
+c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
+      ! temperature, zonal and meridional wind
+      DO l=1,nlayer
+        DO ig=1,ngrid
+          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
+        ENDDO
+      ENDDO
+
+      ! tracers
+      DO iq=1, nq
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
+          ENDDO
+        ENDDO
+      ENDDO
+
+      ! surface pressure
+      DO ig=1,ngrid
+          ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
+      ENDDO
+
+      ! pressure
+      DO l=1,nlayer
+        DO ig=1,ngrid
+             zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
+             zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
+        ENDDO
+      ENDDO
+
+      ! Density 
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
+         ENDDO
+      ENDDO
+
+c    Compute surface stress : (NB: z0 is a common in planete.h)
+c     DO ig=1,ngrid
+c        cd = (0.4/log(zzlay(ig,1)/z0))**2
+c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
+c     ENDDO
+
+c     Sum of fluxes in solar spectral bands (for output only)
+      DO ig=1,ngrid
+	     fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
+	     fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
+      ENDDO
+c ******* TEST ******************************************************
+      ztim1 = 999
+      DO l=1,nlayer
+        DO ig=1,ngrid
+           if (pt(ig,l).lt.ztim1) then
+               ztim1 = pt(ig,l)
+               igmin = ig
+               lmin = l 
+           end if
+        ENDDO
+      ENDDO
+      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then	
+        write(*,*) 'PHYSIQ: stability WARNING :'
+        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
+     &              'ig l =', igmin, lmin
+      end if
+c *******************************************************************
+
+c     ---------------------
+c     Outputs to the screen 
+c     ---------------------
+
+      IF (lwrite) THEN
+         PRINT*,'Global diagnostics for the physics'
+         PRINT*,'Variables and their increments x and dx/dt * dt'
+         WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
+         WRITE(*,'(2f10.5,2f15.5)')
+     s   tsurf(igout),zdtsurf(igout)*ptimestep,
+     s   pplev(igout,1),pdpsrf(igout)*ptimestep
+         WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
+         WRITE(*,'(i4,6f10.5)') (l,
+     s   pu(igout,l),pdu(igout,l)*ptimestep,
+     s   pv(igout,l),pdv(igout,l)*ptimestep,
+     s   pt(igout,l),pdt(igout,l)*ptimestep,
+     s   l=1,nlayer)
+      ENDIF ! of IF (lwrite)
+
+
+      IF (ngrid.NE.1) THEN
+         print*,'Ls =',zls*180./pi,
+     &   ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
+     &   ' tau(Viking1) =',tau(ig_vl1,1)
+
+c        -------------------------------------------------------------------
+c        Writing NetCDF file  "RESTARTFI" at the end of the run
+c        -------------------------------------------------------------------
+c        Note: 'restartfi' is stored just before dynamics are stored
+c              in 'restart'. Between now and the writting of 'restart',
+c              there will have been the itau=itau+1 instruction and
+c              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
+c              thus we store for time=time+dtvr
+
+         IF(lastcall) THEN
+            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) 
+            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
+     .              ptimestep,pday,
+     .              ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
+     .              area,albedodat,inertiedat,zmea,zstd,zsig,
+     .              zgam,zthe)
+         ENDIF
+
+c        -----------------------------------------------------------------
+c        Saving statistics :
+c        -----------------------------------------------------------------
+c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+c        which can later be used to make the statistic files of the run:
+c        "stats")          only possible in 3D runs !
+
+         
+         IF (callstats) THEN
+
+
+            call wstats(ngrid,"ps","Surface pressure","K",2,ps)
+            call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+            call wstats(ngrid,"co2ice","CO2 ice cover",
+     .                  "kg.m-2",2,co2ice)
+c            call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
+c    .                  emis)
+             call wstats(ngrid,"fluxsurf_lw",
+     .                 "Thermal IR radiative flux to surface","W.m-2",2,
+     .                 fluxsurf_lw)
+             call wstats(ngrid,"fluxsurf_sw",
+     .                  "Solar radiative flux to surface","W.m-2",2,
+     .                   fluxsurf_sw_tot)
+             call wstats(ngrid,"fluxtop_lw",
+     .                  "Thermal IR radiative flux to space","W.m-2",2,
+     .                  fluxtop_lw)
+             call wstats(ngrid,"fluxtop_sw",
+     .                  "Solar radiative flux to space","W.m-2",2,
+     .                  fluxtop_sw_tot)
+            call wstats(ngrid,"dod","Dust optical depth"," ",2,tau)
+
+            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
+            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
+            call wstats(ngrid,"v","Meridional (North-South) wind",
+     .                  "m.s-1",3,zv)
+            call wstats(ngrid,"w","Vertical (down-up) wind",
+     .                  "m.s-1",3,pw)
+            call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
+             call wstats(ngrid,"q2",
+     .               "Boundary layer eddy kinetic energy","m2.s-2",3,q2)
+
+           if (tracer) then
+            if (water) then
+               vmr=zq(1:ngridmx,1:nlayermx,nqmx)*mugaz/mmol(nqmx)
+               call wstats(ngrid,"vmr_h2ovapor",
+     .                    "H2O vapor volume mixing ratio","mol/mol",
+     .                    3,vmr)
+               if (iceparty) then
+                  vmr=zq(1:ngridmx,1:nlayermx,nqmx-1)*mugaz/mmol(nqmx-1)
+                  call wstats(ngrid,"vmr_h2oice",
+     .                       "H2O ice volume mixing ratio","mol/mol",
+     .                       3,vmr)
+               endif
+            endif
+
+             if (thermochem.or.photochem) then
+                do iq=1,nq
+                   if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
+     .               (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
+     .               (noms(iq).eq."h2").or.
+     .               (noms(iq).eq."o3")) then
+                        do l=1,nlayer
+                          do ig=1,ngrid
+                            vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
+                          end do
+                        end do
+                        call wstats(ngrid,"vmr_"//trim(noms(iq)),
+     .                     "Volume mixing ratio","mol/mol",3,vmr)
+                   endif
+                enddo
+             endif
+           endif !tracer 
+
+            IF(lastcall) THEN
+              write (*,*) "Writing stats..."
+              call mkstats(ierr)
+            ENDIF
+          ENDIF !if callstats
+
+c        (Store EOF for Mars Climate database software)
+         IF (calleofdump) THEN
+            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
+         ENDIF
+
+
+c        ----------------------------------------------------------------------
+c        output in netcdf file "DIAGFI", containing any variable for diagnostic
+c        (output with  period "ecritphy", set in "run.def")
+c        ----------------------------------------------------------------------
+c        WRITEDIAGFI can ALSO be called from any other subroutines
+c        for any variables !!
+       call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     .                  emis)
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     .                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     .                  co2ice)
+c       call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
+c    .       fluxsurf_lw)
+c       call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
+c    .       fluxsurf_sw_tot)
+c       call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
+c    .       fluxtop_lw)
+c       call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
+c    .       fluxtop_sw_tot)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+c        call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
+c       call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
+c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
+c       call WRITEDIAGFI(ngridm,'Teta','T potentielle','K',3,zh)
+c       call WRITEDIAGFI(ngridm,'Pression','Pression','Pa',3,pplay)
+        call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil)
+
+
+c       OUTPUT of tracer mass mixing ratio and surface value :    
+       if (tracer) then
+c         (for photochemistry, outputs are done in calchim)
+          do iq=1,nqmx
+            write(str2(1:2),'(i2.2)') iq
+            call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq),
+     &                         'kg/kg',3,zq(1,1,iq))
+            call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
+     &                     'kg.m-2',2,qsurf(1,iq))
+          end do
+        end if
+c ***************************************************
+
+c  Outputs for H2O
+       if (tracer) then
+        if (activice) then
+c         call WRITEDIAGFI(ngridmx,'tauice','tau','SI',2,tau(1,2))
+c         call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
+c     &                   'w.m-2',3,zdtsw)
+c         call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
+c     &                   'w.m-2',3,zdtlw)
+        endif  !(activice)
+
+        if (water.and..not.photochem) then
+           iq=nq
+c          write(str2(1:2),'(i2.2)') iq
+c          call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
+c     &                    'kg.m-2',2,zdqscloud(1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
+c     &                    'kg/kg',3,zdqchim(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
+c     &                    'kg/kg',3,zdqdif(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
+c     &                    'kg/kg',3,zdqadj(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
+c     &                    'kg/kg',3,zdqc(1,1,iq))
+        endif  !(water.and..not.photochem)
+
+c       if (iceparty) then
+c          iq=nq-1
+c          write(str2(1:2),'(i2.2)') iq
+c          call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
+c    &                     'kg/kg',3,zq(1,1,iq))
+c       endif  !(iceparty)
+       endif
+
+c  Outputs for dust tracers
+
+         if (tracer.and.(dustbin.ne.0)) then
+            call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
+            if (doubleq) then
+               call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
+     &                         'kg.m-2',2,qsurf(1,1))
+               call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
+     &                         'N.m-2',2,qsurf(1,2))
+               call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
+     &                         'kg.m-2.s-1',2,zdqsdev(1,1))
+               call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
+     &                         'kg.m-2.s-1',2,zdqssed(1,1))
+               do  l=1,nlayer
+                  do ig=1, ngrid
+                     reff(ig,l)= ref_r0 *
+     &               (r3n_q*pq(ig,l,1)/max(pq(ig,l,2),0.01))**(1./3.)
+                     reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6)
+                  end do
+               end do
+               call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff)
+            else
+               do iq=1,dustbin
+                  write(str2(1:2),'(i2.2)') iq
+                  call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
+     &                            'kg/kg',3,zq(1,1,iq))
+                  call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
+     &                            'kg.m-2',2,qsurf(1,iq))
+               end do
+            endif ! (doubleq)
+         end if  ! (tracer.and.(dustbin.ne.0))
+
+
+      ELSE     ! if(ngrid.eq.1)
+
+         print*,'Ls =',zls*180./pi,
+     &  '  tauref(700 Pa) =',tauref,' local tau =',tau(1,1)
+c      ----------------------------------------------------------------------
+c      Output in grads file "g1d" (ONLY when using testphys1d)
+c      (output at every X physical timestep)
+c      ----------------------------------------------------------------------
+c
+c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
+c        CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
+         CALL writeg1d(ngrid,1,ps,'ps','Pa')
+         
+         CALL writeg1d(ngrid,nlayer,zt,'T','K')
+c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
+
+         if(tracer) then
+c           CALL writeg1d(ngrid,1,tau,'tau','SI')
+            do iq=1,nq
+              CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 
+            end do
+         end if
+
+         zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
+
+         do l=2,nlayer-1
+            tmean=zt(1,l)
+            if(zt(1,l).ne.zt(1,l-1))
+     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
+              zlocal(l)= zlocal(l-1)
+     &        -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g
+         enddo
+         zlocal(nlayer)= zlocal(nlayer-1)-
+     &                   log(pplay(1,nlayer)/pplay(1,nlayer-1))*
+     &                   rnew(1,nlayer)*tmean/g
+
+c         if(tracer) then
+c         do l=2,nlayer
+c            do iq=1,nq
+c               hco2(iq)=zq(1,l,iq)/zq(1,l-1,iq)
+c               hco2(iq)=-(zlocal(l)-zlocal(l-1))/log(hco2(iq))/1000.
+c            enddo
+c            write(225,*) l,pt(1,l),(hco2(iq),iq=1,6)
+c            write(226,*) l,(hco2(iq),iq=7,13)
+c         enddo
+c         endif
+
+      END IF       ! if(ngrid.ne.1)
+
+
+      icount=icount+1
+      RETURN
+      END
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/writediagfi.info
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/writediagfi.info	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/writediagfi.info	(revision 214)
@@ -0,0 +1,21 @@
+        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     .                  emis)
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     .                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     .                  co2ice)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil)
+        if (tracer) then
+            do iq=1,nqmx
+              write(str2(1:2),'(i2.2)') iq
+              call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq),
+     &                          'kg/kg',3,zq(1,1,iq))
+              call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
+     &                        'kg.m-2',2,qsurf(1,iq))
+            end do
+          end if
+
Index: trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/calltherm_interface.F90
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/calltherm_interface.F90	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/calltherm_interface.F90	(revision 214)
@@ -0,0 +1,1 @@
+link ../../libf_gcm/phymars/calltherm_interface.F90
Index: trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/calltherm_mars.F90
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/calltherm_mars.F90	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/calltherm_mars.F90	(revision 214)
@@ -0,0 +1,1 @@
+link ../../libf_gcm/phymars/calltherm_mars.F90
Index: trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/thermcell_dqupdown.F90
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/thermcell_dqupdown.F90	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/thermcell_dqupdown.F90	(revision 214)
@@ -0,0 +1,1 @@
+link ../../libf_gcm/phymars/thermcell_dqupdown.F90
Index: trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/thermcell_main_mars.F90
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/thermcell_main_mars.F90	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS.new/libf/phymars/thermcell_main_mars.F90	(revision 214)
@@ -0,0 +1,1 @@
+link ../../libf_gcm/phymars/thermcell_main_mars.F90
Index: trunk/MESOSCALE/LMDZ.MARS/LMDZ.MARS.source.txt
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/LMDZ.MARS.source.txt	(revision 213)
+++ trunk/MESOSCALE/LMDZ.MARS/LMDZ.MARS.source.txt	(revision 214)
@@ -1,1 +1,1 @@
-link ../LMD_MM_MARS/SIMU/in_lmdz_mars/LMDZ.MARS.source.txt
+link in_lmdz_mars/LMDZ.MARS.source.txt
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/LMDZ.MARS.source.txt
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/LMDZ.MARS.source.txt	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/LMDZ.MARS.source.txt	(revision 214)
@@ -0,0 +1,306 @@
+--------------------------------------------------------------------
+Changes in the LMD MARS GCM reference sources 
+(most recent version available in LMDZ.MARS.source.tar)
+--------------------------------------------------------------------
+
+09/2008 dyn3d/physdem1.F
+Correct possible bugs on some computers + english text & comments
+
+09/2008 dyn3d/defrun_new.F
+Corrected small bug (tapedef argument should not be altered)
+
+09/2008 dyn3d/eofdump.F
+Corrected small bug in the sampling time
+
+09/2008 dyn3d/vanleer.F
+removed inapropriate 'EXTERNAL' statement
+
+O9/2008 aeronomars/euvheat.F
+removed print of 'ngrid','nlayer' (undefined and unused variables)
+
+O9/2008 phymars/mkstat.F90
+Cleaned up (removed unused variables)
+
+09/2008 phymars/uniplot.F bibio/uniplot.F are removed because useless
+
+O9/2008 phymars/comcstfi.h
+09/2008 phymars/writediagfi.F
+08/2007 phymars/tabfi.F
+Cosmetic: english text & comments
+
+09/2008 dyn3d/iniwrite.F => phymars/iniwrite.F 
+- moved it from dyn3d to phymars
+- cleanup & english
+
+09/2008 phymars/physdem1.F
+commented unused (wrong) 'external' qualifiers
+
+09/2008 phymars/orosetup.F
+Replaced logical 'xor' with 'neqv' (for standards compliance)
+
+09/2008 phymars/physiq.F
+- Change to ice albedo to fit new TES observation (Franck L.)
+- cleanup + comments in english (Ehouarn)
+
+09/2008 phymars/watercloud.F
+cleaned version with reduced number of nuclei, and TES like output. 
+
+09/2008 phymars/growthrate.F
+updated to improved version provided by F. Montmessin (Feb 2006...)
+ 
+09/2008 phymars/vdifc.F
+Change to remove an old quick fix to take into account condensation in
+vertical diffusion
+
+09/2008  
+- phymars/mkstat.F90 
+- phymars/def_var.F90
+Slight change to improve compilation (include => #include)
+
+O7/2008 MAJOR CHANGE to solve the "warm polar night" problem
+- dyn3d/addfi.F
+- dyn3d/gcm.F
+Recompute mass after incrementation of ps in addfi.F (masse is added in the argument
+list of addfi). 
+Also : the call of groupeun is now controlled by a flag.  
+
+
+07/2008 phymars/physiq.F
+- Change made (by Aymeric) to improve the models when iradia > 1 (iradia
+  is set is callphys.def). 
+(i.e. when one don't want to compute radiative transfer at every timestep):
+The change has no impact when iradia=1
+=> The net radiative surface flux  "fluxrad" is computed outside the
+test on iradia    
+=> Add intermediate variable "fluxrad_sky" (rad. flux from sky absorbed 
+by surface) to be saved when rad. transfer is not computed.
+- Cleaning of old commented lines about radiative transfer in clouds
+
+
+07/2008 Utility concatnc.F90 (not used by the gcm): improvement in order
+to 1) concatenate 1D variable and 2) increase the number of input files
+up to 1000
+
+12/2007 dyn3d/lect_start_archive.F
+some old stuuf is commented to avoid bugs with some start_archive with
+many tracers
+
+12/2007 dyn3d/dump2d.F
+Changed to improve crazy output to screen with newstart...
+
+12/2007 bibio/mxva.F
+Changed to avoid array bounds problems on some compilers
+
+12/2007 
+aeronomars/conduction.F
+aeronomars/molvis.F
+Changed to ensure variable initialisation on some compilers. 
+
+12/2007 phymars/convadj.F
+Improvement regarding convection in the polar night
+
+12/2007 
+-phymars/callsedim.F  
+-phymars/initracer.F  
+-phymars/tracer.h
+-phymars/callsedim2q.F
+Changed to initialize tracer particle densities in initracer.F, and
+transfer it with common in tracer.h
+
+12/2007 script "makegcm" "create_make_gcm"
+Improved to separate object compiled files in LIBO depending on the type
+of machine you are (32bits or 64bits). Also allow to compute with LMDZ.2
+dynamics.  
+
+10/2007 Utility lslin.F90 (not used by the gcm)
+changed evaluation of 'start_var' from hard-coded values to a computed value
+
+09/2007 phymars/wstats.F90 and phymars/readtesassim.F90
+change " include" to "#include" to force recompilation by makegcm
+
+08/2007 dyn3d/newstart.F90
+Translate french questions asked by program in english...
+Improve options to remove topography... (for intercomparison)
+
+08/2007 phymars/dustopacity.F
+add possibility for scenario iaervar=5 (for paleoclimate...)
+
+08/2007 phymars/tabfi.F
+Translate french questions asked by program in english...
+
+11/2006 phymars/convadj.F
+Correction to remove a bug for compiler which did not "save" variable 
+by default (e.g. pgf on opteron, or g95)
+
+11/2006 aeronomars/thermosphere.F
+Correction to remove a bug if tracer=F and callthermos=T 
+
+10/2006 phymars/readtesassim.F90
+Correction to use dust scenario MY24 at high horizontal resolution
+
+10/2006 dyn3d/dynredem.F and phymars/physdem1.F
+Syntax changes to avoid problem when writing start files when compiling
+on 64bits computers.
+
+05/2006 : phymars/writediagfi.F
+Add a test to be sure not to use writediagfi with 1D version of GCM 
+
+09/2005 : phymars/mucorr.F
+Fix a weird bug that occured when running the model with "Diurnal=F"
+on Linux, pgf90
+
+-----------------------------------
+Update in june 2005
+-----------------------------------
+
+06/2005 : deftank/run0
+Improved so that the "run manager" run0 send you an E-mail when a run crash
+(this can easily be commented)
+
+06/2005 : deftank/run_mcd and deftank/callphys.def
+Cosmetic changes to make thing clearer in english
+
+06/2005 : makegcm
+IMPORTANT FOR 1D GCM VERSION ON LINUX : modified to follow a change in
+the Linux syntax of "sed"
+
+06/2005 : aeronomars/calchim.F
+NEW version provided by S. Lebonnois and F. Lefevre on May, 11, 2005
+("Cosmetic changes")+syntax changes by FF (1 line > 72 characters)
+
+06/2005 : phymars/convadj.F
+Modified if one tracer is CO2, take into account the
+Molecular mass variation (e.g. when CO2 condense) to trigger convection
+
+06/2005 : phymars/newcondens.F
+The CO2 deplation induced convection is now skipped, since it is properly
+done in convadj.F (see above)
+
+06/2005 : phymars/watercloud.F
+IMPORTANT (especially for PALEOCLIMATE) : correction if a lot of subliming 
+CO2 fills the 1st layer. The number of input is modified (see below)
+
+06/2005 : phymars/physiq.F
+- Change in the input argument of watercloud.F : "pdpsrf" is added (see
+above)
+- If "caps", the cold trap at the south pole is only active for
+obliquity lower than 27 deg (only useful for paleoclimate...)
+
+06/2005 : dyn3d/groupe.F
+Just cleaned to remove ambiguity of where "ngroup" is set (in groupeun, and
+not in groupe). 
+
+06/2005 : phymars/inistats.F
+remove declaration of "pseudoalt", already in comvert.h 
+(was a problem only for some compiler)
+
+06/2005 : dyn3d/newstart.F
+Cosmetic changes in some format for compatibility with gnu compiler
+
+06/2005 : dyn3d/readhead_NC.F
+Cosmetic changes in some format for compatibility with gnu compiler
+
+06/2005 : dyn3d/sortvarc.F
+Cosmetic changes in some format for compatibility with gnu compiler
+
+06/2005 : phymars/testphys1d.F
+Cleaning and changes in declaration for compatibility with gnu compiler
+
+06/2005 : phymars/lwflux.F
+Initialisation of "ksidb" for compiler that does not unitialized variable.
+
+06/2005 : phymars/lwxd.F
+Initialisation of "ksi_emis" for compiler that does not unitialized variable.
+
+06/2005 : phymars/nirco2abs.F
+Negligible change of one loop for compiler that does not like loop on real.
+
+06/2005 : aeronomars/paramfoto.F
+Negligible change of one loop for compiler that does not like loop on real.
+
+
+--------------------------------------------------------------
+Change in December 2004 - January 2005 
+--------------------------------------------------------------
+
+01/2005 : aeronomars/euvheat.F
+IMPORTANT FOR THERMOSPHERE : correction to avoid EUV flux time variation 
+and keep it constant throughout the run
+
+01/2005 : LMDZ.MARS/create_make_gcm
+add phymars on line 187 to ensure recompilation of aeronomars subroutines
+when any *.h file is updated in phymars
+
+01/2005 : dyn3d/dynetat0.F
+Add an option to reindex the tracer when using newstart if needed
+
+01/2005 : phymars/phyetat0.F
+Fix a bug when using newstart with more tracers than in old startfi
+Add an option to reindex the surface tracer if needed
+
+12/2004 : phymars/watercloud.F
+IMPORTANT: Several Changes to fix a bug that made watercloud
+to create negative mixing ratio in some conditions (past climate)
+
+12/2004 : dyn3d/interp_horiz.F
+IMPORTANT : add a test to fix a bug when using newstart with 
+a new grid size.
+
+12/2004 : phymars/vdifc.F
+FOR compilation on DEC machine :
+Just change order of declaration of pqsurf(ngrid,nq
+
+12/2004 : phymars/phyetat0.F
+FOR compilation on DEC machine :Just change order of declaration 
+
+12/2004 : dyn3d/addfi.F
+FOR compilation on DEC machine : qtestt = 1.0e-30 (before: qtestt = 1.0e-40)
+
+12/2004 : LMDZ.MARS/makefile
+Add new option to allow compilation on DEC machine
+
+--------------------------------------------------------------
+Change in August-september 2004 : before production of MCD V4.0
+--------------------------------------------------------------
+
+09/2004 : phymars/eofdump.F
+Change to store density in EOF to build Mars Climate Database
+
+09/2004 : phymars/physiq.F
+add rho in the input of eofdump:
+CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
+
+09/2004 : phymars/ineofdump.F
+Change to store hybrid coordinates in EOF to build Mars Climate Database
+
+
+08/2004 : dyn3d/inidissip.F
+Form change to simplify the choice of dissipation parameters
+and use the same code for the thermosphere of only atmosphere model
+
+08/2004 : dyn3d/newstart.F
+Add an option "q=x" to initialise tracer manually
+
+08/2004 : aeronomars/chimiedata.h
+remove mmol,noms (now in tracer.h common)
+
+08/2004 : In aeronomars: calchim.F, concentrations.F, euvheat.F,
+inichim_newstart.F, moldiff.F, moldiffcoeff.F, perosat.F
+Add : #include "tracer.h"  (where tracer names are now)
+
+08/2004 : phymars/inifis.F
+Change some test on callphys.def options to allow any tracer transport...
+
+08/2004 : phymars/initracer
+Initialization of tracer name "noms(iq)" usable in many place
+
+08/2004 : phymars/tracer.h
+add noms,mmol to the common tracer
+
+08/2004 : phymars/newcondens.F
+IMPORTANT : major changes to allow proper treatment of non volatile 
+gaz enrichment when CO2 condense and depletion when CO2 sublime
+
+08/2004 : phymars/writediagfi.F
+remove the "abort" command so that the run does not stop due to a NETCDF bug
+(NETCDF cannot write in file larger than 2098184 kb
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/README
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/README	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/README	(revision 214)
@@ -0,0 +1,29 @@
+> tar xvf LMDZ.MARS.source.tar
+ajouter perso
+change resol in compile
+check NETCDF NCDFLIB NCDFINC
+> compile
+> mkdir myGCM
+> cd myGCM
+lier z2sig.def (voir what_I_need)
+lier un start_archive avec bon nombre traceurs
+> ln -sf ../gcm.e
+> ln -sf ../newstart.e
+> ln -sf ../deftank/run0 
+lier un callphys.def
+lier un run.def
+cp num_run
+> newstart.e
+> mv restart.nc start0.nc
+> mv restartfi.nc startfi0.nc
+LANCER AVEC
+> nohup run_mcd_3days > log_run_mcd_3days &
+ATTENDRE LONGTEMPS ...
+> mv tmpsave ../STARTBASE_64_48_32_t2
+> rm -rf start.nc start0.nc startfi.nc startfi0.nc
+
+lier le run.def simu
+modifier physiq.F pour sortir des donnees 
+recompiler
+link diagfi.nc vers un espace temp
+ex: ln -sf $MODELS/TMPDIR/GCMINI/diagfi.nc
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/compile
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/compile	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/compile	(revision 214)
@@ -0,0 +1,47 @@
+#! /bin/bash
+
+## changer eventuellement le callphys.def egalement
+
+##################
+tracers=2
+nx=64
+ny=48
+nz=32
+##################
+#tracers=1
+#nx=64
+#ny=48
+#nz=25
+##################
+
+
+
+# GCM environment variables
+#--------------------------
+export LMDGCM=$PWD
+export LIBOGCM=$PWD/libo
+
+echo $LMDGCM
+echo $LIBOGCM
+echo $NCDFLIB
+echo $NCDFINC 
+
+#rm -rf libo
+mkdir libo
+touch libo/log_compile_newstart
+touch libo/log_compile_gcm
+
+echo bricolo
+\rm libf/grid/dimensions.h 2> /dev/null
+
+echo compile newstart
+echo check in $PWD/libo/log_compile_newstart
+makegcm -t ${tracers} -p mars -d ${nx}x${ny}x${nz} newstart | tee libo/log_compile_newstart | grep 'warnings' >> libo/log_compile_newstart
+
+echo compile gcm
+echo check in $PWD/libo/log_compile_gcm
+makegcm -t ${tracers} -p mars -d ${nx}x${ny}x${nz} gcm | tee libo/log_compile_gcm | grep 'warnings' >> libo/log_compile_gcm
+
+echo well, OK.
+
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/callphys.def
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/callphys.def	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/callphys.def	(revision 214)
@@ -0,0 +1,94 @@
+General options
+~~~~~~~~~~~~~~~
+tracer    (Run with or without tracer transport ?)
+F
+diurnal   (Diurnal cycle ?  if diurnal=F, diurnal averaged solar heating)
+T
+season    (Seasonal cycle ? if season=F, Ls stays constant like in "start")
+T
+lwrite    (want some more output on the screen ?) 
+F
+stats     (Saving statistics in file "stats" ?)
+T
+calleofdump (Saving EOF profiles in file "profiles" for Climate Database ?)
+F
+Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+~~~~~~~~~~~~~
+iaervar  (=1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario)
+4        (=4 Mars Year 24 from TES assimilation) 
+iddist  (Dust vertical distribution: =0: old distrib. (Pollack90) 
+3       (=1: top set by "topdustref"; =2: Viking scenario; =3 MGS scenario )
+topdustref (Dust top altitude (km). Matter only if iddist=1)
+55.
+Physical Parameterizations :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+callrad   (call radiative transfer ?)
+T
+callnlte (call NLTE radiative schemes ?   matter only if callrad=T)
+T
+callnirco2 (call CO2 NIR absorption ?   matter only if callrad=T)
+T
+calldifv  (call turbulent vertical diffusion ?)
+T
+calladj   (call convective adjustment ?)
+T
+callcond  (call CO2 condensation ?)
+T
+callsoil  (call thermal conduction in the soil ?)
+T
+calllott  (call Lott's gravity wave/subgrid topography scheme ?)
+T
+Radiative transfer options :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+iradia    (the rad.transfer is computed every "iradia" physical timestep)
+1
+callg2d   (Output of the exchange coefficient mattrix ? for diagnostic only)
+F
+rayleigh  (Rayleigh scattering : should be =F for now)
+F
+Tracer (dust water, ice and/or chemical species) options (use if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dustbin    (DUST: Transported dust ? (if >0, uses q(1) to q(dustbin))
+0
+active     (DUST: Radiatively active dust ? (uses q(1) to q(dustbin))
+F
+doubleq    (DUST: needs dustbin=1, use mass (q(1)) and number (q(2)) mixing ratio to predict dust size ?)
+F
+lifting    (DUST: lifted by GCM surface winds ?)
+F
+dustdevil  (DUST: lifted by dust devils ?)
+F
+scavenging (DUST: Scavenging by CO2 snowfall ?)
+F
+sedimentation (DUST/WATERICE: Gravitationnal sedimentation ?)
+F              (set to T to simulate water cycle with clouds)                
+iceparty   (WATERICE: includes water ice: q(nqmx-1)). "water" must be T !
+F
+activice   (WATERICE: Radiatively active transported atmospheric water ice ?)
+F
+water      (WATER: Compute water cycle using q(nqmx) )
+F
+caps       (WATER: current permanent caps at both poles. T IS RECOMMENDED
+T          (caps=T means Ncap is a source of water, S pole is a cold trap)
+photochem  (PHOTOCHEMISTRY: chemical species included)
+F
+Thermospheric options (relevant if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+callthermos  (call thermosphere ?)
+F
+thermoswater  (WATER: included without cycle	only if water=F)
+F
+callconduct  (call thermal conduction ?     matter only if callthermos=T)
+F
+calleuv  (call EUV heating ?                matter only if callthermos=T)
+F
+callmolvis  (call molecular viscosity ?     matter only if callthermos=T)
+F
+callmoldiff  (call molecular diffusion ?    matter only if callthermos=T)
+F
+thermochem  (call thermospheric photochemistry ?  matter only if callthermos=T)
+F
+solarcondate (date for solar flux calculation: 1985 < date < 2002))
+1993.4       (Solar min=1996.4 ave=1993.4 max=1990.6)
+
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.32x24x25
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.32x24x25	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.32x24x25	(revision 214)
@@ -0,0 +1,152 @@
+
+----------------------------------------------------------------------- 
+Parametres de controle du run:                                          
+------------------------------                                          
+
+- Nombre de jours d'integration                                         
+     nday
+ 9999 
+
+- nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step
+ 480
+
+- periode pour le pas Matsuno (en pas)                                  
+  iperiod
+ 5
+
+- periode de sortie des variables de controle (en pas)                  
+  iconser
+ 120
+
+- periode d'ecriture du fichier histoire (en jour)                      
+    iecri
+ 100
+
+- periode de stockage fichier histmoy (en jour)                         
+ periodav
+ 60.
+
+- periode de la dissipation (en pas)                                    
+  idissip
+ 1
+
+- choix de l'operateur de dissipation (star ou  non star )              
+ lstardis
+ T
+
+- avec ou sans coordonnee hybrides
+ hybrid
+ T
+
+- nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv
+ 1
+
+- nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot
+ 2
+
+- nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh
+ 2
+
+- temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv
+ 3000.
+
+- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot
+ 9000.
+
+- temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp
+ 9000.
+
+- coefficient pour gamdissip                                            
+  coefdis
+ 0.
+
+- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats
+ F
+
+- avec ou sans physique                                                 
+   physic
+ T
+
+- periode de la physique (en pas)                                       
+  iphysiq
+ 10
+
+- choix d'une grille reguliere                                          
+  grireg
+  T
+
+- frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy
+ 120
+
+- longitude en degres du centre du zoom                                 
+   clon
+ 63.
+
+- latitude en degres du centre du zoom                                  
+   clat
+ 0.
+
+- facteur de grossissement du zoom,selon longitude                      
+  grossismx
+ 1.
+
+- facteur de grossissement du zoom ,selon latitude                      
+ grossismy
+ 1.
+
+-  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb
+ F
+
+- extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx
+ 0.
+
+- extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy
+ 0.
+
+-  raideur du zoom en  X
+    taux
+ 2.
+
+-  raideur du zoom en  Y
+    tauy
+ 2.
+
+-  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus
+  F
+
+- Avec sponge layer
+  callsponge
+  T
+
+- Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge
+  2
+
+- Sponge:  hauteur de sponge (km)
+  hsponge
+  90
+
+- Sponge:  tetasponge (secondes)
+  tetasponge
+  50000
+
+c   *****************************
+c
+c      Autre reglage possible :
+c
+cc  F ) deltay : deplacement ( en degres ) en  y  de la zone du zoom
+cc                                          dans  *** inigeom.F ***
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.32x24x32
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.32x24x32	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.32x24x32	(revision 214)
@@ -0,0 +1,152 @@
+
+----------------------------------------------------------------------- 
+Parametres de controle du run:                                          
+------------------------------                                          
+
+- Nombre de jours d'integration                                         
+     nday
+ 9999
+
+- nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step
+ 480
+
+- periode pour le pas Matsuno (en pas)                                  
+  iperiod
+ 5
+
+- periode de sortie des variables de controle (en pas)                  
+  iconser
+ 120
+
+- periode d'ecriture du fichier histoire (en jour)                      
+    iecri
+ 100
+
+- periode de stockage fichier histmoy (en jour)                         
+ periodav
+ 60.
+
+- periode de la dissipation (en pas)                                    
+  idissip
+ 1
+
+- choix de l'operateur de dissipation (star ou  non star )              
+ lstardis
+ T
+
+- avec ou sans coordonnee hybrides
+ hybrid
+ T
+
+- nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv
+ 1
+
+- nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot
+ 2
+
+- nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh
+ 2
+
+- temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv
+ 3000.
+
+- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot
+ 9000.
+
+- temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp
+ 9000.
+
+- coefficient pour gamdissip                                            
+  coefdis
+ 0.
+
+- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats
+ F
+
+- avec ou sans physique                                                 
+   physic
+ T
+
+- periode de la physique (en pas)                                       
+  iphysiq
+ 10
+
+- choix d'une grille reguliere                                          
+  grireg
+  T
+
+- frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy
+ 120
+
+- longitude en degres du centre du zoom                                 
+   clon
+ 63.
+
+- latitude en degres du centre du zoom                                  
+   clat
+ 0.
+
+- facteur de grossissement du zoom,selon longitude                      
+  grossismx
+ 1.
+
+- facteur de grossissement du zoom ,selon latitude                      
+ grossismy
+ 1.
+
+-  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb
+ F
+
+- extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx
+ 0.
+
+- extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy
+ 0.
+
+-  raideur du zoom en  X
+    taux
+ 2.
+
+-  raideur du zoom en  Y
+    tauy
+ 2.
+
+-  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus
+  F
+
+- Avec sponge layer
+  callsponge
+  T
+
+- Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge
+  2
+
+- Sponge:  hauteur de sponge (km)
+  hsponge
+  130
+
+- Sponge:  tetasponge (secondes)
+  tetasponge
+  50000
+
+c   *****************************
+c
+c      Autre reglage possible :
+c
+cc  F ) deltay : deplacement ( en degres ) en  y  de la zone du zoom
+cc                                          dans  *** inigeom.F ***
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x25
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x25	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x25	(revision 214)
@@ -0,0 +1,152 @@
+
+----------------------------------------------------------------------- 
+Parametres de controle du run:                                          
+------------------------------                                          
+
+- Nombre de jours d'integration                                         
+     nday
+ 9999
+
+- nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step
+ 960
+
+- periode pour le pas Matsuno (en pas)                                  
+  iperiod
+ 5
+
+- periode de sortie des variables de controle (en pas)                  
+  iconser
+ 120
+
+- periode d'ecriture du fichier histoire (en jour)                      
+    iecri
+ 100
+
+- periode de stockage fichier histmoy (en jour)                         
+ periodav
+ 60.
+
+- periode de la dissipation (en pas)                                    
+  idissip
+ 5
+
+- choix de l'operateur de dissipation (star ou  non star )              
+ lstardis
+ T
+
+- avec ou sans coordonnee hybrides
+ hybrid
+ T
+
+- nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv
+ 1
+
+- nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot
+ 2
+
+- nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh
+ 2
+
+- temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv
+ 4000.
+
+- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot
+ 5000.
+
+- temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp
+ 5000.
+
+- coefficient pour gamdissip                                            
+  coefdis
+ 0.
+
+- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats
+ F
+
+- avec ou sans physique                                                 
+   physic
+ T
+
+- periode de la physique (en pas)                                       
+  iphysiq
+ 20
+
+- choix d'une grille reguliere                                          
+  grireg
+  T
+
+- frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy
+ 240
+
+- longitude en degres du centre du zoom                                 
+   clon
+ 63.
+
+- latitude en degres du centre du zoom                                  
+   clat
+ 0.
+
+- facteur de grossissement du zoom,selon longitude                      
+  grossismx
+ 1.
+
+- facteur de grossissement du zoom ,selon latitude                      
+ grossismy
+ 1.
+
+-  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb
+ F
+
+- extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx
+ 0.
+
+- extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy
+ 0.
+
+-  raideur du zoom en  X
+    taux
+ 2.
+
+-  raideur du zoom en  Y
+    tauy
+ 2.
+
+-  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus
+  F
+
+- Avec sponge layer
+  callsponge
+  T
+
+- Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge
+  2
+
+- Sponge:  hauteur de sponge (km)
+  hsponge
+  90
+
+- Sponge:  tetasponge (secondes)
+  tetasponge
+  50000
+
+c   *****************************
+c
+c      Autre reglage possible :
+c
+cc  F ) deltay : deplacement ( en degres ) en  y  de la zone du zoom
+cc                                          dans  *** inigeom.F ***
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x32
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x32	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x32	(revision 214)
@@ -0,0 +1,152 @@
+
+----------------------------------------------------------------------- 
+Parametres de controle du run:                                          
+------------------------------                                          
+
+- Nombre de jours d'integration                                         
+     nday
+ 9999
+
+- nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step
+ 960
+
+- periode pour le pas Matsuno (en pas)                                  
+  iperiod
+ 5
+
+- periode de sortie des variables de controle (en pas)                  
+  iconser
+ 120
+
+- periode d'ecriture du fichier histoire (en jour)                      
+    iecri
+ 100
+
+- periode de stockage fichier histmoy (en jour)                         
+ periodav
+ 60.
+
+- periode de la dissipation (en pas)                                    
+  idissip
+ 5
+
+- choix de l'operateur de dissipation (star ou  non star )              
+ lstardis
+ T
+
+- avec ou sans coordonnee hybrides
+ hybrid
+ T
+
+- nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv
+ 1
+
+- nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot
+ 2
+
+- nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh
+ 2
+
+- temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv
+ 4000.
+
+- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot
+ 5000.
+
+- temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp
+ 5000.
+
+- coefficient pour gamdissip                                            
+  coefdis
+ 0.
+
+- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats
+ F
+
+- avec ou sans physique                                                 
+   physic
+ T
+
+- periode de la physique (en pas) USUALLY = daystep/48               
+  iphysiq
+ 20
+
+- choix d'une grille reguliere                                          
+  grireg
+  T
+
+- frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy
+ 240
+
+- longitude en degres du centre du zoom                                 
+   clon
+ 63.
+
+- latitude en degres du centre du zoom                                  
+   clat
+ 0.
+
+- facteur de grossissement du zoom,selon longitude                      
+  grossismx
+ 1.
+
+- facteur de grossissement du zoom ,selon latitude                      
+ grossismy
+ 1.
+
+-  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb
+ F
+
+- extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx
+ 0.
+
+- extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy
+ 0.
+
+-  raideur du zoom en  X
+    taux
+ 2.
+
+-  raideur du zoom en  Y
+    tauy
+ 2.
+
+-  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus
+  F
+
+- Avec sponge layer
+  callsponge
+  T
+
+- Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge
+  2
+
+- Sponge:  hauteur de sponge (km)
+  hsponge
+  130
+
+- Sponge:  tetasponge (secondes)
+  tetasponge
+  50000
+
+c   *****************************
+c
+c      Autre reglage possible :
+c
+cc  F ) deltay : deplacement ( en degres ) en  y  de la zone du zoom
+cc                                          dans  *** inigeom.F ***
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x50
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x50	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run.def.64x48x50	(revision 214)
@@ -0,0 +1,152 @@
+
+----------------------------------------------------------------------- 
+Parametres de controle du run:                                          
+------------------------------                                          
+
+- Nombre de jours d'integration                                         
+     nday
+ 1
+
+- nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step
+ 960
+
+- periode pour le pas Matsuno (en pas)                                  
+  iperiod
+ 5
+
+- periode de sortie des variables de controle (en pas)                  
+  iconser
+ 120
+
+- periode d'ecriture du fichier histoire (en jour)                      
+    iecri
+ 100
+
+- periode de stockage fichier histmoy (en jour)                         
+ periodav
+ 60.
+
+- periode de la dissipation (en pas)                                    
+  idissip
+ 1
+
+- choix de l'operateur de dissipation (star ou  non star )              
+ lstardis
+ T
+
+- avec ou sans coordonnee hybrides
+ hybrid
+ T
+
+- nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv
+ 1
+
+- nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot
+ 2
+
+- nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh
+ 2
+
+- temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv
+ 3000
+
+- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot
+ 20000
+
+- temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp
+ 20000
+
+- coefficient pour gamdissip                                            
+  coefdis
+ 0.
+
+- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats
+ F
+
+- avec ou sans physique                                                 
+   physic
+ T
+
+- periode de la physique (en pas)                                       
+  iphysiq
+ 20
+
+- choix d'une grille reguliere                                          
+  grireg
+  T
+
+- frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy
+ 240
+
+- longitude en degres du centre du zoom                                 
+   clon
+ 63.
+
+- latitude en degres du centre du zoom                                  
+   clat
+ 0.
+
+- facteur de grossissement du zoom,selon longitude                      
+  grossismx
+ 1.
+
+- facteur de grossissement du zoom ,selon latitude                      
+ grossismy
+ 1.
+
+-  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb
+ F
+
+- extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx
+ 0.
+
+- extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy
+ 0.
+
+-  raideur du zoom en  X
+    taux
+ 2.
+
+-  raideur du zoom en  Y
+    tauy
+ 2.
+
+-  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus
+  F
+
+- Avec sponge layer
+  callsponge
+  T
+
+- Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge
+  2
+
+- Sponge:  hauteur de sponge (km)
+  hsponge
+  230
+
+- Sponge:  tetasponge (secondes)
+  tetasponge
+  30000
+
+c   *****************************
+c
+c      Autre reglage possible :
+c
+cc  F ) deltay : deplacement ( en degres ) en  y  de la zone du zoom
+cc                                          dans  *** inigeom.F ***
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run0
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run0	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run0	(revision 214)
@@ -0,0 +1,102 @@
+#!/bin/csh
+###########################################################################
+# Script to perform several chained LMD Mars GCM simulations
+# SET HERE the maximum total number of simulations
+
+set nummax=1
+
+###########################################################################
+
+
+echo "---------------------------------------------------------"
+echo "starting run0"
+
+set dir=`pwd`
+set machine=`hostname`
+set address=`whoami`
+
+# Look for file "num_run" which should contain 
+# the value of the previously computed season
+# (defaults to 0 if file "num_run" does not exist)
+if ( -f num_run ) then
+ echo "found file num_run"
+  set numold=`cat num_run`
+else
+  set numold=0
+endif
+echo "numold is set to" ${numold}
+
+
+# Set value of current season 
+set numnew=${numold}
+@ numnew = ${numnew} + 1
+echo "numnew is set to" ${numnew}
+
+# Look for initialization data files (exit if none found)
+if ( ( -f start${numold}.nc ) && ( -f startfi${numold}.nc ) ) then
+   \cp -f start${numold}.nc start.nc
+   \cp -f startfi${numold}.nc startfi.nc
+else if ( ${numold} == 99999 ) then
+    echo "No run because previous run crashed ! (99999 in num_run)"
+    exit
+else
+   echo "Where is file start"${numold}".nc??"
+   exit
+endif
+
+
+# Run GCM
+nice gcm.e >! lrun${numnew} 
+
+
+# Check if run ended normaly and copy datafiles
+if ( ( -f restartfi.nc ) && ( -f restart.nc ) ) then
+  echo "Run seems to have ended normaly"
+  \mv -f restartfi.nc startfi${numnew}.nc
+  \mv -f restart.nc start${numnew}.nc
+else
+  if ( -f num_run ) then
+    \mv -f num_run num_run.crash
+  else
+    echo "No file num_run to build num_run.crash from !!"
+    # Impose a default value of 0 for num_run
+    echo 0 >! num_run.crash
+  endif
+ echo 99999 >! num_run
+############## To receive an Email message if the run crash ########
+mail -s "crash run GCM" $address <<ENDMAIL
+The run on $machine in $dir has just crashed.
+ENDMAIL
+############################################"" 
+ exit
+endif
+
+# Copy other datafiles that may have been generated
+if ( -f diagfi.nc ) then
+  \mv -f diagfi.nc diagfi${numnew}.nc
+endif
+if ( -f stats.nc ) then
+  \mv -f stats.nc stats${numnew}.nc
+endif
+if ( -f profiles.dat ) then
+  \mv -f profiles.dat profiles${numnew}.dat
+  \mv -f profiles.hdr profiles${numnew}.hdr
+endif
+if ( -f stats.nc ) then
+  \mv -f stats.nc stats${numnew}.nc
+endif
+
+
+# Prepare things for upcoming runs by writing
+# value of computed season in file num_run
+echo ${numnew} >! num_run
+
+# If we are over nummax : stop
+if ( `expr $numnew + 1 \> $nummax ` == 1) then
+   exit
+else
+   \cp -f run0 exe_mars
+   exe_mars
+endif 
+
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run_farm
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run_farm	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run_farm	(revision 214)
@@ -0,0 +1,28 @@
+# Script to submit simulations with the LMD MARS GCM on 
+# a farm of computers using the run manager named "Torque".
+# To run it, just type 
+# > qsub run_farm
+
+# Shell
+#PBS -S /bin/bash
+# Job Name
+#PBS -N run_mars
+# Merge output and error logs
+#PBS -j oe
+
+####################################
+# SELECT HERE THE RUN DURATION CLASS :
+# Queue std(6h) / day / week /infini
+
+#PBS -q std
+
+# Go to present working directory
+cd  $PBS_O_WORKDIR
+
+# SELECT HERE THE TYPE OF SCRIPT TO USE : "run0" OR "run_mcd"
+
+run0
+
+####################################
+
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run_mcd
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run_mcd	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/run_mcd	(revision 214)
@@ -0,0 +1,78 @@
+#!/bin/csh
+##############################################################################
+# Script to run successive monthly (a month spans 30 deg. in Ls)
+# simulations of the LMD Mars GCM.
+# This script uses the following files:
+# - "run0": script to compute a given month
+# - "num_run": single line file containing a value, which is the value of
+#   the 'initializing' month
+# - "run.def.ref": a "run.def" file, exept that variable "nday" is set to 9999
+# - all files usually neaded to run the GCM. Note that initialization
+#   datafiles "start.nc" and "startfi.nc" must be
+#   "start${num_run}.nc" and "startfi${num_run}.nc"
+##############################################################################
+# DON'T FORGET TO:
+# 1) Set the "number of simulation" to 1 in "run0"
+# 2) Copy "run.def" to "run.def.ref" and SET nday to 9999 in run.def.ref
+# 3) uncomment the various lines below (those which begin with "sed ...")
+# 4) If restarting in the middle of a year: 
+#       -Below, Comment all the month already computed
+#       -Write  "N" in file "num_run" if starting from startN.nc !
+##############################################################################
+
+
+# Remove old "error" file (if any) and create a brand new (blank) one
+\rm -f  error; touch error
+
+# 12 martian months (=30°deg Ls) of unequal length (excentric orbit)
+  sed s/9999/61/ run.def.ref >! run.def ; run0 >>& error     #1
+  sed s/9999/66/ run.def.ref >! run.def ; run0 >>& error     #2
+  sed s/9999/66/ run.def.ref >! run.def ; run0 >>& error     #3
+  sed s/9999/65/ run.def.ref >! run.def ; run0 >>& error     #4
+  sed s/9999/60/ run.def.ref >! run.def ; run0 >>& error     #5
+  sed s/9999/54/ run.def.ref >! run.def ; run0 >>& error     #6
+  sed s/9999/50/ run.def.ref >! run.def ; run0 >>& error     #7
+  sed s/9999/46/ run.def.ref >! run.def ; run0 >>& error     #8
+  sed s/9999/47/ run.def.ref >! run.def ; run0 >>& error     #9
+  sed s/9999/47/ run.def.ref >! run.def ; run0 >>& error     #10
+  sed s/9999/51/ run.def.ref >! run.def ; run0 >>& error     #11 
+  sed s/9999/56/ run.def.ref >! run.def ; run0 >>& error     #12
+#
+# 12 martian months (=30°deg Ls) of unequal length (excentric orbit)
+# sed s/9999/61/ run.def.ref >! run.def ; run0 >>& error     #13
+# sed s/9999/66/ run.def.ref >! run.def ; run0 >>& error     #14
+# sed s/9999/66/ run.def.ref >! run.def ; run0 >>& error     #15
+# sed s/9999/65/ run.def.ref >! run.def ; run0 >>& error     #16
+# sed s/9999/60/ run.def.ref >! run.def ; run0 >>& error     #17
+# sed s/9999/54/ run.def.ref >! run.def ; run0 >>& error     #18
+# sed s/9999/50/ run.def.ref >! run.def ; run0 >>& error     #19
+# sed s/9999/46/ run.def.ref >! run.def ; run0 >>& error     #20
+# sed s/9999/47/ run.def.ref >! run.def ; run0 >>& error     #21
+# sed s/9999/47/ run.def.ref >! run.def ; run0 >>& error     #22
+# sed s/9999/51/ run.def.ref >! run.def ; run0 >>& error     #23
+# sed s/9999/56/ run.def.ref >! run.def ; run0 >>& error     #24
+#
+# 12 martian months (=30°deg Ls) of unequal length (excentric orbit)
+# sed s/9999/61/ run.def.ref >! run.def ; run0 >>& error     #25
+# sed s/9999/66/ run.def.ref >! run.def ; run0 >>& error     #26
+# sed s/9999/66/ run.def.ref >! run.def ; run0 >>& error     #27
+# sed s/9999/65/ run.def.ref >! run.def ; run0 >>& error     #28
+# sed s/9999/60/ run.def.ref >! run.def ; run0 >>& error     #29
+# sed s/9999/54/ run.def.ref >! run.def ; run0 >>& error     #30
+# sed s/9999/50/ run.def.ref >! run.def ; run0 >>& error     #31
+# sed s/9999/46/ run.def.ref >! run.def ; run0 >>& error     #32
+# sed s/9999/47/ run.def.ref >! run.def ; run0 >>& error     #33
+# sed s/9999/47/ run.def.ref >! run.def ; run0 >>& error     #34
+# sed s/9999/51/ run.def.ref >! run.def ; run0 >>& error     #35
+# sed s/9999/56/ run.def.ref >! run.def ; run0 >>& error     #36
+
+# These last lines allow to reconstruct the last num_run after a model crash 
+# (because after crash, num_run was set to 9999 so that no further simulation
+# Could be performed
+   set number=`cat num_run`
+   if ( $number == 99999) then
+      set numold=`cat num_run.crash`
+      echo '-----' >>& error
+      echo 'Seasonal simulation aborted after starting from start'$numold >>& error
+      \mv -f num_run.crash num_run  
+   endif
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/testphys1d.def
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/testphys1d.def	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/testphys1d.def	(revision 214)
@@ -0,0 +1,36 @@
+0      ! Initial day (martian sol; = 0 at Ls=0) 
+0      ! initial local time  (between 0 and 24) 
+48     ! number of time-step per day.
+100      ! Run duration (sols) 
+610.   ! psurf: Surface pressure (Pa) 
+0.2    ! tauref: reference dust opacity at 700 Pa (true tau~tauref*psurf/700)
+0.    ! latitude  (degrees) 
+0.2    ! ground albedo  (under the ice)
+400    ! ground thermal inertia  (SI)
+10.    ! composante vers l'est du vent geostrophique (u) 
+0.     ! composante vers le nord du vent geostrophique (v)
+0.     ! initial CO2 ice mass (kg.m-2)
+T      ! Hybrid vertical coordinate ? T=Hybrid , F=sigma 
+
+Initial temperature profile
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+5       ! ichoice :  Kind of initial profile (see below)
+200     ! tref : reference temperature (K) (used if ichoice=1)
+0       ! isin (if isin=1, add a perturbation to profile)
+26.5220 ! pic  pic perturbation gauss pour ichoice = 6 ou 7
+10      ! largeur de la perturbation gauss pour ichoice = 6 ou 7
+30.     ! hauteur de la perturbation gauss pour ichoice = 6 ou 7 
+
+c=======================================================================
+c   differents profils d'atmospheres. T=f(z)
+c         ichoice=1   Temperature constante T=tref
+c         ichoice=2   profil Savidjari (comme Seiff mais avec dT/dz=cte)
+c         ichoice=3   Lindner (profil polaire)
+c         ichoice=4   inversion
+c         ichoice=5   Seiff  (profile standard base sur les entrees Viking)
+c         ichoice=6    T constante + perturbation gauss (level)
+c         ichoice=7    T constante + perturbation gauss (km)
+c         ichoice=8   Read in an ascii file "profile" 
+
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+CECI EST LE FICHIER .DEF UTILISE EN 1D PAR TESTPHYS1D.F
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/z2sig.def
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/z2sig.def	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/deftank/z2sig.def	(revision 214)
@@ -0,0 +1,68 @@
+10.00000     H: atmospheric scale height (km) (used as a reference only)
+0.0040       Typical pseudo-altitude (m) for 1st layer (z=H*log(sigma))
+0.018       ,, ,,  ,, ,, ,, ,,  ,, ,, ,,  2nd layer, etc...
+0.0400
+0.1000
+0.228200
+0.460400
+0.907000
+1.73630
+3.19040
+5.54010
+8.97780
+13.5138
+18.9666
+25.0626
+31.5527
+38.4369
+45.4369
+52.4369
+59.4369
+66.4369
+73.4369
+80.4369
+87.4369
+94.4369
+101.4369
+108.437
+115.437
+122.437
+129.437
+136.437
+143.437
+150.437
+157.437
+164.437
+171.437
+178.437
+185.437
+192.437
+199.437
+206.437
+213.437
+220.437
+227.437
+234.437
+241.437
+248.437
+255.437
+262.437
+269.437
+276.437
+283.437
+290.437
+297.437
+304.437
+311.437
+318.437
+325.437
+332.437
+339.437
+346.437
+353.437
+360.437
+367.437
+374.437
+381.437
+388.437
+395.437
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/callphys.def
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/callphys.def	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/callphys.def	(revision 214)
@@ -0,0 +1,94 @@
+General options
+~~~~~~~~~~~~~~~
+tracer    (Run with or without tracer transport ?)
+F
+diurnal   (Diurnal cycle ?  if diurnal=F, diurnal averaged solar heating)
+T
+season    (Seasonal cycle ? if season=F, Ls stays constant like in "start")
+T
+lwrite    (want some more output on the screen ?) 
+F
+stats     (Saving statistics in file "stats" ?)
+F
+calleofdump (Saving EOF profiles in file "profiles" for Climate Database ?)
+F
+Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+~~~~~~~~~~~~~
+iaervar  (=1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario)
+4        (=4 Mars Year 24 from TES assimilation) 
+iddist  (Dust vertical distribution: =0: old distrib. (Pollack90) 
+3       (=1: top set by "topdustref"; =2: Viking scenario; =3 MGS scenario )
+topdustref (Dust top altitude (km). Matter only if iddist=1)
+55.
+Physical Parameterizations :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+callrad   (call radiative transfer ?)
+T
+callnlte (call NLTE radiative schemes ?   matter only if callrad=T)
+T
+callnirco2 (call CO2 NIR absorption ?   matter only if callrad=T)
+T
+calldifv  (call turbulent vertical diffusion ?)
+T
+calladj   (call convective adjustment ?)
+T
+callcond  (call CO2 condensation ?)
+T
+callsoil  (call thermal conduction in the soil ?)
+T
+calllott  (call Lott's gravity wave/subgrid topography scheme ?)
+T
+Radiative transfer options :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+iradia    (the rad.transfer is computed every "iradia" physical timestep)
+1
+callg2d   (Output of the exchange coefficient mattrix ? for diagnostic only)
+F
+rayleigh  (Rayleigh scattering : should be =F for now)
+F
+Tracer (dust water, ice and/or chemical species) options (use if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dustbin    (DUST: Transported dust ? (if >0, uses q(1) to q(dustbin))
+0
+active     (DUST: Radiatively active dust ? (uses q(1) to q(dustbin))
+F
+doubleq    (DUST: needs dustbin=1, use mass (q(1)) and number (q(2)) mixing ratio to predict dust size ?)
+F
+lifting    (DUST: lifted by GCM surface winds ?)
+F
+dustdevil  (DUST: lifted by dust devils ?)
+F
+scavenging (DUST: Scavenging by CO2 snowfall ?)
+F
+sedimentation (DUST/WATERICE: Gravitationnal sedimentation ?)
+F              (set to T to simulate water cycle with clouds)                
+iceparty   (WATERICE: includes water ice: q(nqmx-1)). "water" must be T !
+F
+activice   (WATERICE: Radiatively active transported atmospheric water ice ?)
+F
+water      (WATER: Compute water cycle using q(nqmx) )
+F
+caps       (WATER: current permanent caps at both poles. T IS RECOMMENDED
+F          (caps=T means Ncap is a source of water, S pole is a cold trap)
+photochem  (PHOTOCHEMISTRY: chemical species included)
+F
+Thermospheric options (relevant if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+callthermos  (call thermosphere ?)
+F
+thermoswater  (WATER: included without cycle	only if water=F)
+F
+callconduct  (call thermal conduction ?     matter only if callthermos=T)
+F
+calleuv  (call EUV heating ?                matter only if callthermos=T)
+F
+callmolvis  (call molecular viscosity ?     matter only if callthermos=T)
+F
+callmoldiff  (call molecular diffusion ?    matter only if callthermos=T)
+F
+thermochem  (call thermospheric photochemistry ?  matter only if callthermos=T)
+F
+solarcondate (date for solar flux calculation: 1985 < date < 2002))
+1993.4       (Solar min=1996.4 ave=1993.4 max=1990.6)
+
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/callphys.def.water
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/callphys.def.water	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/callphys.def.water	(revision 214)
@@ -0,0 +1,94 @@
+General options
+~~~~~~~~~~~~~~~
+tracer    (Run with or without tracer transport ?)
+T
+diurnal   (Diurnal cycle ?  if diurnal=F, diurnal averaged solar heating)
+T
+season    (Seasonal cycle ? if season=F, Ls stays constant like in "start")
+T
+lwrite    (want some more output on the screen ?) 
+F
+stats     (Saving statistics in file "stats" ?)
+F
+calleofdump (Saving EOF profiles in file "profiles" for Climate Database ?)
+F
+Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+~~~~~~~~~~~~~
+iaervar  (=1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario)
+4        (=4 Mars Year 24 from TES assimilation) 
+iddist  (Dust vertical distribution: =0: old distrib. (Pollack90) 
+3       (=1: top set by "topdustref"; =2: Viking scenario; =3 MGS scenario )
+topdustref (Dust top altitude (km). Matter only if iddist=1)
+55.
+Physical Parameterizations :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+callrad   (call radiative transfer ?)
+T
+callnlte (call NLTE radiative schemes ?   matter only if callrad=T)
+T
+callnirco2 (call CO2 NIR absorption ?   matter only if callrad=T)
+T
+calldifv  (call turbulent vertical diffusion ?)
+T
+calladj   (call convective adjustment ?)
+T
+callcond  (call CO2 condensation ?)
+T
+callsoil  (call thermal conduction in the soil ?)
+T
+calllott  (call Lott's gravity wave/subgrid topography scheme ?)
+T
+Radiative transfer options :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+iradia    (the rad.transfer is computed every "iradia" physical timestep)
+1
+callg2d   (Output of the exchange coefficient mattrix ? for diagnostic only)
+F
+rayleigh  (Rayleigh scattering : should be =F for now)
+F
+Tracer (dust water, ice and/or chemical species) options (use if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dustbin    (DUST: Transported dust ? (if >0, uses q(1) to q(dustbin))
+0
+active     (DUST: Radiatively active dust ? (uses q(1) to q(dustbin))
+F
+doubleq    (DUST: needs dustbin=1, use mass (q(1)) and number (q(2)) mixing ratio to predict dust size ?)
+F
+lifting    (DUST: lifted by GCM surface winds ?)
+F
+dustdevil  (DUST: lifted by dust devils ?)
+F
+scavenging (DUST: Scavenging by CO2 snowfall ?)
+F
+sedimentation (DUST/WATERICE: Gravitationnal sedimentation ?)
+T              (set to T to simulate water cycle with clouds)                
+iceparty   (WATERICE: includes water ice: q(nqmx-1)). "water" must be T !
+T
+activice   (WATERICE: Radiatively active transported atmospheric water ice ?)
+F
+water      (WATER: Compute water cycle using q(nqmx) )
+T
+caps       (WATER: current permanent caps at both poles. T IS RECOMMENDED
+T          (caps=T means Ncap is a source of water, S pole is a cold trap)
+photochem  (PHOTOCHEMISTRY: chemical species included)
+F
+Thermospheric options (relevant if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+callthermos  (call thermosphere ?)
+F
+thermoswater  (WATER: included without cycle	only if water=F)
+F
+callconduct  (call thermal conduction ?     matter only if callthermos=T)
+F
+calleuv  (call EUV heating ?                matter only if callthermos=T)
+F
+callmolvis  (call molecular viscosity ?     matter only if callthermos=T)
+F
+callmoldiff  (call molecular diffusion ?    matter only if callthermos=T)
+F
+thermochem  (call thermospheric photochemistry ?  matter only if callthermos=T)
+F
+solarcondate (date for solar flux calculation: 1985 < date < 2002))
+1993.4       (Solar min=1996.4 ave=1993.4 max=1990.6)
+
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/launch_gcm
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/launch_gcm	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/launch_gcm	(revision 214)
@@ -0,0 +1,68 @@
+#! /bin/bash
+
+##############################
+# use: - launch_gcm          #
+#      - echo 2 | launch_gcm #
+##############################
+
+# Author : A. Spiga - Spring/Winter 2008
+
+################## penser a recompiler si on change ceci
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_32_t2'
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/old_STARTBASE_64_48_25'
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_25'
+startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_32_t2'
+##################
+
+#
+# get simulation start day and compute start/startfi required date (simulation day - 1)
+#
+echo '>>>' choose start sol ? ; read daygcm
+  #if [[ "$daygcm" == 0 ]]
+  #then
+  #daygcm=669 
+  #fi
+daygcm=$(expr ${daygcm} - 1)
+
+#
+# find the correct start/startfi files among 3-day startbase
+#
+ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null
+if [[ "$?" == 0 ]]
+then
+    echo found ! ${daygcm}
+else
+    #echo not found ! ${daygcm}	
+    daygcm=$(expr ${daygcm} - 1)
+    ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null
+    if [[ "$?" == 0 ]]     # probably a 'while' loop would be better if more than 3
+    then
+	echo found ! ${daygcm}
+    else
+	#echo not found ! ${daygcm}    
+	daygcm=$(expr ${daygcm} - 1)
+        ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null	
+        echo found ! ${daygcm}
+    fi	
+fi	
+
+#
+# link start/startfi files
+#
+ln -sf ${startbase}/start${daygcm}.nc start.nc
+ln -sf ${startbase}/startfi${daygcm}.nc startfi.nc
+
+#
+# clean the place
+#
+rm -rf restart.nc 2> /dev/null 
+rm -rf restartfi.nc 2> /dev/null
+rm -rf log_gcm 2> /dev/null
+
+#
+# run !
+#
+echo running the GCM - check in ${PWD}/log_gcm 
+nohup gcm.e > log_gcm
+
+exit
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/num_run.STARTBASE
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/num_run.STARTBASE	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/num_run.STARTBASE	(revision 214)
@@ -0,0 +1,1 @@
+0
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run.def.INIT
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run.def.INIT	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run.def.INIT	(revision 214)
@@ -0,0 +1,152 @@
+
+----------------------------------------------------------------------- 
+Parametres de controle du run:                                          
+------------------------------                                          
+
+- Nombre de jours d'integration                                         
+     nday
+ 7 5 10 30 10 40
+
+- nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step
+ 960
+
+- periode pour le pas Matsuno (en pas)                                  
+  iperiod
+ 5
+
+- periode de sortie des variables de controle (en pas)                  
+  iconser
+ 40
+
+- periode d'ecriture du fichier histoire (en jour)                      
+    iecri
+ 100
+
+- periode de stockage fichier histmoy (en jour)                         
+ periodav
+ 60.
+
+- periode de la dissipation (en pas)                                    
+  idissip
+ 5
+
+- choix de l'operateur de dissipation (star ou  non star )              
+ lstardis
+ T
+
+- avec ou sans coordonnee hybrides
+ hybrid
+ T
+
+- nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv
+ 1
+
+- nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot
+ 2
+
+- nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh
+ 2
+
+- temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv
+ 4000.
+
+- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot
+ 5000.
+
+- temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp
+ 5000.
+
+- coefficient pour gamdissip                                            
+  coefdis
+ 0.
+
+- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats
+ F
+
+- avec ou sans physique                                                 
+   physic
+ T
+
+- periode de la physique (en pas)                                       
+  iphysiq
+ 20
+
+- choix d'une grille reguliere                                          
+  grireg
+  T
+
+- frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy
+ 40 80 40
+
+- longitude en degres du centre du zoom                                 
+   clon
+ 63.
+
+- latitude en degres du centre du zoom                                  
+   clat
+ 0.
+
+- facteur de grossissement du zoom,selon longitude                      
+  grossismx
+ 1.
+
+- facteur de grossissement du zoom ,selon latitude                      
+ grossismy
+ 1.
+
+-  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb
+ F
+
+- extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx
+ 0.
+
+- extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy
+ 0.
+
+-  raideur du zoom en  X
+    taux
+ 2.
+
+-  raideur du zoom en  Y
+    tauy
+ 2.
+
+-  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus
+  F
+
+- Avec sponge layer
+  callsponge
+  T
+
+- Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge
+  2
+
+- Sponge:  hauteur de sponge (km)
+  hsponge
+  90
+
+- Sponge:  tetasponge (secondes)
+  tetasponge
+  50000
+
+c   *****************************
+c
+c      Autre reglage possible :
+c
+cc  F ) deltay : deplacement ( en degres ) en  y  de la zone du zoom
+cc                                          dans  *** inigeom.F ***
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run.def.STARTBASE
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run.def.STARTBASE	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run.def.STARTBASE	(revision 214)
@@ -0,0 +1,152 @@
+
+----------------------------------------------------------------------- 
+Parametres de controle du run:                                          
+------------------------------                                          
+
+- Nombre de jours d'integration                                         
+     nday
+ 3
+
+- nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step
+ 960
+
+- periode pour le pas Matsuno (en pas)                                  
+  iperiod
+ 5
+
+- periode de sortie des variables de controle (en pas)                  
+  iconser
+ 40
+
+- periode d'ecriture du fichier histoire (en jour)                      
+    iecri
+ 100
+
+- periode de stockage fichier histmoy (en jour)                         
+ periodav
+ 60.
+
+- periode de la dissipation (en pas)                                    
+  idissip
+ 5
+
+- choix de l'operateur de dissipation (star ou  non star )              
+ lstardis
+ T
+
+- avec ou sans coordonnee hybrides
+ hybrid
+ T
+
+- nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv
+ 1
+
+- nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot
+ 2
+
+- nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh
+ 2
+
+- temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv
+ 4000.
+
+- temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot
+ 5000.
+
+- temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp
+ 5000.
+
+- coefficient pour gamdissip                                            
+  coefdis
+ 0.
+
+- choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats
+ F
+
+- avec ou sans physique                                                 
+   physic
+ T
+
+- periode de la physique (en pas)                                       
+  iphysiq
+ 20
+
+- choix d'une grille reguliere                                          
+  grireg
+  T
+
+- frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy
+ 100000
+
+- longitude en degres du centre du zoom                                 
+   clon
+ 63.
+
+- latitude en degres du centre du zoom                                  
+   clat
+ 0.
+
+- facteur de grossissement du zoom,selon longitude                      
+  grossismx
+ 1.
+
+- facteur de grossissement du zoom ,selon latitude                      
+ grossismy
+ 1.
+
+-  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb
+ F
+
+- extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx
+ 0.
+
+- extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy
+ 0.
+
+-  raideur du zoom en  X
+    taux
+ 2.
+
+-  raideur du zoom en  Y
+    tauy
+ 2.
+
+-  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus
+  F
+
+- Avec sponge layer
+  callsponge
+  T
+
+- Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge
+  2
+
+- Sponge:  hauteur de sponge (km)
+  hsponge
+  90
+
+- Sponge:  tetasponge (secondes)
+  tetasponge
+  50000
+
+c   *****************************
+c
+c      Autre reglage possible :
+c
+cc  F ) deltay : deplacement ( en degres ) en  y  de la zone du zoom
+cc                                          dans  *** inigeom.F ***
+
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run_mcd_3days
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run_mcd_3days	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/run_mcd_3days	(revision 214)
@@ -0,0 +1,60 @@
+#! /bin/bash
+
+#
+# - please link or copy start0.nc and startfi.nc (or start_archive.nc)
+# - please set 3 in the first parameter of run.def
+# - please set 0 in the first line of num_run 
+# - please set daylim below (669 for the whole year)
+#
+
+##########
+daylim=669
+##########
+
+rm -rf tmpsave
+mkdir tmpsave
+
+runnum=0
+runnum_obsolete=0
+daynum=0
+
+while [ ${daynum} -lt ${daylim} ]  
+do
+	# run GCM during three days
+	echo run starting at day ${daynum}
+	echo '****'
+	time run0 
+
+	# save generated starting files 
+	echo saving files
+	echo '****'
+	runnum=$(expr ${runnum} + 1)
+        save=$(expr ${daynum} + 2)
+	cp start${runnum}.nc tmpsave/start${save}.nc
+        cp startfi${runnum}.nc tmpsave/startfi${save}.nc
+	mv lrun${runnum} tmpsave/lrun_for_start${save}
+
+	# increment day to day+3 for while loop assessement
+	daynum=$(expr ${daynum} + 3)
+
+        # get rid of what was generated two steps before
+	echo cleaning
+	echo '****'
+	if [ ${runnum} -gt 1 ]	#because we want to keep start0.nc and startfi0.nc
+	then
+	runnum_obsolete=$(expr ${runnum} - 1)
+	rm -rf start${runnum_obsolete}.nc
+	rm -rf startfi${runnum_obsolete}.nc
+	rm -rf diagfi${runnum_obsolete}.nc
+	fi
+done
+
+# clean the last files
+rm -rf start${runnum}.nc
+rm -rf startfi${runnum}.nc
+rm -rf diagfi${runnum}.nc
+
+# end
+echo -- end --
+echo please check results in file tmpsave
+echo if you want to keep the results, rename this folder
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/what_I_need.info
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/what_I_need.info	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/myGCM/what_I_need.info	(revision 214)
@@ -0,0 +1,9 @@
+total 24K
+-rw-r--r--  1 aslmd lmdjus 3.4K Mar 11 11:50 callphys.def
+-rw-r--r--  1 aslmd lmdjus 3.4K Mar 11 11:45 callphys.def~
+lrwxrwxrwx  1 aslmd lmdjus   18 Mar 11 11:43 gcm.e -> ../LMDZ.MARS/gcm.e
+lrwxrwxrwx  1 aslmd lmdjus   23 Mar 11 11:43 newstart.e -> ../LMDZ.MARS/newstart.e
+-rw-r--r--  1 aslmd lmdjus 3.2K Mar 11 11:49 run.def
+lrwxrwxrwx  1 aslmd lmdjus   57 Mar 11 11:44 start_archive.nc -> /u/forget/WWW/datagcm/Starts/start_archive_64x48x25_t2.nc
+-rw-r--r--  1 aslmd lmdjus    0 Mar 11 11:53 what_I_need
+lrwxrwxrwx  1 aslmd lmdjus   30 Mar 11 11:52 z2sig.def -> ../LMDZ.MARS/deftank/z2sig.def
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/physiq.F
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/physiq.F	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/physiq.F	(revision 214)
@@ -0,0 +1,1351 @@
+      SUBROUTINE physiq(ngrid,nlayer,nq,
+     $            firstcall,lastcall,
+     $            pday,ptime,ptimestep,
+     $            pplev,pplay,pphi,
+     $            pu,pv,pt,pq,
+     $            pw,
+     $            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
+
+
+      IMPLICIT NONE
+c=======================================================================
+c
+c   subject:
+c   --------
+c
+c   Organisation of the physical parametrisations of the LMD 
+c   martian atmospheric general circulation model.
+c
+c   The GCM can be run without or with tracer transport
+c   depending on the value of Logical "tracer" in file  "callphys.def"
+c   Tracers may be water vapor, ice OR chemical species OR dust particles
+c
+c   SEE comments in initracer.F about numbering of tracer species...
+c
+c   It includes:
+c
+c      1. Initialization:
+c      1.1 First call initializations
+c      1.2 Initialization for every call to physiq
+c      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
+c      2. Compute radiative transfer tendencies
+c         (longwave and shortwave) for CO2 and dust.
+c      3. Gravity wave and subgrid scale topography drag :
+c      4. Vertical diffusion (turbulent mixing):
+c      5. Convective adjustment
+c      6. Condensation and sublimation of carbon dioxide.
+c      7.  TRACERS :
+c       7a. water and water ice
+c       7b. call for photochemistry when tracers are chemical species
+c       7c. other scheme for tracer (dust) transport (lifting, sedimentation)
+c       7d. updates (CO2 pressure variations, surface budget)
+c      8. Contribution to tendencies due to thermosphere
+c      9. Surface and sub-surface temperature calculations
+c     10. Write outputs :
+c           - "startfi", "histfi" (if it's time)
+c           - Saving statistics (if "callstats = .true.")
+c           - Dumping eof (if "calleofdump = .true.")
+c           - Output any needed variables in "diagfi" 
+c      10. Diagnostic: mass conservation of tracers
+c 
+c   author: 
+c   ------- 
+c           Frederic Hourdin	15/10/93
+c           Francois Forget		1994
+c           Christophe Hourdin	02/1997 
+c           Subroutine completly rewritten by F.Forget (01/2000)
+c           Introduction of the photochemical module: S. Lebonnois (11/2002)
+c           Introduction of the thermosphere module: M. Angelats i Coll (2002)
+c           Water ice clouds: Franck Montmessin (update 06/2003)
+c           
+c
+c
+c   arguments:
+c   ----------
+c
+c   input:
+c   ------
+c    ecri                  period (in dynamical timestep) to write output
+c    ngrid                 Size of the horizontal grid.
+c                          All internal loops are performed on that grid.
+c    nlayer                Number of vertical layers.
+c    nq                    Number of advected fields
+c    firstcall             True at the first call
+c    lastcall              True at the last call
+c    pday                  Number of days counted from the North. Spring
+c                          equinoxe.
+c    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
+c    ptimestep             timestep (s)
+c    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
+c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
+c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
+c    pu(ngrid,nlayer)      u component of the wind (ms-1)
+c    pv(ngrid,nlayer)      v component of the wind (ms-1)
+c    pt(ngrid,nlayer)      Temperature (K)
+c    pq(ngrid,nlayer,nq)   Advected fields
+c    pudyn(ngrid,nlayer)    \ 
+c    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
+c    ptdyn(ngrid,nlayer)     / corresponding variables
+c    pqdyn(ngrid,nlayer,nq) /
+c    pw(ngrid,?)           vertical velocity
+c
+c   output:
+c   -------
+c
+c    pdu(ngrid,nlayermx)        \
+c    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
+c    pdt(ngrid,nlayermx)         /  variables due to physical processes.
+c    pdq(ngrid,nlayermx)        /
+c    pdpsrf(ngrid)             /
+c    tracerdyn                 call tracer in dynamical part of GCM ?
+
+c
+c=======================================================================
+c
+c    0.  Declarations :
+c    ------------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comgeomfi.h"
+#include "surfdat.h"
+#include "comdiurn.h"
+#include "callkeys.h"
+#include "comcstfi.h"
+#include "planete.h"
+#include "comsaison.h"
+#include "control.h"
+#include "dimradmars.h"
+#include "comg1d.h"
+#include "tracer.h"
+#include "nlteparams.h"
+
+#include "chimiedata.h"
+#include "watercap.h"
+#include "fisice.h"
+#include "param.h"
+#include "param_v3.h"
+#include "conc.h"
+
+#include "netcdf.inc"
+
+
+
+c Arguments :
+c -----------
+
+c   inputs:
+c   -------
+      INTEGER ngrid,nlayer,nq
+      REAL ptimestep
+      REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
+      REAL pphi(ngridmx,nlayer)
+      REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
+      REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
+      REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
+      REAL zh(ngridmx,nlayermx)      ! potential temperature (K)
+      LOGICAL firstcall,lastcall
+
+      REAL pday
+      REAL ptime 
+      logical tracerdyn
+
+c   outputs:
+c   --------
+c     physical tendencies
+      REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
+      REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
+      REAL pdpsrf(ngridmx) ! surface pressure tendency
+
+
+c Local saved variables:
+c ----------------------
+c     aerosol (dust or ice) extinction optical depth  at reference wavelength 
+c     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
+c      aerosol optical properties  :
+      REAL aerosol(ngridmx,nlayermx,naerkind) 
+
+      INTEGER day_ini  ! Initial date of the run (sol since Ls=0) 
+      INTEGER icount     ! counter of calls to physiq during the run.
+      REAL tsurf(ngridmx)            ! Surface temperature (K)
+      REAL tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
+      REAL co2ice(ngridmx)           ! co2 ice surface layer (kg.m-2)  
+      REAL albedo(ngridmx,2)         ! Surface albedo in each solar band
+      REAL emis(ngridmx)             ! Thermal IR surface emissivity
+      REAL dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
+      REAL fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
+      REAL fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
+      REAL capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
+      REAL fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
+      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
+      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
+      INTEGER ig_vl1                 ! Grid Point near VL1   (for diagnostic) 
+
+      SAVE day_ini, icount
+      SAVE aerosol, tsurf,tsoil
+      SAVE co2ice,albedo,emis, q2
+      SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
+      SAVE ig_vl1
+
+      REAL stephan   
+      DATA stephan/5.67e-08/  ! Stephan Boltzman constant
+      SAVE stephan
+
+c Local variables :
+c -----------------
+
+
+      CHARACTER*80 fichier 
+      INTEGER l,ig,ierr,igout,iq, tapphys
+      INTEGER iqmin                     ! Used if iceparty engaged
+
+      REAL fluxsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2)
+      REAL fluxsurf_sw(ngridmx,2)    !incident SW (solar) surface flux (W.m-2)
+      REAL fluxtop_lw(ngridmx)       !Outgoing LW (IR) flux to space (W.m-2)
+      REAL fluxtop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
+c     for clear area (uncovered by clouds) :
+      REAL clsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2) 
+      REAL clsurf_sw(ngridmx,2)    !incident SW (solar) surface flux (W.m-2)
+      REAL cltop_lw(ngridmx)       !Outgoing LW (IR) flux to space (W.m-2)
+      REAL cltop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
+      REAL tauref(ngridmx)           ! Reference column optical depth at 700 Pa
+                                     ! (used if active=F) 
+      REAL tau(ngridmx,naerkind)     ! Column dust optical depth at each point
+      REAL zls                       !  solar longitude (rad)
+      REAL zday                      ! date (time since Ls=0, in martian days)
+      REAL zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
+      REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
+      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
+
+
+c     Tendancies due to various processes:
+      REAL dqsurf(ngridmx,nqmx)
+      REAL zdtlw(ngridmx,nlayermx)     ! (K/s)
+      REAL zdtsw(ngridmx,nlayermx)     ! (K/s)
+      REAL cldtlw(ngridmx,nlayermx)     ! (K/s) LW heating rate for clear area
+      REAL cldtsw(ngridmx,nlayermx)     ! (K/s) SW heating rate for clear area
+      REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
+      REAL zdtnlte(ngridmx,nlayermx)   ! (K/s)
+      REAL zdtsurf(ngridmx)            ! (K/s)
+      REAL zdtcloud(ngridmx,nlayermx)
+      REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
+      REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhadj(ngridmx,nlayermx)                           ! (K/s)
+      REAL zdtgw(ngridmx,nlayermx)                            ! (K/s)
+      REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
+      REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
+      REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
+
+      REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
+      REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
+      REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
+      REAL zdqadj(ngridmx,nlayermx,nqmx)
+      REAL zdqc(ngridmx,nlayermx,nqmx)
+      REAL zdqcloud(ngridmx,nlayermx,nqmx)
+      REAL zdqscloud(ngridmx,nqmx)
+      REAL zdqchim(ngridmx,nlayermx,nqmx)
+      REAL zdqschim(ngridmx,nqmx)
+
+      REAL zdteuv(ngridmx,nlayermx)    ! (K/s)
+      REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
+      REAL zdumolvis(ngridmx,nlayermx)
+      REAL zdvmolvis(ngridmx,nlayermx)
+      real zdqmoldiff(ngridmx,nlayermx,nqmx)
+
+c     Local variable for local intermediate calcul:
+      REAL zflubid(ngridmx)
+      REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
+      REAL zdum1(ngridmx,nlayermx)
+      REAL zdum2(ngridmx,nlayermx)
+      REAL ztim1,ztim2,ztim3, z1,z2
+      REAL ztime_fin
+      REAL zdh(ngridmx,nlayermx)
+      REAL pclc_min ! minimum of the cloud fraction over the whole domain
+      INTEGER length
+      PARAMETER (length=100)
+
+c local variables only used for diagnostic (output in file "diagfi" or "stats")
+c -----------------------------------------------------------------------------
+      REAL ps(ngridmx), zt(ngridmx,nlayermx)
+      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
+      REAL zq(ngridmx,nlayermx,nqmx)
+      REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
+      character*2 str2
+      character*5 str5
+      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
+      real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
+      real qtot1,qtot2 ! total aerosol mass
+      integer igmin, lmin
+      logical tdiag
+
+
+      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
+      REAL zstress(ngrid), cd
+      real hco2(nqmx),tmean, zlocal(nlayermx)
+      real rho(ngridmx,nlayermx) ! density
+      real vmr(ngridmx,nlayermx) ! volume mixing ratio
+
+      REAL time_phys
+     
+c=======================================================================
+
+c 1. Initialisation:
+c -----------------
+
+c  1.1   Initialisation only at first call
+c  ---------------------------------------
+      IF (firstcall) THEN
+
+c        variables set to 0
+c        ~~~~~~~~~~~~~~~~~~
+         call zerophys(ngrid*nlayer*naerkind,aerosol)
+         call zerophys(ngrid*nlayer,dtrad)
+         call zerophys(ngrid,fluxrad)
+
+c        read startfi 
+c        ~~~~~~~~~~~~
+
+! Read netcdf initial physical parameters.
+         CALL phyetat0 ("startfi.nc",0,0,
+     &         nsoilmx,nq,
+     &         day_ini,time_phys,
+     &         tsurf,tsoil,emis,q2,qsurf,co2ice)
+
+         if (pday.ne.day_ini) then
+           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
+     &                "physics and dynamics"
+           write(*,*) "dynamics day: ",pday
+           write(*,*) "physics day:  ",day_ini
+           stop
+         endif
+
+         write (*,*) 'In physiq day_ini =', day_ini
+
+c        Initialize albedo and orbital calculation
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         CALL surfini(ngrid,co2ice,qsurf,albedo)
+         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+
+c        initialize soil 
+c        ~~~~~~~~~~~~~~~
+         IF (callsoil) THEN
+            CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
+     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+         ELSE
+            PRINT*,'WARNING! Thermal conduction in the soil turned off'
+            DO ig=1,ngrid
+               capcal(ig)=1.e5
+               fluxgrd(ig)=0.
+            ENDDO
+         ENDIF
+         icount=1
+
+
+c        initialize tracers
+c        ~~~~~~~~~~~~~~~~~~
+         tracerdyn=tracer
+         IF (tracer) THEN
+            CALL initracer(qsurf,co2ice)
+         ENDIF  ! end tracer
+
+c        Determining gridpoint near Viking Lander 1 (used for diagnostic only)
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if(ngrid.ne.1) then
+           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.)
+           write(*,*) 'Viking Lander 1 GCM point: lat,lon',
+     &              lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi
+         end if 
+
+c        Initialize thermospheric parameters
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if (callthermos) call param_read
+
+c        Initialize R and Cp as constant
+
+         if (.not.callthermos .and. .not.photochem) then
+                 do l=1,nlayermx
+                  do ig=1,ngridmx
+                   rnew(ig,l)=r
+                   cpnew(ig,l)=cpp
+                   mmean(ig,l)=mugaz
+                   enddo
+                  enddo  
+         endif         
+                   
+      ENDIF        !  (end of "if firstcall")
+
+c ---------------------------------------------------
+c 1.2   Initializations done at every physical timestep:
+c ---------------------------------------------------
+c
+      IF (ngrid.NE.ngridmx) THEN
+         PRINT*,'STOP in PHYSIQ'
+         PRINT*,'Probleme de dimensions :'
+         PRINT*,'ngrid     = ',ngrid
+         PRINT*,'ngridmx   = ',ngridmx
+         STOP
+      ENDIF
+
+c     Initialize various variables
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      call zerophys(ngrid*nlayer, pdv)
+      call zerophys(ngrid*nlayer, pdu)
+      call zerophys(ngrid*nlayer, pdt)
+      call zerophys(ngrid*nlayer*nq, pdq)
+      call zerophys(ngrid, pdpsrf)
+      call zerophys(ngrid, zflubid)
+      call zerophys(ngrid, zdtsurf)
+      call zerophys(ngrid*nq, dqsurf)
+      igout=ngrid/2+1 
+
+
+      zday=pday+ptime ! compute time, in sols (and fraction thereof)
+
+c     Compute Solar Longitude (Ls) :
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (season) then
+         call solarlong(zday,zls)
+      else
+         call solarlong(float(day_ini),zls)
+      end if
+
+c     Compute geopotential at interlayers
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c     ponderation des altitudes au niveau des couches en dp/p
+
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zzlay(ig,l)=pphi(ig,l)/g
+         ENDDO
+      ENDDO
+      DO ig=1,ngrid
+         zzlev(ig,1)=0.
+         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+      ENDDO
+      DO l=2,nlayer
+         DO ig=1,ngrid
+            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
+            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
+            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+         ENDDO
+      ENDDO
+
+
+!     Potential temperature calculation not the same in physiq and dynamic
+
+c     Compute potential temperature
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO l=1,nlayer
+         DO ig=1,ngrid 
+            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
+            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c    1.2.5 Compute mean mass, cp, and R
+c    --------------------------------
+
+      if(photochem.or.callthermos) then
+         call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
+      endif
+c-----------------------------------------------------------------------
+c    2. Compute radiative tendencies :
+c------------------------------------
+
+
+      IF (callrad) THEN
+         IF( MOD(icount-1,iradia).EQ.0) THEN
+
+c          Local Solar zenith angle
+c          ~~~~~~~~~~~~~~~~~~~~~~~~
+           CALL orbite(zls,dist_sol,declin)
+
+           IF(diurnal) THEN
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
+     s         ztim1,ztim2,ztim3, mu0,fract)
+
+           ELSE
+               CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
+           ENDIF
+
+c          NLTE cooling from CO2 emission
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+           IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
+
+c          Find number of layers for LTE radiation calculations
+           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
+     &          CALL nlthermeq(ngrid,nlayer,pplev,pplay)
+
+
+c          Atmospheric dust opacity and aerosol distribution:
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           CALL dustopacity(ngrid,nlayer,nq,zday,pplay,pplev,zls,pq,
+     $      tauref,tau,aerosol)
+         
+c          Call main radiative transfer scheme
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c          Transfer through dust and CO2, except NIR CO2 absorption
+
+           CALL callradite(icount,ngrid,nlayer,aerosol,albedo,
+     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
+     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw)
+
+c          CO2 near infrared absorption
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           call zerophys(ngrid*nlayer,zdtnirco2)
+           if (callnirco2) then
+              call nirco2abs (ngrid,nlayer,pplay,dist_sol,
+     .                       mu0,fract,declin, zdtnirco2)
+           endif
+
+c          Radiative flux from the sky absorbed by the surface (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           DO ig=1,ngrid
+               fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
+     $         +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
+     $         +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
+           ENDDO
+
+
+c          Net atmospheric radiative heating rate (K.s-1)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           IF(callnlte) THEN
+              CALL blendrad(ngrid, nlayer, pplay,
+     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
+           ELSE
+              DO l=1,nlayer
+                 DO ig=1,ngrid
+                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
+     &                          +zdtnirco2(ig,l)
+                  ENDDO
+              ENDDO
+           ENDIF
+
+
+
+        ENDIF ! of if(mod(icount-1,iradia).eq.0)
+
+
+
+c    Transformation of the radiative tendencies:
+c    -------------------------------------------
+
+c          Net radiative surface flux (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+           DO ig=1,ngrid
+               zplanck(ig)=tsurf(ig)*tsurf(ig)
+               zplanck(ig)=emis(ig)*
+     $         stephan*zplanck(ig)*zplanck(ig)
+               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
+           ENDDO
+
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
+            ENDDO
+         ENDDO
+
+      ENDIF ! of IF (callrad)
+
+c-----------------------------------------------------------------------
+c    3. Gravity wave and subgrid scale topography drag :
+c    -------------------------------------------------
+
+
+      IF(calllott)THEN
+
+        CALL calldrag_noro(ngrid,nlayer,ptimestep,
+     &                 pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
+ 
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
+          ENDDO
+        ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c    4. Vertical diffusion (turbulent mixing):
+c    -----------------------------------------
+c
+      IF (calldifv) THEN
+
+
+         DO ig=1,ngrid
+            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
+         ENDDO
+
+         CALL zerophys(ngrid*nlayer,zdum1)
+         CALL zerophys(ngrid*nlayer,zdum2)
+         do l=1,nlayer
+            do ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            enddo
+         enddo
+
+c        Calling vdif (Martian version WITH CO2 condensation)
+         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
+     $        ptimestep,capcal,lwrite,
+     $        pplay,pplev,zzlay,zzlev,z0,
+     $        pu,pv,zh,pq,tsurf,emis,qsurf,
+     $        zdum1,zdum2,zdh,pdq,zflubid,
+     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
+     &        zdqdif,zdqsdif)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
+               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
+
+               zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
+
+            ENDDO
+         ENDDO
+
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
+         ENDDO
+
+         if (tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+           DO iq=1, nq
+              DO ig=1,ngrid
+                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+              ENDDO
+           ENDDO
+         end if ! of if (tracer)
+
+      ELSE    
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+
+     s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+         ENDDO
+      ENDIF ! of IF (calldifv)
+
+
+c-----------------------------------------------------------------------
+c   5. Dry convective adjustment:
+c   -----------------------------
+
+      IF(calladj) THEN
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            ENDDO
+         ENDDO
+         CALL zerophys(ngrid*nlayer,zduadj)
+         CALL zerophys(ngrid*nlayer,zdvadj)
+         CALL zerophys(ngrid*nlayer,zdhadj)
+
+         CALL convadj(ngrid,nlayer,nq,ptimestep,
+     $                pplay,pplev,zpopsk,
+     $                pu,pv,zh,pq,
+     $                pdu,pdv,zdh,pdq,
+     $                zduadj,zdvadj,zdhadj,
+     $                zdqadj)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
+               pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
+
+               zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
+            ENDDO
+         ENDDO
+
+         if(tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         end if
+      ENDIF ! of IF(calladj)
+
+c-----------------------------------------------------------------------
+c   6. Carbon dioxide condensation-sublimation:
+c   -------------------------------------------
+
+      IF (callcond) THEN
+         CALL newcondens(ngrid,nlayer,nq,ptimestep,
+     $              capcal,pplay,pplev,tsurf,pt,
+     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
+     $              co2ice,albedo,emis,
+     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
+     $	            fluxsurf_sw) 
+
+         DO l=1,nlayer
+           DO ig=1,ngrid
+             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
+             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
+             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
+           ENDDO
+         ENDDO
+         DO ig=1,ngrid
+            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
+         ENDDO
+
+         IF (tracer) THEN
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         ENDIF ! of IF (tracer)
+
+      ENDIF  ! of IF (callcond)
+
+c-----------------------------------------------------------------------
+c   7. Specific parameterizations for tracers 
+c:   -----------------------------------------
+
+      if (tracer) then 
+
+c   7a. Water and ice
+c     ---------------
+
+c        ---------------------------------------
+c        Water ice condensation in the atmosphere
+c        ----------------------------------------
+         IF (water) THEN
+           call watercloud(ngrid,nlayer, ptimestep,
+     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
+     &                pq,pdq,zdqcloud,qsurf,zdqscloud,zdtcloud,
+     &                nq,naerkind,tau,icount,zls)
+
+           if (activice) then
+c Temperature variation due to latent heat release
+           DO l=1,nlayer
+             DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
+             ENDDO
+           ENDDO
+           endif
+
+           IF (iceparty) then
+             iqmin=nq-1
+           ELSE
+             iqmin=nq
+           ENDIF
+
+           DO iq=iqmin,nq
+             DO l=1,nlayer
+                DO ig=1,ngrid
+                   pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
+                ENDDO
+             ENDDO
+             DO ig=1,ngrid
+                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
+             ENDDO
+           ENDDO
+
+         END IF  ! of IF (water)
+
+c   7b. Chemical species
+c     ------------------
+
+c        --------------
+c        photochemistry :
+c        --------------
+         IF (photochem .or. thermochem) then
+          call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
+     .      zzlay,zday,pq,pdq,zdqchim,zdqschim,zdqcloud,zdqscloud)
+           
+c Photochemistry includes condensation of H2O2
+
+           do iq=nqchem_min,nq
+            if (noms(iq).eq."h2o2") then
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
+               ENDDO
+             ENDDO
+            else
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
+               ENDDO
+             ENDDO
+            endif
+           ENDDO
+           do iq=nqchem_min,nq
+            if (noms(iq).eq."h2o2") then
+               DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
+               ENDDO
+            else
+               DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
+               ENDDO
+            endif
+           ENDDO
+
+         END IF  ! of IF (photochem.or.thermochem)
+
+c   7c. Aerosol particles
+c     -------------------
+
+c        ----------
+c        Dust devil :
+c        ----------
+         IF(callddevil) then 
+           call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
+     &                zdqdev,zdqsdev)
+ 
+           if (dustbin.ge.1) then
+              do iq=1,nq
+                 DO l=1,nlayer
+                    DO ig=1,ngrid
+                       pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
+                    ENDDO
+                 ENDDO
+              enddo
+              do iq=1,nq
+                 DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
+                 ENDDO
+              enddo
+           endif  ! of if (dustbin.ge.1)
+
+         END IF ! of IF (callddevil)
+
+c        ------------- 
+c        Sedimentation :   acts also on water ice
+c        ------------- 
+         IF (sedimentation) THEN 
+           call zerophys(ngrid*nlayer*nq, zdqsed)
+           call zerophys(ngrid*nq, zdqssed)
+
+           if(doubleq) then
+            call callsedim2q(ngrid,nlayer, ptimestep,
+     &                pplev,zzlev, pt,
+     &                pq, pdq, zdqsed, zdqssed,nq)
+           else
+            call callsedim(ngrid,nlayer, ptimestep,
+     &                pplev,zzlev, pt,
+     &                pq, pdq, zdqsed, zdqssed,nq)
+           end if
+
+
+           DO iq=1, nq
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
+               ENDDO
+             ENDDO
+           ENDDO
+           DO iq=1, nq
+             DO ig=1,ngrid
+                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
+             ENDDO
+           ENDDO
+         END IF   ! of IF (sedimentation)
+
+c   7d. Updates
+c     ---------
+
+        DO iq=1, nq
+          DO ig=1,ngrid
+
+c       ---------------------------------
+c       Updating tracer budget on surface
+c       ---------------------------------
+            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
+
+          ENDDO  ! (ig)
+        ENDDO    ! (iq)
+
+      endif !  of if (tracer) 
+
+
+c-----------------------------------------------------------------------
+c   8. THERMOSPHERE CALCULATION
+c-----------------------------------------------------------------------
+
+      if (callthermos) then
+        call thermosphere(pplev,pplay,dist_sol,
+     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
+     &     pt,pq,pu,pv,pdt,pdq,
+     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
+c           do iq=nqchem_min,nq
+c           write(*,*) 'thermo iq,pq',iq,pq(690,1,iq)
+c           enddo
+
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
+     &                         +zdteuv(ig,l)
+            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
+            DO iq=1, nq
+              pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
+            ENDDO
+          ENDDO
+        ENDDO
+
+
+      endif
+c-----------------------------------------------------------------------
+c   9. Surface  and sub-surface soil temperature
+c-----------------------------------------------------------------------
+c
+c
+c   9.1 Increment Surface temperature:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      DO ig=1,ngrid
+         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
+      ENDDO
+
+c  Prescribe a cold trap at south pole (except at high obliquity !!)
+c  Temperature at the surface is set there to be the temperature
+c  corresponding to equilibrium temperature between phases of CO2
+
+      IF (tracer.AND.water.AND.ngridmx.NE.1) THEN
+         if (caps.and.(obliquit.lt.27.)) then
+           ! NB: Updated surface pressure, at grid point 'ngrid', is
+           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
+           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
+     &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
+         endif
+c       -------------------------------------------------------------
+c       Change of surface albedo (set to 0.4) in case of ground frost
+c       everywhere except on the north permanent cap and in regions
+c       covered by dry ice. 
+c              ALWAYS PLACE these lines after newcondens !!!
+c       -------------------------------------------------------------
+         do ig=1,ngrid
+
+c       -------------- Version temporaire fit TES 2008 ------------
+         if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
+              albedo(ig,1)=0.45                                     
+              albedo(ig,2)=0.45                                     
+          endif
+
+c         if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
+c           if (.not.watercaptag(ig)) then
+c             albedo(ig,1)=0.4
+c             albedo(ig,2)=0.4
+c           endif
+c         endif
+c       -------------- version Francois ---------------
+c          if (co2ice(ig).eq.0.and.
+c    &        ((qsurf(ig,nqmx).gt.0.005).or.(watercaptag(ig)))) then
+c              albedo(ig,1)=max(0.4,albedodat(ig))
+c              albedo(ig,2)=albedo(ig,1)
+c          endif
+c       ---------------------------------------------
+         enddo  ! of do ig=1,ngrid
+      ENDIF  ! of IF (tracer.AND.water.AND.ngridmx.NE.1)
+
+c
+c   9.2 Compute soil temperatures and subsurface heat flux:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (callsoil) THEN
+         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
+     &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+      ENDIF
+
+c-----------------------------------------------------------------------
+c  10. Write output files
+c  ----------------------
+
+c    -------------------------------
+c    Dynamical fields incrementation
+c    -------------------------------
+c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
+      ! temperature, zonal and meridional wind
+      DO l=1,nlayer
+        DO ig=1,ngrid
+          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
+        ENDDO
+      ENDDO
+
+      ! tracers
+      DO iq=1, nq
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
+          ENDDO
+        ENDDO
+      ENDDO
+
+      ! surface pressure
+      DO ig=1,ngrid
+          ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
+      ENDDO
+
+      ! pressure
+      DO l=1,nlayer
+        DO ig=1,ngrid
+             zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
+             zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
+        ENDDO
+      ENDDO
+
+      ! Density 
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
+         ENDDO
+      ENDDO
+
+c    Compute surface stress : (NB: z0 is a common in planete.h)
+c     DO ig=1,ngrid
+c        cd = (0.4/log(zzlay(ig,1)/z0))**2
+c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
+c     ENDDO
+
+c     Sum of fluxes in solar spectral bands (for output only)
+      DO ig=1,ngrid
+	     fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
+	     fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
+      ENDDO
+c ******* TEST ******************************************************
+      ztim1 = 999
+      DO l=1,nlayer
+        DO ig=1,ngrid
+           if (pt(ig,l).lt.ztim1) then
+               ztim1 = pt(ig,l)
+               igmin = ig
+               lmin = l 
+           end if
+        ENDDO
+      ENDDO
+      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then	
+        write(*,*) 'PHYSIQ: stability WARNING :'
+        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
+     &              'ig l =', igmin, lmin
+      end if
+c *******************************************************************
+
+c     ---------------------
+c     Outputs to the screen 
+c     ---------------------
+
+      IF (lwrite) THEN
+         PRINT*,'Global diagnostics for the physics'
+         PRINT*,'Variables and their increments x and dx/dt * dt'
+         WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
+         WRITE(*,'(2f10.5,2f15.5)')
+     s   tsurf(igout),zdtsurf(igout)*ptimestep,
+     s   pplev(igout,1),pdpsrf(igout)*ptimestep
+         WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
+         WRITE(*,'(i4,6f10.5)') (l,
+     s   pu(igout,l),pdu(igout,l)*ptimestep,
+     s   pv(igout,l),pdv(igout,l)*ptimestep,
+     s   pt(igout,l),pdt(igout,l)*ptimestep,
+     s   l=1,nlayer)
+      ENDIF ! of IF (lwrite)
+
+
+      IF (ngrid.NE.1) THEN
+         print*,'Ls =',zls*180./pi,
+     &   ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
+     &   ' tau(Viking1) =',tau(ig_vl1,1)
+
+c        -------------------------------------------------------------------
+c        Writing NetCDF file  "RESTARTFI" at the end of the run
+c        -------------------------------------------------------------------
+c        Note: 'restartfi' is stored just before dynamics are stored
+c              in 'restart'. Between now and the writting of 'restart',
+c              there will have been the itau=itau+1 instruction and
+c              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
+c              thus we store for time=time+dtvr
+
+         IF(lastcall) THEN
+            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) 
+            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
+     .              ptimestep,pday,
+     .              ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
+     .              area,albedodat,inertiedat,zmea,zstd,zsig,
+     .              zgam,zthe)
+         ENDIF
+
+c        -----------------------------------------------------------------
+c        Saving statistics :
+c        -----------------------------------------------------------------
+c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+c        which can later be used to make the statistic files of the run:
+c        "stats")          only possible in 3D runs !
+
+         
+         IF (callstats) THEN
+
+
+            call wstats(ngrid,"ps","Surface pressure","K",2,ps)
+            call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+            call wstats(ngrid,"co2ice","CO2 ice cover",
+     .                  "kg.m-2",2,co2ice)
+c            call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
+c    .                  emis)
+             call wstats(ngrid,"fluxsurf_lw",
+     .                 "Thermal IR radiative flux to surface","W.m-2",2,
+     .                 fluxsurf_lw)
+             call wstats(ngrid,"fluxsurf_sw",
+     .                  "Solar radiative flux to surface","W.m-2",2,
+     .                   fluxsurf_sw_tot)
+             call wstats(ngrid,"fluxtop_lw",
+     .                  "Thermal IR radiative flux to space","W.m-2",2,
+     .                  fluxtop_lw)
+             call wstats(ngrid,"fluxtop_sw",
+     .                  "Solar radiative flux to space","W.m-2",2,
+     .                  fluxtop_sw_tot)
+            call wstats(ngrid,"dod","Dust optical depth"," ",2,tau)
+
+            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
+            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
+            call wstats(ngrid,"v","Meridional (North-South) wind",
+     .                  "m.s-1",3,zv)
+            call wstats(ngrid,"w","Vertical (down-up) wind",
+     .                  "m.s-1",3,pw)
+            call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
+             call wstats(ngrid,"q2",
+     .               "Boundary layer eddy kinetic energy","m2.s-2",3,q2)
+
+           if (tracer) then
+            if (water) then
+               vmr=zq(1:ngridmx,1:nlayermx,nqmx)*mugaz/mmol(nqmx)
+               call wstats(ngrid,"vmr_h2ovapor",
+     .                    "H2O vapor volume mixing ratio","mol/mol",
+     .                    3,vmr)
+               if (iceparty) then
+                  vmr=zq(1:ngridmx,1:nlayermx,nqmx-1)*mugaz/mmol(nqmx-1)
+                  call wstats(ngrid,"vmr_h2oice",
+     .                       "H2O ice volume mixing ratio","mol/mol",
+     .                       3,vmr)
+               endif
+            endif
+
+             if (thermochem.or.photochem) then
+                do iq=1,nq
+                   if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
+     .               (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
+     .               (noms(iq).eq."h2").or.
+     .               (noms(iq).eq."o3")) then
+                        do l=1,nlayer
+                          do ig=1,ngrid
+                            vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
+                          end do
+                        end do
+                        call wstats(ngrid,"vmr_"//trim(noms(iq)),
+     .                     "Volume mixing ratio","mol/mol",3,vmr)
+                   endif
+                enddo
+             endif
+           endif !tracer 
+
+            IF(lastcall) THEN
+              write (*,*) "Writing stats..."
+              call mkstats(ierr)
+            ENDIF
+          ENDIF !if callstats
+
+c        (Store EOF for Mars Climate database software)
+         IF (calleofdump) THEN
+            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
+         ENDIF
+
+
+c        ----------------------------------------------------------------------
+c        output in netcdf file "DIAGFI", containing any variable for diagnostic
+c        (output with  period "ecritphy", set in "run.def")
+c        ----------------------------------------------------------------------
+c        WRITEDIAGFI can ALSO be called from any other subroutines
+c        for any variables !!
+       call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     .                  emis)
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     .                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     .                  co2ice)
+c       call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
+c    .       fluxsurf_lw)
+c       call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
+c    .       fluxsurf_sw_tot)
+c       call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
+c    .       fluxtop_lw)
+c       call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
+c    .       fluxtop_sw_tot)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+c        call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
+c       call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
+c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
+c       call WRITEDIAGFI(ngridm,'Teta','T potentielle','K',3,zh)
+c       call WRITEDIAGFI(ngridm,'Pression','Pression','Pa',3,pplay)
+        call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil)
+
+
+c       OUTPUT of tracer mass mixing ratio and surface value :    
+       if (tracer) then
+c         (for photochemistry, outputs are done in calchim)
+          do iq=1,nqmx
+            write(str2(1:2),'(i2.2)') iq
+            call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq),
+     &                         'kg/kg',3,zq(1,1,iq))
+            call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
+     &                     'kg.m-2',2,qsurf(1,iq))
+          end do
+        end if
+c ***************************************************
+
+c  Outputs for H2O
+       if (tracer) then
+        if (activice) then
+c         call WRITEDIAGFI(ngridmx,'tauice','tau','SI',2,tau(1,2))
+c         call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
+c     &                   'w.m-2',3,zdtsw)
+c         call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
+c     &                   'w.m-2',3,zdtlw)
+        endif  !(activice)
+
+        if (water.and..not.photochem) then
+           iq=nq
+c          write(str2(1:2),'(i2.2)') iq
+c          call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
+c     &                    'kg.m-2',2,zdqscloud(1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
+c     &                    'kg/kg',3,zdqchim(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
+c     &                    'kg/kg',3,zdqdif(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
+c     &                    'kg/kg',3,zdqadj(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
+c     &                    'kg/kg',3,zdqc(1,1,iq))
+        endif  !(water.and..not.photochem)
+
+c       if (iceparty) then
+c          iq=nq-1
+c          write(str2(1:2),'(i2.2)') iq
+c          call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
+c    &                     'kg/kg',3,zq(1,1,iq))
+c       endif  !(iceparty)
+       endif
+
+c  Outputs for dust tracers
+
+         if (tracer.and.(dustbin.ne.0)) then
+            call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
+            if (doubleq) then
+               call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
+     &                         'kg.m-2',2,qsurf(1,1))
+               call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
+     &                         'N.m-2',2,qsurf(1,2))
+               call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
+     &                         'kg.m-2.s-1',2,zdqsdev(1,1))
+               call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
+     &                         'kg.m-2.s-1',2,zdqssed(1,1))
+               do  l=1,nlayer
+                  do ig=1, ngrid
+                     reff(ig,l)= ref_r0 *
+     &               (r3n_q*pq(ig,l,1)/max(pq(ig,l,2),0.01))**(1./3.)
+                     reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6)
+                  end do
+               end do
+               call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff)
+            else
+               do iq=1,dustbin
+                  write(str2(1:2),'(i2.2)') iq
+                  call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
+     &                            'kg/kg',3,zq(1,1,iq))
+                  call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
+     &                            'kg.m-2',2,qsurf(1,iq))
+               end do
+            endif ! (doubleq)
+         end if  ! (tracer.and.(dustbin.ne.0))
+
+
+      ELSE     ! if(ngrid.eq.1)
+
+         print*,'Ls =',zls*180./pi,
+     &  '  tauref(700 Pa) =',tauref,' local tau =',tau(1,1)
+c      ----------------------------------------------------------------------
+c      Output in grads file "g1d" (ONLY when using testphys1d)
+c      (output at every X physical timestep)
+c      ----------------------------------------------------------------------
+c
+c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
+c        CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
+         CALL writeg1d(ngrid,1,ps,'ps','Pa')
+         
+         CALL writeg1d(ngrid,nlayer,zt,'T','K')
+c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
+
+         if(tracer) then
+c           CALL writeg1d(ngrid,1,tau,'tau','SI')
+            do iq=1,nq
+              CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 
+            end do
+         end if
+
+         zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
+
+         do l=2,nlayer-1
+            tmean=zt(1,l)
+            if(zt(1,l).ne.zt(1,l-1))
+     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
+              zlocal(l)= zlocal(l-1)
+     &        -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g
+         enddo
+         zlocal(nlayer)= zlocal(nlayer-1)-
+     &                   log(pplay(1,nlayer)/pplay(1,nlayer-1))*
+     &                   rnew(1,nlayer)*tmean/g
+
+c         if(tracer) then
+c         do l=2,nlayer
+c            do iq=1,nq
+c               hco2(iq)=zq(1,l,iq)/zq(1,l-1,iq)
+c               hco2(iq)=-(zlocal(l)-zlocal(l-1))/log(hco2(iq))/1000.
+c            enddo
+c            write(225,*) l,pt(1,l),(hco2(iq),iq=1,6)
+c            write(226,*) l,(hco2(iq),iq=7,13)
+c         enddo
+c         endif
+
+      END IF       ! if(ngrid.ne.1)
+
+
+      icount=icount+1
+      RETURN
+      END
Index: trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/writediagfi.info
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/writediagfi.info	(revision 214)
+++ trunk/MESOSCALE/LMDZ.MARS/in_lmdz_mars/writediagfi.info	(revision 214)
@@ -0,0 +1,21 @@
+        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     .                  emis)
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     .                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     .                  co2ice)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil)
+        if (tracer) then
+            do iq=1,nqmx
+              write(str2(1:2),'(i2.2)') iq
+              call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq),
+     &                          'kg/kg',3,zq(1,1,iq))
+              call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
+     &                        'kg.m-2',2,qsurf(1,iq))
+            end do
+          end if
+
Index: trunk/MESOSCALE/LMD_MM_MARS/SIMU/RUN/checkrunmpi
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SIMU/RUN/checkrunmpi	(revision 214)
+++ trunk/MESOSCALE/LMD_MM_MARS/SIMU/RUN/checkrunmpi	(revision 214)
@@ -0,0 +1,12 @@
+#! /bin/bash
+
+echo Machines
+$WHERE_MPI/mpiexec -l -n 4 hostname
+echo GCM processes
+$WHERE_MPI/mpiexec -l -n 4 ps -e -o '%u %p %C %x %c' | grep gcm 
+echo Mesoscale processes
+$WHERE_MPI/mpiexec -l -n 4 ps -e -o '%u %p %C %x %c' | grep wrf
+
+
+
+
Index: trunk/MESOSCALE/LMD_MM_MARS/SIMU/checkrunmpi
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SIMU/checkrunmpi	(revision 213)
+++ 	(revision )
@@ -1,12 +1,0 @@
-#! /bin/bash
-
-echo Machines
-$WHERE_MPI/mpiexec -l -n 4 hostname
-echo GCM processes
-$WHERE_MPI/mpiexec -l -n 4 ps -e -o '%u %p %C %x %c' | grep gcm 
-echo Mesoscale processes
-$WHERE_MPI/mpiexec -l -n 4 ps -e -o '%u %p %C %x %c' | grep wrf
-
-
-
-
