Index: trunk/LMDZ.MARS/libf/phymars/aeropacity.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/aeropacity.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/aeropacity.F	(revision 1224)
@@ -6,6 +6,7 @@
       USE ioipsl_getincom, only: getin
       use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
-     &                      igcm_dust_submicron, rho_dust, rho_ice
-      use comgeomfi_h, only: lati ! grid point latitudes (rad)
+     &                      igcm_dust_submicron, rho_dust, rho_ice,
+     &                      nqdust
+      use comgeomfi_h, only: lati, sinlat ! grid point latitudes (rad)
       use yomaer_h, only: tauvis
        IMPLICIT NONE
@@ -111,9 +112,8 @@
                                !   Qabs instead of Qext
                                !   (direct comparison with TES)
+      REAL topdust0(ngrid)
 
 c   local saved variables
 c   ---------------------
-
-      REAL,SAVE,ALLOCATABLE :: topdust0(:) 
 
 c     Level under which the dust mixing ratio is held constant
@@ -126,5 +126,4 @@
 
 ! indexes of water ice and dust tracers:
-      INTEGER,ALLOCATABLE,SAVE :: nqdust(:) ! to store the indexes of dust tracers
       INTEGER,SAVE :: i_ice=0  ! water ice
       real,parameter :: odpref=610. ! DOD reference pressure (Pa)
@@ -140,8 +139,4 @@
 
       IF (firstcall) THEN
-        ! allocate local saved arrays
-        allocate(nqdust(nq))
-        allocate(topdust0(ngrid))
-        
         ! identify scatterers that are dust
         naerdust=0
@@ -167,10 +162,4 @@
           write(*,*) "aeropacity: i_ice=",i_ice
         ENDIF
-
-c       altitude of the top of the aerosol layer (km) at Ls=2.76rad:
-c       in the Viking year scenario
-        DO ig=1,ngrid
-            topdust0(ig)=60. -22.*SIN(lati(ig))**2
-        END DO
 
 c       typical profile of solsir and (1-w)^(-1):
@@ -295,4 +284,7 @@
         else if (iddist.eq.2) then          ! "Viking" scenario
           do ig=1,ngrid
+            ! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
+            ! in the Viking year scenario
+            topdust0(ig)=60. -22.*sinlat(ig)**2
             topdust(ig)=topdust0(ig)+18.*zlsconst
           end do
Index: trunk/LMDZ.MARS/libf/phymars/comdiurn_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/comdiurn_h.F90	(revision 1223)
+++ 	(revision )
@@ -1,11 +1,0 @@
-
-       module comdiurn_h
-
-       implicit none
-
-       ! These arrays are allocated in inifis
-       real, allocatable, save, dimension(:) :: sinlon, coslon, sinlat, coslat
-       logical,save :: ldiurn
-
-       end module comdiurn_h
-
Index: trunk/LMDZ.MARS/libf/phymars/comgeomfi_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/comgeomfi_h.F90	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/comgeomfi_h.F90	(revision 1224)
@@ -4,9 +4,49 @@
        implicit none
 
-       ! These arrays are allocated in inifis
+       ! These arrays are allocated in phys_state_var_init
        REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: long ! longitudes (rad)
        REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: lati ! latitudes (rad)
        REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: area ! mesh area (m2)
        REAL,SAVE :: totarea
+       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: sinlon 
+       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: coslon
+       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: sinlat
+       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: coslat
+
+       contains
+
+         subroutine ini_comgeomfi_h(ngrid)
+
+         implicit none
+         integer,intent(in) :: ngrid ! number of atmospheric columns
+
+         allocate(lati(ngrid))
+         allocate(long(ngrid))
+         allocate(area(ngrid))
+
+         end subroutine ini_comgeomfi_h
+
+         subroutine ini_sincosgeom(ngrid)
+
+         implicit none
+         integer,intent(in) :: ngrid ! number of atmospheric columns
+         integer :: ig
+
+         if (.not.allocated(lati)) then
+           PRINT*,"lati not allocated. did you run ini_comgeomfi_h?"
+           STOP
+         endif
+
+         allocate(sinlat(ngrid))
+         allocate(coslat(ngrid))
+         allocate(sinlon(ngrid))
+         allocate(coslon(ngrid))
+         DO ig=1,ngrid
+            sinlat(ig)=sin(lati(ig))
+            coslat(ig)=cos(lati(ig))
+            sinlon(ig)=sin(long(ig))
+            coslon(ig)=cos(long(ig))
+         ENDDO
+         end subroutine ini_sincosgeom
 
        end module comgeomfi_h
Index: trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 1224)
@@ -14,4 +14,14 @@
        !                 soil_settings.F)
 
+  ! variables
+  REAL,SAVE,ALLOCATABLE :: tsoil(:,:)  ! sub-surface temperatures (K)
+  real,save,allocatable :: mthermdiff(:,:)  ! mid-layer thermal diffusivity
+  real,save,allocatable :: thermdiff(:,:)   ! inter-layer thermal diffusivity
+  real,save,allocatable :: coefq(:)         ! q_{k+1/2} coefficients
+  real,save,allocatable :: coefd(:,:)       ! d_k coefficients
+  real,save,allocatable :: alph(:,:)        ! alpha_k coefficients
+  real,save,allocatable :: beta(:,:)        ! beta_k coefficients
+  real,save :: mu
+
 contains
 
@@ -24,5 +34,15 @@
     allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths
     allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia
