Index: trunk/DOC/chantiers/commit_importants.log
===================================================================
--- trunk/DOC/chantiers/commit_importants.log	(revision 1523)
+++ trunk/DOC/chantiers/commit_importants.log	(revision 1524)
@@ -1645,2 +1645,10 @@
 - turn "iniphysiq.F90" into module "iniphysiq_mod.F90"
   (and of course adapt gcm.F[90] and 1D models accordingly)
+
+**********************
+**** commit_v1523 ****
+**********************
+Ehouarn: More updates to enforce dynamics/physics separation:
+  get rid of references to "temps_mod" from physics packages;
+  make a "time_phylmdz_mod.F90" module to store that
+  information and fill it via "iniphysiq".
Index: trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq_mod.F90	(revision 1523)
+++ trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phytitan/iniphysiq_mod.F90	(revision 1524)
@@ -22,4 +22,6 @@
                         rlond, & ! longitudes
                         rlatd ! latitudes
+  USE temps_mod, ONLY: annee_ref, day_ref, day_ini, day_end
+  USE time_phylmdz_mod, ONLY: init_time
   USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
                                  east, west, north, south, &
@@ -172,4 +174,7 @@
   call suphec
 
+  ! Initialize some "temporal and calendar" related variables
+  CALL init_time(annee_ref,day_ref,day_ini,day_end,ptimestep)
+
 !$OMP END PARALLEL
 
Index: trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq_mod.F90	(revision 1523)
+++ trunk/LMDZ.COMMON/libf/dynlonlat_phylonlat/phyvenus/iniphysiq_mod.F90	(revision 1524)
@@ -22,4 +22,6 @@
                         rlond, & ! longitudes
                         rlatd ! latitudes
+  USE temps_mod, ONLY: annee_ref, day_ref, day_ini, day_end
+  USE time_phylmdz_mod, ONLY: init_time
   USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
                                  east, west, north, south, &
@@ -172,4 +174,7 @@
   call suphec
 
+  ! Initialize some "temporal and calendar" related variables
+  CALL init_time(annee_ref,day_ref,day_ini,day_end,ptimestep)
+
 !$OMP END PARALLEL
 
Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 1523)
+++ trunk/LMDZ.GENERIC/README	(revision 1524)
@@ -1141,2 +1141,9 @@
 - Turned iniphysiq (in dynlonlat_phylonlat/phymars)into module
   "iniphysiq_mod.F90".
+
+== 29/03/2016 == EM
+- Added "time_phylmdz_mod.F90" module to store information on time and
+  calendar in the physics, should be used instead of accessing "temps_mod"
+  from the dynamics. And moved daysec and dtphys from module "comcstfi_mod"
+  to module "time_phylmdz_mod".
+- made "inifis" a module.
Index: trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90	(revision 1524)
@@ -24,4 +24,5 @@
 use planete_mod, only: ini_planete_mod
 USE comvert_mod, ONLY: ap,bp,preff
+use inifis_mod, only: inifis
 use regular_lonlat_mod, only: init_regular_lonlat, &
                               east, west, north, south, &
Index: trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/newstart.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/newstart.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/newstart.F	(revision 1524)
@@ -22,7 +22,5 @@
       USE comgeomfi_h, ONLY: lati, long, area
       use datafile_mod, only: datadir, surfdir
-! to use  'getin'
-!      USE ioipsl_getincom, only: getin
-      USE ioipsl_getincom_p, only: getin_p
+      use ioipsl_getin_p_mod, only: getin_p
       use control_mod, only: day_step, iphysiq, anneeref, planet_type
       use phyredem, only: physdem0, physdem1
@@ -37,4 +35,5 @@
       USE temps_mod, ONLY: day_ini
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      use inifis_mod, only: inifis
       implicit none
 
Index: trunk/LMDZ.GENERIC/libf/phystd/comcstfi_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/comcstfi_mod.F90	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/comcstfi_mod.F90	(revision 1524)
@@ -8,10 +8,8 @@
       REAL,SAVE :: cpp ! Cp of the atmosphere
       REAL,SAVE :: rcp ! r/cpp
-      REAL,SAVE :: dtphys ! physics time step (s)
-      REAL,SAVE :: daysec ! length of day (s)
       REAL,SAVE :: mugaz ! molar mass of the atmosphere (g/mol)
       REAL,SAVE :: omeg ! planet rotation rate (rad/s)
       REAL,SAVE :: avocado ! something like 6.022e23
-!$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg,avocado)
+!$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado)
 
 END MODULE comcstfi_mod
Index: trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F	(revision 1524)
@@ -19,10 +19,12 @@
      &         obliquit,nres,z0,lmixmin,emin_turb,coefvis,coefir,
      &         timeperi,e_elips,p_elips
-      use comcstfi_mod, only: pi, cpp, daysec, dtphys, rad, g, r,
+      use comcstfi_mod, only: pi, cpp, rad, g, r,
      &                        mugaz, rcp, omeg
+      use time_phylmdz_mod, only: daysec, dtphys
       use callkeys_mod, only: tracer,check_cpp_match,rings_shadow,
-     &		specOLR,water,pceil,ok_slab_ocean
+     &	                      specOLR,water,pceil,ok_slab_ocean
       USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff
       USE logic_mod, ONLY: hybrid,autozlevs
+      use inifis_mod, only: inifis
       implicit none
 
@@ -130,5 +132,5 @@
       character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
 
-      real :: latitude, longitude
+      real :: latitude(1), longitude(1)
 
 c=======================================================================
