Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM.newphys
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM.newphys	(revision 1233)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM.newphys	(revision 1234)
@@ -82,34 +82,34 @@
 ####
 ############ name in WRF !!!!! ############# name in NETCDF !!!! ############################################## name in meso_physiq !!!!!
-state  real  TSURF     ij   misc  1  -  rhd  "TSURF"     "SURF TEMPERATURE"                "K"       #SAVEMARS2 tsurf
-state  real  CO2ICE    ij   misc  1  -  rd   "CO2ICE"    "CO2 GROUND ICE"                  "kg m-2"  #SAVEMARS2 co2ice
-state  real  EMISSIV   ij   misc  1  -  rd   "EMISSIV"   "EMISSIVITY"                      ""        #SAVEMARS2 emis
-state  real  THETA     ij   misc  1  -  rd   "THETA"     "SLOPE INCLINATION"               "deg"     #SAVEMARS2 theta_sl
-state  real  PSI       ij   misc  1  -  rd   "PSI"       "SLOPE ORIENTATION"               "deg"     #SAVEMARS2 psi_sl
-state  real  TAU_DUST  ij   misc  1  -  rhd   "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""        #SAVEMARS2 tauref
-state  real  SWDOWNZ   ij   misc  1  -  rhd   "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"   #SAVEMARS2 fluxsurf_sw_tot      
-state  real  LWDOWNZ   ij   misc  1  -  rhd   "LWDOWNZ"   "DOWNWARD LW FLUX AT SURFACE"     "W m-2"   #SAVEMARS2 fluxsurf_lw
-state  real  SWUP      ij   misc  1  -  rd   "SWUP"      "UPWARD SW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_sw_tot    
-state  real  LWUP      ij   misc  1  -  rd   "LWUP"      "UPWARD LW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_lw
-state  real  MTOT      ij   misc  1  -  rhd   "MTOT"      "TOTAL MASS WATER VAPOR in pmic"  "pmic"    #SAVEMARS2 mtot
-state  real  ICETOT    ij   misc  1  -  rhd   "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"  #SAVEMARS2 icetot
-state  real  RAVE      ij   misc  1  -  rhd   "RAVE"      "MEAN ICE RADIUS"                 "m"       #SAVEMARS2 rave
-state  real  RICE      ikj  misc  1  -  rd   "RICE"      "ICE RADIUS"                      "m"       #SAVEMARS3 rice
-state  real  HR_SW     ikj  misc  1  -  rd   "HR_SW"     "HEATING RATE SW"                 "K/s"     #SAVEMARS3 zdtsw
-state  real  HR_LW     ikj  misc  1  -  rd   "HR_LW"     "HEATING RATE LW"                 "K/s"     #SAVEMARS3 zdtlw
-state  real  HR_SH     ikj  misc  1  -  rd   "HR_SH"     "HEATING RATE sens. heat"         "K/s"     #SAVEMARS3 zdtdif
-state  real  QSURFICE  ij   misc  1  -  rhd   "QSURFICE"  "WATER ICE AT SURFACE"            "kg m-2"  #SAVEMARS2 qsurfice
-state  real  QSURFDUST ij   misc  1  -  rd   "QSURFDUST" "DUST MASS AT SURFACE"            "kg m-2"  #SAVEMARS2 qsurfdust
-state  real  RDUST     ikj  misc  1  -  rd   "RDUST"     "DUST RADIUS"                     "m"       #SAVEMARS3 rdust
-state  real  HR_NIR    ikj  misc  1  -  rd   "HR_NIR"    "HEATING RATE nirco2"             "K/s"     #SAVEMARS3 zdtnirco2
-state  real  HR_NLTE   ikj  misc  1  -  rd   "HR_NLTE"   "HEATING RATE nlte"               "K/s"     #SAVEMARS3 zdtnlte
-state  real  ALBBARE   ij   misc  1  -  rhd   "ALBBARE"   "SOIL ALBEDO"                     ""        #SAVEMARS2 albedodat
-state  real  VMR_ICE   ikj  misc  1  -  rhd   "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"     #SAVEMARS3 vmr
-state  real  TAU_ICE   ij   misc  1  -  rhd   "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""        #SAVEMARS2 tauTES
-state  real  PDTZ      ikj  misc  1  -  rd   "PDT"       "TEMP TENDENCY"                   "K s-1"   #SAVEMARS3 pdt
-state  real  ZMAX_TH   ij   misc  1  -  rhd   "ZMAX_TH"   "MAXIMUM LEVEL REACHED IN TH"     "m"       #SAVEMARS2 zmax_th
-state  real  HFMAX_TH  ij   misc  1  -  rhd   "HFMAX_TH"  "MAXIMUM TH HEAT FLUX"            "m.K/s"   #SAVEMARS2 hfmax_th
-state  real  WSTAR     ij   misc  1  -  rhd   "WSTAR"     "FREE CONVECTION VELOCITY FROM TH" "m/s"    #SAVEMARS2 wstar
-state  real  Z0SET     ij   misc  1  -  rd   "Z0SET"     "SET SURFACE ROUGHNESS"           "m"       #SAVEMARS2 z0
+state  real  wTSURF     ij   misc  1  -  rhd  "TSURF"     "SURF TEMPERATURE"                "K"       #SAVEMARS2 tsurf
+state  real  wCO2ICE    ij   misc  1  -  rd   "CO2ICE"    "CO2 GROUND ICE"                  "kg m-2"  #SAVEMARS2 co2ice
+state  real  wEMISSIV   ij   misc  1  -  rd   "EMISSIV"   "EMISSIVITY"                      ""        #SAVEMARS2 emis
+state  real  wTHETA     ij   misc  1  -  rd   "THETA"     "SLOPE INCLINATION"               "deg"     #SAVEMARS2 theta_sl
+state  real  wPSI       ij   misc  1  -  rd   "PSI"       "SLOPE ORIENTATION"               "deg"     #SAVEMARS2 psi_sl
+state  real  wTAU_DUST  ij   misc  1  -  rhd   "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""        #SAVEMARS2 tauref
+state  real  wSWDOWNZ   ij   misc  1  -  rhd   "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"   #SAVEMARS2 fluxsurf_sw_tot      
+state  real  wLWDOWNZ   ij   misc  1  -  rhd   "LWDOWNZ"   "DOWNWARD LW FLUX AT SURFACE"     "W m-2"   #SAVEMARS2 fluxsurf_lw
+state  real  wSWUP      ij   misc  1  -  rd   "SWUP"      "UPWARD SW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_sw_tot    
+state  real  wLWUP      ij   misc  1  -  rd   "LWUP"      "UPWARD LW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_lw
+state  real  wMTOT      ij   misc  1  -  rhd   "MTOT"      "TOTAL MASS WATER VAPOR in pmic"  "pmic"    #SAVEMARS2 mtot
+state  real  wICETOT    ij   misc  1  -  rhd   "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"  #SAVEMARS2 icetot
+state  real  wRAVE      ij   misc  1  -  rhd   "RAVE"      "MEAN ICE RADIUS"                 "m"       #SAVEMARS2 rave
+state  real  wRICE      ikj  misc  1  -  rd   "RICE"      "ICE RADIUS"                      "m"       #SAVEMARS3 rice
+state  real  wHR_SW     ikj  misc  1  -  rd   "HR_SW"     "HEATING RATE SW"                 "K/s"     #SAVEMARS3 zdtsw
+state  real  wHR_LW     ikj  misc  1  -  rd   "HR_LW"     "HEATING RATE LW"                 "K/s"     #SAVEMARS3 zdtlw
+state  real  wHR_SH     ikj  misc  1  -  rd   "HR_SH"     "HEATING RATE sens. heat"         "K/s"     #SAVEMARS3 zdtdif
+state  real  wQSURFICE  ij   misc  1  -  rhd   "QSURFICE"  "WATER ICE AT SURFACE"            "kg m-2"  #SAVEMARS2 qsurfice
+state  real  wQSURFDUST ij   misc  1  -  rd   "QSURFDUST" "DUST MASS AT SURFACE"            "kg m-2"  #SAVEMARS2 qsurfdust
+state  real  wRDUST     ikj  misc  1  -  rd   "RDUST"     "DUST RADIUS"                     "m"       #SAVEMARS3 rdust
+state  real  wHR_NIR    ikj  misc  1  -  rd   "HR_NIR"    "HEATING RATE nirco2"             "K/s"     #SAVEMARS3 zdtnirco2
+state  real  wHR_NLTE   ikj  misc  1  -  rd   "HR_NLTE"   "HEATING RATE nlte"               "K/s"     #SAVEMARS3 zdtnlte
+state  real  wALBBARE   ij   misc  1  -  rhd   "ALBBARE"   "SOIL ALBEDO"                     ""        #SAVEMARS2 albedodat
+state  real  wVMR_ICE   ikj  misc  1  -  rhd   "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"     #SAVEMARS3 vmr
+state  real  wTAU_ICE   ij   misc  1  -  rhd   "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""        #SAVEMARS2 tauTES
+state  real  wPDTZ      ikj  misc  1  -  rd   "PDT"       "TEMP TENDENCY"                   "K s-1"   #SAVEMARS3 pdt
+state  real  wZMAX_TH   ij   misc  1  -  rhd   "ZMAX_TH"   "MAXIMUM LEVEL REACHED IN TH"     "m"       #SAVEMARS2 zmax_th
+state  real  wHFMAX_TH  ij   misc  1  -  rhd   "HFMAX_TH"  "MAXIMUM TH HEAT FLUX"            "m.K/s"   #SAVEMARS2 hfmax_th
+state  real  wWSTAR     ij   misc  1  -  rhd   "WSTAR"     "FREE CONVECTION VELOCITY FROM TH" "m/s"    #SAVEMARS2 wstar
+state  real  wZ0SET     ij   misc  1  -  rd   "Z0SET"     "SET SURFACE ROUGHNESS"           "m"       #SAVEMARS2 z0
 
 ####
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F.new
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F.new	(revision 1234)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/phys/module_lmd_driver.F.new	(revision 1234)
@@ -0,0 +1,1287 @@
+!*******************************************************************************
+! PURPOSE: LMD_driver is the WRF mediation layer routine that provides the
+! interface to LMD physics packages in the WRF model layer. For those familiar
+! with the LMD GCM, the aim of this driver is to do part of the job of the
+! calfis.F routine.
+!*******************************************************************************
+! AUTHOR: A. Spiga - January 2007
+!*******************************************************************************
+! UPDATES:  - included all final updates for the paper - March 2008
+!           - general cleaning of code and comments - October 2008 
+!           - additions for idealized cases - January 2009
+!           - additions for new soil model in physics - January 2010 
+!           - unified module_lmd_driver: old, new phys and LES - February 2011
+!           - a new paradigm to prepare nested simulations - April 2014
+!*******************************************************************************
+MODULE module_lmd_driver
+CONTAINS
+    SUBROUTINE lmd_driver(id,max_dom,DT,ITIMESTEP,XLAT,XLONG, &
+        IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME,KMS,KME, &
+        i_start,i_end,j_start,j_end,kts,kte,num_tiles, &
+        DX,DY, &
+        MSFT,MSFU,MSFV, &
+        GMT,JULYR,JULDAY, &
+        P8W,DZ8W,T8W,Z,HT, &
+        U,V,W,TH,T,P,EXNER,RHO, &
+        PTOP, &
+        RADT,CUDT, &
+        TSK,PSFC, &
+        RTHBLTEN,RUBLTEN,RVBLTEN, &
+        num_3d_s,SCALAR, &
+        MARS_MODE, &
+        MARS_ALB,MARS_TI,MARS_CICE,MARS_EMISS, &
+        MARS_WICE, &
+        MARS_TSOIL, &
+        MARS_Q2, &
+        MARS_TSURF, &
+        MARS_FLUXRAD, &
+        MARS_WSTAR, &
+        MARS_ISOIL, &
+        MARS_DSOIL, &
+        MARS_Z0, &
+        CST_Z0, &
+        MARS_GW, &
+        NUM_SOIL_LAYERS, &
+        CST_AL, &
+        CST_TI, &
+        isfflx, &
+        diff_opt, &
+        km_opt, &
+        UAVE, USTD, &
+        VAVE, VSTD, &
+        WAVE, WSTD, &
+        TAVE, TSTD, &
+        HISTORY_INTERVAL, &
+#ifndef NOPHYS
+#include "module_lmd_driver_output1.inc"
+#endif
+        SLPX,SLPY,RESTART)
+! NB: module_lmd_driver_output1.inc : output arguments generated from Registry
+
+
+
+
+
+!==================================================================
+! USES
+!==================================================================
+   USE module_model_constants
+   USE module_wrf_error
+
+#ifndef NOPHYS
+   !! MODULES FROM LMD PHYSICS
+   use slope_mod, only: theta_sl, psi_sl
+   use surfdat_h, only: phisfi, albedodat,  &
+                        zmea, zstd, zsig, zgam, zthe, z0, &
+                        emissiv,albedice,iceradius, &
+                        emisice,dtemisice, &
+                        z0_default, &
+                        tsurf, co2ice, emis, qsurf
+   use comsoil_h, only: inertiedat,mlayer,layer, &
+                        volcapa, tsoil
+   use planete_h, only: year_day,periheli,aphelie, &
+                        peri_day,obliquit,emin_turb, &
+                        lmixmin
+   use turb_mod, only: q2,wstar,hfmax_th
+   use dimradmars_mod, only: fluxrad
+   use tracer_mod, only: noms
+   use comcstfi_h, only: omeg,mugaz
+   use comgeomfi_h, only: ini_fillgeom
+#endif
+
+!==================================================================
+  IMPLICIT NONE
+!==================================================================
+
+!==================================================================
+! COMMON
+!==================================================================
+
+#ifndef NOPHYS
+!
+! INCLUDE AUTOMATIQUEMENT GENERE A PARTIR DU REGISTRY
+!
+include "../mars_lmd/libf/phymars/wrf_output_2d.h"
+include "../mars_lmd/libf/phymars/wrf_output_3d.h"
+#endif
+
+!==================================================================
+! VARIABLES
+!==================================================================
+
+! WRF Dimensions
+INTEGER, INTENT(IN   )    ::   &
+      ids,ide,jds,jde,kds,kde, &
+      ims,ime,jms,jme,kms,kme, &
+            kts,kte,num_tiles, &
+              NUM_SOIL_LAYERS
+INTEGER, DIMENSION(num_tiles), INTENT(IN) ::  &
+      i_start,i_end,j_start,j_end   
+! Scalars
+INTEGER, INTENT(IN  ) :: JULDAY, itimestep, julyr,id,max_dom,MARS_MODE
+INTEGER, INTENT(IN  ) :: isfflx,diff_opt,km_opt
+REAL, INTENT(IN  ) :: GMT,dt,dx,dy,RADT,CUDT
+REAL, INTENT(IN  ) :: CST_AL, CST_TI 
+REAL, INTENT(IN  ) :: PTOP
+! 2D arrays          
+REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  :: &
+     MSFT,MSFU,MSFV, &
+     XLAT,XLONG,HT,  &
+     MARS_ALB,MARS_TI,MARS_EMISS, &
+     SLPX,SLPY
+REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT   )  :: &
+     MARS_CICE,MARS_WICE, &
+     MARS_TSURF
+! 3D arrays 
+REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
+     dz8w,p8w,p,exner,t,t8w,rho,u,v,w,z,th
+REAL, DIMENSION( ims:ime, kms:kme+1, jms:jme ), INTENT(INOUT ) :: &
+     MARS_Q2
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: &
+     UAVE, USTD, VAVE, VSTD, WAVE, WSTD, TAVE, TSTD
+INTEGER, INTENT(IN   ) :: HISTORY_INTERVAL
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(INOUT )  :: &
+     MARS_TSOIL 
+REAL, INTENT(IN  ) :: CST_Z0
+REAL, DIMENSION( ims:ime, NUM_SOIL_LAYERS, jms:jme ), INTENT(IN   )  :: &
+     MARS_ISOIL, MARS_DSOIL         
+REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  :: &
+     MARS_Z0
+REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT   )  :: &
+     MARS_FLUXRAD,MARS_WSTAR
+REAL, DIMENSION( ims:ime, 5, jms:jme ), INTENT(IN   )  :: &
+     MARS_GW
+! 4D arrays
+INTEGER, INTENT(IN ) :: num_3d_s
+REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:num_3d_s ), INTENT(INOUT ) :: &
+     scalar
+! Logical
+LOGICAL, INTENT(IN ) :: restart
+
+!-------------------------------------------
+! OUTPUT VARIABLES
+!-------------------------------------------
+!
+! Generated from Registry
+!
+! default definitions :
+! 2D : TSK, PSFC
+! 3D : RTHBLTEN,RUBLTEN,RVBLTEN
+#ifndef NOPHYS
+#include "module_lmd_driver_output2.inc"
+   REAL, DIMENSION(:,:), ALLOCATABLE :: output_tab2d
+   REAL, DIMENSION(:,:,:), ALLOCATABLE :: output_tab3d
+#else
+   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)  :: PSFC,TSK
+   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT)  :: RTHBLTEN,RUBLTEN,RVBLTEN
+#endif
+!-------------------------------------------
+! OUTPUT VARIABLES
+!-------------------------------------------
+
+
+!-------------------------------------------   
+! LOCAL  VARIABLES
+!-------------------------------------------
+   INTEGER ::    i,j,k,its,ite,jts,jte,ij,kk
+   INTEGER ::    subs,iii
+   REAL :: lay1,alpha
+   
+
+   ! *** for Mars Mode 20 ***
+   REAL ::    tau_decay
+   ! *** ----------------------- ***
+
+   ! *** for LMD physics
+   ! ------> inputs:
+   INTEGER :: ngrid,nlayer,nq,nsoil 
+   REAL :: pday,ptime,MY  
+   REAL :: aire_val,lat_val,lon_val
+   REAL :: phisfi_val,albedodat_val,inertiedat_val
+   REAL :: tsurf_val,co2ice_val,emis_val
+   REAL :: zmea_val,zstd_val,zsig_val,zgam_val,zthe_val
+   REAL :: theta_val, psi_val
+   LOGICAL :: firstcall,lastcall,tracerdyn
+   REAL,DIMENSION(:),ALLOCATABLE :: plon, plat, parea
+   REAL,DIMENSION(:),ALLOCATABLE :: q2_val, qsurf_val, tsoil_val
+   REAL :: wstar_val,fluxrad_val
+   REAL,DIMENSION(:),ALLOCATABLE :: isoil_val, dsoil_val
+   REAL :: z0_val
+   ! ----------	
+   REAL,DIMENSION(:,:),ALLOCATABLE :: pplev,pplay,pphi,pu,pv,pt,pw
+   REAL,DIMENSION(:,:,:),ALLOCATABLE :: pq  
+
+   ! <------ outputs:
+   !     physical tendencies
+   REAL,DIMENSION(:),ALLOCATABLE :: pdpsrf
+   REAL,DIMENSION(:,:),ALLOCATABLE :: pdu,pdv,pdt
+   REAL,DIMENSION(:,:,:),ALLOCATABLE :: pdq
+   ! ... intermediate arrays
+   REAL, DIMENSION(:), ALLOCATABLE  :: &
+     dz8w_prof,p8w_prof,p_prof,t_prof,t8w_prof, &
+     u_prof,v_prof,z_prof
+   REAL, DIMENSION(:,:), ALLOCATABLE :: q_prof
+
+   ! Additional control variables
+   INTEGER :: sponge_top,relax,ips,ipe,jps,jpe,kps,kpe
+   REAL :: elaps, ptimestep
+   INTEGER :: test
+   REAL :: wappel_phys
+   LOGICAL :: flag_LES
+   LOGICAL, SAVE :: flag_first_restart
+!**************************************************
+! IMPORTANT: pour travailler avec grid nesting
+! --- deux comportements distincts du save
+! ... ex1: ferme planeto avec PGF+MPI: standard
+! ... ex2: iDataPlex avec IFORT+MPI: SPECIAL_NEST_SAVE 
+!**************************************************
+#ifdef SPECIAL_NEST_SAVE
+      !  une dimension supplementaire liee au nest
+      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
+             dp_save
+      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
+             du_save, dv_save, dt_save
+      REAL, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: &
+             dq_save      
+#else
+      REAL, DIMENSION(:), ALLOCATABLE, SAVE :: &
+             dp_save
+      REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: &
+             du_save, dv_save, dt_save
+      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
+             dq_save      
+#endif
+
+!!!IDEALIZED IDEALIZED
+      REAL :: lat_input, lon_input, ls_input, lct_input
+!!!IDEALIZED IDEALIZED
+
+
+!REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: &
+!     UMAX, UMIN, VMAX, VMIN, WMAX, WMIN, TMAX, TMIN
+
+
+
+!==================================================================
+! CODE
+!==================================================================
+
+IF (JULYR .ne. 9999) THEN
+    flag_LES = .false.  ! "True" LES is not available in this version
+    PRINT *, '*** REAL-CASE SIMULATION ***'
+ELSE
+     PRINT *, '*** IDEALIZED SIMULATION ***'
+     IF ((diff_opt .eq. 2) .and. (km_opt .eq. 2)) THEN
+          PRINT *, '*** type: LES ***'
+          PRINT *, '*** diff_opt = 2 *** km_opt = 2'
+          PRINT *, '*** forcing is isfflx = ',isfflx 
+          flag_LES = .true.
+          !! SPECIAL LES
+     ELSE
+          PRINT *, '*** type: not LES ***'
+          PRINT *, '*** diff_opt = ',diff_opt
+          PRINT *, '*** km_opt = ',km_opt
+          flag_LES = .false.
+          !! IDEALIZED, no LES
+          !! cependant, ne veut-on pas pouvoir
+          !! prescrire un flux en idealise ??
+     ENDIF
+ENDIF
+
+
+print *,'** Mars ** DOMAIN',id
+
+!-------------------------!
+! TWEAK on WRF DIMENSIONS !
+!-------------------------!
+its = i_start(1)   ! define here one big tile
+ite = i_end(num_tiles)
+jts = j_start(1)
+jte = j_end(num_tiles)
+!!
+IF (flag_LES .eqv. .false.) THEN
+ relax=0
+ sponge_top=0               ! another value than 0 triggers instabilities  
+ IF (id .gt. 1) relax=2     ! fix to avoid noise in nesting simulations ; 1 >> too much noise ...
+ENDIF
+ips=its
+ipe=ite
+jps=jts
+jpe=jte
+IF (flag_LES .eqv. .false.) THEN
+ IF (ips .eq. ids)   ips=its+relax !! IF tests necesary for parallel runs
+ IF (ipe .eq. ide-1) ipe=ite-relax
+ IF (jps .eq. jds)   jps=jts+relax
+ IF (jpe .eq. jde-1) jpe=jte-relax
+ENDIF
+kps=kts         !! start at surface 
+IF (flag_LES .eqv. .false.) THEN
+ kpe=kte-sponge_top
+ELSE
+ PRINT *, '*** IDEALIZED SIMULATION: LES *** kpe=kte'
+ kpe=kte !-sponge_top
+ENDIF
+
+!----------------------------!
+! DIMENSIONS FOR THE PHYSICS !
+!----------------------------!
+!day_ini = JULDAY - 1      !! GCM convention  !! pday at firstcall is day_ini
+wappel_phys = RADT
+ptimestep = dt*wappel_phys            ! physical timestep (s)
+ngrid=(ipe-ips+1)*(jpe-jps+1)         ! size of the horizontal grid
+nlayer = kpe-kps+1                    ! number of vertical layers: nlayermx
+nsoil = NUM_SOIL_LAYERS               ! number of soil layers: nsoilmx
+if (num_3d_s > 1) then                ! number of advected fields
+        nq = num_3d_s-1               
+else
+        nq = 1
+endif
+! **** needed but hardcoded
+lastcall = .false.
+! **** needed but hardcoded
+
+          !PRINT *, ips, ipe, jps, jpe
+          !PRINT *, ngrid
+
+elaps = (float(itimestep)-1.)*dt  ! elapsed seconds of simulation
+!----------------------------------------------!
+! what is done at the first step of simulation !
+!----------------------------------------------!
+IF (elaps .eq. 0.) THEN
+   flag_first_restart = .false.
+ELSE
+   flag_first_restart=flag_first_restart.OR.(.NOT. ALLOCATED(dp_save))
+ENDIF
+
+IF ((elaps .eq. 0.).or.flag_first_restart) THEN
+firstcall=.true.  !! for continuity with GCM, physics are always called at the first WRF timestep
+  !firstcall=.false. !! just in case you'd want to get rid of the physics 
+test=0
+#ifdef SPECIAL_NEST_SAVE
+PRINT *, 'ALLOCATED dp_save ???', ALLOCATED( dp_save ), id
+IF( .NOT. ALLOCATED( dp_save  ) ) THEN
+   PRINT *, '**** check **** OK I ALLOCATE one save ARRAY for all NESTS ', max_dom, id
+   !! here are the arrays that would be useful to save physics tendencies
+   ALLOCATE(dp_save(ngrid,max_dom))
+   ALLOCATE(du_save(ngrid,nlayer,max_dom))
+   ALLOCATE(dv_save(ngrid,nlayer,max_dom))
+   ALLOCATE(dt_save(ngrid,nlayer,max_dom))
+   ALLOCATE(dq_save(ngrid,nlayer,nq,max_dom))
+   dp_save(:,:)=0.    !! initialize these arrays ...
+   du_save(:,:,:)=0.
+   dv_save(:,:,:)=0.
+   dt_save(:,:,:)=0.
+   dq_save(:,:,:,:)=0.
+ENDIF
+IF (id .lt. max_dom) THEN
+   flag_first_restart=.true.
+ELSE
+   flag_first_restart=.false.
+ENDIF
+#else
+IF( .NOT. ALLOCATED( dp_save  ) ) THEN
+ALLOCATE(dp_save(ngrid)) !! here are the arrays that would be useful to save physics tendencies
+ALLOCATE(du_save(ngrid,nlayer))
+ALLOCATE(dv_save(ngrid,nlayer))
+ALLOCATE(dt_save(ngrid,nlayer))
+ALLOCATE(dq_save(ngrid,nlayer,nq))
+ENDIF
+dp_save(:)=0.    !! initialize these arrays ...
+du_save(:,:)=0.
+dv_save(:,:)=0.
+dt_save(:,:)=0.
+dq_save(:,:,:)=0.
+flag_first_restart=.false.
+#endif
+
+!! put here some general information you'd like to print just once
+    print *, 'TILES: ', i_start,i_end, j_start, j_end  ! numbers for simple runs, arrays for parallel runs
+    print *, 'DOMAIN: ', ids, ide, jds, jde
+    print *, 'MEMORY: ', ims, ime, jms, jme
+    print *, 'ADVECTED TRACERS: ', num_3d_s-1
+    print *, 'PHYSICS IS CALLED EACH...',wappel_phys
+!! put here some general information you'd like to print just once
+ELSE
+!-------------------------------------------------!
+! what is done for the other steps of simulation  !
+!-------------------------------------------------!
+IF (wappel_phys .gt. 0.) THEN
+firstcall = .false.
+test = MODULO(itimestep-1,int(wappel_phys)) ! WRF time is integrated time (itimestep=1 at the end of first step) 
+                                            ! -- same strategy as diagfi in the LMD GCM 
+                                            ! -- arguments of modulo must be of the same kind (here: integers)
+ELSE
+firstcall = .false.
+test = 9999   
+ENDIF
+ENDIF 
+
+!!!******!!
+!!! TIME !!
+!!!******!!
+IF (JULYR .ne. 9999) THEN
+  !
+  ! specified
+  !
+  ptime = (GMT + elaps/3700.) !! universal time (0<ptime<1): ptime=0.5 at 12:00 UT 
+  ptime = MODULO(ptime,24.)   !! the two arguments of MODULO must be of the same type 
+  ptime = ptime / 24.
+  pday = (JULDAY - 1 + INT((3700*GMT+elaps)/88800))
+  pday = MODULO(int(pday),669)
+  MY = (JULYR-2000) + (88800.*(JULDAY - 1)+3700.*GMT+elaps)/59496000.
+  MY = INT(MY)
+ELSE
+  !
+  ! idealized
+  !
+  PRINT *,'** Mars ** IDEALIZED SIMULATION'
+  PRINT *,'** Mars ** BEWARE: input_coord must be here'
+  open(unit=14,file='input_coord',form='formatted',status='old')
+  rewind(14)
+  read(14,*) lon_input
+  read(14,*) lat_input
+  read(14,*) ls_input
+  read(14,*) lct_input
+  close(14)
+  ptime = lct_input - lon_input / 15. + elaps/3700.
+  ptime = MODULO(ptime,24.)
+  ptime = ptime / 24.
+  pday = floor(ls2sol(ls_input)) + INT((3700*(lct_input - lon_input / 15.) + elaps)/88800)
+  pday = MODULO(int(pday),669)
+  MY = 2024
+  !day_ini = floor(ls2sol(ls_input)) !! pday at firstcall is day_ini
+ENDIF
+print *,'** Mars ** TIME IS', pday, ptime*24.
+
+
+!----------!
+! ALLOCATE !
+!----------!
+ALLOCATE(pdpsrf(ngrid))
+ALLOCATE(pdu(ngrid,nlayer))
+ALLOCATE(pdv(ngrid,nlayer))
+ALLOCATE(pdt(ngrid,nlayer))
+ALLOCATE(pdq(ngrid,nlayer,nq))
+!!!
+!!! BIG LOOP : 1. no call for physics, used saved values
+!!!
+IF (test.NE.0) THEN
+print *,'** Mars ** NO CALL FOR PHYSICS, go to next step...',test
+#ifdef SPECIAL_NEST_SAVE
+pdpsrf(:)=dp_save(:,id)
+pdu(:,:)=du_save(:,:,id)
+pdv(:,:)=dv_save(:,:,id)
+pdt(:,:)=dt_save(:,:,id)
+pdq(:,:,:)=dq_save(:,:,:,id)
+#else
+pdpsrf(:)=dp_save(:)
+pdu(:,:)=du_save(:,:)
+pdv(:,:)=dv_save(:,:)
+pdt(:,:)=dt_save(:,:)
+pdq(:,:,:)=dq_save(:,:,:)
+#endif
+!!!
+!!! BIG LOOP : 2. calculate physical tendencies
+!!!
+ELSE
+!----------!
+! ALLOCATE !
+!----------!
+! inputs ...
+#ifndef NOPHYS
+ALLOCATE(output_tab2d(ngrid,n2d))
+ALLOCATE(output_tab3d(ngrid,nlayer,n3d))
+#endif
+ALLOCATE(q2_val(nlayer+1))
+ALLOCATE(qsurf_val(nq))
+ALLOCATE(tsoil_val(nsoil))
+ALLOCATE(isoil_val(nsoil))
+ALLOCATE(dsoil_val(nsoil))
+ALLOCATE(pplev(ngrid,nlayer+1))  !!!!!
+ALLOCATE(pplay(ngrid,nlayer))    !!!!!
+ALLOCATE(pphi(ngrid,nlayer))     !!!!!
+ALLOCATE(pu(ngrid,nlayer))       !!!!!
+ALLOCATE(pv(ngrid,nlayer))       !!!!!
+ALLOCATE(pt(ngrid,nlayer))       !!!!!
+ALLOCATE(pw(ngrid,nlayer))       !!!!!
+ALLOCATE(pq(ngrid,nlayer,nq))    !!!!!
+! interm
+ALLOCATE(dz8w_prof(nlayer))
+ALLOCATE(p8w_prof(nlayer))
+ALLOCATE(p_prof(nlayer))
+ALLOCATE(t_prof(nlayer))
+ALLOCATE(t8w_prof(nlayer))
+ALLOCATE(u_prof(nlayer))
+ALLOCATE(v_prof(nlayer))
+ALLOCATE(z_prof(nlayer))
+ALLOCATE(q_prof(nlayer,nq))
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! PREPARE PHYSICS INPUTS !!  
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!! INITIALIZE AND ALLOCATE EVERYTHING !!
+allocation_firstcall: IF (firstcall .EQV. .true.) THEN
+
+  !! TRACERS
+  ALLOCATE(noms(nq)) !! est fait dans initracer normalement
+                     !! tableau dans tracer_mod.F90
+
+  !!! name of tracers to mimic entries in tracer.def 
+  !!! ----> IT IS IMPORTANT TO KEEP SAME ORDERING AS IN THE REGISTRY !!!!
+  !!!                               SAME NAMING   AS IN THE LMD PHYSICS !!!!
+  SELECT CASE (MARS_MODE)
+    CASE(0,10)  
+      noms(nq) = 'co2'
+    CASE(1)    ! scalar:qh2o,qh2o_ice
+      noms(1)  = 'h2o_vap'
+      noms(2)  = 'h2o_ice'
+    CASE(2)    ! scalar:qdust
+      noms(1)  = 'dust01'
+    CASE(3)    ! scalar:qdust,qdustn
+      noms(1)  = 'dust_mass'
+      noms(2)  = 'dust_number'
+    CASE(11)   ! scalar:qh2o,qh2o_ice,qdust,qdustn
+      noms(1)  = 'h2o_vap'
+      noms(2)  = 'h2o_ice'
+      noms(3)  = 'dust_mass'
+      noms(4)  = 'dust_number'
+    CASE(12)
+      noms(1)  = 'h2o_vap'
+      noms(2)  = 'h2o_ice'
+      noms(3)  = 'dust_mass'
+      noms(4)  = 'dust_number'
+      noms(5)  = 'ccn_mass'
+      noms(6)  = 'ccn_number'
+    CASE(20) 
+      noms(1) = 'qtrac1'
+    CASE(42)
+      noms(1)  = 'co2'
+      noms(2)  = 'co'
+      noms(3)  = 'o'
+      noms(4)  = 'o1d'
+      noms(5)  = 'o2'
+      noms(6)  = 'o3'
+      noms(7)  = 'h'
+      noms(8)  = 'h2'
+      noms(9)  = 'oh'
+      noms(10)  = 'ho2'
+      noms(11)  = 'h2o2'
+      noms(12)  = 'ch4'
+      noms(13)  = 'n2'
+      noms(14)  = 'ar'
+      noms(15)  = 'h2o_ice'
+      noms(16)  = 'h2o_vap'
+      noms(17)  = 'dust_mass'
+      noms(18)  = 'dust_number'
+  END SELECT
+
+!!*******************************************!!
+
+#ifndef NOPHYS
+  !! PHYSICS VARIABLES (cf. iniphysiq in LMD GCM)
+  !! parameters are defined in the module_model_constants.F WRF routine
+  PRINT *,'** Mars ** INITIALIZE ARRAYS FOR PHYSICS'
+  CALL phys_state_var_init(ngrid,nlayer,nq,&
+                           wdaysec,ptimestep,1/reradius,g,r_d,cp)
+
+  !! Fill planetary parameters in modules
+  !! Values defined in the module_model_constants.F WRF routine
+
+  !! comcstfi_h
+  omeg = womeg                
+  mugaz = wmugaz  
+  print*,"check: omeg,mugaz"
+  print*,omeg,mugaz
+  !! planet_h
+  year_day = wyear_day
+  periheli = wperiheli
+  aphelie = waphelie
+  peri_day = wperi_day
+  obliquit = wobliquit
+  emin_turb = wemin_turb
+  lmixmin = wlmixmin
+  print*,"check: year_day,periheli,aphelie,peri_day,obliquit"
+  print*,year_day,periheli,aphelie,peri_day,obliquit
+  print*,"check: emin_turb,lmixmin"
+  print*,emin_turb,lmixmin
+  !! surfdat_h
+  emissiv = wemissiv
+  emisice(1) = wemissiceN
+  emisice(2) = wemissiceS
+  albedice(1) = walbediceN
+  albedice(2) = walbediceS
+  iceradius(1) = wiceradiusN
+  iceradius(2) = wiceradiusS
+  dtemisice(1) = wdtemisiceN
+  dtemisice(2) = wdtemisiceS
+  z0_default = wz0
+  print*,"check: z0def,emissiv,emisice,albedice,iceradius,dtemisice"
+  print*,z0_default,emissiv,emisice,albedice,iceradius,dtemisice
+  !! comsoil_h
+  volcapa = wvolcapa
+  print*,"check: volcapa",volcapa
+
+  !! Read callphys.def !!
+  call conf_phys(nq)
+
+#endif
+
+ENDIF allocation_firstcall
+
+!!*****************************!!
+!! PREPARE "FOOD" FOR PHYSIQ.F !!
+!!*****************************!!
+
+DO j = jps,jpe
+DO i = ips,ipe
+
+!!*******************************************!!
+!! FROM 3D WRF FIELDS TO 1D PHYSICS PROFILES !!
+!!*******************************************!!
+
+!-----------------------------------!
+! 1D subscript for physics "cursor" !
+!-----------------------------------!
+subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
+  
+!--------------------------------------!
+! 1D columns : inputs for the physics  !
+! ... vert. coord., meteor. fields ... !
+!--------------------------------------!
+dz8w_prof(:) = dz8w(i,kps:kpe,j)   ! dz between full levels (m)   
+p8w_prof(:) = p8w(i,kps:kpe,j)     ! pressure full level (Pa) >> pplev
+p_prof(:) = p(i,kps:kpe,j)         ! pressure half level (Pa) >> pplay
+t_prof(:) = t(i,kps:kpe,j)         ! temperature half level (K) >> pt
+t8w_prof(:) = t8w(i,kps:kpe,j)     ! temperature full level (K)
+u_prof(:) = u(i,kps:kpe,j)         ! zonal wind (A-grid: unstaggered) half level (m/s) >> pu  
+v_prof(:) = v(i,kps:kpe,j)         ! meridional wind (A-grid: unstaggered) half level (m/s) >> pv
+z_prof(:) = z(i,kps:kpe,j)         ! geopotential height half level (m) >> pphi/g
+
+!--------------------------------!
+! specific treatment for tracers !
+!--------------------------------!
+IF (MARS_MODE .EQ. 0) THEN
+    q_prof(:,1)=0.95
+ELSE
+    q_prof(:,1:nq) = SCALAR(i,kps:kpe,j,2:nq+1)  !! the names were set above !! one dummy tracer in WRF
+ENDIF
+
+IF (MARS_MODE .EQ. 20) THEN
+  IF (firstcall .EQV. .true. .and. (.not. restart)) THEN
+      q_prof(:,:) = 0.95
+  ENDIF
+ENDIF
+
+IF (firstcall .EQV. .true.) THEN
+  !-----------------!
+  ! Optional output !
+  !-----------------!
+  IF ( (i == ips) .AND. (j == jps) ) THEN
+  PRINT *,'z_prof ',z_prof
+  PRINT *,'dz8w_prof',dz8w_prof
+  PRINT *,'p8w_prof ',p8w_prof
+  PRINT *,'p_prof ',p_prof
+  PRINT *,'t_prof ',t_prof
+  PRINT *,'t8w_prof ',t8w_prof
+  PRINT *,'u_prof ',u_prof
+  PRINT *,'v_prof ',v_prof
+  PRINT *,'q_prof ',q_prof
+  ENDIF
+ENDIF
+
+!---------------------------------------------!
+! in LMD physics, geopotential must be        !
+! expressed with respect to the local surface !
+!---------------------------------------------!
+pphi(subs,:) = g*( z_prof(:)-(z_prof(1)-dz8w_prof(1)/2.) )
+
+!--------------------------------!
+! Dynamic fields for LMD physics !
+!--------------------------------!
+pplev(subs,1:nlayer) = p8w_prof(1:nlayer)  !! NB: last level: no data
+pplay(subs,:) = p_prof(:)
+pt(subs,:) = t_prof(:)
+pu(subs,:) = u_prof(:)
+pv(subs,:) = v_prof(:)
+pw(subs,:) = 0   !! NB: not used in the physics, only diagnostic...
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! for IDEALIZED CASES ONLY
+IF (JULYR .eq. 9999) pplev(subs,nlayer+1)=0.  !! pplev(subs,nlayer+1)=ptop >> NO !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! NOTE:
+! IF ( pplev(subs,nlayer+1) .le. 0 ) pplev(subs,nlayer+1)=ptop
+! cree des diagnostics delirants et aleatoires dans le transfert radiatif 
+
+!---------!
+! Tracers !  
+!---------!
+pq(subs,:,:) = q_prof(:,:)  !! traceurs generiques, seuls noms sont specifiques
+
+ENDDO
+ENDDO
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!! ONE DOMAIN CASE
+!!! --> we simply need to initialize static and physics inputs
+!!! SEVERAL DOMAINS
+!!! --> we update static and physics inputs each time
+!!!     to account for 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+pass_interface: IF ( (firstcall .EQV. .true.) .or. (max_dom .GT. 1) ) THEN
+
+PRINT *,'** Mars ** PASS INTERFACE. EITHER FIRSTCALL or NESTED SIMULATION.'
+
+!!!!!!!!!!!!
+!!! GEOM !!!
+!!!!!!!!!!!!
+
+ALLOCATE(plon(ngrid))
+ALLOCATE(plat(ngrid))
+ALLOCATE(parea(ngrid))
+
+DO j = jps,jpe
+DO i = ips,ipe
+
+  !-----------------------------------!
+  ! 1D subscript for physics "cursor" !
+  !-----------------------------------!
+  subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
+
+  !----------------------------------------!
+  ! Surface of each part of the grid (m^2) !
+  !----------------------------------------!
+  !parea(subs) = dx*dy                           !! 1. idealized cases - computational grid
+  parea(subs) = (dx/msft(i,j))*(dy/msft(i,j))    !! 2. WRF map scale factors - assume that msfx=msfy (msf=covariance)
+  !parea(subs)=dx*dy/msfu(i,j)                   !! 3. special for Mercator GCM-like simulations
+
+  !---------------------------------------------!
+  ! Mass-point latitude and longitude (radians) !
+  !---------------------------------------------!
+  IF (JULYR .ne. 9999) THEN
+   plat(subs) = XLAT(i,j)*DEGRAD
+   plon(subs) = XLONG(i,j)*DEGRAD
+  ELSE
+   !!! IDEALIZED CASE
+   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lat: ',lat_input
+   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION lon: ',lon_input
+   plat(subs) = lat_input*DEGRAD
+   plon(subs) = lon_input*DEGRAD
+  ENDIF 
+
+ENDDO
+ENDDO
+
+!! FILL GEOMETRICAL ARRAYS !!
+call ini_fillgeom(ngrid,plat,plon,parea)
+
+DEALLOCATE(plon)
+DEALLOCATE(plat)
+DEALLOCATE(parea)
+
+!!*******************************************!!
+!!*******************************************!!
+
+DO j = jps,jpe
+DO i = ips,ipe
+
+!-----------------------------------!
+! 1D subscript for physics "cursor" !
+!-----------------------------------!
+subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
+
+!! **************
+!! use slope_mod
+!! **************
+!-------------------!
+! Slope inclination !
+!-------------------!
+IF (JULYR .ne. 9999) THEN
+  theta_val=atan(sqrt( (1000.*SLPX(i,j))**2 + (1000.*SLPY(i,j))**2 ))
+  theta_val=theta_val/DEGRAD
+ELSE
+  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION slope inclination is 0'
+  theta_val=0.
+ENDIF
+!-------------------------------------------!
+! Slope orientation; 0 is north, 90 is east !
+!-------------------------------------------!
+IF (JULYR .ne. 9999) THEN
+  psi_val=-90.*DEGRAD-atan(SLPY(i,j)/SLPX(i,j))
+  if (SLPX(i,j) .ge. 0.) then
+    psi_val=psi_val-180.*DEGRAD
+  endif
+  psi_val=360.*DEGRAD+psi_val
+  psi_val=psi_val/DEGRAD
+  psi_val = MODULO(psi_val+180.,360.)
+ELSE
+  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION slope orientation is 0 (well, whatever)'
+  psi_val=0.
+ENDIF
+theta_sl(subs) = theta_val
+psi_sl(subs) = psi_val
+
+!! **************
+!! use surfdat_h
+!! **************
+!---------------------------------------------------------!
+! Ground geopotential                                     !
+! (=g*HT(i,j), only used in the microphysics: newcondens) !
+!---------------------------------------------------------! 
+phisfi_val=g*(z_prof(1)-dz8w_prof(1)/2.)   !! NB: z_prof are half levels, so z_prof(1) is not the surface
+!---------------------------------!
+! Ground albedo & Thermal Inertia !
+!---------------------------------!
+IF (JULYR .ne. 9999) THEN 
+ IF (CST_AL == 0) THEN
+ albedodat_val=MARS_ALB(i,j)
+ ELSE 
+ albedodat_val=CST_AL
+ IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT ALBEDO ', albedodat_val
+ ENDIF
+ELSE
+ IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION albedo: ', CST_AL
+ albedodat_val=CST_AL
+ENDIF
+!-----------------------------------------!
+! Gravity wave parametrization            !
+! NB: usually 0 in mesoscale applications !
+!-----------------------------------------!
+IF (JULYR .ne. 9999) THEN
+ zmea_val=MARS_GW(i,1,j)
+ zstd_val=MARS_GW(i,2,j)
+ zsig_val=MARS_GW(i,3,j)
+ zgam_val=MARS_GW(i,4,j)
+ zthe_val=MARS_GW(i,5,j)
+ELSE
+ IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION GWdrag OFF'
+ zmea_val=0.
+ zstd_val=0.
+ zsig_val=0.
+ zgam_val=0.
+ zthe_val=0.
+ENDIF
+phisfi(subs) = phisfi_val
+albedodat(subs) = albedodat_val
+zmea(subs) = zmea_val
+zstd(subs) = zstd_val
+zsig(subs) = zsig_val
+zgam(subs) = zgam_val
+zthe(subs) = zthe_val
+!----------------------------!
+! Variable surface roughness !
+!----------------------------!
+IF (JULYR .ne. 9999) THEN
+ IF (CST_Z0 == 0) THEN
+   z0_val = MARS_Z0(i,j)
+ ELSE
+   z0_val = CST_Z0
+   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT SURF ROUGHNESS (m) ',CST_Z0
+ ENDIF
+ELSE
+ IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION z0 (m) ', CST_Z0
+ z0_val=CST_Z0
+ENDIF
+!!!! ADDITIONAL SECURITY. THIS MIGHT HAPPEN WITH OLD INIT FILES.
+IF (z0_val == 0.) THEN
+   IF ( (i == ips) .AND. (j == jps) ) PRINT *, 'WELL, z0 is 0, this is no good. Setting to old defaults value 0.01 m'
+   z0_val = 0.01
+ENDIF
+!!!! ADDITIONAL SECURITY. INTERP+SMOOTH IN GEOGRID MIGHT YIELD NEGATIVE Z0 !!!
+IF (z0_val < 0.) THEN
+   PRINT *, 'WELL, z0 is NEGATIVE, this is impossible. better stop here.'
+   PRINT *, 'advice --> correct interpolation / smoothing of z0 in WPS'
+   PRINT *, '           -- or check the constant value set in namelist.input'
+   STOP
+ENDIF
+z0(subs) = z0_val
+!-----------------------------------------------!
+! Ground temperature, emissivity, CO2 ice cover !
+!-----------------------------------------------!
+IF (.not. restart) THEN
+   tsurf_val=tsk(i,j) ! because tsk is not updated in lmd_driver, the stored value in restart is at t, not t+dt
+                      ! note to aslmd: i dont know why we dont update tsk, so I use MARS_TSURF to be conservative
+ELSE
+   tsurf_val=MARS_TSURF(i,j)
+ENDIF
+tsurf(subs) = tsurf_val
+emis(subs) = MARS_EMISS(i,j)
+co2ice(subs) = MARS_CICE(i,j)
+!-------------------!
+! Tracer at surface !
+!-------------------!
+qsurf_val(:)=0. ! default case
+SELECT CASE (MARS_MODE)
+    CASE(1)
+    qsurf_val(2)=MARS_WICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
+                                   !! ----- retrocompatible ancienne physique
+                                   !! ----- [H2O ice is last tracer in qsurf in LMD physics]
+    CASE(2)   
+    qsurf_val(1)=0.                !! not coupled with lifting for the moment [non remobilise]
+    CASE(3)
+    qsurf_val(1)=q_prof(1,1)                !!! temporaire, a definir       
+    qsurf_val(2)=q_prof(1,2)     
+    CASE(11)
+    qsurf_val(2)=MARS_WICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
+    qsurf_val(3)=0.                !! not coupled with lifting for the moment [non remobilise]
+    CASE(12)
+    qsurf_val(2)=MARS_WICE(i,j)    !! logique avec noms(2) = 'h2o_ice' defini ci-dessus
+    qsurf_val(3)=0.                !! not coupled with lifting for the moment [non remobilise]
+END SELECT
+qsurf(subs,:) = qsurf_val(:)
+
+!! **************
+!! use comsoil_h
+!! **************
+!---------------------------------!
+! Ground albedo & Thermal Inertia !
+!---------------------------------!
+IF (JULYR .ne. 9999) THEN
+ IF (CST_TI == 0) THEN
+ inertiedat_val=MARS_TI(i,j)
+ ELSE
+ inertiedat_val=CST_TI
+ IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** SET CONSTANT THERMAL INERTIA ', inertiedat_val
+ ENDIF
+ELSE
+ IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION inertia: ',CST_TI
+ inertiedat_val=CST_TI
+ENDIF
+!inertiedat(subs) = inertiedat_val
+!--pb de dimensions???!!???
+IF (JULYR .ne. 9999) THEN
+  isoil_val(:)=MARS_ISOIL(i,:,j)
+  dsoil_val(:)=MARS_DSOIL(i,:,j)
+ELSE
+   IF ( nsoil .lt. 18 ) THEN
+       PRINT *,'** Mars ** WRONG NUMBER OF SOIL LAYERS. SHOULD BE 18 AND IT IS ',nsoil
+       STOP
+   ENDIF
+   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION isoil and dsoil standard'
+  do k=1,nsoil
+   isoil_val(k) = inertiedat_val
+!  dsoil_val(k) = sqrt(887.75/3.14)*((2.**(k-0.5))-1.) * inertiedat_val / wvolcapa    ! old setting
+   dsoil_val(k) = 2.E-4 * (2.**(k-0.5-1.))                                            ! new gcm settings
+   IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** ISOIL DSOIL are ',isoil_val(k), dsoil_val(k)
+  enddo
+ENDIF
+inertiedat(subs,:) = isoil_val(:)
+mlayer(0:nsoil-1)=dsoil_val(:)
+    !!!!!!!!!!!!!!!!! DONE in soil_setting.F 
+    ! 1.5 Build layer(); following the same law as mlayer()
+    ! Assuming layer distribution follows mid-layer law:
+    ! layer(k)=lay1*alpha**(k-1)
+    lay1=sqrt(mlayer(0)*mlayer(1))
+    alpha=mlayer(1)/mlayer(0)
+    do iii=1,nsoil
+       layer(iii)=lay1*(alpha**(iii-1))
+    enddo
+!------------------------!
+! Deep soil temperatures ! 
+!------------------------!
+IF (JULYR .ne. 9999) THEN
+  IF (MARS_TSOIL(i,1,j) .gt. 0.) THEN
+   tsoil_val(:)=MARS_TSOIL(i,:,j)
+  ELSE
+   tsoil_val = tsoil_val*0. + tsurf_val
+  ENDIF
+ELSE
+  IF ( (i == ips) .AND. (j == jps) ) PRINT *,'** Mars ** IDEALIZED SIMULATION tsoil is set to tsurf'
+  do k=1,nsoil
+   IF (.not.restart) THEN
+     tsoil_val(k) = tsurf_val
+   ELSE
+     !this is a restart run. We must not set tsoil to tsurf in the init.
+     !tsoil was saved in physiq.F under the name MARS_TSOIL in the restart file
+     !(see Registry)
+     tsoil_val(k)=MARS_TSOIL(i,k,j)
+   ENDIF
+  enddo
+ENDIF
+tsoil(subs,:) = tsoil_val(:)
+
+!! **************
+!! use turb_mod
+!! **************
+IF (.not. restart) THEN
+   q2_val(:) = 1.E-6      !PBL wind variance
+   wstar_val=0.
+ELSE
+   q2_val(:)=MARS_Q2(i,:,j)
+   wstar_val=MARS_WSTAR(i,j)
+ENDIF
+q2(subs,:) = q2_val(:)
+wstar(subs) = wstar_val
+
+!! ********************
+!! use dimradmars_mod
+!! ********************
+IF (.not. restart) THEN
+   fluxrad_val=0.
+ELSE
+   fluxrad_val=MARS_FLUXRAD(i,j)
+ENDIF
+fluxrad(subs) = fluxrad_val
+!! et fluxrad_sky ???!???
+
+ENDDO
+ENDDO
+
+!!---------------------!!
+!! OUTPUT FOR CHECKING !!
+!!---------------------!!
+print*,"check: theta_sl",theta_sl(1),theta_sl(ngrid)
+print*,"check: psi_sl",psi_sl(1),psi_sl(ngrid)
+print*,"check: phisfi",phisfi(1),phisfi(ngrid)
+print*,"check: albedodat",albedodat(1),albedodat(ngrid)
+print*,"check: zmea",zmea(1),zmea(ngrid)
+print*,"check: zstd",zstd(1),zstd(ngrid)
+print*,"check: zsig",zsig(1),zsig(ngrid)
+print*,"check: zgam",zgam(1),zgam(ngrid)
+print*,"check: zthe",zthe(1),zthe(ngrid)
+print*,"check: z0",z0(1),z0(ngrid)
+print*,"check: tsurf",tsurf(1),tsurf(ngrid)
+print*,"check: emis",emis(1),emis(ngrid)
+print*,"check: co2ice",co2ice(1),co2ice(ngrid)
+print*,"check: qsurf",qsurf(1,:),qsurf(ngrid,:)
+print*,"check: inertiedat",inertiedat(1,:),inertiedat(ngrid,:)
+print*,"check: mlayer",mlayer(:)
+print*,"check: layer",layer(:)
+print*,"check: q2",q2(1,:),q2(ngrid,:)
+print*,"check: wstar",wstar(1),wstar(ngrid)
+print*,"check: fluxrad",fluxrad(1),fluxrad(ngrid)
+print*,"check: tsoil",tsoil(1,:),tsoil(ngrid,:)
+print*,"check: noms",noms
+
+ENDIF pass_interface
+
+!!*****************!!
+!! CLEAN THE PLACE !!
+!!*****************!!
+DEALLOCATE(q2_val)
+DEALLOCATE(qsurf_val)
+DEALLOCATE(tsoil_val)
+DEALLOCATE(isoil_val)
+DEALLOCATE(dsoil_val)
+DEALLOCATE(dz8w_prof)
+DEALLOCATE(z_prof)
+DEALLOCATE(p8w_prof)
+DEALLOCATE(p_prof)
+DEALLOCATE(t_prof)
+DEALLOCATE(t8w_prof)
+DEALLOCATE(u_prof)
+DEALLOCATE(v_prof)
+DEALLOCATE(q_prof)
+
+!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!
+!! CALL LMD PHYSICS !!  
+!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!
+
+call_physics : IF (wappel_phys .ne. 0.) THEN
+
+!-------------------------------------------------------------------------------!
+! outputs:									!	
+!    pdu(ngrid,nlayermx)        \						!
+!    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding	!
+!    pdt(ngrid,nlayermx)         /  variables due to physical processes.	!
+!    pdq(ngrid,nlayermx)        /						!
+!    pdpsrf(ngrid)             /						!
+!    tracerdyn                 call tracer in dynamical part of GCM ?		!
+!-------------------------------------------------------------------------------!
+pdpsrf(:)=0.
+pdu(:,:)=0.
+pdv(:,:)=0.
+pdt(:,:)=0.
+pdq(:,:,:)=0.
+#ifndef NOPHYS
+print *, '** Mars ** CALL TO LMD PHYSICS'
+CALL physiq (ngrid,nlayer,nq,                         &
+             firstcall,lastcall,pday,ptime,ptimestep, &
+             pplev,pplay,pphi,pu,pv,pt,pq,pw,         &
+             pdu,pdv,pdt,pdq,pdpsrf,tracerdyn,        &
+             output_tab2d,output_tab3d,flag_LES)
+#endif
+DEALLOCATE(pplev)
+DEALLOCATE(pplay)
+DEALLOCATE(pphi)
+DEALLOCATE(pu)
+DEALLOCATE(pv)
+DEALLOCATE(pt)
+DEALLOCATE(pw)
+DEALLOCATE(pq)
+
+!-------------------------------!
+! PHYSIQ OUTPUT IN THE WRF FILE !
+!-------------------------------!
+#ifndef NOPHYS
+DO j = jps,jpe
+DO i = ips,ipe
+subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
+#include "module_lmd_driver_output3.inc"  
+       !  ^-- generated from Registry
+ENDDO
+ENDDO
+DEALLOCATE(output_tab2d)
+DEALLOCATE(output_tab3d)
+#endif
+
+!---------------------------------------------------------------------------------!
+! PHYSIQ TENDENCIES ARE SAVED TO BE SPLIT WITHIN INTERMEDIATE DYNAMICAL TIMESTEPS !
+!---------------------------------------------------------------------------------!
+#ifdef SPECIAL_NEST_SAVE
+dp_save(:,id)=pdpsrf(:)
+du_save(:,:,id)=pdu(:,:)
+dv_save(:,:,id)=pdv(:,:)
+dt_save(:,:,id)=pdt(:,:)
+dq_save(:,:,:,id)=pdq(:,:,:)
+#else
+dp_save(:)=pdpsrf(:)
+du_save(:,:)=pdu(:,:)
+dv_save(:,:)=pdv(:,:)
+dt_save(:,:)=pdt(:,:)
+dq_save(:,:,:)=pdq(:,:,:)
+#endif
+ENDIF call_physics
+
+ENDIF    ! end of BIG LOOP 2.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!!***************************!!
+!! DEDUCE TENDENCIES FOR WRF !!
+!!***************************!!
+RTHBLTEN(ims:ime,kms:kme,jms:jme)=0.
+RUBLTEN(ims:ime,kms:kme,jms:jme)=0.
+RVBLTEN(ims:ime,kms:kme,jms:jme)=0.
+PSFC(ims:ime,jms:jme)=p8w(ims:ime,kms,jms:jme) ! was done in surface driver in regular WRF 
+!------------------------------------------------------------------!
+! WRF adds the physical and dynamical tendencies                   !
+! thus the tendencies must be passed as 'per second' tendencies    !
+! --when physics is not called, the applied physical tendency      !
+! --is the one calculated during the last call to physics          !
+!------------------------------------------------------------------!
+
+DO j = jps,jpe
+DO i = ips,ipe
+
+  subs = (j-jps)*(ipe-ips+1)+(i-ips+1)
+
+    ! zonal wind
+  RUBLTEN(i,kps:kpe,j) = pdu(subs,kps:kpe) 
+
+    ! meridional wind
+  RVBLTEN(i,kps:kpe,j) = pdv(subs,kps:kpe) 
+
+    ! potential temperature
+    ! (dT = dtheta * exner for isobaric coordinates or if pressure variations are negligible)
+  RTHBLTEN(i,kps:kpe,j) = pdt(subs,kps:kpe) / exner(i,kps:kpe,j)
+
+    ! update surface pressure (cf CO2 cycle in physics)
+    ! here dt is needed
+  PSFC(i,j)=PSFC(i,j)+pdpsrf(subs)*dt 
+
+    ! tracers
+  SCALAR(i,kps:kpe,j,1)=0.
+  SELECT CASE (MARS_MODE)
+    CASE(0)
+      SCALAR(i,kps:kpe,j,:)=0.
+    CASE(20)
+      !! Mars mode 20 : add a passive tracer with radioactive-like decay
+      IF ( (i == ips) .AND. (j == jps) )   print *, 'RADIOACTIVE-LIKE TRACER WITH SOURCE AT SURFACE LAYER.'
+      tau_decay=60.*10. !! why not make it a namelist argument?
+      SCALAR(i,kps:kpe,j,2) = SCALAR(i,kps:kpe,j,2)*exp(-dt/tau_decay)
+      SCALAR(i,1,j,2) = SCALAR(i,1,j,2) + 1. !! this tracer is emitted in the surface layer
+    CASE DEFAULT
+      SCALAR(i,kps:kpe,j,2:nq+1)=SCALAR(i,kps:kpe,j,2:nq+1)+pdq(subs,kps:kpe,1:nq)*dt  !!! here dt is needed
+  END SELECT
+    
+  !------------------------------------!
+  ! Save key variables for restart !  
+  !------------------------------------!
+  MARS_TSOIL(i,:,j) = tsoil(subs,:)
+  MARS_CICE(i,j) = co2ice(subs)
+  MARS_Q2(i,kps:kpe+1,j) = q2(subs,:)
+  MARS_TSURF(i,j) = tsurf(subs)
+  MARS_WSTAR(i,j) = wstar(subs)
+  MARS_FLUXRAD(i,j) = fluxrad(subs)
+  SELECT CASE (MARS_MODE)
+   CASE (1,11,12)
+     MARS_WICE(i,j) = qsurf(subs,2)  !! see above Tracer at surface
+  END SELECT
+
+ENDDO
+ENDDO
+
+DEALLOCATE(pdpsrf)
+DEALLOCATE(pdu)
+DEALLOCATE(pdv)
+DEALLOCATE(pdt)
+DEALLOCATE(pdq)
+
+!!*****!!
+!! END !!
+!!*****!!
+print *,'** Mars ** END LMD PHYSICS'
+!----------------------------------------------------------------!
+! use debug (see solve_em) if you wanna display some message ... !
+!----------------------------------------------------------------!
+END SUBROUTINE lmd_driver
+
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      real function ls2sol(ls)
+
+!c  Returns solar longitude, Ls (in deg.), from day number (in sol),
+!c  where sol=0=Ls=0 at the northern hemisphere spring equinox
+
+      implicit none
+
+!c  Arguments:
+      real ls
+
+!c  Local:
+      double precision xref,zx0,zteta,zz
+!c	xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
+      double precision year_day 
+      double precision peri_day,timeperi,e_elips
+      double precision pi,degrad 
+      parameter (year_day=668.6d0) ! number of sols in a amartian year
+!c      data peri_day /485.0/
+      parameter (peri_day=485.35d0) ! date (in sols) of perihelion
+!c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
+      parameter (timeperi=1.90258341759902d0)
+      parameter (e_elips=0.0934d0)  ! eccentricity of orbit
+      parameter (pi=3.14159265358979d0)
+      parameter (degrad=57.2957795130823d0)
+
+      if (abs(ls).lt.1.0e-5) then
+         if (ls.ge.0.0) then
+            ls2sol = 0.0
+         else
+            ls2sol = year_day
+         end if
+         return
+      end if
+
+      zteta = ls/degrad + timeperi
+      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
+      xref = zx0-e_elips*dsin(zx0)
+      zz = xref/(2.*pi)
+      ls2sol = zz*year_day + peri_day
+      if (ls2sol.lt.0.0) ls2sol = ls2sol + year_day
+      if (ls2sol.ge.year_day) ls2sol = ls2sol - year_day
+
+      return
+      end function ls2sol
+!!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+END MODULE module_lmd_driver
Index: trunk/MESOSCALE/LMD_MM_MARS/makemeso
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/makemeso	(revision 1233)
+++ trunk/MESOSCALE/LMD_MM_MARS/makemeso	(revision 1234)
@@ -162,4 +162,10 @@
 
     conf_wrf="${config}${scenario}${phys}${compilo}_${machine}${single}${testflag}"
+
+    if [ ${debug} -eq 1 ]
+    then
+     conf_wrf="debug_${conf_wrf}"
+    fi
+
     \rm what_folder 2> /dev/null
     echo ${conf_wrf} > what_folder
@@ -242,4 +248,7 @@
               cp ../../SRC/WRFV2/Registry/Registry.EM.newphys Registry/Registry.EM
           fi
+          cd phys
+          ln -sf module_lmd_driver.F.new module_lmd_driver.F
+          cd ..
         else
           ln -sf mars_lmd_new_${scenario} mars_lmd 
@@ -747,4 +756,10 @@
   cd ..
 
+  if [[ "${phys}" == "newphys_" ]]
+  then
+    # include modules for interfacing
+    mv *.mod ../inc/
+  fi
+
   # save a copy 
   #--------------