-  
+ 
+    allocate(tsoil(ngrid,nsoilmx)) ! soil temperatures
+
+    allocate(mthermdiff(ngrid,0:nsoilmx-1))
+    allocate(thermdiff(ngrid,nsoilmx-1))
+    allocate(coefq(0:nsoilmx-1))
+    allocate(coefd(ngrid,nsoilmx-1))
+    allocate(alph(ngrid,nsoilmx-1))
+    allocate(beta(ngrid,nsoilmx-1))
+
+ 
   end subroutine ini_comsoil_h
 
Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1224)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1224)
@@ -0,0 +1,778 @@
+      SUBROUTINE conf_phys(
+     $           nq
+     $           ,day_ini,pdaysec,ptimestep
+     $           ,prad,pg,pr,pcpp
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_inifisinvar.F"
+#endif
+     $           )
+!
+!=======================================================================
+!
+!   purpose:
+!   -------
+!
+!   Initialisation for the physical parametrisations of the LMD 
+!   martian atmospheric general circulation modele.
+!
+!   author: Frederic Hourdin 15 / 10 /93
+!   -------
+!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
+!             Ehouarn Millour (oct. 2008) tracers are now identified
+!              by their names and may not be contiguously
+!              stored in the q(:,:,:,:) array
+!             E.M. (june 2009) use getin routine to load parameters
+!             adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
+!             separated inifis into conf_phys and phys_state_var_init (A. Spiga)
+!
+!
+!   arguments:
+!   ----------
+!
+!   input:
+!   ------
+!
+!    nq                    Number of tracers
+!    pdayref               Day of reference for the simulation
+!    pday                  Number of days counted from the North. Spring
+!                          equinoxe.
+!
+!=======================================================================
+!
+!-----------------------------------------------------------------------
+!   declarations:
+!   -------------
+! to use  'getin'
+      USE ioipsl_getincom, only : getin
+      use tracer_mod, only : nqmx, nuice_sed, ccn_factor
+      use surfdat_h, only: albedo_h2o_ice, inert_h2o_ice,
+     &                     frost_albedo_threshold
+      use yomaer_h,only: tauvis
+      use control_mod, only: ecritphy
+
+#ifdef MESOSCALE
+      !! see meso_inc_inifisini.F 
+      use surfdat_h, only: emissiv,albedice,iceradius,
+     &                     emisice,dtemisice,
+     &                     z0_default,z0,
+     &                     albedodat,phisfi,
+     &                     zmea,zstd,zsig,zgam,zthe
+      use slope_mod, only: theta_sl,psi_sl
+      use comsoil_h, only: volcapa !!MESOSCALE -- needed to fill volcapa
+#endif
+
+      IMPLICIT NONE
+#include "dimensions.h"
+#include "dimphys.h"
+#include "planete.h"
+#include "comcstfi.h"
+!#include "comsaison.h"
+!#include "comdiurn.h"
+!#include "comgeomfi.h"
+#include "callkeys.h"
+!#include "surfdat.h"
+!#include "dimradmars.h"
+!#include "yomaer.h"
+#include "datafile.h"
+!#include "slope.h"
+#include "microphys.h"
+!#include "tracer.h"
+! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
+#include"scatterers.h"
+
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_inifisvar.F"
+#endif
+      REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec
+
+      REAL,INTENT(IN) :: ptimestep 
+      INTEGER,INTENT(IN) :: nq,day_ini
+      INTEGER ig,ierr
+ 
+      CHARACTER ch1*12
+      CHARACTER ch80*80
+
+!      logical chem, h2o
+
+!      chem = .false.
+!      h2o = .false.
+
+
+      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
+
+      rad=prad
+      cpp=pcpp
+      g=pg
+      r=pr
+      rcp=r/cpp
+      daysec=pdaysec
+      dtphys=ptimestep
+
+      nqmx=nq
+
+!! MESOSCALE INITIALIZATIONS 
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_inifisini.F"
+#endif
+
+      ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
+      ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
+      call getin("ecritphy",ecritphy)
+      
+! --------------------------------------------------------------
+!  Reading the "callphys.def" file controlling some key options
+! --------------------------------------------------------------
+     
+      ! check that 'callphys.def' file is around
+      OPEN(99,file='callphys.def',status='old',form='formatted'
+     &     ,iostat=ierr)
+      CLOSE(99)
+      
+      IF(ierr.EQ.0) THEN
+         PRINT*
+         PRINT*
+         PRINT*,'--------------------------------------------'
+         PRINT*,' conf_phys: Parameters for the physics (callphys.def)'
+         PRINT*,'--------------------------------------------'
+
+         write(*,*) "Directory where external input files are:"
+         datafile="/u/forget/WWW/datagcm/datafile"
+         call getin("datadir",datafile) ! default path
+         write(*,*) " datafile = ",trim(datafile)
+
+         write(*,*) "Run with or without tracer transport ?"
+         tracer=.false. ! default value
+         call getin("tracer",tracer)
+         write(*,*) " tracer = ",tracer
+
+         write(*,*) "Diurnal cycle ?"
+         write(*,*) "(if diurnal=False, diurnal averaged solar heating)"
+         diurnal=.true. ! default value
+         call getin("diurnal",diurnal)
+         write(*,*) " diurnal = ",diurnal
+
+         write(*,*) "Seasonal cycle ?"
+         write(*,*) "(if season=False, Ls stays constant, to value ",
+     &   "set in 'start'"
+         season=.true. ! default value
+         call getin("season",season)
+         write(*,*) " season = ",season
+
+         write(*,*) "Write some extra output to the screen ?"
+         lwrite=.false. ! default value
+         call getin("lwrite",lwrite)
+         write(*,*) " lwrite = ",lwrite
+
+         write(*,*) "Save statistics in file stats.nc ?"
+#ifdef MESOSCALE
+         callstats=.false. ! default value
+#else
+         callstats=.true. ! default value
+#endif
+         call getin("callstats",callstats)
+         write(*,*) " callstats = ",callstats
+
+         write(*,*) "Save EOF profiles in file 'profiles' for ",
+     &              "Climate Database?"
+         calleofdump=.false. ! default value
+         call getin("calleofdump",calleofdump)
+         write(*,*) " calleofdump = ",calleofdump
+
+         write(*,*) "Dust scenario: 1=constant dust (read from startfi",
+     &   " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
+     &   "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
+     &   "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
+         iaervar=3 ! default value
+         call getin("iaervar",iaervar)
+         write(*,*) " iaervar = ",iaervar
+
+         write(*,*) "Reference (visible) dust opacity at 610 Pa ",
+     &   "(matters only if iaervar=1)"
+         ! NB: default value of tauvis is set/read in startfi.nc file
+         call getin("tauvis",tauvis)
+         write(*,*) " tauvis = ",tauvis
+
+         write(*,*) "Dust vertical distribution:"
+         write(*,*) "(=1 top set by topdustref parameter;",
+     & " =2 Viking scenario; =3 MGS scenario)"
+         iddist=3 ! default value
+         call getin("iddist",iddist)
+         write(*,*) " iddist = ",iddist
+
+         write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
+         topdustref= 90.0 ! default value
+         call getin("topdustref",topdustref)
+         write(*,*) " topdustref = ",topdustref
+
+         write(*,*) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
+         tke_heat_flux=0. ! default value
+         call getin("tke_heat_flux",tke_heat_flux)
+         write(*,*) " tke_heat_flux = ",tke_heat_flux
+         write(*,*) " 0 means the usual schemes are computing"
+
+         write(*,*) "call radiative transfer ?"
+         callrad=.true. ! default value
+         call getin("callrad",callrad)
+         write(*,*) " callrad = ",callrad
+
+         write(*,*) "call slope insolation scheme ?",
+     &              "(matters only if callrad=T)"
+#ifdef MESOSCALE
+         callslope=.true. ! default value
+#else
+         callslope=.false. ! default value (not supported yet)
+#endif
+         call getin("callslope",callslope)
+         write(*,*) " callslope = ",callslope
+
+         write(*,*) "call NLTE radiative schemes ?",
+     &              "(matters only if callrad=T)"
+         callnlte=.false. ! default value
+         call getin("callnlte",callnlte)
+         write(*,*) " callnlte = ",callnlte
+         
+         nltemodel=0    !default value
+         write(*,*) "NLTE model?"
+         write(*,*) "0 -> old model, static O"
+         write(*,*) "1 -> old model, dynamic O"
+         write(*,*) "2 -> new model"
+         write(*,*) "(matters only if callnlte=T)"
+         call getin("nltemodel",nltemodel)
+         write(*,*) " nltemodel = ",nltemodel
+
+         write(*,*) "call CO2 NIR absorption ?",
+     &              "(matters only if callrad=T)"
+         callnirco2=.false. ! default value
+         call getin("callnirco2",callnirco2)
+         write(*,*) " callnirco2 = ",callnirco2
+
+         write(*,*) "New NIR NLTE correction ?",
+     $              "0-> old model (no correction)",
+     $              "1-> new correction",
+     $              "(matters only if callnirco2=T)"
+#ifdef MESOSCALE
+         nircorr=0      !default value. this is OK below 60 km.
+#else
+         nircorr=0      !default value
+#endif
+         call getin("nircorr",nircorr)
+         write(*,*) " nircorr = ",nircorr
+
+         write(*,*) "call turbulent vertical diffusion ?"
+         calldifv=.true. ! default value
+         call getin("calldifv",calldifv)
+         write(*,*) " calldifv = ",calldifv
+
+         write(*,*) "call thermals ?"
+         calltherm=.false. ! default value
+         call getin("calltherm",calltherm)
+         write(*,*) " calltherm = ",calltherm
+
+         write(*,*) "call convective adjustment ?"
+         calladj=.true. ! default value
+         call getin("calladj",calladj)
+         write(*,*) " calladj = ",calladj
+         
+         if (calltherm .and. (.not. calladj)) then
+          print*,'Convadj has to be activated when using thermals'
+          stop
+         endif
+
+         write(*,*) "call Richardson-based surface layer ?"
+         callrichsl=.false. ! default value
+         call getin("callrichsl",callrichsl)
+         write(*,*) " callrichsl = ",callrichsl
+
+         if (calltherm .and. .not.callrichsl) then
+          print*,'WARNING WARNING WARNING'
+          print*,'if calltherm=T we strongly advise that '
+          print*,'you use the new surface layer scheme '
+          print*,'by setting callrichsl=T '
+         endif
+
+         if (calladj .and. callrichsl .and. (.not. calltherm)) then
+          print*,'You should not be calling the convective adjustment
+     & scheme with the Richardson surface-layer and without the thermals
+     &. This approach is not
+     & physically consistent and can lead to unrealistic friction
+     & values.'
+          print*,'If you want to use the Ri. surface-layer, either
+     & activate thermals OR de-activate the convective adjustment.'
+          stop
+         endif
+
+         write(*,*) "call CO2 condensation ?"
+         callcond=.true. ! default value
+         call getin("callcond",callcond)
+         write(*,*) " callcond = ",callcond
+
+         write(*,*)"call thermal conduction in the soil ?"
+         callsoil=.true. ! default value
+         call getin("callsoil",callsoil)
+         write(*,*) " callsoil = ",callsoil
+         
+
+         write(*,*)"call Lott's gravity wave/subgrid topography ",
+     &             "scheme ?"
+         calllott=.true. ! default value
+         call getin("calllott",calllott)
+         write(*,*)" calllott = ",calllott
+
+
+         write(*,*)"rad.transfer is computed every iradia",
+     &             " physical timestep"
+         iradia=1 ! default value
+         call getin("iradia",iradia)
+         write(*,*)" iradia = ",iradia
+         
+
+         write(*,*)"Output of the exchange coefficient mattrix ?",
+     &             "(for diagnostics only)"
+         callg2d=.false. ! default value
+         call getin("callg2d",callg2d)
+         write(*,*)" callg2d = ",callg2d
+
+         write(*,*)"Rayleigh scattering : (should be .false. for now)"
+         rayleigh=.false.
+         call getin("rayleigh",rayleigh)
+         write(*,*)" rayleigh = ",rayleigh
+
+
+! TRACERS:
+
+! dustbin
+         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
+         dustbin=0 ! default value
+         call getin("dustbin",dustbin)
+         write(*,*)" dustbin = ",dustbin
+! active
+         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
+         active=.false. ! default value
+         call getin("active",active)
+         write(*,*)" active = ",active
+
+! Test of incompatibility:
+! if active is used, then dustbin should be > 0
+
+         if (active.and.(dustbin.lt.1)) then
+           print*,'if active is used, then dustbin should > 0'
+           stop
+         endif
+! doubleq
+         write(*,*)"use mass and number mixing ratios to predict",
+     &             " dust size ?"
+         doubleq=.false. ! default value
+         call getin("doubleq",doubleq)
+         write(*,*)" doubleq = ",doubleq
+! submicron
+         submicron=.false. ! default value
+         call getin("submicron",submicron)
+         write(*,*)" submicron = ",submicron
+
+! Test of incompatibility:
+! if doubleq is used, then dustbin should be 2
+
+         if (doubleq.and.(dustbin.ne.2)) then
+           print*,'if doubleq is used, then dustbin should be 2'
+           stop
+         endif
+         if (doubleq.and.submicron.and.(nq.LT.3)) then
+           print*,'If doubleq is used with a submicron tracer,'
+           print*,' then the number of tracers has to be'
+           print*,' larger than 3.'
+           stop
+         endif
+
+! lifting
+         write(*,*)"dust lifted by GCM surface winds ?"
+         lifting=.false. ! default value
+         call getin("lifting",lifting)
+         write(*,*)" lifting = ",lifting
+
+! Test of incompatibility:
+! if lifting is used, then dustbin should be > 0
+
+         if (lifting.and.(dustbin.lt.1)) then
+           print*,'if lifting is used, then dustbin should > 0'
+           stop
+         endif
+
+! free evolving dust
+! freedust=true just says that there is no lifting and no dust opacity scaling.
+         write(*,*)"dust lifted by GCM surface winds ?"
+         freedust=.false. ! default value
+         call getin("freedust",freedust)
+         write(*,*)" freedust = ",freedust
+         if (freedust.and..not.doubleq) then
+           print*,'freedust should be used with doubleq !'
+           stop
+         endif
+         if (freedust.and.lifting) then
+           print*,'if freedust is used, then lifting should not be used'
+           print*,'lifting forced to false !!'
+           lifting=.false.
+         endif
+
+! callddevil
+         write(*,*)" dust lifted by dust devils ?"
+         callddevil=.false. !default value
+         call getin("callddevil",callddevil)
+         write(*,*)" callddevil = ",callddevil
+
+! Test of incompatibility:
+! if dustdevil is used, then dustbin should be > 0
+
+         if (callddevil.and.(dustbin.lt.1)) then
+           print*,'if dustdevil is used, then dustbin should > 0'
+           stop
+         endif
+! sedimentation
+         write(*,*) "Gravitationnal sedimentation ?"
+         sedimentation=.true. ! default value
+         call getin("sedimentation",sedimentation)
+         write(*,*) " sedimentation = ",sedimentation
+! activice
+         write(*,*) "Radiatively active transported atmospheric ",
+     &              "water ice ?"
+         activice=.false. ! default value
+         call getin("activice",activice)
+         write(*,*) " activice = ",activice
+! water
+         write(*,*) "Compute water cycle ?"
+         water=.false. ! default value
+         call getin("water",water)
+         write(*,*) " water = ",water
+
+! thermal inertia feedback
+         write(*,*) "Activate the thermal inertia feedback ?"
+         tifeedback=.false. ! default value
+         call getin("tifeedback",tifeedback)
+         write(*,*) " tifeedback = ",tifeedback
+
+! Test of incompatibility:
+
+         if (tifeedback.and..not.water) then
+           print*,'if tifeedback is used,'
+           print*,'water should be used too'
+           stop
+         endif
+
+         if (tifeedback.and..not.callsoil) then
+           print*,'if tifeedback is used,'
+           print*,'callsoil should be used too'
+           stop
+         endif
+
+         if (activice.and..not.water) then
+           print*,'if activice is used, water should be used too'
+           stop
+         endif
+
+         if (water.and..not.tracer) then
+           print*,'if water is used, tracer should be used too'
+           stop
+         endif
+         
+! water ice clouds effective variance distribution for sedimentaion       
+        write(*,*) "effective variance for water ice clouds ?"
+        nuice_sed=0.45 
+        call getin("nuice_sed",nuice_sed)
+        write(*,*) "water_param nueff Sedimentation:", nuice_sed
+         
+! ccn factor if no scavenging         
+        write(*,*) "water param CCN reduc. factor ?", ccn_factor
+        ccn_factor = 4.5
+        call getin("ccn_factor",ccn_factor)
+        write(*,*)" ccn_factor = ",ccn_factor
+        write(*,*)"Careful: only used when microphys=F, otherwise"
+        write(*,*)"the contact parameter is used instead;"
+
+! microphys
+         write(*,*)"Microphysical scheme for water-ice clouds?"
+         microphys=.false. ! default value
+         call getin("microphys",microphys)
+         write(*,*)" microphys = ",microphys
+
+! microphysical parameter contact       
+         write(*,*) "water contact parameter ?"
+         mteta  = 0.95
+         call getin("mteta",mteta)
+         write(*,*) "mteta = ", mteta
+
+! scavenging
+         write(*,*)"Dust scavenging by H2O/CO2 snowfall ?"
+         scavenging=.false. ! default value
+         call getin("scavenging",scavenging)
+         write(*,*)" scavenging = ",scavenging
+         
+
+! Test of incompatibility:
+! if scavenging is used, then dustbin should be > 0
+
+         if ((microphys.and..not.doubleq).or.
+     &       (microphys.and..not.water)) then
+             print*,'if microphys is used, then doubleq,'
+             print*,'and water must be used!'
+             stop
+         endif
+         if (microphys.and..not.scavenging) then
+             print*,''
+             print*,'----------------WARNING-----------------'
+             print*,'microphys is used without scavenging !!!'
+             print*,'----------------WARNING-----------------'
+             print*,''
+         endif
+
+         if ((scavenging.and..not.microphys).or.
+     &       (scavenging.and.(dustbin.lt.1))) then
+             print*,'if scavenging is used, then microphys'
+             print*,'must be used!'
+             stop
+         endif
+
+! Test of incompatibility:
+
+         write(*,*) "Permanent water caps at poles ?",
+     &               " .true. is RECOMMENDED"
+         write(*,*) "(with .true., North cap is a source of water ",
+     &   "and South pole is a cold trap)"
+         caps=.true. ! default value
+         call getin("caps",caps)
+         write(*,*) " caps = ",caps
+
+! albedo_h2o_ice
+         write(*,*) "water ice albedo ?"
+         albedo_h2o_ice=0.45
+         call getin("albedo_h2o_ice",albedo_h2o_ice)
+         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
+! inert_h2o_ice
+         write(*,*) "water ice thermal inertia ?"
+         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
+         call getin("inert_h2o_ice",inert_h2o_ice)
+         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
+! frost_albedo_threshold
+         write(*,*) "frost thickness threshold for albedo ?"
+         frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
+         call getin("frost_albedo_threshold",
+     &    frost_albedo_threshold)
+         write(*,*) " frost_albedo_threshold = ",
+     &            frost_albedo_threshold
+
+! call Titus crocus line -- DEFAULT IS NONE
+         write(*,*) "Titus crocus line ?"
+         tituscap=.false.  ! default value
+         call getin("tituscap",tituscap)
+         write(*,*) "tituscap",tituscap
+                     
+
+         write(*,*) "photochemistry: include chemical species"
+         photochem=.false. ! default value
+         call getin("photochem",photochem)
+         write(*,*) " photochem = ",photochem
+
+
+! THERMOSPHERE
+
+         write(*,*) "call thermosphere ?"
+         callthermos=.false. ! default value
+         call getin("callthermos",callthermos)
+         write(*,*) " callthermos = ",callthermos
+         
+
+         write(*,*) " water included without cycle ",
+     &              "(only if water=.false.)"
+         thermoswater=.false. ! default value
+         call getin("thermoswater",thermoswater)
+         write(*,*) " thermoswater = ",thermoswater
+
+         write(*,*) "call thermal conduction ?",
+     &    " (only if callthermos=.true.)"
+         callconduct=.false. ! default value
+         call getin("callconduct",callconduct)
+         write(*,*) " callconduct = ",callconduct
+
+         write(*,*) "call EUV heating ?",
+     &   " (only if callthermos=.true.)"
+         calleuv=.false.  ! default value
+         call getin("calleuv",calleuv)
+         write(*,*) " calleuv = ",calleuv
+
+         write(*,*) "call molecular viscosity ?",
+     &   " (only if callthermos=.true.)"
+         callmolvis=.false. ! default value
+         call getin("callmolvis",callmolvis)
+         write(*,*) " callmolvis = ",callmolvis
+
+         write(*,*) "call molecular diffusion ?",
+     &   " (only if callthermos=.true.)"
+         callmoldiff=.false. ! default value
+         call getin("callmoldiff",callmoldiff)
+         write(*,*) " callmoldiff = ",callmoldiff
+         
+
+         write(*,*) "call thermospheric photochemistry ?",
+     &   " (only if callthermos=.true.)"
+         thermochem=.false. ! default value
+         call getin("thermochem",thermochem)
+         write(*,*) " thermochem = ",thermochem
+
+         write(*,*) "Method to include solar variability"
+         write(*,*) "0-> old method (using solarcondate); ",
+     &                  "1-> variability wit E10.7"
+         solvarmod=1
+         call getin("solvarmod",solvarmod)
+         write(*,*) " solvarmod = ",solvarmod
+
+         write(*,*) "date for solar flux calculation:",
+     &   " (1985 < date < 2002)",
+     $   " (Only used if solvarmod=0)"
+         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
+         solarcondate=1993.4 ! default value
+         call getin("solarcondate",solarcondate)
+         write(*,*) " solarcondate = ",solarcondate
+         
+         write(*,*) "Solar variability as observed for MY: "
+         write(*,*) "Only if solvarmod=1"
+         solvaryear=24
+         call getin("solvaryear",solvaryear)
+         write(*,*) " solvaryear = ",solvaryear
+
+         write(*,*) "UV heating efficiency:",
+     &   "measured values between 0.19 and 0.23 (Fox et al. 1996)",
+     &   "lower values may be used to compensate low 15 um cooling"
+         euveff=0.21 !default value
+         call getin("euveff",euveff)
+         write(*,*) " euveff = ", euveff
+
+         if (.not.callthermos) then
+           if (thermoswater) then
+             print*,'if thermoswater is set, callthermos must be true'
+             stop
+           endif          
+           if (callconduct) then
+             print*,'if callconduct is set, callthermos must be true'
+             stop
+           endif        
+           if (calleuv) then
+             print*,'if calleuv is set, callthermos must be true'
+             stop
+           endif         
+           if (callmolvis) then
+             print*,'if callmolvis is set, callthermos must be true'
+             stop
+           endif        
+           if (callmoldiff) then
+             print*,'if callmoldiff is set, callthermos must be true'
+             stop
+           endif          
+           if (thermochem) then
+             print*,'if thermochem is set, callthermos must be true'
+             stop
+           endif          
+        endif
+
+! Test of incompatibility:
+! if photochem is used, then water should be used too
+
+         if (photochem.and..not.water) then
+           print*,'if photochem is used, water should be used too'
+           stop
+         endif
+
+! if callthermos is used, then thermoswater should be used too 
+! (if water not used already)
+
+         if (callthermos .and. .not.water) then
+           if (callthermos .and. .not.thermoswater) then
+             print*,'if callthermos is used, water or thermoswater 
+     &               should be used too'
+             stop
+           endif
+         endif
+
+         PRINT*,'--------------------------------------------'
+         PRINT*
+         PRINT*
+      ELSE
+         write(*,*)
+         write(*,*) 'Cannot read file callphys.def. Is it here ?'
+         stop
+      ENDIF
+
+8000  FORMAT(t5,a12,l8)
+8001  FORMAT(t5,a12,i8)
+
+      PRINT*
+      PRINT*,'conf_phys: daysec',daysec
+      PRINT*
+      PRINT*,'conf_phys: The radiative transfer is computed:'
+      PRINT*,'           each ',iradia,' physical time-step'
+      PRINT*,'        or each ',iradia*dtphys,' seconds'
+      PRINT*
+! --------------------------------------------------------------
+!  Managing the Longwave radiative transfer
+! --------------------------------------------------------------
+
+!     In most cases, the run just use the following values :
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      callemis=.true.     
+!     ilwd=10*int(daysec/dtphys) ! bug before 22/10/01       
+      ilwd=1
+      ilwn=1 !2
+      ilwb=1 !2
+      linear=.true.        
+      ncouche=3
+      alphan=0.4
+      semi=0
+
+!     BUT people working hard on the LW may want to read them in 'radia.def' 
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      OPEN(99,file='radia.def',status='old',form='formatted'
+     .     ,iostat=ierr)
+      IF(ierr.EQ.0) THEN
+         write(*,*) 'conf_phys: Reading radia.def !!!'
+         READ(99,fmt='(a)') ch1
+         READ(99,*) callemis
+         WRITE(*,8000) ch1,callemis
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) iradia
+         WRITE(*,8001) ch1,iradia
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) ilwd
+         WRITE(*,8001) ch1,ilwd
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) ilwn
+         WRITE(*,8001) ch1,ilwn
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) linear
+         WRITE(*,8000) ch1,linear
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) ncouche
+         WRITE(*,8001) ch1,ncouche
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) alphan
+         WRITE(*,*) ch1,alphan
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) ilwb
+         WRITE(*,8001) ch1,ilwb
+
+
+         READ(99,fmt='(a)') ch1
+         READ(99,'(l1)') callg2d
+         WRITE(*,8000) ch1,callg2d
+
+         READ(99,fmt='(a)') ch1
+         READ(99,*) semi
+         WRITE(*,*) ch1,semi
+      end if
+      CLOSE(99)
+
+      END
Index: trunk/LMDZ.MARS/libf/phymars/datafile.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/datafile.h	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/datafile.h	(revision 1224)
@@ -5,4 +5,4 @@
       COMMON/datadirectory/datafile
       character (len=100) :: datafile