@@ -339,4 +341,160 @@
       longitude=longitude*pi/180.E+0
 
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!! PLANETARY CONSTANTS !!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      g = -99999.
+      PRINT *,'GRAVITY in m s-2 ?'
+      call getin("g",g)
+      IF (g.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED g IN RCM1D.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> g = ",g
+      ENDIF
+
+      rad = -99999.
+      PRINT *,'PLANETARY RADIUS in m ?'
+      call getin("rad",rad)
+      ! Planetary  radius is needed to compute shadow of the rings
+      IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN
+          PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> rad = ",rad
+      ENDIF
+
+      daysec = -99999.
+      PRINT *,'LENGTH OF A DAY in s ?'
+      call getin("daysec",daysec)
+      IF (daysec.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> daysec = ",daysec
+      ENDIF
+      omeg=4.*asin(1.)/(daysec)
+      PRINT *,"OK. FROM THIS I WORKED OUT:"
+      PRINT *,"--> omeg = ",omeg
+
+      year_day = -99999.
+      PRINT *,'LENGTH OF A YEAR in days ?'
+      call getin("year_day",year_day)
+      IF (year_day.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
+          STOP
+      ELSE 
+          PRINT *,"--> year_day = ",year_day
+      ENDIF
+
+      periastr = -99999.
+      PRINT *,'MIN DIST STAR-PLANET in AU ?'
+      call getin("periastr",periastr)
+      IF (periastr.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> periastr = ",periastr
+      ENDIF
+
+      apoastr = -99999.
+      PRINT *,'MAX DIST STAR-PLANET in AU ?'
+      call getin("apoastr",apoastr)
+      IF (apoastr.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> apoastr = ",apoastr
+      ENDIF
+
+      peri_day = -99999.
+      PRINT *,'DATE OF PERIASTRON in days ?'
+      call getin("peri_day",peri_day)
+      IF (peri_day.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."
+          STOP
+      ELSE IF (peri_day.gt.year_day) THEN
+          PRINT *,"STOP. peri_day.gt.year_day"
+          STOP
+      ELSE  
+          PRINT *,"--> peri_day = ", peri_day
+      ENDIF 
+
+      obliquit = -99999. 
+      PRINT *,'OBLIQUITY in deg ?'
+      call getin("obliquit",obliquit)
+      IF (obliquit.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> obliquit = ",obliquit
+      ENDIF 
+
+      psurf = -99999.
+      PRINT *,'SURFACE PRESSURE in Pa ?'
+      call getin("psurf",psurf)
+      IF (psurf.eq.-99999.) THEN
+          PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
+          STOP
+      ELSE
+          PRINT *,"--> psurf = ",psurf
+      ENDIF
+      !! we need reference pressures
+      pa=psurf/30.
+      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!! END PLANETARY CONSTANTS !!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+c  Date et heure locale du debut du run
+c  ------------------------------------
+c    Date (en sols depuis le solstice de printemps) du debut du run
+      day0 = 0 ! default value for day0
+      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
+      call getin("day0",day0)
+      day=float(day0)
+      write(*,*) " day0 = ",day0
+c  Heure de demarrage
+      time=0 ! default value for time
+      write(*,*)'Initial local time (in hours, between 0 and 24)?'
+      call getin("time",time)
+      write(*,*)" time = ",time
+      time=time/24.E+0 ! convert time (hours) to fraction of sol
+
+
+c  Discretisation (Definition de la grille et des pas de temps)
+c  --------------
+c
+      nlayer=llm
+      nlevel=nlayer+1
+      nsoil=nsoilmx
+
+      day_step=48 ! default value for day_step
+      PRINT *,'Number of time steps per sol ?'
+      call getin("day_step",day_step)
+      write(*,*) " day_step = ",day_step
+
+       
+      ecritphy=day_step ! default value for ecritphy
+      PRINT *,'Nunber of steps between writediagfi ?'
+      call getin("ecritphy",ecritphy)
+      write(*,*) " ecritphy = ",ecritphy
+
+      ndt=10 ! default value for ndt
+      PRINT *,'Number of sols to run ?'
+      call getin("ndt",ndt)
+      write(*,*) " ndt = ",ndt
+
+      ndt=ndt*day_step     
+      dtphys=daysec/day_step  
+      write(*,*)"-------------------------------------"
+      write(*,*)"-------------------------------------"
+      write(*,*)"--> Physical timestep is ",dtphys
+      write(*,*)"-------------------------------------"
+      write(*,*)"-------------------------------------"
+
 !!! CALL INIFIS
 !!! - read callphys.def
@@ -355,160 +513,4 @@
       r = 8.314511E+0 * 1000.E+0 / mugaz
       rcp = r / cpp
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!! PLANETARY CONSTANTS !!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      g = -99999.
-      PRINT *,'GRAVITY in m s-2 ?'
-      call getin("g",g)
-      IF (g.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED g IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> g = ",g
-      ENDIF
-
-      rad = -99999.
-      PRINT *,'PLANETARY RADIUS in m ?'
-      call getin("rad",rad)
-      ! Planetary  radius is needed to compute shadow of the rings
-      IF (rad.eq.-99999. .and. rings_shadow .eqv. .true.) THEN
-          PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> rad = ",rad
-      ENDIF
-
-      daysec = -99999.
-      PRINT *,'LENGTH OF A DAY in s ?'
-      call getin("daysec",daysec)
-      IF (daysec.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> daysec = ",daysec
-      ENDIF
-      omeg=4.*asin(1.)/(daysec)
-      PRINT *,"OK. FROM THIS I WORKED OUT:"
-      PRINT *,"--> omeg = ",omeg
-
-      year_day = -99999.
-      PRINT *,'LENGTH OF A YEAR in days ?'
-      call getin("year_day",year_day)
-      IF (year_day.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
-          STOP
-      ELSE 
-          PRINT *,"--> year_day = ",year_day
-      ENDIF
-
-      periastr = -99999.
-      PRINT *,'MIN DIST STAR-PLANET in AU ?'
-      call getin("periastr",periastr)
-      IF (periastr.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> periastr = ",periastr
-      ENDIF
-
-      apoastr = -99999.
-      PRINT *,'MAX DIST STAR-PLANET in AU ?'
-      call getin("apoastr",apoastr)
-      IF (apoastr.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> apoastr = ",apoastr
-      ENDIF
-
-      peri_day = -99999.
-      PRINT *,'DATE OF PERIASTRON in days ?'
-      call getin("peri_day",peri_day)
-      IF (peri_day.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."
-          STOP
-      ELSE IF (peri_day.gt.year_day) THEN
-          PRINT *,"STOP. peri_day.gt.year_day"
-          STOP
-      ELSE  
-          PRINT *,"--> peri_day = ", peri_day
-      ENDIF 
-
-      obliquit = -99999. 
-      PRINT *,'OBLIQUITY in deg ?'
-      call getin("obliquit",obliquit)
-      IF (obliquit.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> obliquit = ",obliquit
-      ENDIF 
-
-      psurf = -99999.
-      PRINT *,'SURFACE PRESSURE in Pa ?'
-      call getin("psurf",psurf)
-      IF (psurf.eq.-99999.) THEN
-          PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
-          STOP
-      ELSE
-          PRINT *,"--> psurf = ",psurf
-      ENDIF
-      !! we need reference pressures
-      pa=psurf/30.
-      preff=psurf ! these values are not needed in 1D  (are you sure JL12)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!! END PLANETARY CONSTANTS !!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-c  Date et heure locale du debut du run
-c  ------------------------------------
-c    Date (en sols depuis le solstice de printemps) du debut du run
-      day0 = 0 ! default value for day0
-      write(*,*) 'Initial date (in martian sols ; =0 at Ls=0)?'
-      call getin("day0",day0)
-      day=float(day0)
-      write(*,*) " day0 = ",day0
-c  Heure de demarrage
-      time=0 ! default value for time
-      write(*,*)'Initial local time (in hours, between 0 and 24)?'
-      call getin("time",time)
-      write(*,*)" time = ",time
-      time=time/24.E+0 ! convert time (hours) to fraction of sol
-
-
-c  Discretisation (Definition de la grille et des pas de temps)
-c  --------------
-c
-      nlayer=llm
-      nlevel=nlayer+1
-      nsoil=nsoilmx
-
-      day_step=48 ! default value for day_step
-      PRINT *,'Number of time steps per sol ?'
-      call getin("day_step",day_step)
-      write(*,*) " day_step = ",day_step
-
-       
-      ecritphy=day_step ! default value for ecritphy
-      PRINT *,'Nunber of steps between writediagfi ?'
-      call getin("ecritphy",ecritphy)
-      write(*,*) " ecritphy = ",ecritphy
-
-      ndt=10 ! default value for ndt
-      PRINT *,'Number of sols to run ?'
-      call getin("ndt",ndt)
-      write(*,*) " ndt = ",ndt
-
-      ndt=ndt*day_step     
-      dtphys=daysec/day_step  
-      write(*,*)"-------------------------------------"
-      write(*,*)"-------------------------------------"
-      write(*,*)"--> Physical timestep is ",dtphys
-      write(*,*)"-------------------------------------"
-      write(*,*)"-------------------------------------"
 
 c output spectrum?
Index: trunk/LMDZ.GENERIC/libf/phystd/dyn1d/temps_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/dyn1d/temps_mod.F90	(revision 1523)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/temps_mod.F90
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 1523)
+++ 	(revision )
@@ -1,720 +1,0 @@
-      SUBROUTINE inifis(ngrid,nlayer,nq,
-     $           day_ini,pdaysec,ptimestep,
-     $           plat,plon,parea,
-     $           prad,pg,pr,pcpp)
-
-      use radinc_h, only : naerkind
-      use datafile_mod, only: datadir
-      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
-      use comgeomfi_h, only: long, lati, area, totarea, totarea_planet
-      use comsoil_h, only: ini_comsoil_h
-      use control_mod, only: ecritphy
-      use comcstfi_mod, only: rad, daysec, dtphys, cpp, g, r, rcp,
-     &                        mugaz, pi, avocado
-      use planete_mod, only: nres
-      use planetwide_mod, only: planetwide_sumval
-      use callkeys_mod
-
-!=======================================================================
-!
-!   purpose:
-!   -------
-!
-!   Initialisation for the physical parametrisations of the LMD 
-!   Generic Model.
-!
-!   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
-!
-!
-!   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:
-!   -------------
-      use datafile_mod, only: datadir
-      use ioipsl_getin_p_mod, only: getin_p
-      IMPLICIT NONE
-
-
-
-      REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
- 
-      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
-      REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
-      integer day_ini
-      INTEGER ig,ierr
- 
-      EXTERNAL iniorbit,orbite
-      EXTERNAL SSUM
-      REAL SSUM
- 
-      CHARACTER ch1*12
-      CHARACTER ch80*80
-
-      logical chem, h2o
-      logical :: parameter, doubleq=.false.
-
-      real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
-
-      rad=prad
-      daysec=pdaysec
-      dtphys=ptimestep
-      cpp=pcpp
-      g=pg
-      r=pr
-      rcp=r/cpp
-
-      avocado = 6.02214179e23   ! added by RW
-
-
-      ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
-      ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
-      call getin_p("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) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
-      
-!!!      IF(ierr.EQ.0) THEN
-      IF(iscallphys) THEN
-         PRINT*
-         PRINT*
-         PRINT*,'--------------------------------------------'
-         PRINT*,' inifis: Parametres pour la physique (callphys.def)'
-         PRINT*,'--------------------------------------------'
-
-         write(*,*) "Directory where external input files are:"
-         ! default 'datadir' is set in "datadir_mod"
-         call getin_p("datadir",datadir) ! default path
-         write(*,*) " datadir = ",trim(datadir)
-
-         write(*,*) "Run with or without tracer transport ?"
-         tracer=.false. ! default value
-         call getin_p("tracer",tracer)
-         write(*,*) " tracer = ",tracer
-
-         write(*,*) "Run with or without atm mass update ",
-     &      " due to tracer evaporation/condensation?"
-         mass_redistrib=.false. ! default value
-         call getin_p("mass_redistrib",mass_redistrib)
-         write(*,*) " mass_redistrib = ",mass_redistrib
-
-         write(*,*) "Diurnal cycle ?"
-         write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
-         diurnal=.true. ! default value
-         call getin_p("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_p("season",season)
-         write(*,*) " season = ",season
-
-         write(*,*) "Tidally resonant rotation ?"
-         tlocked=.false. ! default value
-         call getin_p("tlocked",tlocked)
-         write(*,*) "tlocked = ",tlocked
-
-         write(*,*) "Saturn ring shadowing ?"
-         rings_shadow = .false.
-         call getin_p("rings_shadow", rings_shadow)
-         write(*,*) "rings_shadow = ", rings_shadow
-         
-         write(*,*) "Compute latitude-dependent gravity field?"
-         oblate = .false.
-         call getin_p("oblate", oblate)
-         write(*,*) "oblate = ", oblate
-
-         write(*,*) "Flattening of the planet (a-b)/a "
-         flatten = 0.0
-         call getin_p("flatten", flatten)
-         write(*,*) "flatten = ", flatten
-         
-
-         write(*,*) "Needed if oblate=.true.: J2"
-         J2 = 0.0
-         call getin_p("J2", J2)
-         write(*,*) "J2 = ", J2
-         
-         write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
-         MassPlanet = 0.0
-         call getin_p("MassPlanet", MassPlanet)
-         write(*,*) "MassPlanet = ", MassPlanet         
-
-         write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
-         Rmean = 0.0
-         call getin_p("Rmean", Rmean)
-         write(*,*) "Rmean = ", Rmean
-         
-! Test of incompatibility:
-! if tlocked, then diurnal should be false
-         if (tlocked.and.diurnal) then
-           print*,'If diurnal=true, we should turn off tlocked.'
-           stop
-         endif
-
-         write(*,*) "Tidal resonance ratio ?"
-         nres=0          ! default value
-         call getin_p("nres",nres)
-         write(*,*) "nres = ",nres
-
-         write(*,*) "Write some extra output to the screen ?"
-         lwrite=.false. ! default value
-         call getin_p("lwrite",lwrite)
-         write(*,*) " lwrite = ",lwrite
-
-         write(*,*) "Save statistics in file stats.nc ?"
-         callstats=.true. ! default value
-         call getin_p("callstats",callstats)
-         write(*,*) " callstats = ",callstats
-
-         write(*,*) "Test energy conservation of model physics ?"
-         enertest=.false. ! default value
-         call getin_p("enertest",enertest)
-         write(*,*) " enertest = ",enertest
-
-         write(*,*) "Check to see if cpp values used match gases.def ?"
-         check_cpp_match=.true. ! default value
-         call getin_p("check_cpp_match",check_cpp_match)
-         write(*,*) " check_cpp_match = ",check_cpp_match
-
-         write(*,*) "call radiative transfer ?"
-         callrad=.true. ! default value
-         call getin_p("callrad",callrad)
-         write(*,*) " callrad = ",callrad
-
-         write(*,*) "call correlated-k radiative transfer ?"
-         corrk=.true. ! default value
-         call getin_p("corrk",corrk)
-         write(*,*) " corrk = ",corrk
-
-         write(*,*) "prohibit calculations outside corrk T grid?"
-         strictboundcorrk=.true. ! default value
-         call getin_p("strictboundcorrk",strictboundcorrk)
-         write(*,*) "strictboundcorrk = ",strictboundcorrk
-
-         write(*,*) "call gaseous absorption in the visible bands?",
-     &              "(matters only if callrad=T)"
-         callgasvis=.false. ! default value
-         call getin_p("callgasvis",callgasvis)
-         write(*,*) " callgasvis = ",callgasvis
-        
-         write(*,*) "call continuum opacities in radiative transfer ?",
-     &              "(matters only if callrad=T)"
-         continuum=.true. ! default value
-         call getin_p("continuum",continuum)
-         write(*,*) " continuum = ",continuum
-
-         write(*,*) "use analytic function for H2O continuum ?"
-         H2Ocont_simple=.false. ! default value
-         call getin_p("H2Ocont_simple",H2Ocont_simple)
-         write(*,*) " H2Ocont_simple = ",H2Ocont_simple
- 
-         write(*,*) "call turbulent vertical diffusion ?"
-         calldifv=.true. ! default value
-         call getin_p("calldifv",calldifv)
-         write(*,*) " calldifv = ",calldifv
-
-         write(*,*) "use turbdiff instead of vdifc ?"
-         UseTurbDiff=.true. ! default value
-         call getin_p("UseTurbDiff",UseTurbDiff)
-         write(*,*) " UseTurbDiff = ",UseTurbDiff
-
-         write(*,*) "call convective adjustment ?"
-         calladj=.true. ! default value
-         call getin_p("calladj",calladj)
-         write(*,*) " calladj = ",calladj
-
-         write(*,*) "call CO2 condensation ?"
-         co2cond=.false. ! default value
-         call getin_p("co2cond",co2cond)
-         write(*,*) " co2cond = ",co2cond
-! Test of incompatibility
-         if (co2cond.and.(.not.tracer)) then
-            print*,'We need a CO2 ice tracer to condense CO2'
-            call abort
-         endif 
- 
-         write(*,*) "CO2 supersaturation level ?"
-         co2supsat=1.0 ! default value
-         call getin_p("co2supsat",co2supsat)
-         write(*,*) " co2supsat = ",co2supsat
-
-         write(*,*) "Radiative timescale for Newtonian cooling ?"
-         tau_relax=30. ! default value
-         call getin_p("tau_relax",tau_relax)
-         write(*,*) " tau_relax = ",tau_relax
-         tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
-
-         write(*,*)"call thermal conduction in the soil ?"
-         callsoil=.true. ! default value
-         call getin_p("callsoil",callsoil)
-         write(*,*) " callsoil = ",callsoil
-         
-         write(*,*)"Rad transfer is computed every iradia",
-     &             " physical timestep"
-         iradia=1 ! default value
-         call getin_p("iradia",iradia)
-         write(*,*)" iradia = ",iradia
-       
-         write(*,*)"Rayleigh scattering ?"
-         rayleigh=.false.
-         call getin_p("rayleigh",rayleigh)
-         write(*,*)" rayleigh = ",rayleigh
-
-         write(*,*) "Use blackbody for stellar spectrum ?"
-         stelbbody=.false. ! default value
-         call getin_p("stelbbody",stelbbody)
-         write(*,*) " stelbbody = ",stelbbody
-
-         write(*,*) "Stellar blackbody temperature ?"
-         stelTbb=5800.0 ! default value
-         call getin_p("stelTbb",stelTbb)
-         write(*,*) " stelTbb = ",stelTbb
-
-         write(*,*)"Output mean OLR in 1D?"
-         meanOLR=.false.
-         call getin_p("meanOLR",meanOLR)
-         write(*,*)" meanOLR = ",meanOLR
-
-         write(*,*)"Output spectral OLR in 3D?"
-         specOLR=.false.
-         call getin_p("specOLR",specOLR)
-         write(*,*)" specOLR = ",specOLR
-
-         write(*,*)"Operate in kastprof mode?"
-         kastprof=.false.
-         call getin_p("kastprof",kastprof)
-         write(*,*)" kastprof = ",kastprof
-
-         write(*,*)"Uniform absorption in radiative transfer?"
-         graybody=.false.
-         call getin_p("graybody",graybody)
-         write(*,*)" graybody = ",graybody
-
-! Slab Ocean 
-         write(*,*) "Use slab-ocean ?"
-         ok_slab_ocean=.false.         ! default value
-         call getin_p("ok_slab_ocean",ok_slab_ocean)
-         write(*,*) "ok_slab_ocean = ",ok_slab_ocean
-
-         write(*,*) "Use slab-sea-ice ?"
-         ok_slab_sic=.true.         ! default value
-         call getin_p("ok_slab_sic",ok_slab_sic)
-         write(*,*) "ok_slab_sic = ",ok_slab_sic
-
-         write(*,*) "Use heat transport for the ocean ?"
-         ok_slab_heat_transp=.true.   ! default value
-         call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
-         write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp
-
-
-
-! Test of incompatibility:
-! if kastprof used, we must be in 1D
-         if (kastprof.and.(ngrid.gt.1)) then
-           print*,'kastprof can only be used in 1D!'
-           call abort
-         endif
-
-         write(*,*)"Stratospheric temperature for kastprof mode?"
-         Tstrat=167.0
-         call getin_p("Tstrat",Tstrat)
-         write(*,*)" Tstrat = ",Tstrat
-
-         write(*,*)"Remove lower boundary?"
-         nosurf=.false.
-         call getin_p("nosurf",nosurf)
-         write(*,*)" nosurf = ",nosurf
-
-! Tests of incompatibility:
-         if (nosurf.and.callsoil) then
-           print*,'nosurf not compatible with soil scheme!'
-           print*,'... got to make a choice!'
-           call abort
-         endif
-
-         write(*,*)"Add an internal heat flux?",
-     .             "... matters only if callsoil=F"
-         intheat=0.
-         call getin_p("intheat",intheat)
-         write(*,*)" intheat = ",intheat
-
-         write(*,*)"Use Newtonian cooling for radiative transfer?"
-         newtonian=.false.
-         call getin_p("newtonian",newtonian)
-         write(*,*)" newtonian = ",newtonian
-
-! Tests of incompatibility:
-         if (newtonian.and.corrk) then
-           print*,'newtonian not compatible with correlated-k!'
-           call abort
-         endif
-         if (newtonian.and.calladj) then
-           print*,'newtonian not compatible with adjustment!'
-           call abort
-         endif
-         if (newtonian.and.calldifv) then
-           print*,'newtonian not compatible with a boundary layer!'
-           call abort
-         endif
-
-         write(*,*)"Test physics timescale in 1D?"
-         testradtimes=.false.
-         call getin_p("testradtimes",testradtimes)
-         write(*,*)" testradtimes = ",testradtimes
-
-! Test of incompatibility:
-! if testradtimes used, we must be in 1D
-         if (testradtimes.and.(ngrid.gt.1)) then
-           print*,'testradtimes can only be used in 1D!'
-           call abort
-         endif
-
-         write(*,*)"Default planetary temperature?"
-         tplanet=215.0
-         call getin_p("tplanet",tplanet)
-         write(*,*)" tplanet = ",tplanet
-
-         write(*,*)"Which star?"
-         startype=1 ! default value = Sol
-         call getin_p("startype",startype)
-         write(*,*)" startype = ",startype
-
-         write(*,*)"Value of stellar flux at 1 AU?"
-         Fat1AU=1356.0 ! default value = Sol today
-         call getin_p("Fat1AU",Fat1AU)
-         write(*,*)" Fat1AU = ",Fat1AU
-
-
-! TRACERS:
-
-         write(*,*)"Varying H2O cloud fraction?"
-         CLFvarying=.false.     ! default value
-         call getin_p("CLFvarying",CLFvarying)
-         write(*,*)" CLFvarying = ",CLFvarying
-
-         write(*,*)"Value of fixed H2O cloud fraction?"
-         CLFfixval=1.0                ! default value
-         call getin_p("CLFfixval",CLFfixval)
-         write(*,*)" CLFfixval = ",CLFfixval
-
-         write(*,*)"fixed radii for Cloud particles?"
-         radfixed=.false. ! default value
-         call getin_p("radfixed",radfixed)
-         write(*,*)" radfixed = ",radfixed
-
-         if(kastprof)then
-            radfixed=.true.
-         endif  
-
-	 write(*,*)"Number mixing ratio of CO2 ice particles:"
-         Nmix_co2=1.e6 ! default value
-         call getin_p("Nmix_co2",Nmix_co2)
-         write(*,*)" Nmix_co2 = ",Nmix_co2
-
-!         write(*,*)"Number of radiatively active aerosols:"
-!         naerkind=0. ! default value
-!         call getin_p("naerkind",naerkind)
-!         write(*,*)" naerkind = ",naerkind
-
-         write(*,*)"Opacity of dust (if used):"
-         dusttau=0. ! default value
-         call getin_p("dusttau",dusttau)
-         write(*,*)" dusttau = ",dusttau
-
-         write(*,*)"Radiatively active CO2 aerosols?"
-         aeroco2=.false.     ! default value
-         call getin_p("aeroco2",aeroco2)
-         write(*,*)" aeroco2 = ",aeroco2
-
-         write(*,*)"Fixed CO2 aerosol distribution?"
-         aerofixco2=.false.     ! default value
-         call getin_p("aerofixco2",aerofixco2)
-         write(*,*)" aerofixco2 = ",aerofixco2
-
-         write(*,*)"Radiatively active water ice?"
-         aeroh2o=.false.     ! default value
-         call getin_p("aeroh2o",aeroh2o)
-         write(*,*)" aeroh2o = ",aeroh2o
-
-         write(*,*)"Fixed H2O aerosol distribution?"
-         aerofixh2o=.false.     ! default value
-         call getin_p("aerofixh2o",aerofixh2o)
-         write(*,*)" aerofixh2o = ",aerofixh2o
-
-         write(*,*)"Radiatively active sulfuric acid aersols?"
-         aeroh2so4=.false.     ! default value
-         call getin_p("aeroh2so4",aeroh2so4)
-         write(*,*)" aeroh2so4 = ",aeroh2so4
-	 
-!=================================
-
-	 write(*,*)"Radiatively active two-layer aersols?"
-         aeroback2lay=.false.     ! default value
-         call getin_p("aeroback2lay",aeroback2lay)
-         write(*,*)" aeroback2lay = ",aeroback2lay
-
-         write(*,*)"TWOLAY AEROSOL: total optical depth ",
-     &              "in the tropospheric layer (visible)"
-         obs_tau_col_tropo=8.D0
-         call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
-         write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
-
-         write(*,*)"TWOLAY AEROSOL: total optical depth ",
-     &              "in the stratospheric layer (visible)"
-         obs_tau_col_strato=0.08D0
-         call getin_p("obs_tau_col_strato",obs_tau_col_strato)
-         write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
-
-         write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
-         pres_bottom_tropo=66000.0
-         call getin_p("pres_bottom_tropo",pres_bottom_tropo)
-         write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
-
-         write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
-         pres_top_tropo=18000.0
-         call getin_p("pres_top_tropo",pres_top_tropo)
-         write(*,*)" pres_top_tropo = ",pres_top_tropo
-
-         write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
-         pres_bottom_strato=2000.0
-         call getin_p("pres_bottom_strato",pres_bottom_strato)
-         write(*,*)" pres_bottom_strato = ",pres_bottom_strato
-
-         write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
-         pres_top_strato=100.0
-         call getin_p("pres_top_strato",pres_top_strato)
-         write(*,*)" pres_top_strato = ",pres_top_strato
-
-         write(*,*)"TWOLAY AEROSOL: particle size in the ",
-     &              "tropospheric layer, in meters"
-         size_tropo=2.e-6
-         call getin_p("size_tropo",size_tropo)
-         write(*,*)" size_tropo = ",size_tropo
-
-         write(*,*)"TWOLAY AEROSOL: particle size in the ",
-     &              "stratospheric layer, in meters"
-         size_strato=1.e-7
-         call getin_p("size_strato",size_strato)
-         write(*,*)" size_strato = ",size_strato
-
-!=================================
-
-         write(*,*)"Cloud pressure level (with kastprof only):"
-         cloudlvl=0. ! default value
-         call getin_p("cloudlvl",cloudlvl)
-         write(*,*)" cloudlvl = ",cloudlvl
-
-         write(*,*)"Is the variable gas species radiatively active?"
-         Tstrat=167.0
-         varactive=.false.
-         call getin_p("varactive",varactive)
-         write(*,*)" varactive = ",varactive
-
-         write(*,*)"Is the variable gas species distribution set?"
-         varfixed=.false.
-         call getin_p("varfixed",varfixed)
-         write(*,*)" varfixed = ",varfixed
-
-         write(*,*)"What is the saturation % of the variable species?"
-         satval=0.8
-         call getin_p("satval",satval)
-         write(*,*)" satval = ",satval
-
-
-! Test of incompatibility:
-! if varactive, then varfixed should be false
-         if (varactive.and.varfixed) then
-           print*,'if varactive, varfixed must be OFF!'
-           stop
-         endif
-
-         write(*,*) "Gravitationnal sedimentation ?"
-         sedimentation=.false. ! default value
-         call getin_p("sedimentation",sedimentation)
-         write(*,*) " sedimentation = ",sedimentation
-
-         write(*,*) "Compute water cycle ?"
-         water=.false. ! default value
-         call getin_p("water",water)
-         write(*,*) " water = ",water
-         
-! Test of incompatibility:
-! if water is true, there should be at least a tracer
-         if (water.and.(.not.tracer)) then
-           print*,'if water is ON, tracer must be ON too!'
-           stop
-         endif
-
-         write(*,*) "Include water condensation ?"
-         watercond=.false. ! default value
-         call getin_p("watercond",watercond)
-         write(*,*) " watercond = ",watercond
-
-! Test of incompatibility:
-! if watercond is used, then water should be used too
-         if (watercond.and.(.not.water)) then
-           print*,'if watercond is used, water should be used too'
-           stop
-         endif
-
-         write(*,*) "Include water precipitation ?"
-         waterrain=.false. ! default value
-         call getin_p("waterrain",waterrain)
-         write(*,*) " waterrain = ",waterrain
-
-         write(*,*) "Include surface hydrology ?"
-         hydrology=.false. ! default value
-         call getin_p("hydrology",hydrology)
-         write(*,*) " hydrology = ",hydrology
-
-         write(*,*) "Evolve surface water sources ?"
-         sourceevol=.false. ! default value
-         call getin_p("sourceevol",sourceevol)
-         write(*,*) " sourceevol = ",sourceevol
-
-         write(*,*) "Ice evolution timestep ?"
-         icetstep=100.0 ! default value
-         call getin_p("icetstep",icetstep)
-         write(*,*) " icetstep = ",icetstep
-         
-         write(*,*) "Spectral Dependant albedo ?"
-         albedo_spectral_mode=.false. ! default value
-         call getin_p("albedo_spectral_mode",albedo_spectral_mode)
-         write(*,*) " albedo_spectral_mode = ",albedo_spectral_mode
-
-         write(*,*) "Snow albedo ?"
-         write(*,*) "If albedo_spectral_mode=.true., then this "
-         write(*,*) "corresponds to the 0.5 microns snow albedo."
-         albedosnow=0.5         ! default value
-         call getin_p("albedosnow",albedosnow)
-         write(*,*) " albedosnow = ",albedosnow
-         
-         write(*,*) "CO2 ice albedo ?"
-         albedoco2ice=0.5       ! default value
-         call getin_p("albedoco2ice",albedoco2ice)
-         write(*,*) " albedoco2ice = ",albedoco2ice
-
-         write(*,*) "Maximum ice thickness ?"
-         maxicethick=2.0         ! default value
-         call getin_p("maxicethick",maxicethick)
-         write(*,*) " maxicethick = ",maxicethick
-
-         write(*,*) "Freezing point of seawater ?"
-         Tsaldiff=-1.8          ! default value
-         call getin_p("Tsaldiff",Tsaldiff)
-         write(*,*) " Tsaldiff = ",Tsaldiff
-
-         write(*,*) "Does user want to force cpp and mugaz?"
-         force_cpp=.false. ! default value
-         call getin_p("force_cpp",force_cpp)
-         write(*,*) " force_cpp = ",force_cpp
-
-         if (force_cpp) then
-           mugaz = -99999.
-           PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
-           call getin_p("mugaz",mugaz)
-           IF (mugaz.eq.-99999.) THEN
-               PRINT *, "mugaz must be set if force_cpp = T"
-               STOP
-           ELSE
-               write(*,*) "mugaz=",mugaz
-           ENDIF
-           cpp = -99999.
-           PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
-           call getin_p("cpp",cpp)
-           IF (cpp.eq.-99999.) THEN
-               PRINT *, "cpp must be set if force_cpp = T"
-               STOP
-           ELSE
-               write(*,*) "cpp=",cpp
-           ENDIF
-!         else
-!           mugaz=8.314*1000./pr
-         endif
-         call su_gases
-         call calc_cpp_mugaz
-
-         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*
-
-
-!-----------------------------------------------------------------------
-!     Some more initialization:
-!     ------------------------
-
-      ! ALLOCATE ARRAYS IN comgeomfi_h
-      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid))
-      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid))
-      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid))
-
-      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)
-      call planetwide_sumval(area,totarea_planet)
-
-      !! those are defined in comdiurn_h.F90
-      IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
-      IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
-      IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
-      IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
-
-      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_mod
-
-      ! allocate "comsoil_h" arrays
-      call ini_comsoil_h(ngrid)
-      
-      END
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 1524)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 1524)
@@ -0,0 +1,722 @@
+MODULE inifis_mod
+IMPLICIT NONE
+
+CONTAINS
+
+  SUBROUTINE inifis(ngrid,nlayer,nq, &
+             day_ini,pdaysec,ptimestep, &
+             plat,plon,parea, &
+             prad,pg,pr,pcpp)
+
+  use radinc_h, only : naerkind
+  use datafile_mod, only: datadir
+  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
+  use comgeomfi_h, only: long, lati, area, totarea, totarea_planet
+  use comsoil_h, only: ini_comsoil_h
+  use control_mod, only: ecritphy
+  use comcstfi_mod, only: rad, cpp, g, r, rcp, &
+                          mugaz, pi, avocado
+  use planete_mod, only: nres
+  use planetwide_mod, only: planetwide_sumval
+  use callkeys_mod
+  use time_phylmdz_mod, only: init_time, daysec, dtphys
+
+!=======================================================================
+!
+!   purpose:
+!   -------
+!
+!   Initialisation for the physical parametrisations of the LMD 
+!   Generic Model.
+!
+!   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
+!
+!
+!   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:
+!   -------------
+  use datafile_mod, only: datadir
+  use ioipsl_getin_p_mod, only: getin_p
+  IMPLICIT NONE
+
+
+
+  REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
+ 
+  INTEGER,INTENT(IN) :: ngrid,nlayer,nq
+  REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
+  integer,intent(in) :: day_ini
+  INTEGER ig,ierr
+ 
+  EXTERNAL iniorbit,orbite
+  EXTERNAL SSUM
+  REAL SSUM
+ 
+  CHARACTER ch1*12
+  CHARACTER ch80*80
+
+  logical chem, h2o
+
+  real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
+
+  ! initialize constants in comcstfi_mod
+  rad=prad
+  cpp=pcpp
+  g=pg
+  r=pr
+  rcp=r/cpp
+  pi=2.*asin(1.)
+  avocado = 6.02214179e23   ! added by RW
+
+  ! Initialize some "temporal and calendar" related variables
+  CALL init_time(day_ini,pdaysec,ptimestep)
+
+  ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
+  ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
+  call getin_p("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) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
+      
+!!!      IF(ierr.EQ.0) THEN
+  IF(iscallphys) THEN
+     PRINT*
+     PRINT*
+     PRINT*,'--------------------------------------------'
+     PRINT*,' inifis: Parametres pour la physique (callphys.def)'
+     PRINT*,'--------------------------------------------'
+
+     write(*,*) "Directory where external input files are:"
+     ! default 'datadir' is set in "datadir_mod"
+     call getin_p("datadir",datadir) ! default path
+     write(*,*) " datadir = ",trim(datadir)
+
+     write(*,*) "Run with or without tracer transport ?"
+     tracer=.false. ! default value
+     call getin_p("tracer",tracer)
+     write(*,*) " tracer = ",tracer
+
+     write(*,*) "Run with or without atm mass update ", &
+            " due to tracer evaporation/condensation?"
+     mass_redistrib=.false. ! default value
+     call getin_p("mass_redistrib",mass_redistrib)
+     write(*,*) " mass_redistrib = ",mass_redistrib
+
+     write(*,*) "Diurnal cycle ?"
+     write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
+     diurnal=.true. ! default value
+     call getin_p("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_p("season",season)
+     write(*,*) " season = ",season
+
+     write(*,*) "Tidally resonant rotation ?"
+     tlocked=.false. ! default value
+     call getin_p("tlocked",tlocked)
+     write(*,*) "tlocked = ",tlocked
+
+     write(*,*) "Saturn ring shadowing ?"
+     rings_shadow = .false.
+     call getin_p("rings_shadow", rings_shadow)
+     write(*,*) "rings_shadow = ", rings_shadow
+         
+     write(*,*) "Compute latitude-dependent gravity field?"
+     oblate = .false.
+     call getin_p("oblate", oblate)
+     write(*,*) "oblate = ", oblate
+
+     write(*,*) "Flattening of the planet (a-b)/a "
+     flatten = 0.0
+     call getin_p("flatten", flatten)
+     write(*,*) "flatten = ", flatten
+         
+
+     write(*,*) "Needed if oblate=.true.: J2"
+     J2 = 0.0
+     call getin_p("J2", J2)
+     write(*,*) "J2 = ", J2
+         
+     write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
+     MassPlanet = 0.0
+     call getin_p("MassPlanet", MassPlanet)
+     write(*,*) "MassPlanet = ", MassPlanet         
+
+     write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
+     Rmean = 0.0
+     call getin_p("Rmean", Rmean)
+     write(*,*) "Rmean = ", Rmean
+         
+! Test of incompatibility:
+! if tlocked, then diurnal should be false
+     if (tlocked.and.diurnal) then
+       print*,'If diurnal=true, we should turn off tlocked.'
+       stop
+     endif
+
+     write(*,*) "Tidal resonance ratio ?"
+     nres=0          ! default value
+     call getin_p("nres",nres)
+     write(*,*) "nres = ",nres
+
+     write(*,*) "Write some extra output to the screen ?"
+     lwrite=.false. ! default value
+     call getin_p("lwrite",lwrite)
+     write(*,*) " lwrite = ",lwrite
+
+     write(*,*) "Save statistics in file stats.nc ?"
+     callstats=.true. ! default value
+     call getin_p("callstats",callstats)
+     write(*,*) " callstats = ",callstats
+
+     write(*,*) "Test energy conservation of model physics ?"
+     enertest=.false. ! default value
+     call getin_p("enertest",enertest)
+     write(*,*) " enertest = ",enertest
+
+     write(*,*) "Check to see if cpp values used match gases.def ?"
+     check_cpp_match=.true. ! default value
+     call getin_p("check_cpp_match",check_cpp_match)
+     write(*,*) " check_cpp_match = ",check_cpp_match
+
+     write(*,*) "call radiative transfer ?"
+     callrad=.true. ! default value
+     call getin_p("callrad",callrad)
+     write(*,*) " callrad = ",callrad
+
+     write(*,*) "call correlated-k radiative transfer ?"
+     corrk=.true. ! default value
+     call getin_p("corrk",corrk)
+     write(*,*) " corrk = ",corrk
+
+     write(*,*) "prohibit calculations outside corrk T grid?"
+     strictboundcorrk=.true. ! default value
+     call getin_p("strictboundcorrk",strictboundcorrk)
+     write(*,*) "strictboundcorrk = ",strictboundcorrk
+
+     write(*,*) "call gaseous absorption in the visible bands?", &
+                    "(matters only if callrad=T)"
+     callgasvis=.false. ! default value
+     call getin_p("callgasvis",callgasvis)
+     write(*,*) " callgasvis = ",callgasvis
+        
+     write(*,*) "call continuum opacities in radiative transfer ?", &
+                    "(matters only if callrad=T)"
+     continuum=.true. ! default value
+     call getin_p("continuum",continuum)
+     write(*,*) " continuum = ",continuum
+
+     write(*,*) "use analytic function for H2O continuum ?"
+     H2Ocont_simple=.false. ! default value
+     call getin_p("H2Ocont_simple",H2Ocont_simple)
+     write(*,*) " H2Ocont_simple = ",H2Ocont_simple
+ 
+     write(*,*) "call turbulent vertical diffusion ?"
+     calldifv=.true. ! default value
+     call getin_p("calldifv",calldifv)
+     write(*,*) " calldifv = ",calldifv
+
+     write(*,*) "use turbdiff instead of vdifc ?"
+     UseTurbDiff=.true. ! default value
+     call getin_p("UseTurbDiff",UseTurbDiff)
+     write(*,*) " UseTurbDiff = ",UseTurbDiff
+
+     write(*,*) "call convective adjustment ?"
+     calladj=.true. ! default value
+     call getin_p("calladj",calladj)
+     write(*,*) " calladj = ",calladj
+
+     write(*,*) "call CO2 condensation ?"
+     co2cond=.false. ! default value
+     call getin_p("co2cond",co2cond)
+     write(*,*) " co2cond = ",co2cond
+! Test of incompatibility
+     if (co2cond.and.(.not.tracer)) then
+        print*,'We need a CO2 ice tracer to condense CO2'
+        call abort
+     endif 
+ 
+     write(*,*) "CO2 supersaturation level ?"
+     co2supsat=1.0 ! default value
+     call getin_p("co2supsat",co2supsat)
+     write(*,*) " co2supsat = ",co2supsat
+
+     write(*,*) "Radiative timescale for Newtonian cooling ?"
+     tau_relax=30. ! default value
+     call getin_p("tau_relax",tau_relax)
+     write(*,*) " tau_relax = ",tau_relax
+     tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
+
+     write(*,*)"call thermal conduction in the soil ?"
+     callsoil=.true. ! default value
+     call getin_p("callsoil",callsoil)
+     write(*,*) " callsoil = ",callsoil
+         
+     write(*,*)"Rad transfer is computed every iradia", &
+                   " physical timestep"
+     iradia=1 ! default value
+     call getin_p("iradia",iradia)
+     write(*,*)" iradia = ",iradia
+       
+     write(*,*)"Rayleigh scattering ?"
+     rayleigh=.false.
+     call getin_p("rayleigh",rayleigh)
+     write(*,*)" rayleigh = ",rayleigh
+
+     write(*,*) "Use blackbody for stellar spectrum ?"
+     stelbbody=.false. ! default value
+     call getin_p("stelbbody",stelbbody)
+     write(*,*) " stelbbody = ",stelbbody
+
+     write(*,*) "Stellar blackbody temperature ?"
+     stelTbb=5800.0 ! default value
+     call getin_p("stelTbb",stelTbb)
+     write(*,*) " stelTbb = ",stelTbb
+
+     write(*,*)"Output mean OLR in 1D?"
+     meanOLR=.false.
+     call getin_p("meanOLR",meanOLR)
+     write(*,*)" meanOLR = ",meanOLR
+
+     write(*,*)"Output spectral OLR in 3D?"
+     specOLR=.false.
+     call getin_p("specOLR",specOLR)
+     write(*,*)" specOLR = ",specOLR
+
+     write(*,*)"Operate in kastprof mode?"
+     kastprof=.false.
+     call getin_p("kastprof",kastprof)
+     write(*,*)" kastprof = ",kastprof
+
+     write(*,*)"Uniform absorption in radiative transfer?"
+     graybody=.false.
+     call getin_p("graybody",graybody)
+     write(*,*)" graybody = ",graybody
+
+! Slab Ocean 
+     write(*,*) "Use slab-ocean ?"
+     ok_slab_ocean=.false.         ! default value
+     call getin_p("ok_slab_ocean",ok_slab_ocean)
+     write(*,*) "ok_slab_ocean = ",ok_slab_ocean
+
+     write(*,*) "Use slab-sea-ice ?"
+     ok_slab_sic=.true.         ! default value
+     call getin_p("ok_slab_sic",ok_slab_sic)
+     write(*,*) "ok_slab_sic = ",ok_slab_sic
+
+     write(*,*) "Use heat transport for the ocean ?"
+     ok_slab_heat_transp=.true.   ! default value
+     call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
+     write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp
+
+
+
+! Test of incompatibility:
+! if kastprof used, we must be in 1D
+     if (kastprof.and.(ngrid.gt.1)) then
+       print*,'kastprof can only be used in 1D!'
+       call abort
+     endif
+
+     write(*,*)"Stratospheric temperature for kastprof mode?"
+     Tstrat=167.0
+     call getin_p("Tstrat",Tstrat)
+     write(*,*)" Tstrat = ",Tstrat
+
+     write(*,*)"Remove lower boundary?"
+     nosurf=.false.
+     call getin_p("nosurf",nosurf)
+     write(*,*)" nosurf = ",nosurf
+
+! Tests of incompatibility:
+     if (nosurf.and.callsoil) then
+       print*,'nosurf not compatible with soil scheme!'
+       print*,'... got to make a choice!'
+       call abort
+     endif
+
+     write(*,*)"Add an internal heat flux?", &
+                   "... matters only if callsoil=F"
+     intheat=0.
+     call getin_p("intheat",intheat)
+     write(*,*)" intheat = ",intheat
+
+     write(*,*)"Use Newtonian cooling for radiative transfer?"
+     newtonian=.false.
+     call getin_p("newtonian",newtonian)
+     write(*,*)" newtonian = ",newtonian
+
+! Tests of incompatibility:
+     if (newtonian.and.corrk) then
+       print*,'newtonian not compatible with correlated-k!'
+       call abort
+     endif
+     if (newtonian.and.calladj) then
+       print*,'newtonian not compatible with adjustment!'
+       call abort
+     endif
+     if (newtonian.and.calldifv) then
+       print*,'newtonian not compatible with a boundary layer!'
+       call abort
+     endif
+
+     write(*,*)"Test physics timescale in 1D?"
+     testradtimes=.false.
+     call getin_p("testradtimes",testradtimes)
+     write(*,*)" testradtimes = ",testradtimes
+
+! Test of incompatibility:
+! if testradtimes used, we must be in 1D
+     if (testradtimes.and.(ngrid.gt.1)) then
+       print*,'testradtimes can only be used in 1D!'
+       call abort
+     endif
+
+     write(*,*)"Default planetary temperature?"
+     tplanet=215.0
+     call getin_p("tplanet",tplanet)
+     write(*,*)" tplanet = ",tplanet
+
+     write(*,*)"Which star?"
+     startype=1 ! default value = Sol
+     call getin_p("startype",startype)
+     write(*,*)" startype = ",startype
+
+     write(*,*)"Value of stellar flux at 1 AU?"
+     Fat1AU=1356.0 ! default value = Sol today
+     call getin_p("Fat1AU",Fat1AU)
+     write(*,*)" Fat1AU = ",Fat1AU
+
+
+! TRACERS:
+
+     write(*,*)"Varying H2O cloud fraction?"
+     CLFvarying=.false.     ! default value
+     call getin_p("CLFvarying",CLFvarying)
+     write(*,*)" CLFvarying = ",CLFvarying
+
+     write(*,*)"Value of fixed H2O cloud fraction?"
+     CLFfixval=1.0                ! default value
+     call getin_p("CLFfixval",CLFfixval)
+     write(*,*)" CLFfixval = ",CLFfixval
+
+     write(*,*)"fixed radii for Cloud particles?"
+     radfixed=.false. ! default value
+     call getin_p("radfixed",radfixed)
+     write(*,*)" radfixed = ",radfixed
+
+     if(kastprof)then
+        radfixed=.true.
+     endif  
+
+     write(*,*)"Number mixing ratio of CO2 ice particles:"
+     Nmix_co2=1.e6 ! default value
+     call getin_p("Nmix_co2",Nmix_co2)
+     write(*,*)" Nmix_co2 = ",Nmix_co2
+
+!         write(*,*)"Number of radiatively active aerosols:"
+!         naerkind=0. ! default value
+!         call getin_p("naerkind",naerkind)
+!         write(*,*)" naerkind = ",naerkind
+
+     write(*,*)"Opacity of dust (if used):"
+     dusttau=0. ! default value
+     call getin_p("dusttau",dusttau)
+     write(*,*)" dusttau = ",dusttau
+
+     write(*,*)"Radiatively active CO2 aerosols?"
+     aeroco2=.false.     ! default value
+     call getin_p("aeroco2",aeroco2)
+     write(*,*)" aeroco2 = ",aeroco2
+
+     write(*,*)"Fixed CO2 aerosol distribution?"
+     aerofixco2=.false.     ! default value
+     call getin_p("aerofixco2",aerofixco2)
+     write(*,*)" aerofixco2 = ",aerofixco2
+
+     write(*,*)"Radiatively active water ice?"
+     aeroh2o=.false.     ! default value
+     call getin_p("aeroh2o",aeroh2o)
+     write(*,*)" aeroh2o = ",aeroh2o
+
+     write(*,*)"Fixed H2O aerosol distribution?"
+     aerofixh2o=.false.     ! default value
+     call getin_p("aerofixh2o",aerofixh2o)
+     write(*,*)" aerofixh2o = ",aerofixh2o
+
+     write(*,*)"Radiatively active sulfuric acid aersols?"
+     aeroh2so4=.false.     ! default value
+     call getin_p("aeroh2so4",aeroh2so4)
+     write(*,*)" aeroh2so4 = ",aeroh2so4
+	 
+!=================================
+
+     write(*,*)"Radiatively active two-layer aersols?"
+     aeroback2lay=.false.     ! default value
+     call getin_p("aeroback2lay",aeroback2lay)
+     write(*,*)" aeroback2lay = ",aeroback2lay
+
+     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
+                    "in the tropospheric layer (visible)"
+     obs_tau_col_tropo=8.D0
+     call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
+     write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
+
+     write(*,*)"TWOLAY AEROSOL: total optical depth ", &
+                    "in the stratospheric layer (visible)"
+     obs_tau_col_strato=0.08D0
+     call getin_p("obs_tau_col_strato",obs_tau_col_strato)
+     write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
+
+     write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
+     pres_bottom_tropo=66000.0
+     call getin_p("pres_bottom_tropo",pres_bottom_tropo)
+     write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
+
+     write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
+     pres_top_tropo=18000.0
+     call getin_p("pres_top_tropo",pres_top_tropo)
+     write(*,*)" pres_top_tropo = ",pres_top_tropo
+
+     write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
+     pres_bottom_strato=2000.0
+     call getin_p("pres_bottom_strato",pres_bottom_strato)
+     write(*,*)" pres_bottom_strato = ",pres_bottom_strato
+
+     write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
+     pres_top_strato=100.0
+     call getin_p("pres_top_strato",pres_top_strato)
+     write(*,*)" pres_top_strato = ",pres_top_strato
+
+     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
+                    "tropospheric layer, in meters"
+     size_tropo=2.e-6
+     call getin_p("size_tropo",size_tropo)
+     write(*,*)" size_tropo = ",size_tropo
+
+     write(*,*)"TWOLAY AEROSOL: particle size in the ", &
+                    "stratospheric layer, in meters"
+     size_strato=1.e-7
+     call getin_p("size_strato",size_strato)
+     write(*,*)" size_strato = ",size_strato
+
+!=================================
+
+     write(*,*)"Cloud pressure level (with kastprof only):"
+     cloudlvl=0. ! default value
+     call getin_p("cloudlvl",cloudlvl)
+     write(*,*)" cloudlvl = ",cloudlvl
+
+     write(*,*)"Is the variable gas species radiatively active?"
+     Tstrat=167.0
+     varactive=.false.
+     call getin_p("varactive",varactive)
+     write(*,*)" varactive = ",varactive
+
+     write(*,*)"Is the variable gas species distribution set?"
+     varfixed=.false.
+     call getin_p("varfixed",varfixed)
+     write(*,*)" varfixed = ",varfixed
+
+     write(*,*)"What is the saturation % of the variable species?"
+     satval=0.8
+     call getin_p("satval",satval)
+     write(*,*)" satval = ",satval
+
+
+! Test of incompatibility:
+! if varactive, then varfixed should be false
+     if (varactive.and.varfixed) then
+       print*,'if varactive, varfixed must be OFF!'
+       stop
+     endif
+
+     write(*,*) "Gravitationnal sedimentation ?"
+     sedimentation=.false. ! default value
+     call getin_p("sedimentation",sedimentation)
+     write(*,*) " sedimentation = ",sedimentation
+
+     write(*,*) "Compute water cycle ?"
+     water=.false. ! default value
+     call getin_p("water",water)
+     write(*,*) " water = ",water
+         
+! Test of incompatibility:
+! if water is true, there should be at least a tracer
+     if (water.and.(.not.tracer)) then
+        print*,'if water is ON, tracer must be ON too!'
+        stop
+     endif
+
+     write(*,*) "Include water condensation ?"
+     watercond=.false. ! default value
+     call getin_p("watercond",watercond)
+     write(*,*) " watercond = ",watercond
+
+! Test of incompatibility:
+! if watercond is used, then water should be used too
+     if (watercond.and.(.not.water)) then
+        print*,'if watercond is used, water should be used too'
+        stop
+     endif
+
+     write(*,*) "Include water precipitation ?"
+     waterrain=.false. ! default value
+     call getin_p("waterrain",waterrain)
+     write(*,*) " waterrain = ",waterrain
+
+     write(*,*) "Include surface hydrology ?"
+     hydrology=.false. ! default value
+     call getin_p("hydrology",hydrology)
+     write(*,*) " hydrology = ",hydrology
+
+     write(*,*) "Evolve surface water sources ?"
+     sourceevol=.false. ! default value
+     call getin_p("sourceevol",sourceevol)
+     write(*,*) " sourceevol = ",sourceevol
+
+     write(*,*) "Ice evolution timestep ?"
+     icetstep=100.0 ! default value
+     call getin_p("icetstep",icetstep)
+     write(*,*) " icetstep = ",icetstep
+         
+     write(*,*) "Spectral Dependant albedo ?"
+     albedo_spectral_mode=.false. ! default value
+     call getin_p("albedo_spectral_mode",albedo_spectral_mode)
+     write(*,*) " albedo_spectral_mode = ",albedo_spectral_mode
+
+     write(*,*) "Snow albedo ?"
+     write(*,*) "If albedo_spectral_mode=.true., then this "
+     write(*,*) "corresponds to the 0.5 microns snow albedo."
+     albedosnow=0.5         ! default value
+     call getin_p("albedosnow",albedosnow)
+     write(*,*) " albedosnow = ",albedosnow
+         
+     write(*,*) "CO2 ice albedo ?"
+     albedoco2ice=0.5       ! default value
+     call getin_p("albedoco2ice",albedoco2ice)
+     write(*,*) " albedoco2ice = ",albedoco2ice
+
+     write(*,*) "Maximum ice thickness ?"
+     maxicethick=2.0         ! default value
+     call getin_p("maxicethick",maxicethick)
+     write(*,*) " maxicethick = ",maxicethick
+
+     write(*,*) "Freezing point of seawater ?"
+     Tsaldiff=-1.8          ! default value
+     call getin_p("Tsaldiff",Tsaldiff)
+     write(*,*) " Tsaldiff = ",Tsaldiff
+
+     write(*,*) "Does user want to force cpp and mugaz?"
+     force_cpp=.false. ! default value
+     call getin_p("force_cpp",force_cpp)
+     write(*,*) " force_cpp = ",force_cpp
+
+     if (force_cpp) then
+       mugaz = -99999.
+       PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
+       call getin_p("mugaz",mugaz)
+       IF (mugaz.eq.-99999.) THEN
+           PRINT *, "mugaz must be set if force_cpp = T"
+           STOP
+       ELSE
+           write(*,*) "mugaz=",mugaz
+       ENDIF
+       cpp = -99999.
+       PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
+       call getin_p("cpp",cpp)
+       IF (cpp.eq.-99999.) THEN
+           PRINT *, "cpp must be set if force_cpp = T"
+           STOP
+       ELSE
+           write(*,*) "cpp=",cpp
+       ENDIF
+!         else
+!           mugaz=8.314*1000./pr
+     endif ! of if (force_cpp)
+     call su_gases
+     call calc_cpp_mugaz
+
+     PRINT*,'--------------------------------------------'
+     PRINT*
+     PRINT*
+  ELSE
+     write(*,*)
+     write(*,*) 'Cannot read file callphys.def. Is it here ?'
+     stop
+  ENDIF ! of IF(iscallphys)
+
+  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*
+
+
+!-----------------------------------------------------------------------
+!     Some more initialization:
+!     ------------------------
+
+  ! ALLOCATE ARRAYS IN comgeomfi_h
+  IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid))
+  IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid))
+  IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid))
+
+  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)
+  call planetwide_sumval(area,totarea_planet)
+
+  !! those are defined in comdiurn_h.F90
+  IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
+  IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
+  IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
+  IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
+
+  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
+
+  ! allocate "comsoil_h" arrays
+  call ini_comsoil_h(ngrid)
+      
+  END SUBROUTINE inifis
+
+END MODULE inifis_mod
Index: trunk/LMDZ.GENERIC/libf/phystd/iniwrite.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniwrite.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniwrite.F	(revision 1524)
@@ -2,6 +2,6 @@
 
       use comsoil_h, only: mlayer, nsoilmx
-      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec, dtphys,
-     &                        pi
+      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
+      use time_phylmdz_mod, only: daysec, dtphys
       USE comvert_mod, ONLY: ap,bp,aps,bps,pseudoalt
       USE logic_mod, ONLY: fxyhypb,ysinus
Index: trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 1524)
@@ -3,6 +3,6 @@
       use radinc_h, only: L_NSPECTI
       use radcommon_h, only: WNOI,DWNI
-      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec, dtphys,
-     &                        pi
+      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
+      use time_phylmdz_mod, only: daysec, dtphys
       USE logic_mod, ONLY: fxyhypb,ysinus
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
Index: trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 1524)
@@ -4,6 +4,6 @@
       use radcommon_h, only: WNOV,DWNV
       use comsoil_h
-      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec, dtphys,
-     &                        pi
+      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
+      use time_phylmdz_mod, only: daysec, dtphys
       USE logic_mod, ONLY: fxyhypb,ysinus
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
Index: trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90	(revision 1524)
@@ -19,5 +19,6 @@
   use planete_mod, only: year_day, periastr, apoastr, peri_day, &
                          obliquit, z0, lmixmin, emin_turb
-  use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, daysec
+  use comcstfi_mod, only: rad, omeg, g, mugaz, rcp
+  use time_phylmdz_mod, only: daysec
 
   implicit none
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1524)
@@ -34,5 +34,6 @@
       use planete_mod, only: apoastr, periastr, year_day, peri_day, &
                             obliquit, nres, z0
-      use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp, daysec
+      use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
+      use time_phylmdz_mod, only: daysec
       use callkeys_mod
       implicit none
@@ -127,7 +128,7 @@
 !     Authors
 !     -------
-!           Frederic Hourdin	15/10/93
-!           Francois Forget	1994
-!           Christophe Hourdin	02/1997 
+!           Frederic Hourdin        15/10/93
+!           Francois Forget        1994
+!           Christophe Hourdin        02/1997 
 !           Subroutine completely rewritten by F. Forget (01/2000)
 !           Water ice clouds: Franck Montmessin (update 06/2003)