-! NB: default value for 'datafile' is set in inifis.F
+! NB: default value for 'datafile' is set in conf_phys.F
 !-----------------------------------------------------------------------
Index: trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90	(revision 1224)
@@ -1,5 +1,5 @@
 module dimradmars_mod
 !   Declaration and settings for radiative transfer calculations
-!   Initializations and allocations are done in inifis
+!   Initializations and allocations are done in phys_state_var_init
 implicit none
   ! nflev: number of vertical layer
@@ -72,4 +72,18 @@
   integer,parameter :: nt_pademx=19
 
+!!
+!! variables
+!!
+  REAL,SAVE,ALLOCATABLE :: dtrad(:,:) ! Net atm. radiative heating rate (K.s-1)
+  REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:) ! rad. flux from sky absorbed by surface (W.m-2)
+  REAL,SAVE,ALLOCATABLE :: fluxrad(:) ! Net radiative surface flux (W.m-2)
+  REAL,SAVE,ALLOCATABLE :: albedo(:,:) ! Surface albedo in each solar band
+  REAL,SAVE,ALLOCATABLE :: tauscaling(:)   ! Convertion factor for qdust and Ndust
+! aerosol (dust or ice) extinction optical depth  at reference wavelength 
+! "longrefvis" set in dimradmars_mod , for one of the "naerkind"  kind of
+! aerosol optical properties  :
+  REAL,SAVE,ALLOCATABLE :: aerosol(:,:,:)
+  REAL,SAVE,ALLOCATABLE :: nueffdust(:,:) ! Dust effective variance
+
 contains
   