@@ -249,5 +250,5 @@
       
 !$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&
-	
+        
         !$OMP zdtlw,zdtsw,sensibFlux)
 
@@ -270,5 +271,5 @@
       real dtlscale(ngrid,nlayer)                             ! Largescale routine.
       real zdtc(ngrid,nlayer)                                 ! Condense_co2 routine.
-      real zdtdif(ngrid,nlayer)	                              ! Turbdiff/vdifc routines.
+      real zdtdif(ngrid,nlayer)                                      ! Turbdiff/vdifc routines.
       real zdtmr(ngrid,nlayer)                                ! Mass_redistribution routine.
       real zdtrain(ngrid,nlayer)                              ! Rain routine.
@@ -442,8 +443,8 @@
         ALLOCATE(tsoil(ngrid,nsoilmx))    
         ALLOCATE(albedo(ngrid,L_NSPECTV))
- 	ALLOCATE(albedo_equivalent(ngrid))	
- 	ALLOCATE(albedo_snow_SPECTV(L_NSPECTV))
- 	ALLOCATE(albedo_co2_ice_SPECTV(L_NSPECTV))
- 	ALLOCATE(albedo_bareground(ngrid))           
+         ALLOCATE(albedo_equivalent(ngrid))        
+         ALLOCATE(albedo_snow_SPECTV(L_NSPECTV))
+         ALLOCATE(albedo_co2_ice_SPECTV(L_NSPECTV))
+         ALLOCATE(albedo_bareground(ngrid))           
         ALLOCATE(rnat(ngrid))         
         ALLOCATE(emis(ngrid))   