@@ -80,5 +94,5 @@
   integer,intent(in) :: ngrid ! number of atmospheric columns
   integer,intent(in) :: nlayer ! number of atmospheric layers
-  
+
    nflev=nlayer
 !  ndomainsz=ngrid
@@ -88,4 +102,12 @@
    ndlo2=ndlon
 
+   allocate(albedo(ngrid,2))
+   allocate(dtrad(ngrid,nlayer))
+   allocate(fluxrad_sky(ngrid))
+   allocate(fluxrad(ngrid))
+   allocate(tauscaling(ngrid))
+   allocate(aerosol(ngrid,nlayer,naerkind))
+   allocate(nueffdust(ngrid,nlayer))
+
   end subroutine ini_dimradmars_mod
   
Index: trunk/LMDZ.MARS/libf/phymars/inifis.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/inifis.F	(revision 1223)
+++ 	(revision )
@@ -1,863 +1,0 @@
-      SUBROUTINE inifis(
-     $           ngrid,nlayer,nq
-     $           ,day_ini,pdaysec,ptimestep
-     $           ,plat,plon,parea
-     $           ,prad,pg,pr,pcpp
-#ifdef MESOSCALE
-#include "meso_inc/meso_inc_inifisinvar.F"
-#endif
-     $           )
-!
-!=======================================================================
-!
-!   purpose:
-!   -------
-!
-!   Initialisation for the physical parametrisations of the LMD 
-!   martian atmospheric general circulation modele.
-!
-!   author: Frederic Hourdin 15 / 10 /93
-!   -------
-!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
-!             Ehouarn Millour (oct. 2008) tracers are now identified
-!              by their names and may not be contiguously
-!              stored in the q(:,:,:,:) array
-!             E.M. (june 2009) use getin routine to load parameters
-!             adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
-!
-!
-!   arguments:
-!   ----------
-!
-!   input:
-!   ------
-!
-!    ngrid                 Size of the horizontal grid.
-!                          All internal loops are performed on that grid.
-!    nlayer                Number of vertical layers.
-!    pdayref               Day of reference for the simulation
-!    pday                  Number of days counted from the North. Spring
-!                          equinoxe.
-!
-!=======================================================================
-!
-!-----------------------------------------------------------------------
-!   declarations:
-!   -------------
-! to use  'getin'
-      USE ioipsl_getincom, only : getin
-      use tracer_mod, only : nqmx, nuice_sed, ccn_factor
-      use comsoil_h, only: ini_comsoil_h
-      use comgeomfi_h, only: long, lati, area, totarea
-      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
-      use surfdat_h, only: ini_surfdat_h, albedo_h2o_ice, inert_h2o_ice,
-     &                     frost_albedo_threshold
-      use comsaison_h, only: ini_comsaison_h
-      use slope_mod, only: ini_slope_mod
-      use dimradmars_mod, only: ini_dimradmars_mod
-      use yomaer_h,only: ini_yomaer_h, tauvis
-      use yomlw_h, only: ini_yomlw_h
-      use conc_mod, only: ini_conc_mod
-      use control_mod, only: ecritphy
-
-#ifdef MESOSCALE
-      !! see meso_inc_inifisini.F 
-      use surfdat_h, only: emissiv,albedice,iceradius,
-     &                     emisice,dtemisice,
-     &                     z0_default,z0,
-     &                     albedodat,phisfi,
-     &                     zmea,zstd,zsig,zgam,zthe
-      use slope_mod, only: theta_sl,psi_sl
-      use comsoil_h, only: volcapa !!MESOSCALE -- needed to fill volcapa
-#endif
-
-
-      IMPLICIT NONE
-#include "dimensions.h"
-#include "dimphys.h"
-#include "planete.h"
-#include "comcstfi.h"
-!#include "comsaison.h"
-!#include "comdiurn.h"
-!#include "comgeomfi.h"
-#include "callkeys.h"
-!#include "surfdat.h"
-!#include "dimradmars.h"
-!#include "yomaer.h"
-#include "datafile.h"
-!#include "slope.h"
-#include "microphys.h"
-!#include "tracer.h"
-! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
-#include"scatterers.h"
-
-#ifdef MESOSCALE
-#include "meso_inc/meso_inc_inifisvar.F"
-#endif
-      REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec
-
-      REAL,INTENT(IN) :: ptimestep 
-      INTEGER,INTENT(IN) :: day_ini
-
-      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
-      REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
-      INTEGER ig,ierr
- 
-!      EXTERNAL iniorbit,orbite
-      EXTERNAL SSUM
-      REAL SSUM
- 
-      CHARACTER ch1*12
-      CHARACTER ch80*80
-
-!      logical chem, h2o
-
-!      chem = .false.
-!      h2o = .false.
-
-
-!!! 1. ALLOCATIONS
-!!! --------------
-
-      ! allocate "slope_mod" arrays
-      call ini_slope_mod(ngrid)
-
-      ! allocate "comsaison_h" arrays
-      call ini_comsaison_h(ngrid)
-
-      ! allocate "surfdat_h" arrays
-      call ini_surfdat_h(ngrid)
-
-      ! allocate "comgeomfi_h" arrays
-      allocate(lati(ngrid))
-      allocate(long(ngrid))
-      allocate(area(ngrid))
-
-      ! fill "comgeomfi_h" data
-      CALL SCOPY(ngrid,plon,1,long,1)
-      CALL SCOPY(ngrid,plat,1,lati,1)
-      CALL SCOPY(ngrid,parea,1,area,1)
-      totarea=SSUM(ngrid,area,1)
-
-      ! allocate "comdiurn_h" data
-      allocate(sinlat(ngrid))
-      allocate(coslat(ngrid))
-      allocate(sinlon(ngrid))
-      allocate(coslon(ngrid))
-
-      ! fill "comdiurn_h" data
-      DO ig=1,ngrid 
-         sinlat(ig)=sin(plat(ig))
-         coslat(ig)=cos(plat(ig))
-         sinlon(ig)=sin(plon(ig))
-         coslon(ig)=cos(plon(ig))
-      ENDDO
-
-      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
-
-      ! allocate "comsoil_h" arrays
-      call ini_comsoil_h(ngrid)
-
-      ! set some variables in "dimradmars_mod"
-      call ini_dimradmars_mod(ngrid,nlayer)
-
-      ! allocate arrays in "yomaer_h"
-      call ini_yomaer_h
-
-      ! allocate arrays in "yomlw_h"
-      call ini_yomlw_h(ngrid)
-
-      ! allocate arrays in "conc_mod"
-      call ini_conc_mod(ngrid,nlayer)
-
-
-!!! 2. SETTINGS and CONSTANTS
-!!! -------------------------
-
-      rad=prad
-      cpp=pcpp
-      g=pg
-      r=pr
-      rcp=r/cpp
-      daysec=pdaysec
-      dtphys=ptimestep
-
-      nqmx=nq
-
-!! MESOSCALE INITIALIZATIONS 
-#ifdef MESOSCALE
-#include "meso_inc/meso_inc_inifisini.F"
-#endif
-
-! --------------------------------------------------------
-!     The usual Tests
-!     --------------
-!      IF (nlayer.NE.nlayermx) THEN
-!         PRINT*,'STOP in inifis'
-!         PRINT*,'Probleme de dimensions :'
-!         PRINT*,'nlayer     = ',nlayer
-!         PRINT*,'nlayermx   = ',nlayermx
-!         STOP
-!      ENDIF
-
-      ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
-      ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
-      call getin("ecritphy",ecritphy)
-      
-! --------------------------------------------------------------
-!  Reading the "callphys.def" file controlling some key options
-! --------------------------------------------------------------
-     
-      ! check that 'callphys.def' file is around
-      OPEN(99,file='callphys.def',status='old',form='formatted'
-     &     ,iostat=ierr)
-      CLOSE(99)
-      
-      IF(ierr.EQ.0) THEN
-         PRINT*
-         PRINT*
-         PRINT*,'--------------------------------------------'
-         PRINT*,' inifis: Parameters for the physics (callphys.def)'
-         PRINT*,'--------------------------------------------'
-
-         write(*,*) "Directory where external input files are:"
-         datafile="/u/forget/WWW/datagcm/datafile"
-         call getin("datadir",datafile) ! default path
-         write(*,*) " datafile = ",trim(datafile)
-
-         write(*,*) "Run with or without tracer transport ?"
-         tracer=.false. ! default value
-         call getin("tracer",tracer)
-         write(*,*) " tracer = ",tracer
-
-         write(*,*) "Diurnal cycle ?"
-         write(*,*) "(if diurnal=False, diurnal averaged solar heating)"
-         diurnal=.true. ! default value
-         call getin("diurnal",diurnal)
-         write(*,*) " diurnal = ",diurnal
-
-         write(*,*) "Seasonal cycle ?"
-         write(*,*) "(if season=False, Ls stays constant, to value ",
-     &   "set in 'start'"
-         season=.true. ! default value
-         call getin("season",season)
-         write(*,*) " season = ",season
-
-         write(*,*) "Write some extra output to the screen ?"
-         lwrite=.false. ! default value
-         call getin("lwrite",lwrite)
-         write(*,*) " lwrite = ",lwrite
-
-         write(*,*) "Save statistics in file stats.nc ?"
-#ifdef MESOSCALE
-         callstats=.false. ! default value
-#else
-         callstats=.true. ! default value
-#endif
-         call getin("callstats",callstats)
-         write(*,*) " callstats = ",callstats
-
-         write(*,*) "Save EOF profiles in file 'profiles' for ",
-     &              "Climate Database?"
-         calleofdump=.false. ! default value
-         call getin("calleofdump",calleofdump)
-         write(*,*) " calleofdump = ",calleofdump
-
-         write(*,*) "Dust scenario: 1=constant dust (read from startfi",
-     &   " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
-     &   "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
-     &   "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
-         iaervar=3 ! default value
-         call getin("iaervar",iaervar)
-         write(*,*) " iaervar = ",iaervar
-
-         write(*,*) "Reference (visible) dust opacity at 610 Pa ",
-     &   "(matters only if iaervar=1)"
-         ! NB: default value of tauvis is set/read in startfi.nc file
-         call getin("tauvis",tauvis)
-         write(*,*) " tauvis = ",tauvis
-
-         write(*,*) "Dust vertical distribution:"
-         write(*,*) "(=1 top set by topdustref parameter;",
-     & " =2 Viking scenario; =3 MGS scenario)"
-         iddist=3 ! default value
-         call getin("iddist",iddist)
-         write(*,*) " iddist = ",iddist
-
-         write(*,*) "Dust top altitude (km). (Matters only if iddist=1)"
-         topdustref= 90.0 ! default value
-         call getin("topdustref",topdustref)
-         write(*,*) " topdustref = ",topdustref
-
-         write(*,*) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
-         tke_heat_flux=0. ! default value
-         call getin("tke_heat_flux",tke_heat_flux)
-         write(*,*) " tke_heat_flux = ",tke_heat_flux
-         write(*,*) " 0 means the usual schemes are computing"
-
-         write(*,*) "call radiative transfer ?"
-         callrad=.true. ! default value
-         call getin("callrad",callrad)
-         write(*,*) " callrad = ",callrad
-
-         write(*,*) "call slope insolation scheme ?",
-     &              "(matters only if callrad=T)"
-#ifdef MESOSCALE
-         callslope=.true. ! default value
-#else
-         callslope=.false. ! default value (not supported yet)
-#endif
-         call getin("callslope",callslope)
-         write(*,*) " callslope = ",callslope
-
-         write(*,*) "call NLTE radiative schemes ?",
-     &              "(matters only if callrad=T)"
-         callnlte=.false. ! default value
-         call getin("callnlte",callnlte)
-         write(*,*) " callnlte = ",callnlte
-         
-         nltemodel=0    !default value
-         write(*,*) "NLTE model?"
-         write(*,*) "0 -> old model, static O"
-         write(*,*) "1 -> old model, dynamic O"
-         write(*,*) "2 -> new model"
-         write(*,*) "(matters only if callnlte=T)"
-         call getin("nltemodel",nltemodel)
-         write(*,*) " nltemodel = ",nltemodel
-
-         write(*,*) "call CO2 NIR absorption ?",
-     &              "(matters only if callrad=T)"
-         callnirco2=.false. ! default value
-         call getin("callnirco2",callnirco2)
-         write(*,*) " callnirco2 = ",callnirco2
-
-         write(*,*) "New NIR NLTE correction ?",
-     $              "0-> old model (no correction)",
-     $              "1-> new correction",
-     $              "(matters only if callnirco2=T)"
-#ifdef MESOSCALE
-         nircorr=0      !default value. this is OK below 60 km.
-#else
-         nircorr=0      !default value
-#endif
-         call getin("nircorr",nircorr)
-         write(*,*) " nircorr = ",nircorr
-
-         write(*,*) "call turbulent vertical diffusion ?"
-         calldifv=.true. ! default value
-         call getin("calldifv",calldifv)
-         write(*,*) " calldifv = ",calldifv
-
-         write(*,*) "call thermals ?"
-         calltherm=.false. ! default value
-         call getin("calltherm",calltherm)
-         write(*,*) " calltherm = ",calltherm
-
-         write(*,*) "call convective adjustment ?"
-         calladj=.true. ! default value
-         call getin("calladj",calladj)
-         write(*,*) " calladj = ",calladj
-         
-         if (calltherm .and. (.not. calladj)) then
-          print*,'Convadj has to be activated when using thermals'
-          stop
-         endif
-
-         write(*,*) "call Richardson-based surface layer ?"
-         callrichsl=.false. ! default value
-         call getin("callrichsl",callrichsl)
-         write(*,*) " callrichsl = ",callrichsl
-
-         if (calltherm .and. .not.callrichsl) then
-          print*,'WARNING WARNING WARNING'
-          print*,'if calltherm=T we strongly advise that '
-          print*,'you use the new surface layer scheme '
-          print*,'by setting callrichsl=T '
-         endif
-
-         if (calladj .and. callrichsl .and. (.not. calltherm)) then
-          print*,'You should not be calling the convective adjustment
-     & scheme with the Richardson surface-layer and without the thermals
-     &. This approach is not
-     & physically consistent and can lead to unrealistic friction
-     & values.'
-          print*,'If you want to use the Ri. surface-layer, either
-     & activate thermals OR de-activate the convective adjustment.'
-          stop
-         endif
-
-         write(*,*) "call CO2 condensation ?"
-         callcond=.true. ! default value
-         call getin("callcond",callcond)
-         write(*,*) " callcond = ",callcond
-
-         write(*,*)"call thermal conduction in the soil ?"
-         callsoil=.true. ! default value
-         call getin("callsoil",callsoil)
-         write(*,*) " callsoil = ",callsoil
-         
-
-         write(*,*)"call Lott's gravity wave/subgrid topography ",
-     &             "scheme ?"
-         calllott=.true. ! default value
-         call getin("calllott",calllott)
-         write(*,*)" calllott = ",calllott
-
-
-         write(*,*)"rad.transfer is computed every iradia",
-     &             " physical timestep"
-         iradia=1 ! default value
-         call getin("iradia",iradia)
-         write(*,*)" iradia = ",iradia
-         
-
-         write(*,*)"Output of the exchange coefficient mattrix ?",
-     &             "(for diagnostics only)"
-         callg2d=.false. ! default value
-         call getin("callg2d",callg2d)
-         write(*,*)" callg2d = ",callg2d
-
-         write(*,*)"Rayleigh scattering : (should be .false. for now)"
-         rayleigh=.false.
-         call getin("rayleigh",rayleigh)
-         write(*,*)" rayleigh = ",rayleigh
-
-
-! TRACERS:
-
-! dustbin
-         write(*,*)"Transported dust ? (if >0, use 'dustbin' dust bins)"
-         dustbin=0 ! default value
-         call getin("dustbin",dustbin)
-         write(*,*)" dustbin = ",dustbin
-! active
-         write(*,*)"Radiatively active dust ? (matters if dustbin>0)"
-         active=.false. ! default value
-         call getin("active",active)
-         write(*,*)" active = ",active
-
-! Test of incompatibility:
-! if active is used, then dustbin should be > 0
-
-         if (active.and.(dustbin.lt.1)) then
-           print*,'if active is used, then dustbin should > 0'
-           stop
-         endif
-! doubleq
-         write(*,*)"use mass and number mixing ratios to predict",
-     &             " dust size ?"
-         doubleq=.false. ! default value
-         call getin("doubleq",doubleq)
-         write(*,*)" doubleq = ",doubleq
-! submicron
-         submicron=.false. ! default value
-         call getin("submicron",submicron)
-         write(*,*)" submicron = ",submicron
-
-! Test of incompatibility:
-! if doubleq is used, then dustbin should be 2
-
-         if (doubleq.and.(dustbin.ne.2)) then
-           print*,'if doubleq is used, then dustbin should be 2'
-           stop
-         endif
-         if (doubleq.and.submicron.and.(nq.LT.3)) then
-           print*,'If doubleq is used with a submicron tracer,'
-           print*,' then the number of tracers has to be'
-           print*,' larger than 3.'
-           stop
-         endif
-
-! lifting
-         write(*,*)"dust lifted by GCM surface winds ?"
-         lifting=.false. ! default value
-         call getin("lifting",lifting)
-         write(*,*)" lifting = ",lifting
-
-! Test of incompatibility:
-! if lifting is used, then dustbin should be > 0
-
-         if (lifting.and.(dustbin.lt.1)) then
-           print*,'if lifting is used, then dustbin should > 0'
-           stop
-         endif
-
-! free evolving dust
-! freedust=true just says that there is no lifting and no dust opacity scaling.
-         write(*,*)"dust lifted by GCM surface winds ?"
-         freedust=.false. ! default value
-         call getin("freedust",freedust)
-         write(*,*)" freedust = ",freedust
-         if (freedust.and..not.doubleq) then
-           print*,'freedust should be used with doubleq !'
-           stop
-         endif
-         if (freedust.and.lifting) then
-           print*,'if freedust is used, then lifting should not be used'
-           print*,'lifting forced to false !!'
-           lifting=.false.
-         endif
-
-! callddevil
-         write(*,*)" dust lifted by dust devils ?"
-         callddevil=.false. !default value
-         call getin("callddevil",callddevil)
-         write(*,*)" callddevil = ",callddevil
-
-! Test of incompatibility:
-! if dustdevil is used, then dustbin should be > 0
-
-         if (callddevil.and.(dustbin.lt.1)) then
-           print*,'if dustdevil is used, then dustbin should > 0'
-           stop
-         endif
-! sedimentation
-         write(*,*) "Gravitationnal sedimentation ?"
-         sedimentation=.true. ! default value
-         call getin("sedimentation",sedimentation)
-         write(*,*) " sedimentation = ",sedimentation
-! activice
-         write(*,*) "Radiatively active transported atmospheric ",
-     &              "water ice ?"
-         activice=.false. ! default value
-         call getin("activice",activice)
-         write(*,*) " activice = ",activice
-! water
-         write(*,*) "Compute water cycle ?"
-         water=.false. ! default value
-         call getin("water",water)
-         write(*,*) " water = ",water
-
-! thermal inertia feedback
-         write(*,*) "Activate the thermal inertia feedback ?"
-         tifeedback=.false. ! default value
-         call getin("tifeedback",tifeedback)
-         write(*,*) " tifeedback = ",tifeedback
-
-! Test of incompatibility:
-
-         if (tifeedback.and..not.water) then
-           print*,'if tifeedback is used,'
-           print*,'water should be used too'
-           stop
-         endif
-
-         if (tifeedback.and..not.callsoil) then
-           print*,'if tifeedback is used,'
-           print*,'callsoil should be used too'
-           stop
-         endif
-
-         if (activice.and..not.water) then
-           print*,'if activice is used, water should be used too'
-           stop
-         endif
-
-         if (water.and..not.tracer) then
-           print*,'if water is used, tracer should be used too'
-           stop
-         endif
-         
-! water ice clouds effective variance distribution for sedimentaion       
-        write(*,*) "effective variance for water ice clouds ?"
-        nuice_sed=0.45 
-        call getin("nuice_sed",nuice_sed)
-        write(*,*) "water_param nueff Sedimentation:", nuice_sed
-         
-! ccn factor if no scavenging         
-        write(*,*) "water param CCN reduc. factor ?", ccn_factor
-        ccn_factor = 4.5
-        call getin("ccn_factor",ccn_factor)
-        write(*,*)" ccn_factor = ",ccn_factor
-        write(*,*)"Careful: only used when microphys=F, otherwise"
-        write(*,*)"the contact parameter is used instead;"
-
-! microphys
-         write(*,*)"Microphysical scheme for water-ice clouds?"
-         microphys=.false. ! default value
-         call getin("microphys",microphys)
-         write(*,*)" microphys = ",microphys
-
-! microphysical parameter contact       
-         write(*,*) "water contact parameter ?"
-         mteta  = 0.95
-         call getin("mteta",mteta)
-         write(*,*) "mteta = ", mteta
-
-! scavenging
-         write(*,*)"Dust scavenging by H2O/CO2 snowfall ?"
-         scavenging=.false. ! default value
-         call getin("scavenging",scavenging)
-         write(*,*)" scavenging = ",scavenging
-         
-
-! Test of incompatibility:
-! if scavenging is used, then dustbin should be > 0
-
-         if ((microphys.and..not.doubleq).or.
-     &       (microphys.and..not.water)) then
-             print*,'if microphys is used, then doubleq,'
-             print*,'and water must be used!'
-             stop
-         endif
-         if (microphys.and..not.scavenging) then
-             print*,''
-             print*,'----------------WARNING-----------------'
-             print*,'microphys is used without scavenging !!!'
-             print*,'----------------WARNING-----------------'
-             print*,''
-         endif
-
-         if ((scavenging.and..not.microphys).or.
-     &       (scavenging.and.(dustbin.lt.1))) then
-             print*,'if scavenging is used, then microphys'
-             print*,'must be used!'
-             stop
-         endif
-
-! Test of incompatibility:
-
-         write(*,*) "Permanent water caps at poles ?",
-     &               " .true. is RECOMMENDED"
-         write(*,*) "(with .true., North cap is a source of water ",
-     &   "and South pole is a cold trap)"
-         caps=.true. ! default value
-         call getin("caps",caps)
-         write(*,*) " caps = ",caps
-
-! albedo_h2o_ice
-         write(*,*) "water ice albedo ?"
-         albedo_h2o_ice=0.45
-         call getin("albedo_h2o_ice",albedo_h2o_ice)
-         write(*,*) " albedo_h2o_ice = ",albedo_h2o_ice
-! inert_h2o_ice
-         write(*,*) "water ice thermal inertia ?"
-         inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
-         call getin("inert_h2o_ice",inert_h2o_ice)
-         write(*,*) " inert_h2o_ice = ",inert_h2o_ice
-! frost_albedo_threshold
-         write(*,*) "frost thickness threshold for albedo ?"
-         frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
-         call getin("frost_albedo_threshold",
-     &    frost_albedo_threshold)
-         write(*,*) " frost_albedo_threshold = ",
-     &            frost_albedo_threshold
-
-! call Titus crocus line -- DEFAULT IS NONE
-         write(*,*) "Titus crocus line ?"
-         tituscap=.false.  ! default value
-         call getin("tituscap",tituscap)
-         write(*,*) "tituscap",tituscap
-                     
-
-         write(*,*) "photochemistry: include chemical species"
-         photochem=.false. ! default value
-         call getin("photochem",photochem)
-         write(*,*) " photochem = ",photochem
-
-
-! THERMOSPHERE
-
-         write(*,*) "call thermosphere ?"
-         callthermos=.false. ! default value
-         call getin("callthermos",callthermos)
-         write(*,*) " callthermos = ",callthermos
-         
-
-         write(*,*) " water included without cycle ",
-     &              "(only if water=.false.)"
-         thermoswater=.false. ! default value
-         call getin("thermoswater",thermoswater)
-         write(*,*) " thermoswater = ",thermoswater
-
-         write(*,*) "call thermal conduction ?",
-     &    " (only if callthermos=.true.)"
-         callconduct=.false. ! default value
-         call getin("callconduct",callconduct)
-         write(*,*) " callconduct = ",callconduct
-
-         write(*,*) "call EUV heating ?",
-     &   " (only if callthermos=.true.)"
-         calleuv=.false.  ! default value
-         call getin("calleuv",calleuv)
-         write(*,*) " calleuv = ",calleuv
-
-         write(*,*) "call molecular viscosity ?",
-     &   " (only if callthermos=.true.)"
-         callmolvis=.false. ! default value
-         call getin("callmolvis",callmolvis)
-         write(*,*) " callmolvis = ",callmolvis
-
-         write(*,*) "call molecular diffusion ?",
-     &   " (only if callthermos=.true.)"
-         callmoldiff=.false. ! default value
-         call getin("callmoldiff",callmoldiff)
-         write(*,*) " callmoldiff = ",callmoldiff
-         
-
-         write(*,*) "call thermospheric photochemistry ?",
-     &   " (only if callthermos=.true.)"
-         thermochem=.false. ! default value
-         call getin("thermochem",thermochem)
-         write(*,*) " thermochem = ",thermochem
-
-         write(*,*) "Method to include solar variability"
-         write(*,*) "0-> old method (using solarcondate); ",
-     &                  "1-> variability wit E10.7"
-         solvarmod=1
-         call getin("solvarmod",solvarmod)
-         write(*,*) " solvarmod = ",solvarmod
-
-         write(*,*) "date for solar flux calculation:",
-     &   " (1985 < date < 2002)",
-     $   " (Only used if solvarmod=0)"
-         write(*,*) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
-         solarcondate=1993.4 ! default value
-         call getin("solarcondate",solarcondate)
-         write(*,*) " solarcondate = ",solarcondate
-         
-         write(*,*) "Solar variability as observed for MY: "
-         write(*,*) "Only if solvarmod=1"
-         solvaryear=24
-         call getin("solvaryear",solvaryear)
-         write(*,*) " solvaryear = ",solvaryear
-
-         write(*,*) "UV heating efficiency:",
-     &   "measured values between 0.19 and 0.23 (Fox et al. 1996)",
-     &   "lower values may be used to compensate low 15 um cooling"
-         euveff=0.21 !default value
-         call getin("euveff",euveff)
-         write(*,*) " euveff = ", euveff
-
-         if (.not.callthermos) then
-           if (thermoswater) then
-             print*,'if thermoswater is set, callthermos must be true'
-             stop
-           endif          
-           if (callconduct) then
-             print*,'if callconduct is set, callthermos must be true'
-             stop
-           endif        
-           if (calleuv) then
-             print*,'if calleuv is set, callthermos must be true'
-             stop
-           endif         
-           if (callmolvis) then
-             print*,'if callmolvis is set, callthermos must be true'
-             stop
-           endif        
-           if (callmoldiff) then
-             print*,'if callmoldiff is set, callthermos must be true'
-             stop
-           endif          
-           if (thermochem) then
-             print*,'if thermochem is set, callthermos must be true'
-             stop
-           endif          
-        endif
-
-! Test of incompatibility:
-! if photochem is used, then water should be used too
-
-         if (photochem.and..not.water) then
-           print*,'if photochem is used, water should be used too'
-           stop
-         endif
-
-! if callthermos is used, then thermoswater should be used too 
-! (if water not used already)
-
-         if (callthermos .and. .not.water) then
-           if (callthermos .and. .not.thermoswater) then
-             print*,'if callthermos is used, water or thermoswater 
-     &               should be used too'
-             stop
-           endif
-         endif
-
-         PRINT*,'--------------------------------------------'
-         PRINT*
-         PRINT*
-      ELSE
-         write(*,*)
-         write(*,*) 'Cannot read file callphys.def. Is it here ?'
-         stop
-      ENDIF
-
-8000  FORMAT(t5,a12,l8)
-8001  FORMAT(t5,a12,i8)
-
-      PRINT*
-      PRINT*,'inifis: daysec',daysec
-      PRINT*
-      PRINT*,'inifis: The radiative transfer is computed:'
-      PRINT*,'           each ',iradia,' physical time-step'
-      PRINT*,'        or each ',iradia*dtphys,' seconds'
-      PRINT*
-! --------------------------------------------------------------
-!  Managing the Longwave radiative transfer
-! --------------------------------------------------------------
-
-!     In most cases, the run just use the following values :
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      callemis=.true.     
-!     ilwd=10*int(daysec/dtphys) ! bug before 22/10/01       
-      ilwd=1
-      ilwn=1 !2
-      ilwb=1 !2
-      linear=.true.        
-      ncouche=3
-      alphan=0.4
-      semi=0
-
-!     BUT people working hard on the LW may want to read them in 'radia.def' 
-!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      OPEN(99,file='radia.def',status='old',form='formatted'
-     .     ,iostat=ierr)
-      IF(ierr.EQ.0) THEN
-         write(*,*) 'inifis: Reading radia.def !!!'
-         READ(99,fmt='(a)') ch1
-         READ(99,*) callemis
-         WRITE(*,8000) ch1,callemis
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) iradia
-         WRITE(*,8001) ch1,iradia
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) ilwd
-         WRITE(*,8001) ch1,ilwd
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) ilwn
-         WRITE(*,8001) ch1,ilwn
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) linear
-         WRITE(*,8000) ch1,linear
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) ncouche
-         WRITE(*,8001) ch1,ncouche
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) alphan
-         WRITE(*,*) ch1,alphan
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) ilwb
-         WRITE(*,8001) ch1,ilwb
-
-
-         READ(99,fmt='(a)') ch1
-         READ(99,'(l1)') callg2d
-         WRITE(*,8000) ch1,callg2d
-
-         READ(99,fmt='(a)') ch1
-         READ(99,*) semi
-         WRITE(*,*) ch1,semi
-      end if
-      CLOSE(99)
-
-      END
Index: trunk/LMDZ.MARS/libf/phymars/iniphysiq.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/iniphysiq.F90	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/iniphysiq.F90	(revision 1224)
@@ -69,6 +69,6 @@
 ! copy some fundamental parameters to physics 
 ! and do some initializations 