@@ -534,7 +535,7 @@
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          albedo(:,:)=0.0
- 	 albedo_bareground(:)=0.0
- 	 albedo_snow_SPECTV(:)=0.0
- 	 albedo_co2_ice_SPECTV(:)=0.0
+          albedo_bareground(:)=0.0
+          albedo_snow_SPECTV(:)=0.0
+          albedo_co2_ice_SPECTV(:)=0.0
          call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 
          
@@ -736,5 +737,5 @@
               zlss=Mod(-(2.*pi*(zday/year_day)*nres - right_ascen),2.*pi)
       elseif (diurnal) then
-	      zlss=-2.*pi*(zday-.5)
+              zlss=-2.*pi*(zday-.5)
       else if(diurnal .eqv. .false.) then
               zlss=9999.
@@ -782,5 +783,5 @@
             zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
             mass(ig,l)  = (pplev(ig,l) - pplev(ig,l+1))/glat(ig)
-	    massarea(ig,l)=mass(ig,l)*area(ig)
+            massarea(ig,l)=mass(ig,l)*area(ig)
          enddo
       enddo
@@ -845,11 +846,11 @@
                   call abort
                endif
-	       if(water) then
+               if(water) then
                   muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap)) 
-		  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) 
-		  ! take into account water effect on mean molecular weight
-	       else
+                  muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap)) 
+                  ! take into account water effect on mean molecular weight
+               else
                   muvar(1:ngrid,1:nlayer+1)=mugaz
-	       endif 
+               endif 
       
 
@@ -904,6 +905,6 @@
                      zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) + tf*zdtsw(ig,1:nlayer)
 
-                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
-                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)		       
+                     OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)                       
+                     OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)                       
                   enddo                               
 
@@ -978,12 +979,12 @@
          !----------------------------
          if(enertest)then
-	    call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
-	    call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW)
-	    !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
+            call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
+            call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW)
+            !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
             call planetwide_sumval(fluxsurfabs_sw(:)*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
-	    call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:)/totarea_planet,dEtotsLW)
-	    dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
-	    dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
-	    if (is_master) then
+            call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:)/totarea_planet,dEtotsLW)
+            dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
+            dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
+            if (is_master) then
                print*,'---------------------------------------------------------------'
                print*,'In corrk SW atmospheric heating       =',dEtotSW,' W m-2'
@@ -993,5 +994,5 @@
                print*,'In corrk LW surface heating           =',dEtotsLW,' W m-2'
                print*,'surface net rad heating (SW+LW)       =',dEtotsLW+dEtotsSW,' W m-2'
-	    endif
+            endif
          endif ! end of 'enertest'
 
@@ -1009,6 +1010,6 @@
 
          ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
-	 if (UseTurbDiff) then
-	 
+         if (UseTurbDiff) then
+         
             call turbdiff(ngrid,nlayer,nq,rnat,                  &
                           ptimestep,capcal,lwrite,               &
@@ -1017,9 +1018,9 @@
                           pdt,pdq,zflubid,                       &
                           zdudif,zdvdif,zdtdif,zdtsdif,          &
-	                  sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
+                          sensibFlux,q2,zdqdif,zdqevap,zdqsdif,  &
                           taux,tauy,lastcall)
 
-	 else
-	 
+         else
+         
             zdh(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)/zpopsk(1:ngrid,1:nlayer)
  
@@ -1030,9 +1031,9 @@
                        zdh,pdq,zflubid,                       &
                        zdudif,zdvdif,zdhdif,zdtsdif,          &
-	               sensibFlux,q2,zdqdif,zdqsdif,          &
+                       sensibFlux,q2,zdqdif,zdqsdif,          &
                        taux,tauy,lastcall)
 
             zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
-	    zdqevap(1:ngrid,1:nlayer)=0.
+            zdqevap(1:ngrid,1:nlayer)=0.
 
          end if !end of 'UseTurbDiff'
@@ -1060,27 +1061,27 @@
          if(enertest)then
          
-	    dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
+            dEzdiff(:,:)=cpp*mass(:,:)*zdtdif(:,:)
             do ig = 1, ngrid
-	       dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
-	       dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
+               dEdiff(ig)=SUM(dEzdiff (ig,:))+ sensibFlux(ig)! subtract flux to the ground
+               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
             enddo
             
-	    call planetwide_sumval(dEdiff(:)*area(:)/totarea_planet,dEtot)
-	    dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
-	    call planetwide_sumval(dEdiffs(:)*area(:)/totarea_planet,dEtots)
-	    call planetwide_sumval(sensibFlux(:)*area(:)/totarea_planet,AtmToSurf_TurbFlux)
+            call planetwide_sumval(dEdiff(:)*area(:)/totarea_planet,dEtot)
+            dEdiffs(:)=capcal(:)*zdtsdif(:)-zflubid(:)-sensibFlux(:)
+            call planetwide_sumval(dEdiffs(:)*area(:)/totarea_planet,dEtots)
+            call planetwide_sumval(sensibFlux(:)*area(:)/totarea_planet,AtmToSurf_TurbFlux)
             
-	    if (is_master) then
+            if (is_master) then
             
                if (UseTurbDiff) then
-	       	  print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
-               	  print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
+                         print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
+                         print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
                   print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
-	       else
-	       	  print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
-               	  print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
-               	  print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
-	       end if
-	    endif ! end of 'is_master'
+               else
+                         print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
+                         print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
+                         print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+               end if
+            endif ! end of 'is_master'
             
          ! JL12 : note that the black body radiative flux emitted by the surface has been updated by the implicit scheme but not given back elsewhere.
@@ -1091,19 +1092,19 @@
          if(watertest.and.water)then
          
-	    call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
-	    call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots_tmp)
+            call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
+            call planetwide_sumval(zdqsdif(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots_tmp)
             do ig = 1, ngrid
-	       vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
-	    enddo
-	    call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
-	    call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
-	    dWtot = dWtot + dWtot_tmp
-	    dWtots = dWtots + dWtots_tmp
+               vdifcncons(ig)=SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_vap))
+            enddo
+            call planetwide_sumval(massarea(:,:)*zdqdif(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
+            call planetwide_sumval(zdqsdif(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
+            dWtot = dWtot + dWtot_tmp
+            dWtots = dWtots + dWtots_tmp
             do ig = 1, ngrid
-	       vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
-	    enddo	    
-	    call planetwide_maxval(vdifcncons(:),nconsMAX)
-
-	    if (is_master) then
+               vdifcncons(ig)=vdifcncons(ig) + SUM(mass(ig,:)*zdqdif(ig,:,igcm_h2o_ice))
+            enddo            
+            call planetwide_maxval(vdifcncons(:),nconsMAX)
+
+            if (is_master) then
                print*,'---------------------------------------------------------------'
                print*,'In difv atmospheric water change        =',dWtot,' kg m-2'
@@ -1111,5 +1112,5 @@
                print*,'In difv non-cons factor                 =',dWtot+dWtots,' kg m-2'
                print*,'In difv MAX non-cons factor             =',nconsMAX,' kg m-2 s-1'
-	    endif
+            endif
 
          endif ! end of 'watertest'
@@ -1157,5 +1158,5 @@
          ! Test energy conservation
          if(enertest)then
-	    call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
+            call planetwide_sumval(cpp*massarea(:,:)*zdtadj(:,:)/totarea_planet,dEtot)
             if (is_master) print*,'In convadj atmospheric energy change  =',dEtot,' W m-2'
          endif
@@ -1163,19 +1164,19 @@
          ! Test water conservation
          if(watertest)then
-	    call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
+            call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dWtot_tmp)
             do ig = 1, ngrid
-	       cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
-	    enddo
-	    call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
-	    dWtot = dWtot + dWtot_tmp
+               cadjncons(ig)=SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_vap))
+            enddo
+            call planetwide_sumval(massarea(:,:)*zdqadj(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
+            dWtot = dWtot + dWtot_tmp
             do ig = 1, ngrid
-	       cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
-	    enddo	    
-	    call planetwide_maxval(cadjncons(:),nconsMAX)
+               cadjncons(ig)=cadjncons(ig) + SUM(mass(ig,:)*zdqadj(ig,:,igcm_h2o_ice))
+            enddo            
+            call planetwide_maxval(cadjncons(:),nconsMAX)
 
             if (is_master) then
-	       print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
+               print*,'In convadj atmospheric water change     =',dWtot,' kg m-2'
                print*,'In convadj MAX non-cons factor          =',nconsMAX,' kg m-2 s-1'
-	    endif
+            endif
             
          endif ! end of 'watertest'
@@ -1208,10 +1209,10 @@
          ! test energy conservation
          if(enertest)then
-	    call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot)
-	    call planetwide_sumval(capcal(:)*zdtsurfc(:)*area(:)/totarea_planet,dEtots)
-	    if (is_master) then
-	       print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
+            call planetwide_sumval(cpp*massarea(:,:)*zdtc(:,:)/totarea_planet,dEtot)
+            call planetwide_sumval(capcal(:)*zdtsurfc(:)*area(:)/totarea_planet,dEtots)
+            if (is_master) then
+               print*,'In co2cloud atmospheric energy change   =',dEtot,' W m-2'
                print*,'In co2cloud surface energy change       =',dEtots,' W m-2'
-	    endif
+            endif
          endif
 
@@ -1233,13 +1234,13 @@
             if(watercond.and.(RLVTT.gt.1.e-8))then
 
-	       dqmoist(1:ngrid,1:nlayer,1:nq)=0.
-	       dtmoist(1:ngrid,1:nlayer)=0.
+               dqmoist(1:ngrid,1:nlayer,1:nq)=0.
+               dtmoist(1:ngrid,1:nlayer)=0.
                
                ! Moist Convective Adjustment.
                ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-	       call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
+               call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
 
                pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap)     &
-	                                          + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
+                                                  + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)
                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice)     &
                                                   + dqmoist(1:ngrid,1:nlayer,igcm_h2o_ice)