-call inifis(klon_omp,nlayer,nqtot,pdayref,punjours,ptimestep, &
-            rlatd,rlond,airephy,prad,pg,pr,pcpp)
+call phys_state_var_init(klon_omp,nlayer,nqtot,rlatd,rlond,airephy)
+call conf_phys(nqtot,pdayref,punjours,ptimestep,prad,pg,pr,pcpp)
 
 !$OMP END PARALLEL
Index: trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 1224)
@@ -1,3 +1,3 @@
-      SUBROUTINE initracer(ngrid,nq,qsurf,co2ice)
+      SUBROUTINE initracer(ngrid,nq,qsurf)
 
 #ifndef MESOSCALE
@@ -39,5 +39,4 @@
       integer,intent(in) :: nq ! number of tracers
       real,intent(out) :: qsurf(ngrid,nq) ! tracer on surface (e.g.  kg.m-2)
-      real,intent(out) :: co2ice(ngrid) ! co2 ice mass on surface (e.g.  kg.m-2)
 
       integer iq,ig,count
@@ -68,6 +67,6 @@
       allocate(alpha_lift(nq))
       allocate(alpha_devil(nq))
-      allocate(dryness(ngrid))
       allocate(igcm_dustbin(nq))
+      allocate(nqdust(nq))
 
 #ifndef MESOSCALE