@@ -1249,24 +1250,24 @@
                if(enertest)then
                
-		  call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
-		  call planetwide_maxval(dtmoist(:,:),dtmoist_max)
-		  call planetwide_minval(dtmoist(:,:),dtmoist_min)
-		  madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:)
+                  call planetwide_sumval(cpp*massarea(:,:)*dtmoist(:,:)/totarea_planet,dEtot)
+                  call planetwide_maxval(dtmoist(:,:),dtmoist_max)
+                  call planetwide_minval(dtmoist(:,:),dtmoist_min)
+                  madjdEz(:,:)=cpp*mass(:,:)*dtmoist(:,:)
                   do ig=1,ngrid
                      madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
                   enddo
                   
-		  if (is_master) then
+                  if (is_master) then
                      print*,'In moistadj atmospheric energy change   =',dEtot,' W m-2'
                      print*,'In moistadj MAX atmospheric energy change   =',dtmoist_max*ptimestep,'K/step'
                      print*,'In moistadj MIN atmospheric energy change   =',dtmoist_min*ptimestep,'K/step'
-		  endif
-		  
-		  call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+	&
-		  	                 massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
-		  if (is_master) print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
+                  endif
+                  
+                  call planetwide_sumval(massarea(:,:)*dqmoist(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+        &
+                                           massarea(:,:)*dqmoist(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
+                  if (is_master) print*,'In moistadj atmospheric water change    =',dWtot,' kg m-2'
                   
                endif ! end of 'enertest'
-	       
+               
 
                ! Large scale condensation/evaporation.
@@ -1284,13 +1285,13 @@
                      lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
                   enddo
-		  call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot)
-
-		  if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
+                  call planetwide_sumval(cpp*massarea(:,:)*dtlscale(:,:)/totarea_planet,dEtot)
+
+                  if (is_master) print*,'In largescale atmospheric energy change =',dEtot,' W m-2'
 
                   ! Test water conservation.
-		  call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+	&
-		  	                 SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
+                  call planetwide_sumval(massarea(:,:)*dqvaplscale(:,:)*ptimestep/totarea_planet+        &
+                                           SUM(massarea(:,:)*dqcldlscale(:,:))*ptimestep/totarea_planet,dWtot)
                         
-		  if (is_master) print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
+                  if (is_master) print*,'In largescale atmospheric water change  =',dWtot,' kg m-2'
                endif ! end of 'enertest'
 
@@ -1316,7 +1317,7 @@
 
                pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) &
-	                                          + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
+                                                  + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)
                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = pdq(1:ngrid,1:nlayer,igcm_h2o_ice) &
-	                                          + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
+                                                  + zdqrain(1:ngrid,1:nlayer,igcm_h2o_ice)
                pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtrain(1:ngrid,1:nlayer)
                
@@ -1327,30 +1328,30 @@
                if(enertest)then
                
-		  call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
-		  if (is_master) print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
-		  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp,dItot_tmp)
-		  call planetwide_sumval(area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp,dItot)
-		  dItot = dItot + dItot_tmp
-		  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dVtot_tmp)
-		  call planetwide_sumval(area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp,dVtot)
-		  dVtot = dVtot + dVtot_tmp
-		  dEtot = dItot + dVtot
+                  call planetwide_sumval(cpp*massarea(:,:)*zdtrain(:,:)/totarea_planet,dEtot)
+                  if (is_master) print*,'In rain atmospheric T energy change       =',dEtot,' W m-2'
+                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)/totarea_planet*RLVTT/cpp,dItot_tmp)
+                  call planetwide_sumval(area(:)*zdqssnow(:)/totarea_planet*RLVTT/cpp,dItot)
+                  dItot = dItot + dItot_tmp
+                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet,dVtot_tmp)
+                  call planetwide_sumval(area(:)*zdqsrain(:)/totarea_planet*RLVTT/cpp,dVtot)
+                  dVtot = dVtot + dVtot_tmp
+                  dEtot = dItot + dVtot
                   
-		  if (is_master) then
+                  if (is_master) then
                      print*,'In rain dItot =',dItot,' W m-2'
                      print*,'In rain dVtot =',dVtot,' W m-2'
                      print*,'In rain atmospheric L energy change       =',dEtot,' W m-2'
-		  endif
+                  endif
 
                   ! Test water conservation
-		  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+	&
-		  	massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
-		  call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*area(:)*ptimestep/totarea_planet,dWtots)
+                  call planetwide_sumval(massarea(:,:)*zdqrain(:,:,igcm_h2o_vap)*ptimestep/totarea_planet+        &
+                          massarea(:,:)*zdqrain(:,:,igcm_h2o_ice)*ptimestep/totarea_planet,dWtot)
+                  call planetwide_sumval((zdqsrain(:)+zdqssnow(:))*area(:)*ptimestep/totarea_planet,dWtots)
                   
-		  if (is_master) then
-                  	print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
-                  	print*,'In rain surface water change            =',dWtots,' kg m-2'
-                  	print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
-		  endif
+                  if (is_master) then
+                          print*,'In rain atmospheric water change        =',dWtot,' kg m-2'
+                          print*,'In rain surface water change            =',dWtots,' kg m-2'
+                          print*,'In rain non-cons factor                 =',dWtot+dWtots,' kg m-2'
+                  endif
                   
                endif ! end of 'enertest'
@@ -1373,10 +1374,10 @@
            
                iq=igcm_h2o_ice
-	       call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
-	       call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
-	       if (is_master) then
-	      	  print*,'Before sedim pq  =',dWtot,' kg m-2'
+               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
+               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
+               if (is_master) then
+                        print*,'Before sedim pq  =',dWtot,' kg m-2'
                   print*,'Before sedim pdq =',dWtots,' kg m-2'
-	       endif
+               endif
             endif
            
@@ -1387,10 +1388,10 @@
             if(watertest)then
                iq=igcm_h2o_ice
-	       call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
-	       call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
-	       if (is_master) then
-              	  print*,'After sedim pq  =',dWtot,' kg m-2'
-              	  print*,'After sedim pdq =',dWtots,' kg m-2'
-  	       endif
+               call planetwide_sumval(massarea(:,:)*pq(:,:,iq)*ptimestep/totarea_planet,dWtot)
+               call planetwide_sumval(massarea(:,:)*pdq(:,:,iq)*ptimestep/totarea_planet,dWtots)
+               if (is_master) then
+                        print*,'After sedim pq  =',dWtot,' kg m-2'
+                        print*,'After sedim pdq =',dWtots,' kg m-2'
+                 endif
             endif
 
@@ -1401,11 +1402,11 @@
             ! Test water conservation
             if(watertest)then
-	       call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot)
-	       call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:)*ptimestep/totarea_planet,dWtots)
-	       if (is_master) then
-              	  print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
-              	  print*,'In sedim surface ice change             =',dWtots,' kg m-2'
-              	  print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
-	       endif
+               call planetwide_sumval(massarea(:,:)*(zdqsed(:,:,igcm_h2o_vap)+zdqsed(:,:,igcm_h2o_ice))*ptimestep/totarea_planet,dWtot)
+               call planetwide_sumval((zdqssed(:,igcm_h2o_vap)+zdqssed(:,igcm_h2o_ice))*area(:)*ptimestep/totarea_planet,dWtots)
+               if (is_master) then
+                        print*,'In sedim atmospheric ice change         =',dWtot,' kg m-2'
+                        print*,'In sedim surface ice change             =',dWtots,' kg m-2'
+                        print*,'In sedim non-cons factor                =',dWtot+dWtots,' kg m-2'
+               endif
             endif
 
@@ -1421,22 +1422,22 @@
 
             zdmassmr(1:ngrid,1:nlayer) = mass(1:ngrid,1:nlayer) *     &
-	       (   zdqevap(1:ngrid,1:nlayer)                          &
-	         + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
-		 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)             &
-		 + dqvaplscale(1:ngrid,1:nlayer) )
+               (   zdqevap(1:ngrid,1:nlayer)                          &
+                 + zdqrain(1:ngrid,1:nlayer,igcm_h2o_vap)             &
+                 + dqmoist(1:ngrid,1:nlayer,igcm_h2o_vap)             &
+                 + dqvaplscale(1:ngrid,1:nlayer) )
 
             do ig = 1, ngrid
-	       zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
+               zdmassmr_col(ig)=SUM(zdmassmr(ig,1:nlayer))
             enddo
             
-	    call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
-	    call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
-	    call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
-
-	    call mass_redistribution(ngrid,nlayer,nq,ptimestep,                     &
+            call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
+            call writediagfi(ngrid,"mass_evap_col","mass gain col"," ",2,zdmassmr_col)
+            call writediagfi(ngrid,"mass","mass","kg/m2",3,mass)
+
+            call mass_redistribution(ngrid,nlayer,nq,ptimestep,                     &
                                      rnat,capcal,pplay,pplev,pt,tsurf,pq,qsurf,     &
-		                     pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
-		                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
-	 
+                                     pu,pv,pdt,zdtsurf,pdq,pdu,pdv,zdmassmr,        &
+                                     zdtmr,zdtsurfmr,zdpsrfmr,zdumr,zdvmr,zdqmr,zdqsurfmr)
+         
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmr(1:ngrid,1:nlayer,1:nq)
             dqsurf(1:ngrid,1:nq)       = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
@@ -1444,8 +1445,8 @@
             pdu(1:ngrid,1:nlayer)      = pdu(1:ngrid,1:nlayer) + zdumr(1:ngrid,1:nlayer)
             pdv(1:ngrid,1:nlayer)      = pdv(1:ngrid,1:nlayer) + zdvmr(1:ngrid,1:nlayer)
-	    pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
+            pdpsrf(1:ngrid)            = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
             zdtsurf(1:ngrid)           = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
-	    
-	 endif
+            
+         endif
 
   ! ------------------
@@ -1493,5 +1494,5 @@
             call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
                         capcal,albedo,albedo_bareground,                    &
-	                albedo_snow_SPECTV,albedo_co2_ice_SPECTV,           &
+                        albedo_snow_SPECTV,albedo_co2_ice_SPECTV,           &
                         mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,            &
                         sea_ice)
@@ -1504,17 +1505,17 @@
             ! Test energy conservation
             if(enertest)then
-	       call planetwide_sumval(area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots)
-	       if (is_master) print*,'In hydrol surface energy change     =',dEtots,' W m-2'
+               call planetwide_sumval(area(:)*capcal(:)*zdtsurf_hyd(:)/totarea_planet,dEtots)
+               if (is_master) print*,'In hydrol surface energy change     =',dEtots,' W m-2'
             endif
 
             ! test water conservation
             if(watertest)then
-	       call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
-	       if (is_master) print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
-	          call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots)
-	       if (is_master) then
+               call planetwide_sumval(dqs_hyd(:,igcm_h2o_ice)*area(:)*ptimestep/totarea_planet,dWtots)
+               if (is_master) print*,'In hydrol surface ice change            =',dWtots,' kg m-2'
+                  call planetwide_sumval(dqs_hyd(:,igcm_h2o_vap)*area(:)*ptimestep/totarea_planet,dWtots)
+               if (is_master) then
                   print*,'In hydrol surface water change          =',dWtots,' kg m-2'
                   print*,'---------------------------------------------------------------'
-	       endif
+               endif
             endif
 
@@ -1590,6 +1591,6 @@
       ! Test energy conservation
       if(enertest)then
-	 call planetwide_sumval(area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)	 
-	 if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
+         call planetwide_sumval(area(:)*capcal(:)*zdtsurf(:)/totarea_planet,dEtots)         
+         if (is_master) print*,'Surface energy change                 =',dEtots,' W m-2'
       endif
 
@@ -1638,8 +1639,8 @@
          print*,Ts1,Ts2,Ts3,TsS
       else
-      	 if (is_master) then
+               if (is_master) then
             print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
             print*,Ts1,Ts2,Ts3
-	 endif
+         endif
       end if
 
@@ -1648,9 +1649,9 @@
       if(corrk)then
 
-	 call planetwide_sumval(area(:)*fluxtop_dn(:)/totarea_planet,ISR)
-	 call planetwide_sumval(area(:)*fluxabs_sw(:)/totarea_planet,ASR)
-	 call planetwide_sumval(area(:)*fluxtop_lw(:)/totarea_planet,OLR)
-	 call planetwide_sumval(area(:)*fluxgrd(:)/totarea_planet,GND)
-	 call planetwide_sumval(area(:)*fluxdyn(:)/totarea_planet,DYN)
+         call planetwide_sumval(area(:)*fluxtop_dn(:)/totarea_planet,ISR)
+         call planetwide_sumval(area(:)*fluxabs_sw(:)/totarea_planet,ASR)
+         call planetwide_sumval(area(:)*fluxtop_lw(:)/totarea_planet,OLR)
+         call planetwide_sumval(area(:)*fluxgrd(:)/totarea_planet,GND)
+         call planetwide_sumval(area(:)*fluxdyn(:)/totarea_planet,DYN)
          do ig=1,ngrid
             if(fluxtop_dn(ig).lt.0.0)then
@@ -1668,9 +1669,9 @@
             DYN=0.0
          endif
-	 
-	 if (is_master) then
+         
+         if (is_master) then
             print*,'  ISR            ASR            OLR            GND            DYN [W m^-2]'
             print*, ISR,ASR,OLR,GND,DYN
-	 endif
+         endif
 
          if(enertest .and. is_master)then