Index: trunk/LMDZ.MARS/libf/phymars/newcondens.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/newcondens.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/newcondens.F	(revision 1224)
@@ -145,4 +145,6 @@
       LOGICAL condsub(ngrid)
 
+      real :: emisref(ngrid)
+
 c variable speciale diagnostique
       real tconda1(ngrid,nlayer)
@@ -164,5 +166,4 @@
       integer,save :: ico2 ! index of CO2 tracer
       real,save :: qco2min,qco2,mmean
-      real,allocatable,save :: emisref(:)
       real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
       real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
@@ -182,6 +183,4 @@
 c
       IF (firstcall) THEN
-         
-         allocate(emisref(ngrid))
          
          bcond=1./tcond1mb
Index: trunk/LMDZ.MARS/libf/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/newstart.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/newstart.F	(revision 1224)
@@ -383,6 +383,6 @@
 !    (for instance initracer needs to know about some flags, and/or
 !      'datafile' path may be changed by user)
-      call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
-     &                latfi,lonfi,airefi,rad,g,r,cpp)
+      call phys_state_var_init(ngridmx,llm,nqtot,latfi,lonfi,airefi)
+      call conf_phys(nqtot,day_ini,daysec,dtphys,rad,g,r,cpp)
 
 c=======================================================================
@@ -484,5 +484,5 @@
       
       ! initialize tracer names noms(:) and indexes (igcm_co2, igcm_h2o_vap, ...)