@@ -1688,8 +1689,8 @@
                open(93,file="tem_bal.out",form='formatted',position='append')
                if(callsoil)then
-               	  write(93,*) zday,Ts1,Ts2,Ts3,TsS
-	       else
-	          write(93,*) zday,Ts1,Ts2,Ts3
-	       endif
+                         write(93,*) zday,Ts1,Ts2,Ts3,TsS
+               else
+                  write(93,*) zday,Ts1,Ts2,Ts3
+               endif
                close(93)
             endif
@@ -1744,16 +1745,16 @@
       if(water)then
 
-	 call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf)
-	 call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf)
-	 call planetwide_sumval(area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol)
-	 call planetwide_sumval(area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol)
+         call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_ice)/totarea_planet,icesrf)
+         call planetwide_sumval(area(:)*qsurf_hist(:,igcm_h2o_vap)/totarea_planet,liqsrf)
+         call planetwide_sumval(area(:)*qcol(:,igcm_h2o_ice)/totarea_planet,icecol)
+         call planetwide_sumval(area(:)*qcol(:,igcm_h2o_vap)/totarea_planet,vapcol)
 
          h2otot = icesrf + liqsrf + icecol + vapcol
-	 
-	 if (is_master) then
+         
+         if (is_master) then
             print*,' Total water amount [kg m^-2]: ',h2otot
             print*,' Surface ice    Surface liq.   Atmos. con.     Atmos. vap. [kg m^-2] '
             print*, icesrf,liqsrf,icecol,vapcol
-	 endif
+         endif
 
          if(meanOLR .and. is_master)then
@@ -2013,32 +2014,32 @@
          
       endif ! end of 'callrad'
-	
+        
       if(enertest) then
       
-	 if (calldifv) then
-         
-	    call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
+         if (calldifv) then
+         
+            call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
             call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
             
-!	     call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
-!	     call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
-!	     call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
-	     
-	 endif
+!             call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
+!             call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
+!             call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
+             
+         endif
           
-	 if (corrk) then
-	    call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
-	    call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
-	 endif
+         if (corrk) then
+            call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
+            call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
+         endif
           
          if(watercond) then
 
-	    call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
-	    call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
-	    call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
+            call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
+            call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
+            call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
              
-!	     call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 
-!	     call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)             
-!	     call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
+!             call writediagfi(ngrid,"lscaledEz","heat from largescale","W m-2",3,lscaledEz) 
+!             call writediagfi(ngrid,"madjdEz","heat from moistadj","W m-2",3,madjdEz)             
+!             call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
 
          endif
@@ -2060,11 +2061,11 @@
       ! Output aerosols.
       if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-	 call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
+         call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
       if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-	 call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
+         call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
       if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
-	 call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
+         call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
       if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
-	 call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
+         call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
 
       ! Output tracers.
@@ -2086,5 +2087,5 @@
                call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
                call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
-	       call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
+               call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
             endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/soil.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1524)
@@ -5,5 +5,6 @@
 
       use comsoil_h, only: layer, mlayer, volcapa, inertiedat
-      use comcstfi_mod, only: pi, daysec
+      use comcstfi_mod, only: pi
+      use time_phylmdz_mod, only: daysec
       use planete_mod, only: year_day
       use comgeomfi_h, only: long, lati
Index: trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 1524)
@@ -53,6 +53,6 @@
       use planete_mod, only: year_day, periastr, apoastr, peri_day,
      &                       obliquit, z0, lmixmin, emin_turb
-      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, dtphys,
-     &                        daysec, r
+      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, cpp, r
+      use time_phylmdz_mod, only: dtphys, daysec
       use callkeys_mod, only: check_cpp_match,force_cpp
       implicit none
@@ -507,5 +507,5 @@
       write(*,*) 
 
-      ENDIF 			!	of if (Lmodif == 1)
+      ENDIF ! of if (Lmodif == 1)
 
 c-----------------------------------------------------------------------
Index: trunk/LMDZ.GENERIC/libf/phystd/time_phylmdz_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/time_phylmdz_mod.F90	(revision 1524)
+++ trunk/LMDZ.GENERIC/libf/phystd/time_phylmdz_mod.F90	(revision 1524)
@@ -0,0 +1,25 @@
+MODULE time_phylmdz_mod
+
+    IMPLICIT NONE
+    REAL,SAVE    :: dtphys     ! physics time step (s)
+!$OMP THREADPRIVATE(dtphys)
+    REAL,SAVE    :: daysec     ! length of day (s)
+!$OMP THREADPRIVATE(daysec)
+    INTEGER,SAVE :: day_ini     ! initial day of the run
+!$OMP THREADPRIVATE(day_ini)
+
+CONTAINS
+
+  SUBROUTINE init_time(day_ini_, daysec_, dtphys_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: day_ini_
+    REAL,INTENT(IN) :: daysec_
+    REAL,INTENT(IN) :: dtphys_
+    
+    day_ini=day_ini_
+    daysec=daysec_
+    dtphys=dtphys_
+
+  END SUBROUTINE init_time
+
+END MODULE time_phylmdz_mod      
Index: trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F	(revision 1524)
@@ -44,5 +44,5 @@
      &                               is_master, gather
       USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-      USE temps_mod, ONLY: day_ini
+      USE time_phylmdz_mod, ONLY: day_ini
       implicit none
 
Index: trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 1524)
@@ -50,5 +50,5 @@
       use control_mod, only: ecritphy, iphysiq, day_step
       use callkeys_mod, only: iradia
-      USE temps_mod, ONLY: day_ini
+      use time_phylmdz_mod, ONLY: day_ini
 
       implicit none
Index: trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 1524)
@@ -50,5 +50,5 @@
       use control_mod, only: ecritphy, iphysiq, day_step
       use callkeys_mod, only: iradia
-      USE temps_mod, ONLY: day_ini
+      use time_phylmdz_mod, ONLY: day_ini
 
       implicit none
Index: trunk/LMDZ.GENERIC/libf/phystd/writeg1d.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/writeg1d.F	(revision 1523)
+++ trunk/LMDZ.GENERIC/libf/phystd/writeg1d.F	(revision 1524)
@@ -129,5 +129,5 @@
 
       SUBROUTINE endg1d(ngrid,nlayer,zlayer,ndt)
-      USE comcstfi_mod, ONLY: dtphys, daysec
+      USE time_phylmdz_mod, ONLY: dtphys, daysec
       USE comg1d_mod, ONLY: g1d_nomfich,g1d_unitfich,g1d_nvar,
      &	g1d_nomvar,g1d_titrevar,g1d_dimvar,g1d_nlayer,g1d_unitctl,
Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 1523)
+++ trunk/LMDZ.MARS/README	(revision 1524)
@@ -2225,2 +2225,9 @@
 - Turned iniphysiq (in dynlonlat_phylonlat/phymars)into module
   "iniphysiq_mod.F90".
+
+== 29/03/2016 == EM
+- Added "time_phylmdz_mod.F90" module to store information on time and
+  calendar in the physics, should be used instead of accessing "temps_mod"
+  from the dynamics. And moved daysec and dtphys from module "comcstfi_h"
+  to module "time_phylmdz_mod".
+- Made "phys_state_var_init" a module
Index: trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq_mod.F90	(revision 1523)
+++ trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/iniphysiq_mod.F90	(revision 1524)
@@ -23,4 +23,6 @@
 use infotrac, only : nqtot ! number of advected tracers
 use comgeomfi_h, only: ini_fillgeom
+use temps_mod, only: day_ini, hour_ini
+use phys_state_var_init_mod, only: phys_state_var_init
 use regular_lonlat_mod, only: init_regular_lonlat, &
                               east, west, north, south, &
@@ -171,7 +173,11 @@
 ! and do some initializations 
 call phys_state_var_init(klon_omp,nlayer,nqtot, &
-                         punjours,ptimestep,prad,pg,pr,pcpp)
+                         day_ini,hour_ini,punjours,ptimestep, &
+                         prad,pg,pr,pcpp)
 call ini_fillgeom(klon_omp,rlatd,rlond,airephy)
 call conf_phys(klon_omp,nlayer,nqtot)
+
+! Initialize some "temporal and calendar" related variables
+!CALL init_time(day_ini,hour_ini,punjours,ptimestep)
 
 !$OMP END PARALLEL
Index: trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/newstart.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/newstart.F	(revision 1524)
@@ -41,4 +41,5 @@
       USE temps_mod, ONLY: day_ini,hour_ini
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE phys_state_var_init_mod, ONLY: phys_state_var_init
 
       implicit none
@@ -373,5 +374,6 @@
 !      'datafile' path may be changed by user)
       call phys_state_var_init(ngridmx,llm,nqtot,
-     .                         daysec,dtphys,rad,g,r,cpp)
+     &                         day_ini,hour_ini,daysec,dtphys,
+     &                         rad,g,r,cpp)
       call ini_fillgeom(ngridmx,latfi,lonfi,airefi)
       call conf_phys(ngridmx,llm,nqtot)
Index: trunk/LMDZ.MARS/libf/phymars/callradite.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/callradite.F	(revision 1524)
@@ -10,5 +10,6 @@
      &            iaer_dust_submicron,iaer_h2o_ice
       use yomlw_h, only: gcp, nlaylte
-      USE comcstfi_h
+      use comcstfi_h, only: g,cpp
+      use time_phylmdz_mod, only: daysec
       IMPLICIT NONE
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/phymars/co2snow.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2snow.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/co2snow.F	(revision 1524)
@@ -4,5 +4,5 @@
       use surfdat_h, only: iceradius, dtemisice
       use comgeomfi_h, only: lati ! grid point latitudes (rad)
-      USE comcstfi_h
+      use time_phylmdz_mod, only: daysec
       IMPLICIT NONE
 
Index: trunk/LMDZ.MARS/libf/phymars/comcstfi_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/comcstfi_h.F90	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/comcstfi_h.F90	(revision 1524)
@@ -2,14 +2,12 @@
       IMPLICIT NONE
 
-      REAL pi
-      REAL rad
-      REAL g
-      REAL r
-      REAL cpp
-      REAL rcp
-      REAL dtphys
-      REAL daysec
-      REAL mugaz
-      REAL omeg
+      REAL pi ! something like 3.14159
+      REAL rad ! radius of the planet (m)
+      REAL g ! gravity (m/s2)
+      REAL r ! reduced gas constant (r=8.314511/(mugaz/1000.0))
+      REAL cpp ! Cp of the atmosphere
+      REAL rcp ! r/cpp
+      REAL mugaz ! molar mass of the atmosphere (g/mol)
+      REAL omeg ! planet rotation rate (rad/s)
 
       END MODULE comcstfi_h
Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1524)
@@ -40,5 +40,5 @@
       use control_mod, only: ecritphy
       use planete_h
-      USE comcstfi_h, only: daysec,dtphys
+      use time_phylmdz_mod, only: daysec,dtphys
       use dimradmars_mod, only: naerkind, name_iaer,
      &                      ini_scatterers,tauvis
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/temps_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/temps_mod.F90	(revision 1523)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/temps_mod.F90
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1524)
@@ -16,9 +16,10 @@
       use planete_h, only: year_day, periheli, aphelie, peri_day,
      &                     obliquit, emin_turb, lmixmin
-      use comcstfi_h, only: pi, rad, daysec, omeg, g, mugaz, rcp, r,
-     &                      cpp, dtphys
+      use comcstfi_h, only: pi, rad, omeg, g, mugaz, rcp, r, cpp
+      use time_phylmdz_mod, only: daysec, dtphys
       use dimradmars_mod, only: tauscaling,tauvis
       USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig
       USE logic_mod, ONLY: hybrid
+      USE phys_state_var_init_mod, ONLY: phys_state_var_init
       IMPLICIT NONE
 
@@ -478,5 +479,5 @@
 !Mars possible matter with dtphys in input and include!!!
       call phys_state_var_init(1,llm,nq,
-     .          daysec,dtphys,rad,g,r,cpp)
+     .          day0,time,daysec,dtphys,rad,g,r,cpp)
       call ini_fillgeom(1,latitude,longitude,(/1.0/))
       call conf_phys(1,llm,nq)
Index: trunk/LMDZ.MARS/libf/phymars/eofdump_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/eofdump_mod.F90	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/eofdump_mod.F90	(revision 1524)
@@ -79,5 +79,6 @@
 
       use comgeomfi_h, only: long, lati
-      use comcstfi_h
+      use comcstfi_h, only: pi
+      use time_phylmdz_mod, only: daysec, dtphys
       USE comvert_mod, ONLY: aps,bps
       implicit none
Index: trunk/LMDZ.MARS/libf/phymars/inistats.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/inistats.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/inistats.F	(revision 1524)
@@ -3,12 +3,13 @@
       use mod_phys_lmdz_para, only : is_master
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
-      USE comcstfi_h, ONLY: daysec,dtphys,pi
+      USE comcstfi_h, ONLY: pi
+      USE time_phylmdz_mod, ONLY: daysec,dtphys
       implicit none
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "statto.h"
-#include "netcdf.inc"
+      include "dimensions.h"
+      include "paramet.h"
+      include "comgeom.h"
+      include "statto.h"
+      include "netcdf.inc"
 
       integer,intent(out) :: ierr
Index: trunk/LMDZ.MARS/libf/phymars/iniwrite.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/iniwrite.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/iniwrite.F	(revision 1524)
@@ -2,9 +2,9 @@
 
       use comsoil_h, only: mlayer, nsoilmx
-      USE comcstfi_h
+      USE comcstfi_h, only: g, mugaz, omeg, rad, rcp, pi 
       USE comvert_mod, ONLY: ap,bp,aps,bps,pseudoalt
       USE logic_mod, ONLY: fxyhypb,ysinus
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
-      USE temps_mod, ONLY: hour_ini
+      USE time_phylmdz_mod, ONLY: hour_ini, daysec, dtphys
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
       IMPLICIT NONE
Index: trunk/LMDZ.MARS/libf/phymars/lwi.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/lwi.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/lwi.F	(revision 1524)
@@ -6,5 +6,6 @@
       use dimradmars_mod, only: ndlo2, ndlon, nflev, nir
       use yomlw_h, only: gcp, nlaylte, xi
-      USE comcstfi_h
+      USE comcstfi_h, ONLY: g, cpp
+      USE time_phylmdz_mod, ONLY: dtphys
       implicit none
 
Index: trunk/LMDZ.MARS/libf/phymars/nirco2abs.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nirco2abs.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/nirco2abs.F	(revision 1524)
@@ -4,5 +4,6 @@
        use tracer_mod, only: igcm_co2, igcm_o
        use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
-       USE comcstfi_h
+       USE comcstfi_h, ONLY: pi
+       USE time_phylmdz_mod, ONLY: daysec
        IMPLICIT NONE
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/phymars/nlte_tcool.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/nlte_tcool.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/nlte_tcool.F	(revision 1524)
@@ -1357,5 +1357,5 @@
       function hrkday_convert                        
      @     ( mmean_nlte,cpmean_nlte )         
-      USE comcstfi_h
+      use time_phylmdz_mod, only: daysec
       use param_v4_h, only: n_avog
       implicit none                           
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1524)
@@ -22,6 +22,8 @@
                       put_var, put_field, length
   use mod_grid_phy_lmdz, only : klon_glo
-  use planete_h
-  use comcstfi_h
+  use planete_h, only: aphelie, emin_turb, lmixmin, obliquit, &
+                       peri_day, periheli, year_day
+  use comcstfi_h, only: g, mugaz, omeg, rad, rcp
+  use time_phylmdz_mod, only: daysec
 
   implicit none
Index: trunk/LMDZ.MARS/libf/phymars/phys_state_var_init.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phys_state_var_init.F	(revision 1523)
+++ 	(revision )
@@ -1,90 +1,0 @@
-      SUBROUTINE phys_state_var_init(ngrid,nlayer,nq
-     .               ,pdaysec,ptimestep,prad,pg,pr,pcpp)
-
-!=======================================================================
-!
-!   purpose:
-!   -------
-!
-!   Allocate arrays in modules
-!   Fill geometrical arrays
-!   Fill a first set of physical constants
-!   -- 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
-      use comsoil_h, only: ini_comsoil_h
-      use dimradmars_mod, only: ini_dimradmars_mod
-      use yomlw_h, only: ini_yomlw_h
-      use conc_mod, only: ini_conc_mod
-      use turb_mod, only: ini_turb_mod
-      use comcstfi_h, only: pi,rad,cpp,g,r,rcp,daysec,dtphys
-      use tracer_mod, only: nqmx
-
-      IMPLICIT NONE
-      
-      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
-      REAL,INTENT(IN) :: pdaysec,ptimestep,prad,pg,pr,pcpp
-      EXTERNAL SSUM
-      REAL SSUM
-
-      ! set dimension in tracer_mod
-      nqmx=nq
-
-      ! set parameters in comcstfi_h
-      pi=2.*asin(1.) 
-      rad=prad
-      cpp=pcpp
-      g=pg
-      r=pr
-      rcp=r/cpp
-      daysec=pdaysec
-      dtphys=ptimestep
-
-      ! 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)
-
-      ! allocate "comsoil_h" arrays
-      call ini_comsoil_h(ngrid)
-
-      ! set some variables in "dimradmars_mod"
-      call ini_dimradmars_mod(ngrid,nlayer)
-
-      ! 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/phys_state_var_init_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 1524)
+++ trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 1524)
@@ -0,0 +1,99 @@
+MODULE phys_state_var_init_mod
+
+CONTAINS
+
+      SUBROUTINE phys_state_var_init(ngrid,nlayer,nq, &
+                                     day_ini,hour_ini,pdaysec,ptimestep, &
+                                     prad,pg,pr,pcpp)
+
+!=======================================================================
+!
+!   purpose:
+!   -------
+!
+!   Allocate arrays in modules
+!   Fill geometrical arrays
+!   Fill a first set of physical constants
+!   -- 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
+      use comsoil_h, only: ini_comsoil_h
+      use dimradmars_mod, only: ini_dimradmars_mod
+      use yomlw_h, only: ini_yomlw_h
+      use conc_mod, only: ini_conc_mod
+      use turb_mod, only: ini_turb_mod
+      use comcstfi_h, only: pi,rad,cpp,g,r,rcp
+      use tracer_mod, only: nqmx
+      use time_phylmdz_mod, only: init_time
+
+      IMPLICIT NONE
+      
+      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
+      INTEGER,INTENT(IN) :: day_ini
+      REAL,INTENT(IN) :: hour_ini
+      REAL,INTENT(IN) :: pdaysec,ptimestep,prad,pg,pr,pcpp
+
+      ! set dimension in tracer_mod
+      nqmx=nq
+
+      ! set parameters in comcstfi_h
+      pi=2.*asin(1.) 
+      rad=prad
+      cpp=pcpp
+      g=pg
+      r=pr
+      rcp=r/cpp
+
+      ! Initialize some "temporal and calendar" related variables
+      call init_time(day_ini,hour_ini,pdaysec,ptimestep)
+
+      ! 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)
+
+      ! allocate "comsoil_h" arrays
+      call ini_comsoil_h(ngrid)
+
+      ! set some variables in "dimradmars_mod"
+      call ini_dimradmars_mod(ngrid,nlayer)
+
+      ! 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
+
+END MODULE phys_state_var_init_mod
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1524)
@@ -33,5 +33,6 @@
       use planete_h, only: aphelie, periheli, year_day, peri_day,
      &                     obliquit
-      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad, daysec
+      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
+      use time_phylmdz_mod, only: daysec
       use param_v4_h, only: nreact,n_avog, 
      &                      fill_data_thermos, allocate_param_thermos
Index: trunk/LMDZ.MARS/libf/phymars/tabfi.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tabfi.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/tabfi.F	(revision 1524)
@@ -50,6 +50,8 @@
       use iostart, only: get_var
       use mod_phys_lmdz_para, only: is_parallel
-      use comcstfi_h
-      use planete_h
+      use comcstfi_h, only: g, mugaz, omeg, rad, rcp
+      use time_phylmdz_mod, only: daysec, dtphys
+      use planete_h, only: aphelie, emin_turb, lmixmin, obliquit,
+     &                     peri_day, periheli, year_day
       implicit none
  
@@ -541,5 +543,5 @@
       write(*,*) 
 
-      ENDIF 			!	of if (Lmodif == 1)
+      ENDIF ! of if (Lmodif == 1)
 
 c-----------------------------------------------------------------------
Index: trunk/LMDZ.MARS/libf/phymars/time_phylmdz_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/time_phylmdz_mod.F90	(revision 1524)
+++ trunk/LMDZ.MARS/libf/phymars/time_phylmdz_mod.F90	(revision 1524)
@@ -0,0 +1,30 @@
+MODULE time_phylmdz_mod
+
+    IMPLICIT NONE
+    REAL,SAVE    :: dtphys      ! physics time step (s)
+!$OMP THREADPRIVATE(pdtphys)
+    REAL,SAVE    :: daysec     ! length of day (s)
+!$OMP THREADPRIVATE(daysec)
+    INTEGER,SAVE :: day_ini     ! initial day of the run
+!$OMP THREADPRIVATE(day_ini)
+    REAL,SAVE :: hour_ini       ! start time (fraction of day) of the run
+                                ! 0=<hour_ini<1
+!$OMP THREADPRIVATE(hour_ini)
+
+CONTAINS
+
+  SUBROUTINE init_time(day_ini_, hour_ini_, daysec_, dtphys_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: day_ini_
+    REAL,INTENT(IN) :: hour_ini_
+    REAL,INTENT(IN) :: daysec_
+    REAL,INTENT(IN) :: dtphys_
+    
+    day_ini=day_ini_
+    hour_ini=hour_ini_
+    daysec=daysec_
+    dtphys=dtphys_
+
+  END SUBROUTINE init_time
+
+END MODULE time_phylmdz_mod      
Index: trunk/LMDZ.MARS/libf/phymars/writediagfi.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 1524)
@@ -44,5 +44,5 @@
      &                               is_master, gather
       USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-      USE temps_mod, ONLY: day_ini
+      USE time_phylmdz_mod, ONLY: day_ini
       implicit none
 
Index: trunk/LMDZ.MARS/libf/phymars/writeg1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/writeg1d.F	(revision 1523)
+++ trunk/LMDZ.MARS/libf/phymars/writeg1d.F	(revision 1524)
@@ -124,5 +124,5 @@
 
       SUBROUTINE endg1d(ngrid,nlayer,zlayer,ndt)
-      USE comcstfi_h
+      USE time_phylmdz_mod, ONLY: dtphys, daysec
       IMPLICIT NONE
 c.......................................................................
Index: trunk/LMDZ.TITAN/libf/phy_common
===================================================================
--- trunk/LMDZ.TITAN/libf/phy_common	(revision 1524)
+++ trunk/LMDZ.TITAN/libf/phy_common	(revision 1524)
@@ -0,0 +1,1 @@
+link ../../LMDZ.COMMON/libf/phy_common
Index: trunk/LMDZ.TITAN/libf/phytitan/dyn1d/temps_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/dyn1d/temps_mod.F90	(revision 1523)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/temps_mod.F90
Index: trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F90	(revision 1523)
+++ trunk/LMDZ.TITAN/libf/phytitan/phyetat0.F90	(revision 1524)
@@ -15,5 +15,5 @@
       USE comgeomphy,  only: rlatd,rlond
       USE control_mod, only: raz_date
-      USE temps_mod, only: itau_phy
+      USE time_phylmdz_mod, only: itau_phy
 
 implicit none
Index: trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90	(revision 1523)
+++ trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90	(revision 1524)
@@ -14,5 +14,5 @@
       USE comgeomphy,  only: rlatd,rlond
       USE control_mod, only: raz_date
-      USE temps_mod, only: day_end, annee_ref, itau_phy
+      USE time_phylmdz_mod, only: day_end, annee_ref, itau_phy
 
       implicit none
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 1523)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 1524)
@@ -72,5 +72,5 @@
       USE moyzon_mod
       USE write_field_phy
-      USE temps_mod, only: itau_phy,day_ref,annee_ref
+      USE time_phylmdz_mod, only: itau_phy,day_ref,annee_ref
       USE logic_mod, only: iflag_trac,moyzon_ch,moyzon_mu
       IMPLICIT none
Index: trunk/LMDZ.TITAN/libf/phytitan/time_phylmdz_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/time_phylmdz_mod.F90	(revision 1524)
+++ trunk/LMDZ.TITAN/libf/phytitan/time_phylmdz_mod.F90	(revision 1524)
@@ -0,0 +1,35 @@
+MODULE time_phylmdz_mod
+
+    IMPLICIT NONE
+    REAL,SAVE    :: pdtphys     ! physics time step (s)
+!$OMP THREADPRIVATE(pdtphys)
+    INTEGER,SAVE :: annee_ref   ! reference year from the origin
+!$OMP THREADPRIVATE(annee_ref)
+    INTEGER,SAVE :: day_ref     ! reference day of the origin
+!$OMP THREADPRIVATE(day_ref)
+    INTEGER,SAVE :: day_ini     ! initial day of the run since first day of annee_ref
+!$OMP THREADPRIVATE(day_ini)
+    INTEGER,SAVE :: day_end     ! final day of the run since first day of annee_ref
+!$OMP THREADPRIVATE(day_end)
+    INTEGER,SAVE :: itau_phy     ! number of physics iterations
+!$OMP THREADPRIVATE(itau_phy)
+
+CONTAINS
+
+  SUBROUTINE init_time(annee_ref_, day_ref_, day_ini_, day_end_, pdtphys_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: annee_ref_
+    INTEGER,INTENT(IN) :: day_ref_
+    INTEGER,INTENT(IN) :: day_ini_
+    INTEGER,INTENT(IN) :: day_end_
+    REAL,INTENT(IN) :: pdtphys_
+    
+    annee_ref=annee_ref_
+    day_ref=day_ref_
+    day_ini=day_ini_
+    day_end=day_end_
+    pdtphys=pdtphys_
+
+  END SUBROUTINE init_time
+
+END MODULE time_phylmdz_mod      
Index: trunk/LMDZ.VENUS/libf/phy_common
===================================================================
--- trunk/LMDZ.VENUS/libf/phy_common	(revision 1524)
+++ trunk/LMDZ.VENUS/libf/phy_common	(revision 1524)
@@ -0,0 +1,1 @@
+link ../../LMDZ.COMMON/libf/phy_common
Index: trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/temps_mod.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/dyn1d/temps_mod.F90	(revision 1523)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/temps_mod.F90
Index: trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90	(revision 1523)
+++ trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90	(revision 1524)
@@ -15,5 +15,5 @@
       USE comgeomphy,  only: rlatd,rlond
       USE control_mod, only: raz_date
-      USE temps_mod, only: itau_phy
+      USE time_phylmdz_mod, only: itau_phy
 
 implicit none
Index: trunk/LMDZ.VENUS/libf/phyvenus/phyredem.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phyredem.F90	(revision 1523)
+++ trunk/LMDZ.VENUS/libf/phyvenus/phyredem.F90	(revision 1524)
@@ -14,5 +14,5 @@
       USE comgeomphy,  only: rlatd,rlond
       USE control_mod, only: raz_date
-      USE temps_mod, only: day_end, annee_ref, itau_phy
+      USE time_phylmdz_mod, only: day_end, annee_ref, itau_phy
 
       implicit none
Index: trunk/LMDZ.VENUS/libf/phyvenus/physiq.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1523)
+++ trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1524)
@@ -70,5 +70,5 @@
       use moyzon_mod, only: tmoy
 !      use ieee_arithmetic
-      use temps_mod, only: itau_phy,day_ref,annee_ref
+      use time_phylmdz_mod, only: annee_ref, day_ref, itau_phy
       use logic_mod, only: iflag_trac
       IMPLICIT none
Index: trunk/LMDZ.VENUS/libf/phyvenus/time_phylmdz_mod.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/time_phylmdz_mod.F90	(revision 1524)
+++ trunk/LMDZ.VENUS/libf/phyvenus/time_phylmdz_mod.F90	(revision 1524)
@@ -0,0 +1,35 @@
+MODULE time_phylmdz_mod
+
+    IMPLICIT NONE
+    REAL,SAVE    :: pdtphys     ! physics time step (s)
+!$OMP THREADPRIVATE(pdtphys)
+    INTEGER,SAVE :: annee_ref   ! reference year from the origin
+!$OMP THREADPRIVATE(annee_ref)
+    INTEGER,SAVE :: day_ref     ! reference day of the origin
+!$OMP THREADPRIVATE(day_ref)
+    INTEGER,SAVE :: day_ini     ! initial day of the run since first day of annee_ref
+!$OMP THREADPRIVATE(day_ini)
+    INTEGER,SAVE :: day_end     ! final day of the run since first day of annee_ref
+!$OMP THREADPRIVATE(day_end)
+    INTEGER,SAVE :: itau_phy     ! number of physics iterations
+!$OMP THREADPRIVATE(itau_phy)
+
+CONTAINS
+
+  SUBROUTINE init_time(annee_ref_, day_ref_, day_ini_, day_end_, pdtphys_)
+    IMPLICIT NONE
+    INTEGER,INTENT(IN) :: annee_ref_
+    INTEGER,INTENT(IN) :: day_ref_
+    INTEGER,INTENT(IN) :: day_ini_
+    INTEGER,INTENT(IN) :: day_end_
+    REAL,INTENT(IN) :: pdtphys_
+    
+    annee_ref=annee_ref_
+    day_ref=day_ref_
+    day_ini=day_ini_
+    day_end=day_end_
+    pdtphys=pdtphys_
+
+  END SUBROUTINE init_time
+
+END MODULE time_phylmdz_mod      