-      call initracer(ngridmx,nqtot,qsurf,co2ice)
+      call initracer(ngridmx,nqtot,qsurf)
       
       if (count.eq.nqtot) then
Index: trunk/LMDZ.MARS/libf/phymars/nirco2abs.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nirco2abs.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/nirco2abs.F	(revision 1224)
@@ -3,5 +3,5 @@
                                                    
        use tracer_mod, only: igcm_co2, igcm_o
-       use comdiurn_h, only: sinlon, coslon, sinlat, coslat
+       use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
        IMPLICIT NONE
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/phymars/phys_state_var_init.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phys_state_var_init.F90	(revision 1224)
+++ trunk/LMDZ.MARS/libf/phymars/phys_state_var_init.F90	(revision 1224)
@@ -0,0 +1,85 @@
+      SUBROUTINE phys_state_var_init(ngrid,nlayer,nq,plat,plon,parea)
+
+!=======================================================================
+!
+!   purpose:
+!   -------
+!
+!   Allocate arrays in modules
+!   Fill geometrical arrays
+!   -- was done previously in inifis
+!
+!=======================================================================
+!   
+!   authors: Ehouarn Millour and Aymeric Spiga
+!            14/04/2014
+!
+!   arguments:
+!   ----------
+!
+!   input:
+!   ------
+!
+!    ngrid                 Size of the horizontal grid.
+!    nlayer                Number of vertical layers.
+!    nq                    Number of tracers.
+!
+!=======================================================================
+
+      use slope_mod, only: ini_slope_mod
+      use comsaison_h, only: ini_comsaison_h
+      use surfdat_h, only: ini_surfdat_h
+      use comgeomfi_h, only: ini_comgeomfi_h,long,lati,area,totarea,ini_sincosgeom
+      use comsoil_h, only: ini_comsoil_h
+      use dimradmars_mod, only: ini_dimradmars_mod
+      use yomaer_h,only: ini_yomaer_h
+      use yomlw_h, only: ini_yomlw_h
+      use conc_mod, only: ini_conc_mod
+      use turb_mod, only: ini_turb_mod
+
+      IMPLICIT NONE
+      
+      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
+      REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
+      EXTERNAL SSUM
+      REAL SSUM
+
+      ! allocate "slope_mod" arrays
+      call ini_slope_mod(ngrid)
+
+      ! allocate "comsaison_h" arrays
+      call ini_comsaison_h(ngrid)
+
+      ! allocate "surfdat_h" arrays
+      call ini_surfdat_h(ngrid,nq)
+
+      ! allocate "comgeomfi_h" arrays
+      call ini_comgeomfi_h(ngrid)
+      ! fill "comgeomfi_h" data
+      call SCOPY(ngrid,plon,1,long,1)
+      call SCOPY(ngrid,plat,1,lati,1)
+      call SCOPY(ngrid,parea,1,area,1)
+      totarea=SSUM(ngrid,area,1)
+      ! create sin cos lat lon arrays
+      ! -- previously in comdiurn
+      CALL ini_sincosgeom(ngrid)
+
+      ! allocate "comsoil_h" arrays
+      call ini_comsoil_h(ngrid)
+
+      ! set some variables in "dimradmars_mod"
+      call ini_dimradmars_mod(ngrid,nlayer)
+
+      ! allocate arrays in "yomaer_h"
+      call ini_yomaer_h
+
+      ! allocate arrays in "yomlw_h"
+      call ini_yomlw_h(ngrid)
+
+      ! allocate arrays in "conc_mod" (aeronomars)
+      call ini_conc_mod(ngrid,nlayer)
+
+      ! allocate arrays in "turb_mod"
+      call ini_turb_mod(ngrid,nlayer)
+
+      END SUBROUTINE phys_state_var_init
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1224)
@@ -18,14 +18,19 @@
      &                      nuice_ref, rho_ice, rho_dust, ref_r0
       use comsoil_h, only: inertiedat, ! soil thermal inertia
-     &                     nsoilmx ! number of subsurface layers
-      use comgeomfi_h, only: long, lati, area
-      use comdiurn_h, only: sinlon, coslon, sinlat, coslat
+     &                     tsoil, nsoilmx ! number of subsurface layers
+      use comgeomfi_h, only: long, lati, area,
+     &                       sinlon, coslon, sinlat, coslat
       use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
      &                     zthe, z0, albedo_h2o_ice,
-     &                     frost_albedo_threshold
+     &                     frost_albedo_threshold,
+     &                     tsurf, co2ice, emis,
+     &                     capcal, fluxgrd, qsurf
       use comsaison_h, only: dist_sol, declin, mu0, fract
       use slope_mod, only: theta_sl, psi_sl
       use conc_mod, only: rnew, cpnew, mmean
       use control_mod, only: iphysiq, day_step, ecritstart
+      use dimradmars_mod, only: tauscaling, aerosol,
+     &                          dtrad, fluxrad_sky, fluxrad, albedo
+      use turb_mod, only: q2, wstar, hfmax_th
 #ifdef MESOSCALE
       use comsoil_h, only: mlayer,layer
@@ -200,24 +205,7 @@
 c Local saved variables:
 c ----------------------
-c     aerosol (dust or ice) extinction optical depth  at reference wavelength 
-c     "longrefvis" set in dimradmars_mod , for one of the "naerkind"  kind of
-c      aerosol optical properties  :
-      REAL,SAVE,ALLOCATABLE :: aerosol(:,:,:)
-
       INTEGER,SAVE :: day_ini  ! Initial date of the run (sol since Ls=0) 
       INTEGER,SAVE :: icount     ! counter of calls to physiq during the run.
-      REAL,SAVE,ALLOCATABLE :: tsurf(:)   ! Surface temperature (K)
-      REAL,SAVE,ALLOCATABLE :: tsoil(:,:) ! sub-surface temperatures (K)
-      REAL,SAVE,ALLOCATABLE :: co2ice(:)  ! co2 ice surface layer (kg.m-2)  
-      REAL,SAVE,ALLOCATABLE :: albedo(:,:) ! Surface albedo in each solar band
-      REAL,SAVE,ALLOCATABLE :: emis(:)    ! Thermal IR surface emissivity
-      REAL,SAVE,ALLOCATABLE :: dtrad(:,:) ! Net atm. radiative heating rate (K.s-1)
-      REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:) ! rad. flux from sky absorbed by surface (W.m-2)
-      REAL,SAVE,ALLOCATABLE :: fluxrad(:) ! Net radiative surface flux (W.m-2)
-      REAL,SAVE,ALLOCATABLE :: capcal(:) ! surface heat capacity (J m-2 K-1)
-      REAL,SAVE,ALLOCATABLE :: fluxgrd(:) ! surface conduction flux (W.m-2)
-      REAL,ALLOCATABLE,SAVE :: qsurf(:,:) ! tracer on surface (e.g. kg.m-2)
-      REAL,SAVE,ALLOCATABLE :: q2(:,:)    ! Turbulent Kinetic Energy
-      
+     
 c     Variables used by the water ice microphysical scheme:
       REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
@@ -318,5 +306,4 @@
 !      character*5 str5
       real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
-      REAL,SAVE,ALLOCATABLE :: tauscaling(:)   ! Convertion factor for qdust and Ndust
       real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
       integer igmin, lmin
@@ -372,6 +359,4 @@
       REAL zz1(ngrid)
       REAL lmax_th_out(ngrid),zmax_th(ngrid)
-      REAL,SAVE,ALLOCATABLE :: wstar(:)
-      REAL,SAVE,ALLOCATABLE :: hfmax_th(:)
       REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
       REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
@@ -399,22 +384,4 @@
       IF (firstcall) THEN
 
-        ! allocate local (saved) arrays:
-        allocate(qsurf(ngrid,nq))
-        allocate(tsoil(ngrid,nsoilmx))
-        allocate(tsurf(ngrid))
-        allocate(aerosol(ngrid,nlayer,naerkind))
-        allocate(co2ice(ngrid))
-        allocate(albedo(ngrid,2))
-        allocate(emis(ngrid))
-        allocate(dtrad(ngrid,nlayer))
-        allocate(fluxrad_sky(ngrid))
-        allocate(fluxrad(ngrid))
-        allocate(capcal(ngrid))
-        allocate(fluxgrd(ngrid))
-        allocate(q2(ngrid,nlayer+1))
-        allocate(tauscaling(ngrid))
-        allocate(wstar(ngrid))
-        allocate(hfmax_th(ngrid))
-        
 c        variables set to 0
 c        ~~~~~~~~~~~~~~~~~~
@@ -455,5 +422,5 @@
          tracerdyn=tracer
          IF (tracer) THEN
-            CALL initracer(ngrid,nq,qsurf,co2ice)
+            CALL initracer(ngrid,nq,qsurf)
          ENDIF  ! end tracer
 
Index: trunk/LMDZ.MARS/libf/phymars/soil.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 1224)
@@ -3,6 +3,10 @@
      &          timestep,tsurf,tsoil,
      &          capcal,fluxgrd)
-      use comsoil_h, only: layer, mlayer, volcapa
+
+      use comsoil_h, only: layer, mlayer, volcapa,
+     &                     mthermdiff, thermdiff, coefq,
+     &                     coefd, alph, beta, mu
       use surfdat_h, only: watercaptag, inert_h2o_ice
+
       implicit none
 
@@ -39,27 +43,8 @@
       real fluxgrd(ngrid) ! surface diffusive heat flux
 
-! local saved variables:
-      real,save,allocatable :: mthermdiff(:,:)  ! mid-layer thermal diffusivity
-      real,save,allocatable :: thermdiff(:,:)   ! inter-layer thermal diffusivity
-      real,save,allocatable :: coefq(:)		! q_{k+1/2} coefficients
-      real,save,allocatable :: coefd(:,:)	! d_k coefficients
-      real,save,allocatable :: alph(:,:)	! alpha_k coefficients
-      real,save,allocatable :: beta(:,:)	! beta_k coefficients
-      real,save :: mu
-      
 ! local variables:
       integer ig,ik
 
 ! 0. Initialisations and preprocessing step
-      if (firstcall) then
-        ! allocate local saved arrays:
-        allocate(mthermdiff(ngrid,0:nsoil-1))
-        allocate(thermdiff(ngrid,nsoil-1))
-        allocate(coefq(0:nsoil-1))
-        allocate(coefd(ngrid,nsoil-1))
-        allocate(alph(ngrid,nsoil-1))
-        allocate(beta(ngrid,nsoil-1))
-      endif
-      
       if (firstcall.or.tifeedback) then
       ! note: firstcall is set to .true. or .false. by the caller
Index: trunk/LMDZ.MARS/libf/phymars/soil_tifeedback.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/soil_tifeedback.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/soil_tifeedback.F	(revision 1224)
@@ -14,5 +14,5 @@
 c   - One layer of surface water ice (the thickness is given
 c     by the variable icecover (in kg of ice per m2) and the thermal
-c     inertia is prescribed by inert_h2o_ice (see surfdat_h and inifis)); 
+c     inertia is prescribed by inert_h2o_ice (see surfdat_h)); 
 c   - A transitional layer of mixed thermal inertia;
 c   - A last layer of regolith below the ice cover whose thermal inertia
Index: trunk/LMDZ.MARS/libf/phymars/start2archive.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/start2archive.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/start2archive.F	(revision 1224)
@@ -134,7 +134,6 @@
 ! allocate arrays:
       allocate(q(ip1jmp1,llm,nqtot))
-      allocate(qsurf(ngridmx,nqtot))
       allocate(qsurfS(ip1jmp1,nqtot))
-      call ini_surfdat_h(ngridmx)
+      call ini_surfdat_h(ngridmx,nqtot)
       call ini_comsoil_h(ngridmx)
       
Index: trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 1224)
@@ -1,5 +1,5 @@
 module surfdat_h
 
-  ! arrays are allocated in inifis
+  ! arrays are allocated in conf_phys
   real,save,allocatable :: albedodat(:) ! albedo of bare ground
   real,save,allocatable :: phisfi(:) ! geopotential at ground level
@@ -9,4 +9,5 @@
   logical,save :: TESicealbedo ! use TES ice cap albedoes (if set to .true.)
   logical,save,allocatable :: watercaptag(:) ! flag for water ice surface
+  real, save, allocatable :: dryness(:)
       
   logical,save :: temptag !temp tag for water caps
@@ -22,14 +23,24 @@
   real,save :: z0_default ! default (constant over planet) surface roughness (m)
 
+  !! variables
+  REAL,SAVE,ALLOCATABLE :: tsurf(:)   ! Surface temperature (K)
+  REAL,SAVE,ALLOCATABLE :: co2ice(:)  ! co2 ice surface layer (kg.m-2)  
+  REAL,SAVE,ALLOCATABLE :: emis(:)    ! Thermal IR surface emissivity
+  REAL,SAVE,ALLOCATABLE :: capcal(:) ! surface heat capacity (J m-2 K-1)
+  REAL,SAVE,ALLOCATABLE :: fluxgrd(:) ! surface conduction flux (W.m-2)
+  REAL,ALLOCATABLE,SAVE :: qsurf(:,:) ! tracer on surface (e.g. kg.m-2)
+
 contains
 
-  subroutine ini_surfdat_h(ngrid)
+  subroutine ini_surfdat_h(ngrid,nq)
   
   implicit none
   integer,intent(in) :: ngrid ! number of atmospheric columns
-  
+  integer,intent(in) :: nq ! number of tracers  
+
     allocate(albedodat(ngrid))
     allocate(phisfi(ngrid))
     allocate(watercaptag(ngrid))
+    allocate(dryness(ngrid))
     allocate(zmea(ngrid))
     allocate(zstd(ngrid))
@@ -38,5 +49,11 @@
     allocate(zthe(ngrid))
     allocate(z0(ngrid))
-  
+    allocate(qsurf(ngrid,nq))
+    allocate(tsurf(ngrid))
+    allocate(co2ice(ngrid))
+    allocate(emis(ngrid))
+    allocate(capcal(ngrid))
+    allocate(fluxgrd(ngrid))
+ 
   end subroutine ini_surfdat_h
 
Index: trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 1224)
@@ -3,9 +3,9 @@
       USE ioipsl_getincom
       use netcdf
-      use tracer_mod, only: nqmx, noms, dryness
+      use tracer_mod, only: nqmx, noms
       use comgeomfi_h, only: long, lati
       use surfdat_h, only: watercaptag, frost_albedo_threshold,
      &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
-     &                     albedice
+     &                     albedice, dryness
 #ifndef MESOSCALE
       use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
Index: trunk/LMDZ.MARS/libf/phymars/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/testphys1d.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/testphys1d.F	(revision 1224)
@@ -5,6 +5,5 @@
       use infotrac, only: nqtot, tname
       use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx
-      use comgeomfi_h, only: lati, long, area
-      use comdiurn_h, only: sinlat
+      use comgeomfi_h, only: lati, long, area, sinlat
       use surfdat_h, only: albedodat, z0_default, emissiv, emisice,
      &                     albedice, iceradius, dtemisice, z0,
@@ -473,10 +472,10 @@
       longitude=longitude*pi/180.E+0
 
-!  "inifis" does some initializations (some of which have already been
+!  some initializations (some of which have already been
 !  done above!) and loads parameters set in callphys.def
 !  and allocates some arrays 
 !Mars possible matter with dtphys in input and include!!!
-      CALL inifis(1,llm,nq,day0,daysec,dtphys,
-     &            latitude,longitude,1.0,rad,g,r,cpp)
+      call phys_state_var_init(1,llm,nq,latitude,longitude,1.0)
+      call conf_phys(nq,day0,daysec,dtphys,rad,g,r,cpp)
 
 
Index: trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 1224)
@@ -4,5 +4,5 @@
  
       ! number of tracers:
-      integer,save :: nqmx ! initialized in inifis
+      integer,save :: nqmx ! initialized in conf_phys
       
       character*20,allocatable,save ::  noms(:)  ! name of the tracer
@@ -23,6 +23,6 @@
       real,save :: ccn_factor  ! ratio of nuclei for water ice particles
 
-      real,allocatable,save :: dryness(:)!"Dryness coefficient" for grnd water ice sublimation
-      
+      INTEGER,ALLOCATABLE,SAVE :: nqdust(:) ! to store the indexes of dust tracers (cf aeropacity)
+
 ! tracer indexes: these are initialized in initracer and should be 0 if the
 !                 corresponding tracer does not exist
Index: trunk/LMDZ.MARS/libf/phymars/turb_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/turb_mod.F90	(revision 1224)
+++ trunk/LMDZ.MARS/libf/phymars/turb_mod.F90	(revision 1224)
@@ -0,0 +1,24 @@
+module turb_mod
+
+  !! variables
+  REAL,SAVE,ALLOCATABLE :: q2(:,:)    ! Turbulent Kinetic Energy
+  REAL,allocatable,SAVE :: l0(:)
+  REAL,SAVE,ALLOCATABLE :: wstar(:)
+  REAL,SAVE,ALLOCATABLE :: hfmax_th(:)
+
+contains
+
+  subroutine ini_turb_mod(ngrid,nlayer)
+  
+  implicit none
+  integer,intent(in) :: ngrid ! number of atmospheric columns
+  integer,intent(in) :: nlayer ! number of atmospheric layers
+
+    allocate(q2(ngrid,nlayer+1))
+    allocate(l0(ngrid))
+    allocate(wstar(ngrid))
+    allocate(hfmax_th(ngrid))
+ 
+  end subroutine ini_turb_mod
+
+end module turb_mod
Index: trunk/LMDZ.MARS/libf/phymars/updatereffrad.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/updatereffrad.F	(revision 1224)
@@ -8,4 +8,5 @@
      &                       igcm_ccn_number, nuice_ref, varian,
      &                       ref_r0, igcm_dust_submicron
+       USE dimradmars_mod, only: nueffdust
        IMPLICIT NONE
 c=======================================================================
@@ -84,19 +85,8 @@
       REAL tauscaling(ngrid)         ! Convertion factor for qccn and Nccn
 
-      LOGICAL,SAVE :: firstcall=.true. , firstcall_alloc=.true.
+      LOGICAL,SAVE :: firstcall=.true.
 
       REAL CBRT
       EXTERNAL CBRT
-
-      REAL,SAVE,ALLOCATABLE :: nueffdust(:,:) ! Dust effective variance
-
-c     Local saved variables:
-c     ---------------------
-
-        if (firstcall_alloc) then
-          !allocate local saved arrays:
-          allocate(nueffdust(ngrid,nlayer))
-          firstcall_alloc=.false.
-        endif
 
 c==================================================================
Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 1224)
@@ -13,6 +13,6 @@
       use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
      &                      igcm_dust_submicron, igcm_h2o_vap,
-     &                      igcm_h2o_ice, dryness, alpha_lift
-      use surfdat_h, only: watercaptag, frost_albedo_threshold
+     &                      igcm_h2o_ice, alpha_lift
+      use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
       IMPLICIT NONE
 
Index: trunk/LMDZ.MARS/libf/phymars/yamada4.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/yamada4.F	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/yamada4.F	(revision 1224)
@@ -12,4 +12,5 @@
      s   ,iflag_pbl)
       use tracer_mod, only: noms
+      use turb_mod, only: l0
       IMPLICIT NONE
 !.......................................................................
@@ -91,5 +92,4 @@
       REAL m2cstat,mcstat,kmcstat
       REAL l(ngrid,nlay+1)
-      REAL,allocatable,SAVE :: l0(:)
       REAL sq(ngrid),sqz(ngrid),zz(ngrid,nlay+1)
       INTEGER iter
@@ -151,5 +151,4 @@
            enddo
         endif
-      allocate(l0(ngrid))
       ndt=ceiling(3840./(3699.*24./dt))
       firstcall=.false.
Index: trunk/LMDZ.MARS/libf/phymars/yomlw_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/yomlw_h.F90	(revision 1223)
+++ trunk/LMDZ.MARS/libf/phymars/yomlw_h.F90	(revision 1224)
@@ -24,5 +24,5 @@
 contains
 
-! Allocate array (subroutine ini_yomlw_h is called by inifis)
+! Allocate array (subroutine ini_yomlw_h is called by phys_state_var_init)
   subroutine ini_yomlw_h(ngrid)
     
