Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/gnu_64.configure.wrf
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/gnu_64.configure.wrf	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/gnu_64.configure.wrf	(revision 31)
@@ -0,0 +1,256 @@
+# configure.wrf
+#
+# This file was automatically generated by the configure script in the
+# top level directory. You may make changes to the settings in this
+# file but be aware they will be overwritten each time you run configure.
+# Ordinarily, it is necessary to run configure once, when the code is
+# first installed.
+#
+# To permanently change options, change the settings for your platform
+# in the file arch/configure.defaults then rerun configure.
+#
+SHELL           =       /bin/sh
+DEVTOP          =       `pwd`
+LIBINCLUDE      =       .
+.SUFFIXES: .F .i .o .f90 .c
+
+#### Get core settings from environment (set in compile script)
+#### Note to add a core, this has to be added to.
+
+COREDEFS = -DEM_CORE=$(WRF_EM_CORE) \
+           -DNMM_CORE=$(WRF_NMM_CORE) -DNMM_MAX_DIM=2600 \
+	   -DCOAMPS_CORE=$(WRF_COAMPS_CORE) \
+           -DEXP_CORE=$(WRF_EXP_CORE)
+
+#### Single location for defining total number of domains.  You need
+#### at least 1 + 2*(number of total nests).  For example, 1 coarse
+#### grid + three fine grids = 1 + 2(3) = 7, so MAX_DOMAINS=7.
+
+MAX_DOMAINS	=	21
+
+#### DM buffer length for the configuration flags.
+
+CONFIG_BUF_LEN	=	16384
+
+
+##############################################################################
+#### The settings in this section are defaults that may be overridden by the 
+#### architecture-specific settings in the next section.  
+##############################################################################
+
+##############################################################################
+#### NOTE:  Do not modify these default values here.  To override these 
+####        values, make changes after "Architecture specific settings".  
+##############################################################################
+
+#### Native size (in bytes) of Fortran REAL data type on this architecture ####
+#### Note:  to change real wordsize (for example, to promote REALs from 
+####        4-byte to 8-byte), modify the definition of RWORDSIZE in the 
+####        section following "Architecture specific settings".  Do not 
+####        change NATIVE_RWORDSIZE as is it architecture-specific.  
+NATIVE_RWORDSIZE = 4
+
+#### Default sed command and script for Fortran source files ####
+SED_FTN = sed -f ../arch/standard.sed
+
+# Hack to work around $(PWD) not changing during OSF1 build.  
+# $(IO_GRIB_SHARE_DIR) is reset during the OSF1 build only.  
+IO_GRIB_SHARE_DIR = 
+
+#### ESMF switches                 ####
+#### These are set up by Config.pl ####
+# switch to use separately installed ESMF library for coupling:  1==true
+ESMF_COUPLING       = 0
+# select dependences on module_utility.o
+ESMF_MOD_DEPENDENCE = ../external/esmf_time_f90/module_utility.o
+# select -I options for separately installed ESMF library, if present
+ESMF_MOD_INC        = 
+# select -I options for external/io_esmf vs. external/esmf_time_f90
+ESMF_IO_INC         = -I../external/esmf_time_f90
+# select cpp token for external/io_esmf vs. external/esmf_time_f90
+ESMF_IO_DEFS        = 
+# select build target for external/io_esmf vs. external/esmf_time_f90
+ESMF_TARGET         = esmf_time
+
+##############################################################################
+
+
+#### Architecture specific settings ####
+
+# Settings for PC Linux x86_64 (IA64 and Opteron), GNU Fortran compiler   (Single-threaded, RSL, Allows nesting)
+#
+DMPARALLEL      =       1
+MAX_PROC        =       1024
+FC              =       gfortran
+LD              =       gfortran
+CC              =       gcc -DFSEEKO64_OK 
+SCC		=	$(CC)
+SFC             =       gfortran
+RWORDSIZE       =       $(NATIVE_RWORDSIZE)
+FCOPTIM		=       #-O3 -funroll-loops #-ftree-vectorize -ftree-loop-linear
+FCDEBUG		=	#-g
+FCBASEOPTS      =       -fno-second-underscore -w $(FCDEBUG) -ffree-form -ffree-line-length-none -fconvert=big-endian -frecord-marker=4
+FCFLAGS         =       $(FCOPTIM) $(FCBASEOPTS)
+ARCHFLAGS       =       -DGRIB1 -DINTIO -DWRF_RSL_IO -DRSL -DDM_PARALLEL \
+                        -DIWORDSIZE=4 -DDWORDSIZE=8 -DRWORDSIZE=$(RWORDSIZE) -DLWORDSIZE=4 -DNETCDF \
+			-DLIMIT_ARGS -DSTUBMPI
+PROMOTION	=
+CFLAGS          =       -DLANDREADSTUB -DDM_PARALLEL -DWRF_RSL_IO \
+                        -DMAXDOM_MAKE=$(MAX_DOMAINS) -DMAXPROC_MAKE=$(MAX_PROC) -I../external/RSL/RSL $(ARCHFLAGS)
+# machine-specific flags needed to link in ESMF library (C++ run-time-library, etc.)
+ESMF_LIB_FLAGS  =       
+ESMF_IO_LIB     =       -L../external/esmf_time_f90 -lesmf_time
+ESMF_IO_LIB_EXT =       -L../../external/esmf_time_f90 -lesmf_time
+INCLUDE_MODULES =       -J../main -I../main -I../external/io_netcdf -I../external/io_int -I../external/esmf_time_f90 \
+                        -I../external -I../frame -I../share -I../phys -I../chem -I../inc
+PERL            =       perl
+REGISTRY        =       Registry
+LIB             =       -L../external/io_netcdf -lwrfio_nf -L/donnees/aslmd/MODELES/NETCDF/gfortran_64/netcdf-3.6.1//lib -lnetcdf -L../external/RSL/RSL -lrsl \
+                        -L../external/io_grib1 -lio_grib1 \
+			-L../external/io_int -lwrfio_int \
+			-L../external/io_grib_share -lio_grib_share \
+			../frame/module_internal_header_util.o ../frame/pack_utils.o -L../external/esmf_time_f90 -lesmf_time -L../mars_lmd/libo -llmd /donnees/aslmd/MODELES/NETCDF/gfortran_64/netcdf-3.6.1//lib/libnetcdf.a 
+LDFLAGS         =       $(FCFLAGS)  -fconvert=big-endian
+ENVCOMPDEFS	=	
+WRF_CHEM	=	0 
+CPP             =       /lib/cpp -C -P -traditional
+POUND_DEF	=	-DNO_RRTM_PHYSICS  $(COREDEFS) -DNONSTANDARD_SYSTEM -DCONFIG_BUF_LEN=$(CONFIG_BUF_LEN) -DMAX_DOMAINS_F=$(MAX_DOMAINS)
+CPPFLAGS        =       -I$(LIBINCLUDE) -C -P $(ARCHFLAGS) -I../external/RSL/RSL -C -P `cat ../inc/dm_comm_cpp_flags` $(ENVCOMPDEFS) $(POUND_DEF)
+AR              =       ar ru
+M4              =       m4 -B 14000
+RANLIB          =       ranlib
+NETCDFPATH	=	/donnees/aslmd/MODELES/NETCDF/gfortran_64/netcdf-3.6.1/
+CC_TOOLS        =       cc
+
+externals : wrf_ioapi_includes wrfio_nf wrfio_grib_share wrfio_grib1 ../external/RSL/RSL/librsl.a wrfio_int gen_comms.c module_dm.F esmf_time
+
+gen_comms.c : ../external/RSL/gen_comms.c
+	( /bin/cp ../tools/gen_comms_warning ../tools/gen_comms.c ; \
+          cat ../external/RSL/gen_comms.c >> ../tools/gen_comms.c )
+
+module_dm.F : ../external/RSL/module_dm.F
+	( /bin/cp module_dm_warning module_dm.F ; cat ../external/RSL/module_dm.F >> module_dm.F )
+
+wrfio_nf : 
+	( cd ../external/io_netcdf ; \
+	  make NETCDFPATH=/donnees/aslmd/MODELES/NETCDF/gfortran_64/netcdf-3.6.1/ RANLIB="$(RANLIB)" CPP="$(CPP)" FC="$(FC) $(PROMOTION) $(FCFLAGS) -w" TRADFLAG="-traditional" )
+
+wrfio_grib_share :
+	( cd ../external/io_grib_share ; \
+          make CC="$(SCC)" CFLAGS="$(CFLAGS)" RM="$(RM)" RANLIB="$(RANLIB)" CPP="$(CPP)" FC="$(SFC) $(PROMOTION) -I. $(FCDEBUG) $(FCBASEOPTS) -w" TRADFLAG="-traditional" archive)
+
+wrfio_grib1 :
+	( cd ../external/io_grib1 ; \
+          make CC="$(SCC)" CFLAGS="$(CFLAGS)" RANLIB="$(RANLIB)" CPP="$(CPP)" FC="$(SFC) $(PROMOTION) -I. $(FCDEBUG) $(FCBASEOPTS) -w" TRADFLAG="-traditional" archive)
+
+wrfio_grib2 :
+	( cd ../external/io_grib2 ; \
+          make CC="$(SCC)" CFLAGS="$(CFLAGS) " RM="$(RM)" RANLIB="$(RANLIB)" CPP="$(CPP)" FC="$(SFC) $(PROMOTION) -I. $(FCDEBUG) $(FCBASEOPTS) -w" TRADFLAG="-traditional" FIXED="-Mfixed" archive)
+
+wrfio_int : 
+	( cd ../external/io_int ; \
+          make CC="$(CC)" RANLIB="$(RANLIB)" CPP="$(CPP)" SFC="$(SFC) $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS)" FC="$(SFC) $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS) -w" TRADFLAG="-traditional" all )
+
+esmf_time : 
+	( cd ../external/esmf_time_f90 ; \
+	  make FC="$(FC) $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS)" RANLIB="$(RANLIB)" CPP="$(CPP) -I../../inc -I. $(POUND_DEF)" )
+
+../external/RSL/RSL/librsl.a :
+	( cd ../external/RSL/RSL ; make CC="$(CC)" CFLAGS='-DSTUBS' FC="$(FC) $(PROMOTION)" FFLAGS='' MAX_DOMAINS=$(MAX_DOMAINS) MAX_PROC=$(MAX_PROC) LEARN_BCAST=-DLEARN_BCAST stub )
+
+# compile these without high optimization to speed compile
+solve_interface.o : solve_interface.F
+shift_domain_em.o : shift_domain_em.F
+module_io_mm5.o : module_io_mm5.F
+module_si_io.o : module_si_io.F
+module_io_wrf.o : module_io_wrf.F
+module_domain.o : module_domain.F
+module_start.o : module_start.F
+module_initialize.o : module_initialize.F
+module_initialize_b_wave.o : module_initialize_b_wave.F
+module_initialize_hill2d_x.o : module_initialize_hill2d_x.F
+module_initialize_quarter_ss.o : module_initialize_quarter_ss.F
+module_initialize_squall2d_x.o : module_initialize_squall2d_x.F
+module_initialize_squall2d_y.o : module_initialize_squall2d_y.F
+module_initialize_real.o : module_initialize_real.F
+module_dm.o : module_dm.F
+start_domain.o : start_domain.F
+mediation_integrate.o : mediation_integrate.F
+module_configure.o : module_configure.F
+
+solve_interface.o                                              \
+shift_domain_em.o                                              \
+module_io_mm5.o module_si_io.o module_io_wrf.o module_domain.o \
+module_start.o module_initialize.o module_initialize_b_wave.o  \
+module_initialize_hill2d_x.o module_initialize_quarter_ss.o    \
+module_initialize_squall2d_x.o module_initialize_squall2d_y.o  \
+module_configure.o                                             \
+mediation_integrate.o \
+module_initialize_real.o module_dm.o start_domain.o :
+	$(RM) $@
+	$(SED_FTN) $*.F > $*.b 
+	$(CPP) -I../inc $(CPPFLAGS) $*.b  > $*.f90
+	$(RM) $*.b
+	$(FC) -c $(PROMOTION) $(FCBASEOPTS) $(MODULE_DIRS) $*.f90
+
+
+###########################################################
+
+# These sub-directory builds are identical across all architectures
+wrfio_esmf :
+	( cd ../external/io_esmf ; \
+	make FC="$(FC) $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS) $(ESMF_MOD_INC)" RANLIB="$(RANLIB)" CPP="$(CPP) $(POUND_DEF) " )
+
+
+#
+#	Macros, these should be generic for all machines
+
+LN	=	ln -sf
+MAKE	=	make -i -r
+RM	= 	rm -f
+
+#	There is probably no reason to modify these rules
+
+wrf_ioapi_includes :
+	( cd ../external/ioapi_share ; \
+	$(MAKE) NATIVE_RWORDSIZE="$(NATIVE_RWORDSIZE)" RWORDSIZE="$(RWORDSIZE)" )
+
+.F.i:
+	$(RM) $@
+	$(CPP) -I../inc $(CPPFLAGS) $*.F > $@
+	mv $*.i $(DEVTOP)/pick/$*.f90
+	cp $*.F $(DEVTOP)/pick
+
+.F.o:
+	$(RM) $@
+	$(SED_FTN) $*.F > $*.b 
+	$(CPP) -I../inc $(CPPFLAGS) $*.b  > $*.f90
+	$(RM) $*.b
+	$(FC) -c $(FCFLAGS) $(MODULE_DIRS) $(PROMOTION) $(FCSUFFIX) $*.f90
+#PGI	if [ ! -e $@ ] ; then \
+#PGI	sleep 10 ; $(FC) -c $(FCFLAGS) $(MODULE_DIRS) $(PROMOTION) $(FCSUFFIX) $*.f90 ; \
+#PGI	fi
+#PGI	if [ ! -e $@ ] ; then \
+#PGI	sleep 10 ; $(FC) -c $(FCFLAGS) $(MODULE_DIRS) $(PROMOTION) $(FCSUFFIX) $*.f90 ; \
+#PGI	fi
+
+.F.f90:
+	$(RM) $@
+	$(SED_FTN) $*.F > $*.b 
+	$(CPP) -I../inc $(CPPFLAGS) $*.b  > $@
+	$(RM) $*.b
+
+.f90.o:
+	$(RM) $@
+	$(FC) -c $(FCFLAGS) $(PROMOTION) $(FCSUFFIX) $*.f90
+#PGI	if [ ! -e $@ ] ; then \
+#PGI	sleep 10 ; $(FC) -c $(FCFLAGS) $(PROMOTION) $(FCSUFFIX) $*.f90 ; \
+#PGI	fi
+#PGI	if [ ! -e $@ ] ; then \
+#PGI	sleep 10 ; $(FC) -c $(FCFLAGS) $(PROMOTION) $(FCSUFFIX) $*.f90 ; \
+#PGI	fi
+
+.c.o:
+	$(RM) $@
+	$(CC) -c $(CFLAGS) $*.c
+
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/README
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/README	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/README	(revision 31)
@@ -0,0 +1,29 @@
+> tar xvf LMDZ.MARS.source.tar
+ajouter perso
+change resol in compile
+check NETCDF NCDFLIB NCDFINC
+> compile
+> mkdir myGCM
+> cd myGCM
+lier z2sig.def (voir what_I_need)
+lier un start_archive avec bon nombre traceurs
+> ln -sf ../gcm.e
+> ln -sf ../newstart.e
+> ln -sf ../deftank/run0 
+lier un callphys.def
+lier un run.def
+cp num_run
+> newstart.e
+> mv restart.nc start0.nc
+> mv restartfi.nc startfi0.nc
+LANCER AVEC
+> nohup run_mcd_3days > log_run_mcd_3days &
+ATTENDRE LONGTEMPS ...
+> mv tmpsave ../STARTBASE_64_48_32_t2
+> rm -rf start.nc start0.nc startfi.nc startfi0.nc
+
+lier le run.def simu
+modifier physiq.F pour sortir des donnees 
+recompiler
+link diagfi.nc vers un espace temp
+ex: ln -sf $MODELS/TMPDIR/GCMINI/diagfi.nc
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/compile
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/compile	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/compile	(revision 31)
@@ -0,0 +1,56 @@
+#! /bin/bash
+
+## changer eventuellement le callphys.def egalement
+
+##################
+tracers=2
+nx=64
+ny=48
+nz=32
+##################
+tracers=1
+nx=64
+ny=48
+nz=25
+##################
+tracers=1
+nx=64
+ny=48
+nz=18
+##################
+#tracers=1
+#nx=64
+#ny=48
+#nz=32
+##################
+
+
+# GCM environment variables
+#--------------------------
+export LMDGCM=$PWD
+export LIBOGCM=$PWD/libo
+
+echo $LMDGCM
+echo $LIBOGCM
+echo $NCDFLIB
+echo $NCDFINC 
+
+#rm -rf libo
+mkdir libo
+touch libo/log_compile_newstart
+touch libo/log_compile_gcm
+
+echo bricolo
+\rm libf/grid/dimensions.h 2> /dev/null
+
+echo compile newstart
+echo check in $PWD/libo/log_compile_newstart
+makegcm -t ${tracers} -p mars -d ${nx}x${ny}x${nz} newstart | tee libo/log_compile_newstart | grep 'warnings' >> libo/log_compile_newstart
+
+echo compile gcm
+echo check in $PWD/libo/log_compile_gcm
+makegcm -t ${tracers} -p mars -d ${nx}x${ny}x${nz} gcm | tee libo/log_compile_gcm | grep 'warnings' >> libo/log_compile_gcm
+
+echo well, OK.
+
+
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/diff_avec_premiere_version_new
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/diff_avec_premiere_version_new	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/diff_avec_premiere_version_new	(revision 31)
@@ -0,0 +1,63 @@
+21,25d20
+< #tracers=1
+< #nx=64
+< #ny=48
+< #nz=32
+< ##################
+Only in /donnees/aslmd/SVN/trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_new/myGCM: .svn
+diff myGCM/callphys.def /donnees/aslmd/SVN/trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_new/myGCM/callphys.def
+5d4
+< ##true???
+17c16
+< stats =.false.
+---
+> stats =.true.
+116,122d114
+< ## albedoes of CO2 ice caps from TES ?
+< TESicealbedo=.true.
+< ## Coefficient for Northern cap albedoes
+< TESice_Ncoef=1.4
+< ## Coefficient for Southern cap albedoes
+< TESice_Scoef=1.4
+< 
+diff myGCM/launch_gcm /donnees/aslmd/SVN/trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_new/myGCM/launch_gcm
+15,19c15
+< #startbase='../../STARTS/'
+< #startbase='/d2/emlmd/run_LMDZ.MARS.BETA.TASI/tau1/'
+< #startbase='/tmp7/aslmd/flush/'
+< #startbase='/tmp7/aslmd/flush2/'
+< startbase='/d6/vblmd/MERIDIANI_EXOMARS/starts_en_cours/tau05/'
+---
+> startbase='../../STARTS/'
+26,29d21
+< if [[ "$daygcm" == 0 ]]
+< then
+< daygcm=669 
+< fi
+5d4
+< ##true???
+17c16
+< stats =.false.
+---
+> stats =.true.
+116,122d114
+< ## albedoes of CO2 ice caps from TES ?
+< TESicealbedo=.true.
+< ## Coefficient for Northern cap albedoes
+< TESice_Ncoef=1.4
+< ## Coefficient for Southern cap albedoes
+< TESice_Scoef=1.4
+< 
+15,19c15
+< #startbase='../../STARTS/'
+< #startbase='/d2/emlmd/run_LMDZ.MARS.BETA.TASI/tau1/'
+< #startbase='/tmp7/aslmd/flush/'
+< #startbase='/tmp7/aslmd/flush2/'
+< startbase='/d6/vblmd/MERIDIANI_EXOMARS/starts_en_cours/tau05/'
+---
+> startbase='../../STARTS/'
+26,29d21
+< if [[ "$daygcm" == 0 ]]
+< then
+< daygcm=669 
+< fi
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def	(revision 31)
@@ -0,0 +1,123 @@
+##General options
+##~~~~~~~~~~~~~~~
+#Run with or without tracer transport ?
+tracer=.false.
+##true???
+
+#Diurnal cycle ?  if diurnal=False, diurnal averaged solar heating
+diurnal=.true.
+
+#Seasonal cycle ? if season=False, Ls stays constant, to value set in "start"
+season = .true.   
+
+#write some more output on the screen ?
+lwrite = .false.    
+
+#Save statistics in file "stats.nc" ?
+stats =.false.
+
+#Save EOF profiles in file "profiles" for Climate Database?
+calleofdump = .false.
+
+## Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+## ~~~~~~~~~~~~~
+#  =1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario,
+#  =4 Mars Year 24 from TES assimilation (same as =24 for now)
+#  =24 Mars Year 24 from TES assimilation (ie: MCD reference case)
+#  =25 Mars Year 25 from TES assimilation (ie: a year with a global dust storm)
+#  =26 Mars Year 26 from TES assimilation
+iaervar = 24 
+# Dust opacity at 700 Pa (when constant, i.e. for the iaervar=1 case)
+tauvis=0.2
+# Dust vertical distribution: 
+# (=0: old distrib. (Pollack90), =1: top set by "topdustref",
+#  =2: Viking scenario; =3 MGS scenario)
+iddist  = 3
+# Dust top altitude (km). (Matters only if iddist=1)
+topdustref = 55.
+
+## Physical Parameterizations :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call radiative transfer ?
+callrad   = .true.
+# call NLTE radiative schemes ?   matters only if callrad=T
+callnlte = .true.
+# call CO2 NIR absorption ?   matters only if callrad=T
+callnirco2 = .true.
+# call turbulent vertical diffusion ?
+calldifv  = .true.
+# call convective adjustment ?
+calladj  = .true.
+# call CO2 condensation ?
+callcond =.true.
+# call thermal conduction in the soil ?
+callsoil = .true.
+# call Lott's gravity wave/subgrid topography scheme ?
+calllott  = .true.
+
+## Radiative transfer options :
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~
+# the rad.transfer is computed every "iradia" physical timestep
+iradia    = 1
+# Output of the exchange coefficient mattrix ? for diagnostic only
+callg2d  = .false.
+# Rayleigh scattering : (should be .false. for now)
+rayleigh  = .false.
+
+## Tracer (dust water, ice and/or chemical species) options (used if tracer=T):
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# DUST: Transported dust ? (if >0, use 'dustbin' dust bins)
+dustbin    = 0
+# DUST: Radiatively active dust ? (matters if dustbin>0)
+active  = .false.
+# DUST: use mass and number mixing ratios to predict dust size ?
+# (must also have dustbin=1)
+doubleq   = .false.
+# DUST: lifted by GCM surface winds ?
+lifting = .false.
+# DUST: lifted by dust devils ?
+callddevil = .false.
+# DUST: Scavenging by CO2 snowfall ?
+scavenging = .false.
+# DUST/WATERICE: Gravitationnal sedimentation ?
+sedimentation = .false.
+# WATERICE: Radiatively active transported atmospheric water ice ?
+activice   = .false.
+# WATER: Compute water cycle
+water = .false.
+# WATER: current permanent caps at both poles. True IS RECOMMENDED
+#        (with .true., North cap is a source of water and South pole
+#         is a cold trap)
+caps  = .true.
+# PHOTOCHEMISTRY: include chemical species 
+photochem  = .false.
+
+## Thermospheric options (relevant if tracer=T) :
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# call thermosphere ?
+callthermos = .false.
+# WATER: included without cycle (only if water=.false.)
+thermoswater = .false.
+# call thermal conduction ? (only if callthermos=.true.)
+callconduct = .false.
+# call EUV heating ? (only if callthermos=.true.)
+calleuv=.false.
+# call molecular viscosity ? (only if callthermos=.true.)
+callmolvis = .false.
+# call molecular diffusion ? (only if callthermos=.true.)
+callmoldiff = .false.
+# call thermospheric photochemistry ? (only if callthermos=.true.)
+thermochem = .false.
+# date for solar flux calculation: (1985 < date < 2002)
+## (Solar min=1996.4 ave=1993.4 max=1990.6)
+solarcondate = 1993.4
+       
+
+## albedoes of CO2 ice caps from TES ?
+TESicealbedo=.true.
+## Coefficient for Northern cap albedoes
+TESice_Ncoef=1.4
+## Coefficient for Southern cap albedoes
+TESice_Scoef=1.4
+
+
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def.water
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def.water	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/callphys.def.water	(revision 31)
@@ -0,0 +1,94 @@
+General options
+~~~~~~~~~~~~~~~
+tracer    (Run with or without tracer transport ?)
+T
+diurnal   (Diurnal cycle ?  if diurnal=F, diurnal averaged solar heating)
+T
+season    (Seasonal cycle ? if season=F, Ls stays constant like in "start")
+T
+lwrite    (want some more output on the screen ?) 
+F
+stats     (Saving statistics in file "stats" ?)
+F
+calleofdump (Saving EOF profiles in file "profiles" for Climate Database ?)
+F
+Dust scenario. Used if the dust is prescribed (i.e. if tracer=F or active=F)
+~~~~~~~~~~~~~
+iaervar  (=1 Dust opt.deph read in startfi; =2 Viking scenario; =3 MGS scenario)
+4        (=4 Mars Year 24 from TES assimilation) 
+iddist  (Dust vertical distribution: =0: old distrib. (Pollack90) 
+3       (=1: top set by "topdustref"; =2: Viking scenario; =3 MGS scenario )
+topdustref (Dust top altitude (km). Matter only if iddist=1)
+55.
+Physical Parameterizations :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+callrad   (call radiative transfer ?)
+T
+callnlte (call NLTE radiative schemes ?   matter only if callrad=T)
+T
+callnirco2 (call CO2 NIR absorption ?   matter only if callrad=T)
+T
+calldifv  (call turbulent vertical diffusion ?)
+T
+calladj   (call convective adjustment ?)
+T
+callcond  (call CO2 condensation ?)
+T
+callsoil  (call thermal conduction in the soil ?)
+T
+calllott  (call Lott's gravity wave/subgrid topography scheme ?)
+T
+Radiative transfer options :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+iradia    (the rad.transfer is computed every "iradia" physical timestep)
+1
+callg2d   (Output of the exchange coefficient mattrix ? for diagnostic only)
+F
+rayleigh  (Rayleigh scattering : should be =F for now)
+F
+Tracer (dust water, ice and/or chemical species) options (use if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dustbin    (DUST: Transported dust ? (if >0, uses q(1) to q(dustbin))
+0
+active     (DUST: Radiatively active dust ? (uses q(1) to q(dustbin))
+F
+doubleq    (DUST: needs dustbin=1, use mass (q(1)) and number (q(2)) mixing ratio to predict dust size ?)
+F
+lifting    (DUST: lifted by GCM surface winds ?)
+F
+dustdevil  (DUST: lifted by dust devils ?)
+F
+scavenging (DUST: Scavenging by CO2 snowfall ?)
+F
+sedimentation (DUST/WATERICE: Gravitationnal sedimentation ?)
+T              (set to T to simulate water cycle with clouds)                
+iceparty   (WATERICE: includes water ice: q(nqmx-1)). "water" must be T !
+T
+activice   (WATERICE: Radiatively active transported atmospheric water ice ?)
+F
+water      (WATER: Compute water cycle using q(nqmx) )
+T
+caps       (WATER: current permanent caps at both poles. T IS RECOMMENDED
+T          (caps=T means Ncap is a source of water, S pole is a cold trap)
+photochem  (PHOTOCHEMISTRY: chemical species included)
+F
+Thermospheric options (relevant if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+callthermos  (call thermosphere ?)
+F
+thermoswater  (WATER: included without cycle	only if water=F)
+F
+callconduct  (call thermal conduction ?     matter only if callthermos=T)
+F
+calleuv  (call EUV heating ?                matter only if callthermos=T)
+F
+callmolvis  (call molecular viscosity ?     matter only if callthermos=T)
+F
+callmoldiff  (call molecular diffusion ?    matter only if callthermos=T)
+F
+thermochem  (call thermospheric photochemistry ?  matter only if callthermos=T)
+F
+solarcondate (date for solar flux calculation: 1985 < date < 2002))
+1993.4       (Solar min=1996.4 ave=1993.4 max=1990.6)
+
+
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/launch_gcm
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/launch_gcm	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/launch_gcm	(revision 31)
@@ -0,0 +1,73 @@
+#! /bin/bash
+
+##############################
+# use: - launch_gcm          #
+#      - echo 2 | launch_gcm #
+##############################
+
+# Author : A. Spiga - Spring/Winter 2008
+
+################## penser a recompiler si on change ceci
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_32_t2'
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/old_STARTBASE_64_48_25'
+#startbase='/d5/aslmd/LMD_MM_MARS_DATA/STARTBASE_64_48_25'
+#startbase='/tmp7/emlmd/run_LMDZ.MARS.BETA.fit.widor/RUN/'
+#startbase='../../STARTS/'
+#startbase='/d2/emlmd/run_LMDZ.MARS.BETA.TASI/tau1/'
+#startbase='/tmp7/aslmd/flush/'
+#startbase='/tmp7/aslmd/flush2/'
+startbase='/d6/vblmd/MERIDIANI_EXOMARS/starts_en_cours/tau05/'
+##################
+
+#
+# get simulation start day and compute start/startfi required date (simulation day - 1)
+#
+echo '>>>' choose start sol ? ; read daygcm
+if [[ "$daygcm" == 0 ]]
+then
+daygcm=669 
+fi
+daygcm=$(expr ${daygcm} - 1)
+
+#
+# find the correct start/startfi files among 3-day startbase
+#
+ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null
+if [[ "$?" == 0 ]]
+then
+    echo found ! ${daygcm}
+else
+    #echo not found ! ${daygcm}	
+    daygcm=$(expr ${daygcm} - 1)
+    ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null
+    if [[ "$?" == 0 ]]     # probably a 'while' loop would be better if more than 3
+    then
+	echo found ! ${daygcm}
+    else
+	#echo not found ! ${daygcm}    
+	daygcm=$(expr ${daygcm} - 1)
+        ls ${startbase}/start${daygcm}.nc > /dev/null 2> /dev/null	
+        echo found ! ${daygcm}
+    fi	
+fi	
+
+#
+# link start/startfi files
+#
+ln -sf ${startbase}/start${daygcm}.nc start.nc
+ln -sf ${startbase}/startfi${daygcm}.nc startfi.nc
+
+#
+# clean the place
+#
+rm -rf restart.nc 2> /dev/null 
+rm -rf restartfi.nc 2> /dev/null
+rm -rf log_gcm 2> /dev/null
+
+#
+# run !
+#
+echo running the GCM - check in ${PWD}/log_gcm 
+nohup gcm.e > log_gcm
+
+exit
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/num_run.STARTBASE
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/num_run.STARTBASE	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/num_run.STARTBASE	(revision 31)
@@ -0,0 +1,1 @@
+0
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/run.def.64x48x25
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/run.def.64x48x25	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/run.def.64x48x25	(revision 31)
@@ -0,0 +1,112 @@
+#
+#----------------------------------------------------------------------- 
+#Parametres de controle du run:                                          
+#------------------------------                                          
+
+# Nombre de jours d'integration                                         
+     nday=10
+
+# nombre de pas par jour (multiple de iperiod) ( ici pour  dt = 1 min )      
+ day_step = 960
+
+# periode pour le pas Matsuno (en pas)                                  
+  iperiod=5
+
+# periode de sortie des variables de controle (en pas)                  
+  iconser=120
+
+# periode d'ecriture du fichier histoire (en jour)                      
+    iecri=100
+
+# periode de stockage fichier histmoy (en jour)                         
+ periodav=60.
+
+# periode de la dissipation (en pas)                                    
+  idissip=5
+
+# choix de l'operateur de dissipation (star ou  non star )              
+ lstardis=.true.
+
+# avec ou sans coordonnee hybrides
+ hybrid=.true.
+
+# nombre d'iterations de l'operateur de dissipation   gradiv            
+nitergdiv=1
+
+# nombre d'iterations de l'operateur de dissipation  nxgradrot          
+nitergrot=2
+
+# nombre d'iterations de l'operateur de dissipation  divgrad            
+   niterh=2
+
+# temps de dissipation des plus petites long.d ondes pour u,v (gradiv)  
+ tetagdiv= 4000.
+
+# temps de dissipation des plus petites long.d ondes pour u,v(nxgradrot)
+ tetagrot=5000.
+
+# temps de dissipation des plus petites long.d ondes pour  h ( divgrad) 
+ tetatemp=5000.
+
+# coefficient pour gamdissip                                            
+  coefdis=0.
+
+# choix du shema d'integration temporelle (Matsuno ou Matsuno-leapfrog) 
+  purmats=.false.
+
+# avec ou sans physique                                                 
+   physic=.true.
+
+# periode de la physique (en pas)                                       
+  iphysiq=20
+
+# choix d'une grille reguliere                                          
+  grireg=.true.
+
+# frequence (en pas) de l'ecriture du fichier diagfi               
+ ecritphy=40
+
+# longitude en degres du centre du zoom                                 
+   clon=63.
+
+# latitude en degres du centre du zoom                                  
+   clat=0.
+
+# facteur de grossissement du zoom,selon longitude                      
+  grossismx=1.
+
+# facteur de grossissement du zoom ,selon latitude                      
+ grossismy=1.
+
+#  Fonction  f(y)  hyperbolique  si = .true.  , sinon  sinusoidale         
+  fxyhypb=.false.
+
+# extension en longitude  de la zone du zoom  ( fraction de la zone totale)
+   dzoomx= 0.
+
+# extension en latitude de la zone  du zoom  ( fraction de la zone totale)
+   dzoomy=0.
+
+#  raideur du zoom en  X
+    taux=2.
+
+#  raideur du zoom en  Y
+    tauy=2.
+
+#  Fonction  f(y) avec y = Sin(latit.) si = .TRUE. ,  Sinon  y = latit.
+  ysinus= .false.
+
+# Avec sponge layer
+  callsponge  = .true.
+ 
+# Sponge:  mode0(u=v=0), mode1(u=umoy,v=0), mode2(u=umoy,v=vmoy)
+  mode_sponge= 2
+
+# Sponge:  hauteur de sponge (km)
+  hsponge= 90
+
+# Sponge:  tetasponge (secondes)
+  tetasponge = 50000
+
+# some definitions for the physics, in file 'callphys.def'
+INCLUDEDEF=callphys.def
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/run_mcd_3days
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/run_mcd_3days	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/run_mcd_3days	(revision 31)
@@ -0,0 +1,60 @@
+#! /bin/bash
+
+#
+# - please link or copy start0.nc and startfi.nc (or start_archive.nc)
+# - please set 3 in the first parameter of run.def
+# - please set 0 in the first line of num_run 
+# - please set daylim below (669 for the whole year)
+#
+
+##########
+daylim=669
+##########
+
+rm -rf tmpsave
+mkdir tmpsave
+
+runnum=0
+runnum_obsolete=0
+daynum=0
+
+while [ ${daynum} -lt ${daylim} ]  
+do
+	# run GCM during three days
+	echo run starting at day ${daynum}
+	echo '****'
+	time run0 
+
+	# save generated starting files 
+	echo saving files
+	echo '****'
+	runnum=$(expr ${runnum} + 1)
+        save=$(expr ${daynum} + 2)
+	cp start${runnum}.nc tmpsave/start${save}.nc
+        cp startfi${runnum}.nc tmpsave/startfi${save}.nc
+	mv lrun${runnum} tmpsave/lrun_for_start${save}
+
+	# increment day to day+3 for while loop assessement
+	daynum=$(expr ${daynum} + 3)
+
+        # get rid of what was generated two steps before
+	echo cleaning
+	echo '****'
+	if [ ${runnum} -gt 1 ]	#because we want to keep start0.nc and startfi0.nc
+	then
+	runnum_obsolete=$(expr ${runnum} - 1)
+	rm -rf start${runnum_obsolete}.nc
+	rm -rf startfi${runnum_obsolete}.nc
+	rm -rf diagfi${runnum_obsolete}.nc
+	fi
+done
+
+# clean the last files
+rm -rf start${runnum}.nc
+rm -rf startfi${runnum}.nc
+rm -rf diagfi${runnum}.nc
+
+# end
+echo -- end --
+echo please check results in file tmpsave
+echo if you want to keep the results, rename this folder
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/what_I_need.info
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/what_I_need.info	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/myGCM/what_I_need.info	(revision 31)
@@ -0,0 +1,9 @@
+total 24K
+-rw-r--r--  1 aslmd lmdjus 3.4K Mar 11 11:50 callphys.def
+-rw-r--r--  1 aslmd lmdjus 3.4K Mar 11 11:45 callphys.def~
+lrwxrwxrwx  1 aslmd lmdjus   18 Mar 11 11:43 gcm.e -> ../LMDZ.MARS/gcm.e
+lrwxrwxrwx  1 aslmd lmdjus   23 Mar 11 11:43 newstart.e -> ../LMDZ.MARS/newstart.e
+-rw-r--r--  1 aslmd lmdjus 3.2K Mar 11 11:49 run.def
+lrwxrwxrwx  1 aslmd lmdjus   57 Mar 11 11:44 start_archive.nc -> /u/forget/WWW/datagcm/Starts/start_archive_64x48x25_t2.nc
+-rw-r--r--  1 aslmd lmdjus    0 Mar 11 11:53 what_I_need
+lrwxrwxrwx  1 aslmd lmdjus   30 Mar 11 11:52 z2sig.def -> ../LMDZ.MARS/deftank/z2sig.def
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/physiq.F
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/physiq.F	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/physiq.F	(revision 31)
@@ -0,0 +1,1351 @@
+      SUBROUTINE physiq(ngrid,nlayer,nq,
+     $            firstcall,lastcall,
+     $            pday,ptime,ptimestep,
+     $            pplev,pplay,pphi,
+     $            pu,pv,pt,pq,
+     $            pw,
+     $            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
+
+
+      IMPLICIT NONE
+c=======================================================================
+c
+c   subject:
+c   --------
+c
+c   Organisation of the physical parametrisations of the LMD 
+c   martian atmospheric general circulation model.
+c
+c   The GCM can be run without or with tracer transport
+c   depending on the value of Logical "tracer" in file  "callphys.def"
+c   Tracers may be water vapor, ice OR chemical species OR dust particles
+c
+c   SEE comments in initracer.F about numbering of tracer species...
+c
+c   It includes:
+c
+c      1. Initialization:
+c      1.1 First call initializations
+c      1.2 Initialization for every call to physiq
+c      1.2.5 Compute mean mass and cp, R and thermal conduction coeff.
+c      2. Compute radiative transfer tendencies
+c         (longwave and shortwave) for CO2 and dust.
+c      3. Gravity wave and subgrid scale topography drag :
+c      4. Vertical diffusion (turbulent mixing):
+c      5. Convective adjustment
+c      6. Condensation and sublimation of carbon dioxide.
+c      7.  TRACERS :
+c       7a. water and water ice
+c       7b. call for photochemistry when tracers are chemical species
+c       7c. other scheme for tracer (dust) transport (lifting, sedimentation)
+c       7d. updates (CO2 pressure variations, surface budget)
+c      8. Contribution to tendencies due to thermosphere
+c      9. Surface and sub-surface temperature calculations
+c     10. Write outputs :
+c           - "startfi", "histfi" (if it's time)
+c           - Saving statistics (if "callstats = .true.")
+c           - Dumping eof (if "calleofdump = .true.")
+c           - Output any needed variables in "diagfi" 
+c      10. Diagnostic: mass conservation of tracers
+c 
+c   author: 
+c   ------- 
+c           Frederic Hourdin	15/10/93
+c           Francois Forget		1994
+c           Christophe Hourdin	02/1997 
+c           Subroutine completly rewritten by F.Forget (01/2000)
+c           Introduction of the photochemical module: S. Lebonnois (11/2002)
+c           Introduction of the thermosphere module: M. Angelats i Coll (2002)
+c           Water ice clouds: Franck Montmessin (update 06/2003)
+c           
+c
+c
+c   arguments:
+c   ----------
+c
+c   input:
+c   ------
+c    ecri                  period (in dynamical timestep) to write output
+c    ngrid                 Size of the horizontal grid.
+c                          All internal loops are performed on that grid.
+c    nlayer                Number of vertical layers.
+c    nq                    Number of advected fields
+c    firstcall             True at the first call
+c    lastcall              True at the last call
+c    pday                  Number of days counted from the North. Spring
+c                          equinoxe.
+c    ptime                 Universal time (0<ptime<1): ptime=0.5 at 12:00 UT
+c    ptimestep             timestep (s)
+c    pplay(ngrid,nlayer)   Pressure at the middle of the layers (Pa)
+c    pplev(ngrid,nlayer+1) intermediate pressure levels (pa)
+c    pphi(ngrid,nlayer)    Geopotential at the middle of the layers (m2s-2)
+c    pu(ngrid,nlayer)      u component of the wind (ms-1)
+c    pv(ngrid,nlayer)      v component of the wind (ms-1)
+c    pt(ngrid,nlayer)      Temperature (K)
+c    pq(ngrid,nlayer,nq)   Advected fields
+c    pudyn(ngrid,nlayer)    \ 
+c    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
+c    ptdyn(ngrid,nlayer)     / corresponding variables
+c    pqdyn(ngrid,nlayer,nq) /
+c    pw(ngrid,?)           vertical velocity
+c
+c   output:
+c   -------
+c
+c    pdu(ngrid,nlayermx)        \
+c    pdv(ngrid,nlayermx)         \  Temporal derivative of the corresponding
+c    pdt(ngrid,nlayermx)         /  variables due to physical processes.
+c    pdq(ngrid,nlayermx)        /
+c    pdpsrf(ngrid)             /
+c    tracerdyn                 call tracer in dynamical part of GCM ?
+
+c
+c=======================================================================
+c
+c    0.  Declarations :
+c    ------------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comgeomfi.h"
+#include "surfdat.h"
+#include "comdiurn.h"
+#include "callkeys.h"
+#include "comcstfi.h"
+#include "planete.h"
+#include "comsaison.h"
+#include "control.h"
+#include "dimradmars.h"
+#include "comg1d.h"
+#include "tracer.h"
+#include "nlteparams.h"
+
+#include "chimiedata.h"
+#include "watercap.h"
+#include "fisice.h"
+#include "param.h"
+#include "param_v3.h"
+#include "conc.h"
+
+#include "netcdf.inc"
+
+
+
+c Arguments :
+c -----------
+
+c   inputs:
+c   -------
+      INTEGER ngrid,nlayer,nq
+      REAL ptimestep
+      REAL pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
+      REAL pphi(ngridmx,nlayer)
+      REAL pu(ngridmx,nlayer),pv(ngridmx,nlayer)
+      REAL pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
+      REAL pw(ngridmx,nlayer) !Mars pvervel transmit par dyn3d
+      REAL zh(ngridmx,nlayermx)      ! potential temperature (K)
+      LOGICAL firstcall,lastcall
+
+      REAL pday
+      REAL ptime 
+      logical tracerdyn
+
+c   outputs:
+c   --------
+c     physical tendencies
+      REAL pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
+      REAL pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
+      REAL pdpsrf(ngridmx) ! surface pressure tendency
+
+
+c Local saved variables:
+c ----------------------
+c     aerosol (dust or ice) extinction optical depth  at reference wavelength 
+c     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
+c      aerosol optical properties  :
+      REAL aerosol(ngridmx,nlayermx,naerkind) 
+
+      INTEGER day_ini  ! Initial date of the run (sol since Ls=0) 
+      INTEGER icount     ! counter of calls to physiq during the run.
+      REAL tsurf(ngridmx)            ! Surface temperature (K)
+      REAL tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
+      REAL co2ice(ngridmx)           ! co2 ice surface layer (kg.m-2)  
+      REAL albedo(ngridmx,2)         ! Surface albedo in each solar band
+      REAL emis(ngridmx)             ! Thermal IR surface emissivity
+      REAL dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
+      REAL fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
+      REAL fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
+      REAL capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
+      REAL fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
+      REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
+      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy 
+      INTEGER ig_vl1                 ! Grid Point near VL1   (for diagnostic) 
+
+      SAVE day_ini, icount
+      SAVE aerosol, tsurf,tsoil
+      SAVE co2ice,albedo,emis, q2
+      SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
+      SAVE ig_vl1
+
+      REAL stephan   
+      DATA stephan/5.67e-08/  ! Stephan Boltzman constant
+      SAVE stephan
+
+c Local variables :
+c -----------------
+
+
+      CHARACTER*80 fichier 
+      INTEGER l,ig,ierr,igout,iq, tapphys
+      INTEGER iqmin                     ! Used if iceparty engaged
+
+      REAL fluxsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2)
+      REAL fluxsurf_sw(ngridmx,2)    !incident SW (solar) surface flux (W.m-2)
+      REAL fluxtop_lw(ngridmx)       !Outgoing LW (IR) flux to space (W.m-2)
+      REAL fluxtop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
+c     for clear area (uncovered by clouds) :
+      REAL clsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2) 
+      REAL clsurf_sw(ngridmx,2)    !incident SW (solar) surface flux (W.m-2)
+      REAL cltop_lw(ngridmx)       !Outgoing LW (IR) flux to space (W.m-2)
+      REAL cltop_sw(ngridmx,2)     !Outgoing SW (solar) flux to space (W.m-2)
+      REAL tauref(ngridmx)           ! Reference column optical depth at 700 Pa
+                                     ! (used if active=F) 
+      REAL tau(ngridmx,naerkind)     ! Column dust optical depth at each point
+      REAL zls                       !  solar longitude (rad)
+      REAL zday                      ! date (time since Ls=0, in martian days)
+      REAL zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
+      REAL zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
+      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
+
+
+c     Tendancies due to various processes:
+      REAL dqsurf(ngridmx,nqmx)
+      REAL zdtlw(ngridmx,nlayermx)     ! (K/s)
+      REAL zdtsw(ngridmx,nlayermx)     ! (K/s)
+      REAL cldtlw(ngridmx,nlayermx)     ! (K/s) LW heating rate for clear area
+      REAL cldtsw(ngridmx,nlayermx)     ! (K/s) SW heating rate for clear area
+      REAL zdtnirco2(ngridmx,nlayermx) ! (K/s)
+      REAL zdtnlte(ngridmx,nlayermx)   ! (K/s)
+      REAL zdtsurf(ngridmx)            ! (K/s)
+      REAL zdtcloud(ngridmx,nlayermx)
+      REAL zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
+      REAL zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
+      REAL zdhadj(ngridmx,nlayermx)                           ! (K/s)
+      REAL zdtgw(ngridmx,nlayermx)                            ! (K/s)
+      REAL zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
+      REAL zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
+      REAL zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
+
+      REAL zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
+      REAL zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
+      REAL zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
+      REAL zdqadj(ngridmx,nlayermx,nqmx)
+      REAL zdqc(ngridmx,nlayermx,nqmx)
+      REAL zdqcloud(ngridmx,nlayermx,nqmx)
+      REAL zdqscloud(ngridmx,nqmx)
+      REAL zdqchim(ngridmx,nlayermx,nqmx)
+      REAL zdqschim(ngridmx,nqmx)
+
+      REAL zdteuv(ngridmx,nlayermx)    ! (K/s)
+      REAL zdtconduc(ngridmx,nlayermx) ! (K/s)
+      REAL zdumolvis(ngridmx,nlayermx)
+      REAL zdvmolvis(ngridmx,nlayermx)
+      real zdqmoldiff(ngridmx,nlayermx,nqmx)
+
+c     Local variable for local intermediate calcul:
+      REAL zflubid(ngridmx)
+      REAL zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
+      REAL zdum1(ngridmx,nlayermx)
+      REAL zdum2(ngridmx,nlayermx)
+      REAL ztim1,ztim2,ztim3, z1,z2
+      REAL ztime_fin
+      REAL zdh(ngridmx,nlayermx)
+      REAL pclc_min ! minimum of the cloud fraction over the whole domain
+      INTEGER length
+      PARAMETER (length=100)
+
+c local variables only used for diagnostic (output in file "diagfi" or "stats")
+c -----------------------------------------------------------------------------
+      REAL ps(ngridmx), zt(ngridmx,nlayermx)
+      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
+      REAL zq(ngridmx,nlayermx,nqmx)
+      REAL fluxtop_sw_tot(ngridmx), fluxsurf_sw_tot(ngridmx)
+      character*2 str2
+      character*5 str5
+      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
+      real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
+      real qtot1,qtot2 ! total aerosol mass
+      integer igmin, lmin
+      logical tdiag
+
+
+      REAL zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
+      REAL zstress(ngrid), cd
+      real hco2(nqmx),tmean, zlocal(nlayermx)
+      real rho(ngridmx,nlayermx) ! density
+      real vmr(ngridmx,nlayermx) ! volume mixing ratio
+
+      REAL time_phys
+     
+c=======================================================================
+
+c 1. Initialisation:
+c -----------------
+
+c  1.1   Initialisation only at first call
+c  ---------------------------------------
+      IF (firstcall) THEN
+
+c        variables set to 0
+c        ~~~~~~~~~~~~~~~~~~
+         call zerophys(ngrid*nlayer*naerkind,aerosol)
+         call zerophys(ngrid*nlayer,dtrad)
+         call zerophys(ngrid,fluxrad)
+
+c        read startfi 
+c        ~~~~~~~~~~~~
+
+! Read netcdf initial physical parameters.
+         CALL phyetat0 ("startfi.nc",0,0,
+     &         nsoilmx,nq,
+     &         day_ini,time_phys,
+     &         tsurf,tsoil,emis,q2,qsurf,co2ice)
+
+         if (pday.ne.day_ini) then
+           write(*,*) "PHYSIQ: ERROR: bad synchronization between ",
+     &                "physics and dynamics"
+           write(*,*) "dynamics day: ",pday
+           write(*,*) "physics day:  ",day_ini
+           stop
+         endif
+
+         write (*,*) 'In physiq day_ini =', day_ini
+
+c        Initialize albedo and orbital calculation
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         CALL surfini(ngrid,co2ice,qsurf,albedo)
+         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
+
+c        initialize soil 
+c        ~~~~~~~~~~~~~~~
+         IF (callsoil) THEN
+            CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
+     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+         ELSE
+            PRINT*,'WARNING! Thermal conduction in the soil turned off'
+            DO ig=1,ngrid
+               capcal(ig)=1.e5
+               fluxgrd(ig)=0.
+            ENDDO
+         ENDIF
+         icount=1
+
+
+c        initialize tracers
+c        ~~~~~~~~~~~~~~~~~~
+         tracerdyn=tracer
+         IF (tracer) THEN
+            CALL initracer(qsurf,co2ice)
+         ENDIF  ! end tracer
+
+c        Determining gridpoint near Viking Lander 1 (used for diagnostic only)
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if(ngrid.ne.1) then
+           latvl1= 22.27 
+           lonvl1= -47.94 
+           ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
+     &              int(1.5+(lonvl1+180)*iim/360.)
+           write(*,*) 'Viking Lander 1 GCM point: lat,lon',
+     &              lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi
+         end if 
+
+c        Initialize thermospheric parameters
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if (callthermos) call param_read
+
+c        Initialize R and Cp as constant
+
+         if (.not.callthermos .and. .not.photochem) then
+                 do l=1,nlayermx
+                  do ig=1,ngridmx
+                   rnew(ig,l)=r
+                   cpnew(ig,l)=cpp
+                   mmean(ig,l)=mugaz
+                   enddo
+                  enddo  
+         endif         
+                   
+      ENDIF        !  (end of "if firstcall")
+
+c ---------------------------------------------------
+c 1.2   Initializations done at every physical timestep:
+c ---------------------------------------------------
+c
+      IF (ngrid.NE.ngridmx) THEN
+         PRINT*,'STOP in PHYSIQ'
+         PRINT*,'Probleme de dimensions :'
+         PRINT*,'ngrid     = ',ngrid
+         PRINT*,'ngridmx   = ',ngridmx
+         STOP
+      ENDIF
+
+c     Initialize various variables
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      call zerophys(ngrid*nlayer, pdv)
+      call zerophys(ngrid*nlayer, pdu)
+      call zerophys(ngrid*nlayer, pdt)
+      call zerophys(ngrid*nlayer*nq, pdq)
+      call zerophys(ngrid, pdpsrf)
+      call zerophys(ngrid, zflubid)
+      call zerophys(ngrid, zdtsurf)
+      call zerophys(ngrid*nq, dqsurf)
+      igout=ngrid/2+1 
+
+
+      zday=pday+ptime ! compute time, in sols (and fraction thereof)
+
+c     Compute Solar Longitude (Ls) :
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (season) then
+         call solarlong(zday,zls)
+      else
+         call solarlong(float(day_ini),zls)
+      end if
+
+c     Compute geopotential at interlayers
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c     ponderation des altitudes au niveau des couches en dp/p
+
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            zzlay(ig,l)=pphi(ig,l)/g
+         ENDDO
+      ENDDO
+      DO ig=1,ngrid
+         zzlev(ig,1)=0.
+         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+      ENDDO
+      DO l=2,nlayer
+         DO ig=1,ngrid
+            z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l))
+            z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l))
+            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+         ENDDO
+      ENDDO
+
+
+!     Potential temperature calculation not the same in physiq and dynamic
+
+c     Compute potential temperature
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO l=1,nlayer
+         DO ig=1,ngrid 
+            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
+            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
+         ENDDO
+      ENDDO
+
+c-----------------------------------------------------------------------
+c    1.2.5 Compute mean mass, cp, and R
+c    --------------------------------
+
+      if(photochem.or.callthermos) then
+         call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
+      endif
+c-----------------------------------------------------------------------
+c    2. Compute radiative tendencies :
+c------------------------------------
+
+
+      IF (callrad) THEN
+         IF( MOD(icount-1,iradia).EQ.0) THEN
+
+c          Local Solar zenith angle
+c          ~~~~~~~~~~~~~~~~~~~~~~~~
+           CALL orbite(zls,dist_sol,declin)
+
+           IF(diurnal) THEN
+               ztim1=SIN(declin)
+               ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+               ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+
+               CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
+     s         ztim1,ztim2,ztim3, mu0,fract)
+
+           ELSE
+               CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad)
+           ENDIF
+
+c          NLTE cooling from CO2 emission
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+           IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
+
+c          Find number of layers for LTE radiation calculations
+           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
+     &          CALL nlthermeq(ngrid,nlayer,pplev,pplay)
+
+
+c          Atmospheric dust opacity and aerosol distribution:
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           CALL dustopacity(ngrid,nlayer,nq,zday,pplay,pplev,zls,pq,
+     $      tauref,tau,aerosol)
+         
+c          Call main radiative transfer scheme
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c          Transfer through dust and CO2, except NIR CO2 absorption
+
+           CALL callradite(icount,ngrid,nlayer,aerosol,albedo,
+     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
+     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw)
+
+c          CO2 near infrared absorption
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           call zerophys(ngrid*nlayer,zdtnirco2)
+           if (callnirco2) then
+              call nirco2abs (ngrid,nlayer,pplay,dist_sol,
+     .                       mu0,fract,declin, zdtnirco2)
+           endif
+
+c          Radiative flux from the sky absorbed by the surface (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           DO ig=1,ngrid
+               fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
+     $         +fluxsurf_sw(ig,1)*(1.-albedo(ig,1))
+     $         +fluxsurf_sw(ig,2)*(1.-albedo(ig,2))
+           ENDDO
+
+
+c          Net atmospheric radiative heating rate (K.s-1)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           IF(callnlte) THEN
+              CALL blendrad(ngrid, nlayer, pplay,
+     &             zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad)
+           ELSE
+              DO l=1,nlayer
+                 DO ig=1,ngrid
+                    dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
+     &                          +zdtnirco2(ig,l)
+                  ENDDO
+              ENDDO
+           ENDIF
+
+
+
+        ENDIF ! of if(mod(icount-1,iradia).eq.0)
+
+
+
+c    Transformation of the radiative tendencies:
+c    -------------------------------------------
+
+c          Net radiative surface flux (W.m-2)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~
+c
+           DO ig=1,ngrid
+               zplanck(ig)=tsurf(ig)*tsurf(ig)
+               zplanck(ig)=emis(ig)*
+     $         stephan*zplanck(ig)*zplanck(ig)
+               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
+           ENDDO
+
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
+            ENDDO
+         ENDDO
+
+      ENDIF ! of IF (callrad)
+
+c-----------------------------------------------------------------------
+c    3. Gravity wave and subgrid scale topography drag :
+c    -------------------------------------------------
+
+
+      IF(calllott)THEN
+
+        CALL calldrag_noro(ngrid,nlayer,ptimestep,
+     &                 pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw)
+ 
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdugw(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l)
+          ENDDO
+        ENDDO
+      ENDIF
+
+c-----------------------------------------------------------------------
+c    4. Vertical diffusion (turbulent mixing):
+c    -----------------------------------------
+c
+      IF (calldifv) THEN
+
+
+         DO ig=1,ngrid
+            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
+         ENDDO
+
+         CALL zerophys(ngrid*nlayer,zdum1)
+         CALL zerophys(ngrid*nlayer,zdum2)
+         do l=1,nlayer
+            do ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            enddo
+         enddo
+
+c        Calling vdif (Martian version WITH CO2 condensation)
+         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
+     $        ptimestep,capcal,lwrite,
+     $        pplay,pplev,zzlay,zzlev,z0,
+     $        pu,pv,zh,pq,tsurf,emis,qsurf,
+     $        zdum1,zdum2,zdh,pdq,zflubid,
+     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
+     &        zdqdif,zdqsdif)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l)
+               pdu(ig,l)=pdu(ig,l)+zdudif(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l)
+
+               zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only
+
+            ENDDO
+         ENDDO
+
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
+         ENDDO
+
+         if (tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+           DO iq=1, nq
+              DO ig=1,ngrid
+                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+              ENDDO
+           ENDDO
+         end if ! of if (tracer)
+
+      ELSE    
+         DO ig=1,ngrid
+            zdtsurf(ig)=zdtsurf(ig)+
+     s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+         ENDDO
+      ENDIF ! of IF (calldifv)
+
+
+c-----------------------------------------------------------------------
+c   5. Dry convective adjustment:
+c   -----------------------------
+
+      IF(calladj) THEN
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
+            ENDDO
+         ENDDO
+         CALL zerophys(ngrid*nlayer,zduadj)
+         CALL zerophys(ngrid*nlayer,zdvadj)
+         CALL zerophys(ngrid*nlayer,zdhadj)
+
+         CALL convadj(ngrid,nlayer,nq,ptimestep,
+     $                pplay,pplev,zpopsk,
+     $                pu,pv,zh,pq,
+     $                pdu,pdv,zdh,pdq,
+     $                zduadj,zdvadj,zdhadj,
+     $                zdqadj)
+
+         DO l=1,nlayer
+            DO ig=1,ngrid
+               pdu(ig,l)=pdu(ig,l)+zduadj(ig,l)
+               pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l)
+               pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l)
+
+               zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only
+            ENDDO
+         ENDDO
+
+         if(tracer) then 
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                 pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         end if
+      ENDIF ! of IF(calladj)
+
+c-----------------------------------------------------------------------
+c   6. Carbon dioxide condensation-sublimation:
+c   -------------------------------------------
+
+      IF (callcond) THEN
+         CALL newcondens(ngrid,nlayer,nq,ptimestep,
+     $              capcal,pplay,pplev,tsurf,pt,
+     $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
+     $              co2ice,albedo,emis,
+     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
+     $	            fluxsurf_sw) 
+
+         DO l=1,nlayer
+           DO ig=1,ngrid
+             pdt(ig,l)=pdt(ig,l)+zdtc(ig,l)
+             pdv(ig,l)=pdv(ig,l)+zdvc(ig,l)
+             pdu(ig,l)=pdu(ig,l)+zduc(ig,l)
+           ENDDO
+         ENDDO
+         DO ig=1,ngrid
+            zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
+         ENDDO
+
+         IF (tracer) THEN
+           DO iq=1, nq
+            DO l=1,nlayer
+              DO ig=1,ngrid
+                pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) 
+              ENDDO
+            ENDDO
+           ENDDO
+         ENDIF ! of IF (tracer)
+
+      ENDIF  ! of IF (callcond)
+
+c-----------------------------------------------------------------------
+c   7. Specific parameterizations for tracers 
+c:   -----------------------------------------
+
+      if (tracer) then 
+
+c   7a. Water and ice
+c     ---------------
+
+c        ---------------------------------------
+c        Water ice condensation in the atmosphere
+c        ----------------------------------------
+         IF (water) THEN
+           call watercloud(ngrid,nlayer, ptimestep,
+     &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
+     &                pq,pdq,zdqcloud,qsurf,zdqscloud,zdtcloud,
+     &                nq,naerkind,tau,icount,zls)
+
+           if (activice) then
+c Temperature variation due to latent heat release
+           DO l=1,nlayer
+             DO ig=1,ngrid
+               pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
+             ENDDO
+           ENDDO
+           endif
+
+           IF (iceparty) then
+             iqmin=nq-1
+           ELSE
+             iqmin=nq
+           ENDIF
+
+           DO iq=iqmin,nq
+             DO l=1,nlayer
+                DO ig=1,ngrid
+                   pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
+                ENDDO
+             ENDDO
+             DO ig=1,ngrid
+                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
+             ENDDO
+           ENDDO
+
+         END IF  ! of IF (water)
+
+c   7b. Chemical species
+c     ------------------
+
+c        --------------
+c        photochemistry :
+c        --------------
+         IF (photochem .or. thermochem) then
+          call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
+     .      zzlay,zday,pq,pdq,zdqchim,zdqschim,zdqcloud,zdqscloud)
+           
+c Photochemistry includes condensation of H2O2
+
+           do iq=nqchem_min,nq
+            if (noms(iq).eq."h2o2") then
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqcloud(ig,l,iq)
+               ENDDO
+             ENDDO
+            else
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqchim(ig,l,iq)
+               ENDDO
+             ENDDO
+            endif
+           ENDDO
+           do iq=nqchem_min,nq
+            if (noms(iq).eq."h2o2") then
+               DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqscloud(ig,iq)
+               ENDDO
+            else
+               DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqschim(ig,iq)
+               ENDDO
+            endif
+           ENDDO
+
+         END IF  ! of IF (photochem.or.thermochem)
+
+c   7c. Aerosol particles
+c     -------------------
+
+c        ----------
+c        Dust devil :
+c        ----------
+         IF(callddevil) then 
+           call dustdevil(ngrid,nlayer,nq, pplev,pu,pv,pt, tsurf,q2,
+     &                zdqdev,zdqsdev)
+ 
+           if (dustbin.ge.1) then
+              do iq=1,nq
+                 DO l=1,nlayer
+                    DO ig=1,ngrid
+                       pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq)
+                    ENDDO
+                 ENDDO
+              enddo
+              do iq=1,nq
+                 DO ig=1,ngrid
+                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
+                 ENDDO
+              enddo
+           endif  ! of if (dustbin.ge.1)
+
+         END IF ! of IF (callddevil)
+
+c        ------------- 
+c        Sedimentation :   acts also on water ice
+c        ------------- 
+         IF (sedimentation) THEN 
+           call zerophys(ngrid*nlayer*nq, zdqsed)
+           call zerophys(ngrid*nq, zdqssed)
+
+           if(doubleq) then
+            call callsedim2q(ngrid,nlayer, ptimestep,
+     &                pplev,zzlev, pt,
+     &                pq, pdq, zdqsed, zdqssed,nq)
+           else
+            call callsedim(ngrid,nlayer, ptimestep,
+     &                pplev,zzlev, pt,
+     &                pq, pdq, zdqsed, zdqssed,nq)
+           end if
+
+
+           DO iq=1, nq
+             DO l=1,nlayer
+               DO ig=1,ngrid
+                    pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq)
+               ENDDO
+             ENDDO
+           ENDDO
+           DO iq=1, nq
+             DO ig=1,ngrid
+                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
+             ENDDO
+           ENDDO
+         END IF   ! of IF (sedimentation)
+
+c   7d. Updates
+c     ---------
+
+        DO iq=1, nq
+          DO ig=1,ngrid
+
+c       ---------------------------------
+c       Updating tracer budget on surface
+c       ---------------------------------
+            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
+
+          ENDDO  ! (ig)
+        ENDDO    ! (iq)
+
+      endif !  of if (tracer) 
+
+
+c-----------------------------------------------------------------------
+c   8. THERMOSPHERE CALCULATION
+c-----------------------------------------------------------------------
+
+      if (callthermos) then
+        call thermosphere(pplev,pplay,dist_sol,
+     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
+     &     pt,pq,pu,pv,pdt,pdq,
+     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
+c           do iq=nqchem_min,nq
+c           write(*,*) 'thermo iq,pq',iq,pq(690,1,iq)
+c           enddo
+
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l)
+            pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)
+     &                         +zdteuv(ig,l)
+            pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l)
+            pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l)
+            DO iq=1, nq
+              pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq)
+            ENDDO
+          ENDDO
+        ENDDO
+
+
+      endif
+c-----------------------------------------------------------------------
+c   9. Surface  and sub-surface soil temperature
+c-----------------------------------------------------------------------
+c
+c
+c   9.1 Increment Surface temperature:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      DO ig=1,ngrid
+         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
+      ENDDO
+
+c  Prescribe a cold trap at south pole (except at high obliquity !!)
+c  Temperature at the surface is set there to be the temperature
+c  corresponding to equilibrium temperature between phases of CO2
+
+      IF (tracer.AND.water.AND.ngridmx.NE.1) THEN
+         if (caps.and.(obliquit.lt.27.)) then
+           ! NB: Updated surface pressure, at grid point 'ngrid', is
+           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
+           tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*
+     &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
+         endif
+c       -------------------------------------------------------------
+c       Change of surface albedo (set to 0.4) in case of ground frost
+c       everywhere except on the north permanent cap and in regions
+c       covered by dry ice. 
+c              ALWAYS PLACE these lines after newcondens !!!
+c       -------------------------------------------------------------
+         do ig=1,ngrid
+
+c       -------------- Version temporaire fit TES 2008 ------------
+         if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
+              albedo(ig,1)=0.45                                     
+              albedo(ig,2)=0.45                                     
+          endif
+
+c         if (co2ice(ig).eq.0.and.qsurf(ig,nqmx).gt.0.005) then
+c           if (.not.watercaptag(ig)) then
+c             albedo(ig,1)=0.4
+c             albedo(ig,2)=0.4
+c           endif
+c         endif
+c       -------------- version Francois ---------------
+c          if (co2ice(ig).eq.0.and.
+c    &        ((qsurf(ig,nqmx).gt.0.005).or.(watercaptag(ig)))) then
+c              albedo(ig,1)=max(0.4,albedodat(ig))
+c              albedo(ig,2)=albedo(ig,1)
+c          endif
+c       ---------------------------------------------
+         enddo  ! of do ig=1,ngrid
+      ENDIF  ! of IF (tracer.AND.water.AND.ngridmx.NE.1)
+
+c
+c   9.2 Compute soil temperatures and subsurface heat flux:
+c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (callsoil) THEN
+         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
+     &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+      ENDIF
+
+c-----------------------------------------------------------------------
+c  10. Write output files
+c  ----------------------
+
+c    -------------------------------
+c    Dynamical fields incrementation
+c    -------------------------------
+c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics)
+      ! temperature, zonal and meridional wind
+      DO l=1,nlayer
+        DO ig=1,ngrid
+          zt(ig,l)=pt(ig,l)  + pdt(ig,l)*ptimestep
+          zu(ig,l)=pu(ig,l)  + pdu(ig,l)*ptimestep
+          zv(ig,l)=pv(ig,l)  + pdv(ig,l)*ptimestep
+        ENDDO
+      ENDDO
+
+      ! tracers
+      DO iq=1, nq
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep
+          ENDDO
+        ENDDO
+      ENDDO
+
+      ! surface pressure
+      DO ig=1,ngrid
+          ps(ig)=pplev(ig,1) + pdpsrf(ig)*ptimestep
+      ENDDO
+
+      ! pressure
+      DO l=1,nlayer
+        DO ig=1,ngrid
+             zplev(ig,l)=pplev(ig,l)/pplev(ig,1)*ps(ig)
+             zplay(ig,l)=pplay(ig,l)/pplev(ig,1)*ps(ig)
+        ENDDO
+      ENDDO
+
+      ! Density 
+      DO l=1,nlayer
+         DO ig=1,ngrid
+            rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l))
+         ENDDO
+      ENDDO
+
+c    Compute surface stress : (NB: z0 is a common in planete.h)
+c     DO ig=1,ngrid
+c        cd = (0.4/log(zzlay(ig,1)/z0))**2
+c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
+c     ENDDO
+
+c     Sum of fluxes in solar spectral bands (for output only)
+      DO ig=1,ngrid
+	     fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2)
+	     fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2)
+      ENDDO
+c ******* TEST ******************************************************
+      ztim1 = 999
+      DO l=1,nlayer
+        DO ig=1,ngrid
+           if (pt(ig,l).lt.ztim1) then
+               ztim1 = pt(ig,l)
+               igmin = ig
+               lmin = l 
+           end if
+        ENDDO
+      ENDDO
+      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then	
+        write(*,*) 'PHYSIQ: stability WARNING :'
+        write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin),
+     &              'ig l =', igmin, lmin
+      end if
+c *******************************************************************
+
+c     ---------------------
+c     Outputs to the screen 
+c     ---------------------
+
+      IF (lwrite) THEN
+         PRINT*,'Global diagnostics for the physics'
+         PRINT*,'Variables and their increments x and dx/dt * dt'
+         WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
+         WRITE(*,'(2f10.5,2f15.5)')
+     s   tsurf(igout),zdtsurf(igout)*ptimestep,
+     s   pplev(igout,1),pdpsrf(igout)*ptimestep
+         WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
+         WRITE(*,'(i4,6f10.5)') (l,
+     s   pu(igout,l),pdu(igout,l)*ptimestep,
+     s   pv(igout,l),pdv(igout,l)*ptimestep,
+     s   pt(igout,l),pdt(igout,l)*ptimestep,
+     s   l=1,nlayer)
+      ENDIF ! of IF (lwrite)
+
+
+      IF (ngrid.NE.1) THEN
+         print*,'Ls =',zls*180./pi,
+     &   ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
+     &   ' tau(Viking1) =',tau(ig_vl1,1)
+
+c        -------------------------------------------------------------------
+c        Writing NetCDF file  "RESTARTFI" at the end of the run
+c        -------------------------------------------------------------------
+c        Note: 'restartfi' is stored just before dynamics are stored
+c              in 'restart'. Between now and the writting of 'restart',
+c              there will have been the itau=itau+1 instruction and
+c              a reset of 'time' (lastacll = .true. when itau+1= itaufin)
+c              thus we store for time=time+dtvr
+
+         IF(lastcall) THEN
+            ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) 
+            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+            call physdem1("restartfi.nc",long,lati,nsoilmx,nq,
+     .              ptimestep,pday,
+     .              ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf,
+     .              area,albedodat,inertiedat,zmea,zstd,zsig,
+     .              zgam,zthe)
+         ENDIF
+
+c        -----------------------------------------------------------------
+c        Saving statistics :
+c        -----------------------------------------------------------------
+c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+c        which can later be used to make the statistic files of the run:
+c        "stats")          only possible in 3D runs !
+
+         
+         IF (callstats) THEN
+
+
+            call wstats(ngrid,"ps","Surface pressure","K",2,ps)
+            call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+            call wstats(ngrid,"co2ice","CO2 ice cover",
+     .                  "kg.m-2",2,co2ice)
+c            call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
+c    .                  emis)
+             call wstats(ngrid,"fluxsurf_lw",
+     .                 "Thermal IR radiative flux to surface","W.m-2",2,
+     .                 fluxsurf_lw)
+             call wstats(ngrid,"fluxsurf_sw",
+     .                  "Solar radiative flux to surface","W.m-2",2,
+     .                   fluxsurf_sw_tot)
+             call wstats(ngrid,"fluxtop_lw",
+     .                  "Thermal IR radiative flux to space","W.m-2",2,
+     .                  fluxtop_lw)
+             call wstats(ngrid,"fluxtop_sw",
+     .                  "Solar radiative flux to space","W.m-2",2,
+     .                  fluxtop_sw_tot)
+            call wstats(ngrid,"dod","Dust optical depth"," ",2,tau)
+
+            call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
+            call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
+            call wstats(ngrid,"v","Meridional (North-South) wind",
+     .                  "m.s-1",3,zv)
+            call wstats(ngrid,"w","Vertical (down-up) wind",
+     .                  "m.s-1",3,pw)
+            call wstats(ngrid,"rho","Atmospheric density","none",3,rho)
+             call wstats(ngrid,"q2",
+     .               "Boundary layer eddy kinetic energy","m2.s-2",3,q2)
+
+           if (tracer) then
+            if (water) then
+               vmr=zq(1:ngridmx,1:nlayermx,nqmx)*mugaz/mmol(nqmx)
+               call wstats(ngrid,"vmr_h2ovapor",
+     .                    "H2O vapor volume mixing ratio","mol/mol",
+     .                    3,vmr)
+               if (iceparty) then
+                  vmr=zq(1:ngridmx,1:nlayermx,nqmx-1)*mugaz/mmol(nqmx-1)
+                  call wstats(ngrid,"vmr_h2oice",
+     .                       "H2O ice volume mixing ratio","mol/mol",
+     .                       3,vmr)
+               endif
+            endif
+
+             if (thermochem.or.photochem) then
+                do iq=1,nq
+                   if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or.
+     .               (noms(iq).eq."co").or.(noms(iq).eq."n2").or.
+     .               (noms(iq).eq."h2").or.
+     .               (noms(iq).eq."o3")) then
+                        do l=1,nlayer
+                          do ig=1,ngrid
+                            vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
+                          end do
+                        end do
+                        call wstats(ngrid,"vmr_"//trim(noms(iq)),
+     .                     "Volume mixing ratio","mol/mol",3,vmr)
+                   endif
+                enddo
+             endif
+           endif !tracer 
+
+            IF(lastcall) THEN
+              write (*,*) "Writing stats..."
+              call mkstats(ierr)
+            ENDIF
+          ENDIF !if callstats
+
+c        (Store EOF for Mars Climate database software)
+         IF (calleofdump) THEN
+            CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
+         ENDIF
+
+
+c        ----------------------------------------------------------------------
+c        output in netcdf file "DIAGFI", containing any variable for diagnostic
+c        (output with  period "ecritphy", set in "run.def")
+c        ----------------------------------------------------------------------
+c        WRITEDIAGFI can ALSO be called from any other subroutines
+c        for any variables !!
+       call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     .                  emis)
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     .                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     .                  co2ice)
+c       call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
+c    .       fluxsurf_lw)
+c       call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
+c    .       fluxsurf_sw_tot)
+c       call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
+c    .       fluxtop_lw)
+c       call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
+c    .       fluxtop_sw_tot)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+c        call WRITEDIAGFI(ngrid,"tau","tau"," ",2,tau)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
+c       call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
+c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
+c       call WRITEDIAGFI(ngridm,'Teta','T potentielle','K',3,zh)
+c       call WRITEDIAGFI(ngridm,'Pression','Pression','Pa',3,pplay)
+        call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil)
+
+
+c       OUTPUT of tracer mass mixing ratio and surface value :    
+       if (tracer) then
+c         (for photochemistry, outputs are done in calchim)
+          do iq=1,nqmx
+            write(str2(1:2),'(i2.2)') iq
+            call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq),
+     &                         'kg/kg',3,zq(1,1,iq))
+            call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
+     &                     'kg.m-2',2,qsurf(1,iq))
+          end do
+        end if
+c ***************************************************
+
+c  Outputs for H2O
+       if (tracer) then
+        if (activice) then
+c         call WRITEDIAGFI(ngridmx,'tauice','tau','SI',2,tau(1,2))
+c         call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate',
+c     &                   'w.m-2',3,zdtsw)
+c         call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
+c     &                   'w.m-2',3,zdtlw)
+        endif  !(activice)
+
+        if (water.and..not.photochem) then
+           iq=nq
+c          write(str2(1:2),'(i2.2)') iq
+c          call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud',
+c     &                    'kg.m-2',2,zdqscloud(1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim',
+c     &                    'kg/kg',3,zdqchim(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif',
+c     &                    'kg/kg',3,zdqdif(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj',
+c     &                    'kg/kg',3,zdqadj(1,1,iq))
+c          call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c',
+c     &                    'kg/kg',3,zdqc(1,1,iq))
+        endif  !(water.and..not.photochem)
+
+c       if (iceparty) then
+c          iq=nq-1
+c          write(str2(1:2),'(i2.2)') iq
+c          call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
+c    &                     'kg/kg',3,zq(1,1,iq))
+c       endif  !(iceparty)
+       endif
+
+c  Outputs for dust tracers
+
+         if (tracer.and.(dustbin.ne.0)) then
+            call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
+            if (doubleq) then
+               call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
+     &                         'kg.m-2',2,qsurf(1,1))
+               call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
+     &                         'N.m-2',2,qsurf(1,2))
+               call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
+     &                         'kg.m-2.s-1',2,zdqsdev(1,1))
+               call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
+     &                         'kg.m-2.s-1',2,zdqssed(1,1))
+               do  l=1,nlayer
+                  do ig=1, ngrid
+                     reff(ig,l)= ref_r0 *
+     &               (r3n_q*pq(ig,l,1)/max(pq(ig,l,2),0.01))**(1./3.)
+                     reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6)
+                  end do
+               end do
+               call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff)
+            else
+               do iq=1,dustbin
+                  write(str2(1:2),'(i2.2)') iq
+                  call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio',
+     &                            'kg/kg',3,zq(1,1,iq))
+                  call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf',
+     &                            'kg.m-2',2,qsurf(1,iq))
+               end do
+            endif ! (doubleq)
+         end if  ! (tracer.and.(dustbin.ne.0))
+
+
+      ELSE     ! if(ngrid.eq.1)
+
+         print*,'Ls =',zls*180./pi,
+     &  '  tauref(700 Pa) =',tauref,' local tau =',tau(1,1)
+c      ----------------------------------------------------------------------
+c      Output in grads file "g1d" (ONLY when using testphys1d)
+c      (output at every X physical timestep)
+c      ----------------------------------------------------------------------
+c
+c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
+c        CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
+         CALL writeg1d(ngrid,1,ps,'ps','Pa')
+         
+         CALL writeg1d(ngrid,nlayer,zt,'T','K')
+c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
+c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
+
+         if(tracer) then
+c           CALL writeg1d(ngrid,1,tau,'tau','SI')
+            do iq=1,nq
+              CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 
+            end do
+         end if
+
+         zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g
+
+         do l=2,nlayer-1
+            tmean=zt(1,l)
+            if(zt(1,l).ne.zt(1,l-1))
+     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
+              zlocal(l)= zlocal(l-1)
+     &        -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g
+         enddo
+         zlocal(nlayer)= zlocal(nlayer-1)-
+     &                   log(pplay(1,nlayer)/pplay(1,nlayer-1))*
+     &                   rnew(1,nlayer)*tmean/g
+
+c         if(tracer) then
+c         do l=2,nlayer
+c            do iq=1,nq
+c               hco2(iq)=zq(1,l,iq)/zq(1,l-1,iq)
+c               hco2(iq)=-(zlocal(l)-zlocal(l-1))/log(hco2(iq))/1000.
+c            enddo
+c            write(225,*) l,pt(1,l),(hco2(iq),iq=1,6)
+c            write(226,*) l,(hco2(iq),iq=7,13)
+c         enddo
+c         endif
+
+      END IF       ! if(ngrid.ne.1)
+
+
+      icount=icount+1
+      RETURN
+      END
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/writediagfi.info
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/writediagfi.info	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/in_lmdz_mars_newphys/writediagfi.info	(revision 31)
@@ -0,0 +1,21 @@
+        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     .                  emis)
+        call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
+     .                  tsurf)
+        call WRITEDIAGFI(ngrid,"ps","surface pressure","K",2,ps)
+        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+     .                  co2ice)
+        call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+        call WRITEDIAGFI(ngrid,"tsoil","soil temperature","K",3,tsoil)
+        if (tracer) then
+            do iq=1,nqmx
+              write(str2(1:2),'(i2.2)') iq
+              call WRITEDIAGFI(ngridmx,'q'//str2,noms(iq),
+     &                          'kg/kg',3,zq(1,1,iq))
+              call WRITEDIAGFI(ngridmx,'qsurf'//str2,noms(iq),
+     &                        'kg.m-2',2,qsurf(1,iq))
+            end do
+          end if
+
Index: /trunk/mesoscale/LMD_MM_MARS/SIMU/run.def
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SIMU/run.def	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SIMU/run.def	(revision 31)
@@ -0,0 +1,2 @@
+# some definitions for the physics, in file 'callphys.def'
+INCLUDEDEF=callphys.def
Index: /trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo_newphys.F90
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo_newphys.F90	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SRC/PREP_MARS/readmeteo_newphys.F90	(revision 31)
@@ -0,0 +1,1380 @@
+program readmeteo
+
+implicit none
+include "netcdf.inc"
+
+!------------------------------------------------------------------------!
+! Readmeteo prepares an initial state for the WPS pre-processor of WRF   !
+!                                                                        !
+! Input is a diagfi.nc NETCDF file from a LMD/GCM simulation             !
+!                                                                        !
+! Outputs are WRFSI intermediate format files ready for metgrid.exe      !
+!       http://www.mmm.ucar.edu/wrf/OnLineTutorial/WPS/IM_si.htm         !
+!                                                                        !
+! **** Please create a WPSFEED folder (or a symbolic link) ****          !       
+!                                                                        !
+! A. Spiga - 16/04/2007                                                  !
+!	     22/05/2007 : need to run zrecast as a preliminary           !
+!            07/06/2007 : external parameter file 'readmeteo.def         !
+!	     30/07/2007 : manage additional variables (tsoil, emiss,...) !	
+!            19/10/2007 : no need to run zrecast at all (eta levels)     !   
+!		12/2007 : include co2ice and emissivity			 ! 
+!               02/2008 : include water vapor and ice                    !
+!               01/2010 : possible use of diagsoil for new physics       ! 
+!									 !
+! spiga@lmd.jussieu.fr							 !	
+!------------------------------------------------------------------------!
+
+
+
+!***************************************************************************
+!***************************************************************************
+REAL, PARAMETER :: emiss_prescribed=0.95
+REAL, PARAMETER :: co2ice_prescribed=0.
+CHARACTER (LEN=3), PARAMETER :: ident='LMD'
+
+
+
+!***************************************************************************
+!***************************************************************************
+
+REAL :: ptop
+REAL, PARAMETER :: grav=3.72
+LOGICAL, PARAMETER :: blank=.false.
+
+
+
+!! Variables to be written in the output file
+!! *** NB: conformity with metgrid/src/read_met_module.F90
+CHARACTER (LEN=33) :: output
+INTEGER :: IFV
+CHARACTER (LEN=24) :: HDATE
+REAL :: XFCST
+CHARACTER (LEN=32) :: SOURCE
+CHARACTER (LEN=9) :: FIELD
+CHARACTER (LEN=25) :: UNITS
+CHARACTER (LEN=46) :: DESC
+REAL :: XLVL
+INTEGER :: NX,NY,IPROJ
+CHARACTER (LEN=8) :: STARTLOC
+REAL :: STARTLAT,STARTLON,DELTALAT,DELTALON
+REAL, POINTER, DIMENSION(:,:) :: SLAB
+
+!! Variables related to NETCDF format
+integer :: nid,nid2,nid3,nvarid,ierr,ierr2
+integer :: timedim,londim,latdim,altdim,altdim2
+integer :: rlatvdim,rlonudim
+integer :: timelen,lonlen,latlen,altlen,altlen2
+integer :: rlatvlen,rlonulen
+integer :: soillen,soildim
+integer :: nphys
+integer, dimension(4) :: corner,edges
+
+!! Intermediate data arrays
+integer :: k,l,m,n,p
+real, dimension(:), allocatable :: lat,lon,time,alt,aps,bps,levels
+real, dimension(:,:), allocatable :: vide,ones,ghtsfile
+real, dimension(:,:), allocatable :: interm
+real, dimension(:,:,:), allocatable :: gwparam
+real, dimension(:,:,:), allocatable :: psfile,tsfile
+real, dimension(:,:,:), allocatable :: emissfile,co2icefile
+real, dimension(:,:,:), allocatable :: tnsfile,unsfile,vnsfile
+real, dimension(:,:,:,:), allocatable :: tfile,ufile,vfile,rfile,hfile
+real, dimension(:,:,:,:), allocatable :: eta_gcm
+real, dimension(:,:,:,:), allocatable :: tfileorig,ufileorig,vfileorig
+real, dimension(:,:,:,:), allocatable :: tsoilfile, dsoilfile, isoilfile
+real, dimension(:,:,:,:), allocatable :: waterfile, watericefile
+
+!! Reading the parameter file 
+integer :: tmp,increment,FILES
+integer :: tmp2,tmp3,tmp4
+character*1 :: answer
+character*4 :: str
+character*2 :: str2, str3, str4
+integer, dimension(:), allocatable :: time_out
+character*13, dimension(:), allocatable :: date_out
+character*19, dimension(:), allocatable :: date_out2
+!***************************************************************************
+!***************************************************************************
+
+
+!!---------------------
+!! Open the datafile
+!!---------------------
+ierr=NF_OPEN ("input_diagfi.nc",NF_NOWRITE,nid)
+IF (ierr.NE.NF_NOERR) THEN
+   write(*,*)'**** Please create a symbolic link called input_diagfi.nc'
+   CALL ABORT
+ENDIF
+
+
+!!--------------------------
+!! Ask for data references
+!!--------------------------
+!write(*,*) "Create n files. What is n ?"
+read(*,*) FILES
+allocate(time_out(FILES))
+allocate(date_out(FILES))
+allocate(date_out2(FILES))
+
+!write(*,*) ""
+!write(*,*) "INPUT: Diagfi time"
+!write(*,*) "list of n subscripts, separated by <Return>s"
+increment=1
+do while (increment.ne.FILES+1)
+    read(*,*) tmp
+    time_out(increment)=tmp
+    increment=increment+1
+enddo
+
+!write(*,*) ""
+!write(*,*) "OUTPUT: WRF time"
+!write(*,*) "fill 4 lines indicating"
+!write(*,*) "-year (4 digits)"
+!write(*,*) "-month (2 digits)" 
+!write(*,*) "-day (2 digits)"
+!write(*,*) "-hour (2 digits)"
+increment=1
+do while (increment.ne.FILES+1)
+    read(*,*) str
+    read(*,*) str2
+    read(*,*) str3
+    read(*,*) str4
+    date_out(increment)=str//'-'//str2//'-'//str3//'_'//str4
+    date_out2(increment)=str//'-'//str2//'-'//str3//'_'//str4//':00:00'
+    !print *,date_out(increment)
+    !write(*,*) "ok? (y/n)"
+    read(*,*) answer
+    if (answer.eq.'n') then
+        !write(*,*) "please write again"
+    else 
+        !write(*,*) "next one, please"    
+        increment=increment+1
+    endif
+enddo
+
+
+!!-------------------
+!! Get array sizes
+!!-------------------
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"Time",timedim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"time",timedim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
+  write(*,*) "timelen: ",timelen
+
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"latitude",latdim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"lat",latdim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
+  write(*,*) "latlen: ",latlen
+
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"longitude",londim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"lon",londim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
+  write(*,*) "lonlen: ",lonlen
+
+SELECT CASE(ident)
+CASE('LMD')
+ierr=NF_INQ_DIMID(nid,"altitude",altdim)
+CASE('OXF')
+ierr=NF_INQ_DIMID(nid,"sigma",altdim)
+END SELECT
+ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
+  write(*,*) "altlen: ",altlen
+
+
+
+!!-------------------------
+!! Allocate local arrays
+!!-------------------------
+allocate(eta_gcm(lonlen,latlen,altlen,timelen))
+allocate(tfile(lonlen,latlen,altlen,timelen))
+allocate(tsoilfile(lonlen,latlen,altlen,timelen))
+allocate(dsoilfile(lonlen,latlen,altlen,timelen))
+allocate(isoilfile(lonlen,latlen,altlen,timelen))
+allocate(tfileorig(lonlen,latlen,altlen,timelen))
+allocate(ufile(lonlen,latlen,altlen,timelen))
+allocate(ufileorig(lonlen,latlen,altlen,timelen))
+allocate(vfile(lonlen,latlen,altlen,timelen))
+allocate(vfileorig(lonlen,latlen,altlen,timelen))
+allocate(rfile(lonlen,latlen,altlen,timelen))
+allocate(hfile(lonlen,latlen,altlen,timelen))  
+allocate(waterfile(lonlen,latlen,altlen,timelen))
+allocate(watericefile(lonlen,latlen,altlen,timelen))
+allocate(psfile(lonlen,latlen,timelen)) 
+allocate(tsfile(lonlen,latlen,timelen))
+allocate(tnsfile(lonlen,latlen,timelen))
+allocate(unsfile(lonlen,latlen,timelen))
+allocate(vnsfile(lonlen,latlen,timelen))
+allocate(emissfile(lonlen,latlen,timelen))
+allocate(co2icefile(lonlen,latlen,timelen))
+allocate(interm(lonlen,latlen))
+allocate(gwparam(lonlen,latlen,5))
+allocate(ghtsfile(lonlen,latlen))  !! no scan axis
+allocate(vide(lonlen,latlen))
+allocate(ones(lonlen,latlen))
+allocate(lat(latlen), lon(lonlen), alt(altlen), time(timelen))
+allocate(aps(altlen),bps(altlen),levels(altlen))
+
+tfile(:,:,:,:)=0
+tsoilfile(:,:,:,:)=0 
+isoilfile(:,:,:,:)=0 
+dsoilfile(:,:,:,:)=0 
+tfileorig(:,:,:,:)=0
+ufileorig(:,:,:,:)=0
+vfileorig(:,:,:,:)=0
+ufile(:,:,:,:)=0 
+vfile(:,:,:,:)=0 
+rfile(:,:,:,:)=0 
+hfile(:,:,:,:)=0
+waterfile(:,:,:,:)=0
+watericefile(:,:,:,:)=0
+psfile(:,:,:)=0 
+tsfile(:,:,:)=0
+emissfile(:,:,:)=0
+co2icefile(:,:,:)=0
+tnsfile(:,:,:)=0
+unsfile(:,:,:)=0
+vnsfile(:,:,:)=0
+interm(:,:)=0
+gwparam(:,:,:)=0
+ghtsfile(:,:)=0
+vide(:,:)=0
+ones(:,:)=0
+
+
+!!-------------------
+!! Read dimensions
+!!-------------------
+
+print *,'>>> Read dimensions !'
+
+print *,'Time'
+SELECT CASE(ident)
+CASE('LMD')
+   ierr = NF_INQ_VARID (nid, "Time",nvarid)
+CASE('OXF')
+   ierr = NF_INQ_VARID (nid, "time",nvarid)
+END SELECT
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <Time> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
+#endif
+   print *,time(1),' ... to ... ',time(timelen)
+
+print *,'Latitude'
+SELECT CASE(ident)
+CASE('LMD')
+   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
+CASE('OXF')
+   ierr = NF_INQ_VARID (nid, "lat",nvarid)
+END SELECT
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <latitude> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
+#endif
+   print *,lat(1),' ... to ... ',lat(latlen),' ... step: ',lat(latlen)-lat(latlen-1)
+
+print *,'Longitude'
+SELECT CASE(ident)
+CASE('LMD')
+   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
+CASE('OXF')
+   ierr = NF_INQ_VARID (nid, "lon",nvarid)
+END SELECT
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <longitude> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
+#endif
+   print *,lon(1),' ... to ... ',lon(lonlen),' ... step: ',lon(lonlen)-lon(lonlen-1)
+
+SELECT CASE(ident)
+CASE('LMD')
+print *, 'Hybrid coordinates'
+   ierr = NF_INQ_VARID (nid, "aps", nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <aps> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, aps)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, aps)
+#endif
+   ierr = NF_INQ_VARID (nid, "bps", nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <bps> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
+#endif
+   print *,aps(1),' ... to ... ',aps(altlen)
+   print *,bps(1),' ... to ... ',bps(altlen)
+CASE('OXF')
+   ierr = NF_INQ_VARID (nid, "sigma", nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <sigma> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, bps)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, bps)
+#endif
+   print *,bps(1),' ... to ... ',bps(altlen)
+END SELECT
+
+
+!!-------------------------------------------
+!! Reading 2D and 3D meteorological fields
+!!-------------------------------------------
+
+IF (blank .EQV. .false.) THEN
+
+print *,'>>> Read 2D optional fields !'
+
+print *,'Emissivity'
+   ierr = NF_INQ_VARID (nid,"emis",nvarid)
+IF (ierr .NE. NF_NOERR) THEN
+        PRINT *, '...warning: not found in diagfi !'
+        PRINT *, '...will be filled with a prescribed value', emiss_prescribed
+        emissfile(:,:,:)=emiss_prescribed
+ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emissfile)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, emissfile)
+#endif
+ENDIF   
+
+print *,'CO2 ice'
+   ierr = NF_INQ_VARID (nid,"co2ice",nvarid)
+IF (ierr .NE. NF_NOERR) THEN
+        PRINT *, '...warning: not found in diagfi !'
+        PRINT *, '...will be filled with a prescribed value', co2ice_prescribed
+        co2icefile(:,:,:)=co2ice_prescribed
+ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2icefile)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, co2icefile)
+#endif
+ENDIF
+
+print *,'>>> Read 2D surface fields !'
+
+print *,'Surface Pressure'
+   ierr = NF_INQ_VARID (nid,"ps",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <ps> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, psfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, psfile)
+#endif
+
+print *,'Ground Temperature'
+   ierr = NF_INQ_VARID (nid,"tsurf",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <tsurf> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tsfile)
+#endif
+
+
+!!"atmospheric" surface temperature is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Temperature'
+   ierr = NF_INQ_VARID (nid,"temp",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      ierr = NF_INQ_VARID (nid,"t",nvarid)
+        IF (ierr .NE. NF_NOERR) THEN
+           PRINT *, "Error: Readmeteo <temp> not found"
+           stop
+        ENDIF
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tfileorig)
+#endif
+   tnsfile=tfileorig(:,:,1,:)
+
+!!"atmospheric" surface u is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Zonal Wind'
+   ierr = NF_INQ_VARID (nid,"u",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+     PRINT *, "Error: Readmeteo <u> not found"
+     stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ufileorig)
+#endif
+   unsfile=ufileorig(:,:,1,:)
+
+!!"atmospheric" surface v is taken
+!!from original diagfi.nc first level
+!!... level is ~3-5 meters
+print *,'Near-Surface Meridional Wind'
+   ierr = NF_INQ_VARID (nid,"v",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+     PRINT *, "Error: Readmeteo <v> not found"
+     stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfileorig)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, vfileorig)
+#endif
+   vnsfile=vfileorig(:,:,1,:)
+
+SELECT CASE(ident)
+CASE('LMD')
+   print *,'Geopotential height at the ground'
+   ierr = NF_INQ_VARID (nid,"phisinit",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <phisinit> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ghtsfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ghtsfile)
+#endif
+!**** from geopotential to geopotential height
+ghtsfile=ghtsfile/grav
+!****
+CASE('OXF')
+!
+! geopotential height ~ altimetry
+!
+print *,'Geopotential height at the ground from file mountain_new.nc'
+ierr=NF_OPEN("mountain_new.nc",NF_NOWRITE,nid3)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "Error: Could not open that file either"
+      stop "Might as well stop here"
+    endif
+ierr=NF_INQ_VARID(nid3,"orography",nvarid)
+ierr=NF_GET_VAR_REAL(nid3,nvarid,ghtsfile)
+if (ierr.ne.NF_NOERR) then
+    stop "Error: Failed reading phisinit"
+endif
+ierr=NF_CLOSE(nid3)
+END SELECT
+
+
+print *,'>>> Read 3D meteorological fields ! - This may take some time ...'
+
+print *,'Temperature'
+   ierr = NF_INQ_VARID (nid,"temp",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      ierr = NF_INQ_VARID (nid,"t",nvarid)
+        IF (ierr .NE. NF_NOERR) THEN
+          PRINT *, "Error: Readmeteo <t> not found"
+          stop
+        ENDIF
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tfile)
+#endif
+
+print *,'Zonal wind'   
+   ierr = NF_INQ_VARID (nid,"u",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <u> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ufile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, ufile)
+#endif
+
+print *,'Meridional wind'
+   ierr = NF_INQ_VARID (nid,"v",nvarid)
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT *, "Error: Readmeteo <v> not found"
+      stop
+   ENDIF
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, vfile)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, vfile)
+#endif
+
+
+!!------------------------
+!! special water stuff
+!!------------------------
+    print *,'Water vapor'
+    ierr=NF_INQ_VARID(nid,"q02",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No q02 - Water vapor set to 0"
+      waterfile(:,:,:,:)=0.  
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,waterfile)
+
+    print *,'Water ice'
+    ierr=NF_INQ_VARID(nid,"q01",nvarid)
+    if (ierr.ne.NF_NOERR) then
+      write(*,*) "...No q01 - Water ice set to 0"  
+      watericefile(:,:,:,:)=0.
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,watericefile)
+!!------------------------
+!! special water stuff
+!!------------------------
+
+
+!SELECT CASE(ident)
+!CASE('LMD')
+
+    print *,'Soil temperatures'
+    ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No tsoil - Soil temperatures set isothermal with surface temperature"
+        DO l=1,altlen
+                tsoilfile(:,:,l,:)=tsfile(:,:,:)
+        ENDDO
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,tsoilfile)
+  
+!!!!!!!!
+!!!!!!!! new physics (but still compatible with old physics)
+
+    print *,'Soil depths'
+    ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No soildepth - Set to default"  !!! see soil_settings in LMD physics
+        DO l=1,altlen
+                dsoilfile(1,1,l,:)=-999.
+        ENDDO
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,dsoilfile(1,1,:,1))
+    print *, dsoilfile(1,1,:,1)
+    DO m=1,lonlen
+    DO n=1,latlen
+    DO p=1,timelen
+      dsoilfile(m,n,:,p) = dsoilfile(1,1,:,1) 
+    ENDDO
+    ENDDO
+    ENDDO
+
+    print *,'Soil thermal inertia'
+    ierr=NF_INQ_VARID(nid,"inertiedat",nvarid)
+    if (ierr.ne.NF_NOERR) then
+        write(*,*) "...No soil therm. inert. - Set to -999"
+        DO l=1,altlen
+                isoilfile(:,:,l,:)=-999.
+        ENDDO
+    endif
+    ierr=NF_GET_VAR_REAL(nid,nvarid,isoilfile)
+
+!!!!!!!!
+!!!!!!!! new physics
+
+
+!END SELECT
+
+ierr=NF_CLOSE(nid)
+
+!!!!lott stuff
+!print *,'get lott'
+!ierr=NF_OPEN("init_lott.nc",NF_NOWRITE,nid3)
+!ierr=NF_INQ_VARID(nid3,"ZMEA",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,1)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZSTD",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,2)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZSIG",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,3)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZGAM",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,4)=interm(:,:)
+!ierr=NF_INQ_VARID(nid3,"ZTHE",nvarid)
+!ierr=NF_GET_VAR_REAL(nid3,nvarid,interm)
+!gwparam(:,:,5)=interm(:,:)
+!ierr=NF_CLOSE(nid3)
+
+ENDIF
+
+
+!!-----------------------------
+!! Loop on the written files
+!! >>> what follows is pretty repetitive ... 
+!!        gotta do something (one more loop?)
+!!-----------------------------
+
+!!! Equivalent eta levels for WRF interpolation in initialize_real
+!print *,'Computing equivalent eta levels'
+!	
+!	!ptop will be passed through RH(surface)
+!	ptop=MINVAL(aps(altlen)+bps(altlen)*psfile(:,:,:))
+!	print *, 'ptop', ptop
+!
+!print *, 'sample: equivalent eta levels at i=1,j=1'
+!DO k = 1,altlen
+!        levels(k)=-k
+!        eta_gcm(:,:,k,:)=(aps(k)+bps(k)*psfile(:,:,:)-ptop)/(psfile(:,:,:)-ptop)
+!        print *, eta_gcm(1,1,k,1) 
+!END DO
+
+!!---better to pass pressure values 
+!!---the eta treatment is done in module_initialize_real
+DO k = 1,altlen
+        levels(k)=-k
+        !!dummy decreasing levels
+END DO
+
+
+
+
+!! Dummy surface level is XLVL=200100.
+
+
+DO l=1,FILES
+
+
+!!---------------------------------------------
+!! Write file in intermediate format for WPS
+!! 1. Surface data
+!!---------------------------------------------
+
+!
+! These variables remain the same in the "loop"
+!
+HDATE=date_out(l)
+IFV=4
+XFCST=0.
+SOURCE=ident
+NX=lonlen
+NY=latlen
+ALLOCATE(SLAB(NX,NY))
+IPROJ=0
+STARTLOC='SWCORNER'
+STARTLAT=lat(1) 
+STARTLON=lon(1)
+DELTALAT=lat(latlen)-lat(latlen-1)
+DELTALON=lon(lonlen)-lon(lonlen-1)
+!
+! Open the data destination file
+!
+output="./WPSFEED/"//ident//":"//date_out2(l)       
+open(UNIT=1,			&
+	FILE=output,		&
+	STATUS='REPLACE',	&
+	FORM='UNFORMATTED',	&
+	ACCESS='SEQUENTIAL')
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='T'
+UNITS='K'
+DESC='Atmospheric temperature'
+XLVL=200100.
+!SLAB=tsfile(:,:,time_out(l))
+!SLAB=tfileorig(:,:,1,time_out(l))
+SLAB=tnsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='U'
+UNITS='m s{-1}'
+DESC='Zonal wind'
+XLVL=200100.
+!SLAB=ufile(:,:,1,time_out(l))
+!SLAB=ufileorig(:,:,1,time_out(l))
+SLAB=unsfile(:,:,time_out(l))
+        ! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='V'
+UNITS='m s{-1}'
+DESC='Meridional wind'
+XLVL=200100.
+!SLAB=vfile(:,:,1,time_out(l))
+!SLAB=vfileorig(:,:,1,time_out(l))
+SLAB=vnsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='RH'
+UNITS=''
+DESC='Customized 2D static field'
+XLVL=200100.
+!SLAB=co2icefile(:,:,time_out(l))
+SLAB=vide(:,:)
+!SLAB=vide(:,:)+ptop
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SOILHGT'
+UNITS='m' 
+DESC='Terrain field of source analysis'  !Ground geopotential height
+XLVL=200100.
+SLAB=ghtsfile(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='PSFC'
+UNITS='Pa'
+DESC='Surface Pressure'
+XLVL=200100.
+SLAB=psfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST000010'
+UNITS=''
+DESC='Emissivity'
+XLVL=200100.
+SLAB=emissfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST010040'
+UNITS=''
+DESC='CO2 ice'
+XLVL=200100.
+SLAB=co2icefile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST040100'
+UNITS=''
+DESC='ZMEA'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,1)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='ST100200'
+UNITS=''
+DESC='ZSTD'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,2)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='LANDSEA'
+UNITS='0/1 Flag'
+DESC='Land/Sea flag'
+XLVL=200100.
+SLAB=ones(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SKINTEMP'
+UNITS='K'
+DESC='Ground temperature'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=tsfile(:,:,time_out(l))
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM000010'
+UNITS=''
+DESC='ZSIG'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,3)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM010040'
+UNITS=''
+DESC='ZGAM'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,4)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM040100'
+UNITS=''
+DESC='ZTHE'
+XLVL=200100.
+!SLAB=vide(:,:)
+SLAB=gwparam(:,:,5)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='SM100200'
+UNITS='fraction'
+DESC='Relative humidity'
+XLVL=200100.
+SLAB=vide(:,:)
+	! And now put everything in the destination file
+	! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON 
+	! ... Data
+	write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HV'
+UNITS='kg/kg'
+DESC='Water vapor'
+XLVL=200100.
+SLAB=waterfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HI'
+UNITS='kg/kg'
+DESC='Water ice'
+XLVL=200100.
+SLAB=watericefile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='TSOIL'
+UNITS='K'
+DESC='Soil temperature'
+XLVL=200100.
+SLAB=tsoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DSOIL'
+UNITS='m'
+DESC='Soil depths'
+XLVL=200100.
+SLAB=dsoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='ISOIL'
+UNITS='tiu'
+DESC='Soil thermal inertia'
+XLVL=200100.
+SLAB=isoilfile(:,:,1,time_out(l))
+        ! And now put everything in the destination file
+        ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+        ! ... Data
+        write(1) SLAB
+!print *,'The field '//DESC//' was written to '//output
+
+
+!!----------------------------------------------------
+!! Write file in intermediate format for WPS
+!! 2. 3D meteorological data
+!! >>> same stuff, but handle vertical profile
+!! >>> !! XLVL must be decreasing (pressure levels) 
+!!----------------------------------------------------
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='T'
+UNITS='K'
+DESC='Atmospheric temperature'
+DO k = 1,altlen
+	XLVL=levels(k)
+	SLAB=tfile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='U'
+UNITS='m s{-1}'
+DESC='Zonal wind'
+DO k = 1,altlen
+        XLVL=levels(k)
+	SLAB=ufile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='V'
+UNITS='m s{-1}'
+DESC='Meridional wind'
+DO k = 1,altlen
+        XLVL=levels(k)
+	SLAB=vfile(:,:,k,time_out(l))
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='RH'
+UNITS=''
+DESC='Customized 2D static field'
+DESC='Eta-levels from the GCM'
+DESC='Pressure values from the GCM'
+DO k = 1,altlen
+        XLVL=levels(k)
+SELECT CASE(ident)
+CASE('LMD')
+        SLAB=aps(k)+bps(k)*psfile(:,:,time_out(l))
+CASE('OXF')
+	SLAB=bps(k)*psfile(:,:,time_out(l))
+END SELECT
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!    
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------! 
+FIELD='HGT'
+UNITS='m'
+DESC='Height'
+DO k = 1,altlen
+        XLVL=levels(k)
+!!*******
+!! PROVISOIRE: normalement, il faudrait la hauteur geopotentielle
+!!*******
+!!however, not used by initialize_real
+SELECT CASE(ident)
+CASE('LMD')
+        SLAB=10000.*alog(610./(aps(k)+bps(k)*psfile(:,:,time_out(l))))
+CASE('OXF')
+	SLAB=10000.*alog(610./(bps(k)*psfile(:,:,time_out(l))))
+END SELECT
+		! And now put everything in the destination file
+		! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON  
+		! ... Data
+		write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HV'
+UNITS='kg/kg'
+DESC='Water vapor'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=waterfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='HI'
+UNITS='kg/kg'
+DESC='Water ice'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=watericefile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='TSOIL'
+UNITS='K'
+DESC='Soil temperature'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=tsoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='DSOIL'
+UNITS='m'
+DESC='Soil depths'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=dsoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+!------------------------!
+! >>> Write a variable   !
+!    ... Copy&Paste part !
+!------------------------!
+FIELD='ISOIL'
+UNITS='tiu'
+DESC='Soil thermal inertia'
+DO k = 1,altlen
+        XLVL=levels(k)
+        SLAB=isoilfile(:,:,k,time_out(l))
+                ! And now put everything in the destination file
+                ! ... Header
+        write(1) IFV
+        write(1) HDATE,XFCST,SOURCE,FIELD,UNITS,DESC,XLVL,NX,NY,IPROJ
+        write(1) STARTLOC,STARTLAT,STARTLON,DELTALAT,DELTALON
+                ! ... Data
+                write(1) SLAB
+END DO
+!print *,'The field '//DESC//' was written to '//output
+
+print *,'****done file '//output, int(100.*float(l)/float(FILES)), ' % '
+close(1)
+
+DEALLOCATE(SLAB)
+
+END DO !! end of the files loop
+
+
+!!-------------------------
+!! DeAllocate local arrays
+!!-------------------------
+deallocate(eta_gcm)
+deallocate(tfile)
+deallocate(tsoilfile)
+deallocate(isoilfile)
+deallocate(dsoilfile)
+deallocate(tfileorig)
+deallocate(ufile)
+deallocate(ufileorig)
+deallocate(vfile)
+deallocate(vfileorig)
+deallocate(rfile)
+deallocate(hfile)  
+deallocate(waterfile)
+deallocate(watericefile)
+deallocate(psfile) 
+deallocate(tsfile)
+deallocate(tnsfile)
+deallocate(unsfile)
+deallocate(vnsfile)
+deallocate(emissfile)
+deallocate(co2icefile)
+deallocate(ghtsfile)	!! no scan axis
+deallocate(vide)
+deallocate(ones)
+deallocate(lat, lon, alt, time)
+deallocate(aps,bps,levels)
+
+
+print *, '------------------------'
+print *, 'End of pre-WPS process !'
+print *, '------------------------'
+print *, 'Here is what you should set in namelist.wps:' 
+print *, " start_date = '"//date_out2(1)//"'"
+print *, " end_date = '"//date_out2(FILES)//"'"
+print *, '------------------------'
+print *, 'Any date between those is ok'
+print *, 'for example, second available date is ... '//date_out2(2)
+print *, '------------------------'
+
+end
Index: /trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SRC/WPS/wps_mars/metgrid/METGRID.TBL.ARW_MarsBase_newphys	(revision 31)
@@ -0,0 +1,562 @@
+========================================
+name=ST
+	z_dim_name=num_st_levels
+        derived=yes
+# IF
+        fill_lev = 1 : ST000010(200100)
+        fill_lev = 2 : ST010040(200100)
+        fill_lev = 3 : ST040100(200100)
+        fill_lev = 4 : ST100200(200100)
+# ELSE IF
+        fill_lev = 1 : ST000010(200100)
+        fill_lev = 2 : ST010200(200100)
+# ELSE
+        fill_lev = 1 : ST000007(200100)
+        fill_lev = 2 : ST007028(200100)
+        fill_lev = 3 : ST028100(200100)
+        fill_lev = 4 : ST100255(200100)
+========================================
+name=SM
+        z_dim_name=num_sm_levels
+        derived=yes
+# IF
+        fill_lev = 1 : SM000010(200100)
+        fill_lev = 2 : SM010040(200100)
+        fill_lev = 3 : SM040100(200100)
+        fill_lev = 4 : SM100200(200100)
+# ELSE IF
+        fill_lev = 1 : SM000010(200100)
+        fill_lev = 2 : SM010200(200100)
+# ELSE
+        fill_lev = 1 : SM000007(200100)
+        fill_lev = 2 : SM007028(200100)
+        fill_lev = 3 : SM028100(200100)
+        fill_lev = 4 : SM100255(200100)
+========================================
+name=SW
+        z_dim_name=num_sw_levels
+        derived=yes
+# IF
+        fill_lev = 1 : SW000010(200100)
+        fill_lev = 2 : SW010040(200100)
+        fill_lev = 3 : SW040100(200100)
+        fill_lev = 4 : SW100200(200100)
+# ELSE IF
+        fill_lev = 1 : SW000010(200100)
+        fill_lev = 2 : SW010200(200100)
+========================================
+name=SOILM
+        z_dim_name=num_soilm_levels
+        derived=yes
+        fill_lev = 1 : SOILM000(200100)      
+        fill_lev = 2 : SOILM005(200100)      
+        fill_lev = 3 : SOILM020(200100)      
+        fill_lev = 4 : SOILM040(200100)      
+        fill_lev = 5 : SOILM160(200100)      
+	fill_lev = 6 : SOILM300(200100)      
+========================================
+name=SOILT
+        z_dim_name=num_soilt_levels
+        derived=yes
+        fill_lev = 1 : SOILT000(200100)      
+        fill_lev = 2 : SOILT005(200100)      
+        fill_lev = 3 : SOILT020(200100)      
+        fill_lev = 4 : SOILT040(200100)      
+        fill_lev = 5 : SOILT160(200100)      
+        fill_lev = 6 : SOILT300(200100)      
+========================================
+name=PRES
+        z_dim_name=num_metgrid_levels
+        derived=yes
+        mandatory=yes    # MUST HAVE THIS FIELD
+        fill_lev=all:PRESSURE
+        fill_lev=200100:PSFC(200100)
+        fill_lev=all:vertical_index; level_template=TT
+========================================
+name=LANDSEA 
+        interp_option=nearest_neighbor
+        fill_missing=-1.
+        fill_lev=200100:LANDMASK(1)
+========================================
+name=SEAICE 
+        interp_option=four_pt+average_4pt
+        interp_mask=LANDSEA(1) 
+        masked=land
+        fill_missing=0.
+========================================
+name=H0ML 
+        interp_option=four_pt+average_4pt
+        interp_mask=LANDSEA(1) 
+        masked=land
+        fill_missing=0.
+========================================
+name=SNOW
+        interp_option=four_pt+average_4pt
+        masked=water 
+        fill_missing=0.
+========================================
+name=SKINTEMP
+#        interp_option=sixteen_pt+four_pt+wt_average_4pt+wt_average_16pt+search
+#        interp_option=four_pt
+#        mandatory=yes    # MUST HAVE THIS FIELD
+## commenter les trois prochaines lignes ?
+#        masked=both
+#        interp_land_mask  = LANDSEA(1)
+#        interp_water_mask = LANDSEA(0)
+#        fill_missing=0.
+        mandatory=yes    # MUST HAVE THIS FIELD
+        interp_option=four_pt+average_4pt
+        #careful: inputs from the GCMs are in a "polygonal" world !
+        #        interp_option=sixteen_pt+four_pt+average_4pt
+        interp_option=sixteen_pt+four_pt
+        fill_missing=-1.0e+30
+        missing_value=-9.99e+33
+========================================
+name=PSFC
+#        interp_option=four_pt+average_4pt
+#        interp_option=sixteen_pt+average_16pt
+#careful: inputs from the GCMs are in a "polygonal" world !
+        interp_option=sixteen_pt+four_pt+average_4pt
+        fill_lev=200100:const(200100.)
+        flag_in_output=FLAG_PSFC
+========================================
+name=VEGCAT
+        interp_option=nearest_neighbor
+        fill_missing=0.
+========================================
+name=CANWAT
+        interp_option=four_pt
+        fill_missing=0.
+========================================
+name=SOIL_CAT
+        interp_option=nearest_neighbor
+        fill_missing=0.
+========================================
+name=SW000010
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SW000010
+========================================
+name=SW010040
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SW010040
+========================================
+name=SW040100
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SW040100
+========================================
+name=SW100200
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SW100200
+========================================
+name=SW010200
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SW010200
+========================================
+name=SM000010
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM000010
+========================================
+name=SM010040
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM010040
+========================================
+name=SM040100
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM040100
+========================================
+name=SM100200
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM100200
+========================================
+name=SM010200
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM010200
+========================================
+# EMISSIVITY
+name=ST000010
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=sixteen_pt+four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=0.
+        flag_in_output=FLAG_ST000010
+========================================
+# CO2 ICE
+name=ST010040
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+#        interp_option=sixteen_pt+four_pt
+        interp_option=four_pt
+###ATTENTION sixteen_pt donne des valeurs negatives
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=0.
+        flag_in_output=FLAG_ST010040
+========================================
+name=ST040100
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_ST040100
+========================================
+name=ST100200
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_ST100200
+========================================
+name=ST010200
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+#        masked=water
+#        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_ST010200
+========================================
+name=SM000007
+#        interp_option=sixteen_pt+four_pt+average_4pt+search
+        interp_option=four_pt
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM000007
+========================================
+name=SM007028
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM007028
+========================================
+name=SM028100
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM028100
+========================================
+name=SM100255
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SM100255
+========================================
+name=ST000007
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_ST000007
+========================================
+name=ST007028
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_ST007028
+========================================
+name=ST028100
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_ST028100
+========================================
+name=ST100255
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_ST100255
+========================================
+name=SOILM000
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SOILM000
+========================================
+name=SOILM005
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SOILM005
+========================================
+name=SOILM020
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SOILM020
+========================================
+name=SOILM040
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SOILM040
+========================================
+name=SOILM160
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SOILM160
+========================================
+name=SOILM300
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=1.
+        flag_in_output=FLAG_SOILM300
+========================================
+name=SOILT000
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_SOILT000
+========================================
+name=SOILT005
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_SOILT005
+========================================
+name=SOILT020
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_SOILT020
+========================================
+name=SOILT040
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_SOILT040
+========================================
+name=SOILT160
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_SOILT160
+========================================
+name=SOILT300
+        interp_option=sixteen_pt+four_pt+average_4pt+search
+        masked=water
+        interp_mask=LANDSEA(0)
+        fill_missing=285.
+        flag_in_output=FLAG_SOILT300
+========================================
+name=PMSL
+        interp_option=sixteen_pt+four_pt+average_4pt
+========================================
+name=RH
+        interp_option=sixteen_pt+four_pt+average_4pt
+        fill_missing=100.
+        fill_lev=200100:const(0.)        
+========================================
+name=TAVGSFC
+        interp_option=four_pt
+        fill_missing=0.
+        fill_lev=200100:TT(200100)
+        flag_in_output=FLAG_TAVGSFC
+========================================
+name=T ; output_name=TT  # If we get T, use entry from TT and
+                         #   write the field out as TT
+========================================
+#average_4pt, average_16pt, wt_average_4pt, wt_average_16pt, 
+#nearest_neighbor, four_pt, sixteen_pt, search, average_gcell(r)
+#
+#wt_average: bof
+#average_gcell: bof aussi
+name=TT
+        mandatory=yes    # MUST HAVE THIS FIELD
+#        interp_option=four_pt+average_4pt
+#careful: inputs from the GCMs are in a "polygonal" world !
+        interp_option=sixteen_pt+four_pt+average_4pt
+#        interp_option=sixteen_pt+four_pt
+        fill_missing=-1.0e+30
+        missing_value=-9.99e+33
+========================================
+name=U ; output_name=UU  # If we get U, use entry from UU and
+                         #   write the field out as UU
+========================================
+name=UU
+        mandatory=yes    # MUST HAVE THIS FIELD
+        interp_option=sixteen_pt+four_pt+average_4pt
+#        interp_option=sixteen_pt+four_pt
+        is_u_field=yes
+        output_stagger=U
+        fill_missing=-1.0e+30
+	missing_value=-9.99e+33
+========================================
+name=V ; output_name=VV  # If we get V, use entry from VV and
+                         #   write the field out as VV
+========================================
+name=VV
+        mandatory=yes    # MUST HAVE THIS FIELD
+        interp_option=sixteen_pt+four_pt+average_4pt
+#        interp_option=sixteen_pt+four_pt
+        is_v_field=yes
+        output_stagger=V
+        fill_missing=-1.0e+30
+        missing_value=-9.99e+33
+========================================
+name=SST
+        interp_option=sixteen_pt+four_pt
+        fill_missing=0. 
+        flag_in_output=FLAG_SST
+========================================
+name=QV
+        interp_option=four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:const(0.)
+        flag_in_output=FLAG_QV
+========================================
+name=QR
+        interp_option=four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:const(0.)
+        flag_in_output=FLAG_QR
+========================================
+name=QC
+        interp_option=four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:const(0.)
+        flag_in_output=FLAG_QC
+========================================
+name=QI
+        interp_option=four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:const(0.)
+        flag_in_output=FLAG_QI
+========================================
+name=QS
+        interp_option=four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:const(0.)
+========================================
+name=QG
+        interp_option=four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:const(0.)
+        flag_in_output=FLAG_QG
+========================================
+name=VPTMP
+        interp_option=sixteen_pt+four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:const(0.)
+========================================
+name=PRESSURE
+        interp_option=sixteen_pt+four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:PSFC(200100)
+========================================
+name=GHT
+        interp_option=sixteen_pt+four_pt+average_4pt
+        fill_missing=0.
+        fill_lev=200100:SOILHGT(200100)
+        fill_lev=200100:HGT_M(1)
+========================================
+name=HGTT
+        output=no
+        interp_option=nearest_neighbor
+========================================
+name=WEASD
+        interp_option=four_pt+average_4pt
+        masked=water 
+        fill_missing=0.
+========================================
+name=SOILHGT
+#        interp_option=four_pt+average_4pt
+#careful: inputs from the GCMs are in a "polygonal" world !
+        interp_option=sixteen_pt+four_pt+average_4pt
+        flag_in_output=FLAG_SOILHGT
+========================================
+name=HV 
+#        mandatory=yes    # MUST HAVE THIS FIELD
+#        interp_option=sixteen_pt+four_pt+average_4pt
+#        interp_option=sixteen_pt+four_pt
+        interp_option=four_pt
+#        fill_lev=200100:const(0.)
+        fill_missing=-1.0e+30
+	missing_value=-9.99e+33
+========================================
+name=HI 
+#        mandatory=yes    # MUST HAVE THIS FIELD
+#        interp_option=sixteen_pt+four_pt+average_4pt
+#        interp_option=sixteen_pt+four_pt
+        interp_option=four_pt        
+#        fill_lev=200100:const(0.)
+        fill_missing=-1.0e+30
+	missing_value=-9.99e+33
+========================================
+name=TSOIL
+#        mandatory=yes    # MUST HAVE THIS FIELD
+#        interp_option=sixteen_pt+four_pt+average_4pt
+         interp_option=sixteen_pt+four_pt
+#        fill_lev=200100:const(0.)
+         fill_missing=-1.0e+30
+         missing_value=-9.99e+33
+========================================
+name=DSOIL
+#        mandatory=yes    # MUST HAVE THIS FIELD
+#        interp_option=sixteen_pt+four_pt+average_4pt
+         interp_option=sixteen_pt+four_pt
+#        fill_lev=200100:const(0.)
+         fill_missing=-1.0e+30
+         missing_value=-9.99e+33
+========================================
+name=ISOIL
+#        mandatory=yes    # MUST HAVE THIS FIELD
+#        interp_option=sixteen_pt+four_pt+average_4pt
+         interp_option=sixteen_pt+four_pt
+#        fill_lev=200100:const(0.)
+         fill_missing=-1.0e+30
+         missing_value=-9.99e+33
+========================================
Index: /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F_gcm
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F_gcm	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F_gcm	(revision 31)
@@ -0,0 +1,1056 @@
+!IDEAL:MODEL_LAYER:INITIALIZATION
+!
+
+!  This MODULE holds the routines which are used to perform various initializations
+!  for the individual domains.  
+
+!  This MODULE CONTAINS the following routines:
+
+!  initialize_field_test - 1. Set different fields to different constant
+!                             values.  This is only a test.  If the correct
+!                             domain is not found (based upon the "id")
+!                             then a fatal error is issued.               
+
+!-----------------------------------------------------------------------
+
+MODULE module_initialize
+
+   USE module_domain
+   USE module_io_domain
+   USE module_state_description
+   USE module_model_constants
+   USE module_bc
+   USE module_timing
+   USE module_configure
+   USE module_init_utilities
+#ifdef DM_PARALLEL
+   USE module_dm
+#endif
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------
+! this is a wrapper for the solver-specific init_domain routines.
+! Also dereferences the grid variables and passes them down as arguments.
+! This is crucial, since the lower level routines may do message passing
+! and this will get fouled up on machines that insist on passing down
+! copies of assumed-shape arrays (by passing down as arguments, the 
+! data are treated as assumed-size -- ie. f77 -- arrays and the copying
+! business is avoided).  Fie on the F90 designers.  Fie and a pox.
+
+   SUBROUTINE init_domain ( grid )
+
+   IMPLICIT NONE
+
+   !  Input data.
+   TYPE (domain), POINTER :: grid 
+   !  Local data.
+   INTEGER                :: dyn_opt 
+   INTEGER :: idum1, idum2
+
+   CALL nl_get_dyn_opt( 1,dyn_opt )
+   
+   CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
+
+   IF (      dyn_opt .eq. 1 &
+        .or. dyn_opt .eq. 2 &
+        .or. dyn_opt .eq. 3 &
+                                       ) THEN
+     CALL init_domain_rk( grid &
+!
+#include <em_actual_new_args.inc>
+!
+                        )
+
+   ELSE
+     WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt
+     CALL wrf_error_fatal ( ' init_domain: unknown or unimplemented dyn_opt ' )
+   ENDIF
+
+   END SUBROUTINE init_domain
+
+!-------------------------------------------------------------------
+
+   SUBROUTINE init_domain_rk ( grid &
+!
+# include <em_dummy_new_args.inc>
+!
+)
+   IMPLICIT NONE
+
+   !  Input data.
+   TYPE (domain), POINTER :: grid
+
+# include <em_dummy_new_decl.inc>
+
+   TYPE (grid_config_rec_type)              :: config_flags
+
+   !  Local data
+   INTEGER                             ::                       &
+                                  ids, ide, jds, jde, kds, kde, &
+                                  ims, ime, jms, jme, kms, kme, &
+                                  its, ite, jts, jte, kts, kte, &
+                                  i, j, k
+
+   ! Local data
+!****Mars
+    REAL    :: x_param,y_param,rho_param,dilat
+    REAL    :: mulu, mulv, addu, addv
+!****Mars
+   INTEGER, PARAMETER :: nl_max = 1000
+   REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
+   INTEGER :: nl_in
+
+
+   INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc
+   REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
+   REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
+!   REAL, EXTERNAL :: interp_0
+   REAL    :: hm, xa
+   REAL    :: pi
+
+!  stuff from original initialization that has been dropped from the Registry 
+   REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
+   REAL    :: qvf1, qvf2, pd_surf
+   INTEGER :: it
+   real :: thtmp, ptmp, temp(3)
+
+   LOGICAL :: moisture_init
+   LOGICAL :: stretch_grid, dry_sounding
+
+  INTEGER :: xs , xe , ys , ye
+  REAL :: mtn_ht
+   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
+
+!!MARS
+ REAL :: lon_input, lat_input, alt_input, tsurf_input
+!!MARS
+
+
+#ifdef DM_PARALLEL
+#    include <em_data_calls.inc>
+#endif
+
+
+   SELECT CASE ( model_data_order )
+         CASE ( DATA_ORDER_ZXY )
+   kds = grid%sd31 ; kde = grid%ed31 ;
+   ids = grid%sd32 ; ide = grid%ed32 ;
+   jds = grid%sd33 ; jde = grid%ed33 ;
+
+   kms = grid%sm31 ; kme = grid%em31 ;
+   ims = grid%sm32 ; ime = grid%em32 ;
+   jms = grid%sm33 ; jme = grid%em33 ;
+
+   kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
+   its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
+   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
+         CASE ( DATA_ORDER_XYZ )
+   ids = grid%sd31 ; ide = grid%ed31 ;
+   jds = grid%sd32 ; jde = grid%ed32 ;
+   kds = grid%sd33 ; kde = grid%ed33 ;
+
+   ims = grid%sm31 ; ime = grid%em31 ;
+   jms = grid%sm32 ; jme = grid%em32 ;
+   kms = grid%sm33 ; kme = grid%em33 ;
+
+   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
+   jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
+   kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
+         CASE ( DATA_ORDER_XZY )
+   ids = grid%sd31 ; ide = grid%ed31 ;
+   kds = grid%sd32 ; kde = grid%ed32 ;
+   jds = grid%sd33 ; jde = grid%ed33 ;
+
+   ims = grid%sm31 ; ime = grid%em31 ;
+   kms = grid%sm32 ; kme = grid%em32 ;
+   jms = grid%sm33 ; jme = grid%em33 ;
+
+   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
+   kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
+   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
+
+   END SELECT
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!MARS : mountain
+!!MARS : mountain ex. hm = 2000. xa = 6.0
+  open(unit=22,file='ze_hill',form='formatted',status='old')
+  rewind(22) 
+  read(22,*) hm, xa
+  write(6,*) 'height, width ', hm, xa
+  close(22)
+!!MARS
+!!MARS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+   icm = ide/2
+!****Mars
+   jcm = jde/2
+!
+!   xa1  = 5000./500.
+!   xal1 = 4000./500.
+!   pii  = 2.*asin(1.0)
+!   hm1  = 250.
+!!   hm1  = 1000.
+!****Mars
+
+   delt = 3.
+!   delt = 10.
+
+!****Mars
+   stretch_grid = .true.
+!   stretch_grid = .false.
+!****Mars
+!   z_scale = .50
+   z_scale = .40
+   pi = 2.*asin(1.0)
+   write(6,*) ' pi is ',pi
+   nxc = (ide-ids)/2
+   nyc = (jde-jds)/2
+
+!!!MARS
+!!!MARS
+!  open(unit=16,file='input_vert',form='formatted',status='old')
+!  rewind(16) 
+!  read(16,*) delt, z_scale
+!  write(6,*) 'delt, z_scale are  ', delt, z_scale
+!  close(16)
+!!!MARS
+!!!MARS
+
+   CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
+
+! here we check to see if the boundary conditions are set properly
+
+   CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
+
+   moisture_init = .true.
+
+    grid%itimestep=0
+
+#ifdef DM_PARALLEL
+   CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
+   CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
+#endif
+
+    CALL nl_set_mminlu(1,'    ')
+    CALL nl_set_iswater(1,0)
+    CALL nl_set_cen_lat(1,40.)
+    CALL nl_set_cen_lon(1,-105.)
+    CALL nl_set_truelat1(1,0.)
+    CALL nl_set_truelat2(1,0.)
+    CALL nl_set_moad_cen_lat (1,0.)
+    CALL nl_set_stand_lon (1,0.)
+    CALL nl_set_map_proj(1,0)
+
+
+!  here we initialize data we currently is not initialized 
+!  in the input data
+
+    DO j = jts, jte
+      DO i = its, ite
+         grid%msft(i,j)     = 1.
+         grid%msfu(i,j)     = 1.
+         grid%msfv(i,j)     = 1.
+         grid%sina(i,j)     = 0.
+         grid%cosa(i,j)     = 1.
+         grid%e(i,j)        = 0.
+         grid%f(i,j)        = 0.         !! MARS: put coriolis here if needed
+
+      END DO
+   END DO
+
+    DO j = jts, jte
+    DO k = kts, kte
+      DO i = its, ite
+         grid%em_ww(i,k,j)     = 0.
+      END DO
+   END DO
+   END DO
+
+   grid%step_number = 0
+
+!! set up the grid
+!
+!   IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
+!     DO k=1, kde
+!      grid%em_znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
+!                                (1.-exp(-1./z_scale))
+!     ENDDO
+!   ELSE
+!     DO k=1, kde
+!      grid%em_znw(k) = 1. - float(k-1)/float(kde-1)
+!     ENDDO
+!   ENDIF
+
+!!MARS
+!!MARS
+  open(unit=12,file='levels',form='formatted',status='old')
+  rewind(12)
+  DO k=1, kde
+  read(12,*) grid%em_znw(k)
+  write(6,*) 'read level ', k,grid%em_znw(k)
+  ENDDO
+  close(12)
+!!MARS
+!!MARS
+   
+
+   DO k=1, kde-1
+    grid%em_dnw(k) = grid%em_znw(k+1) - grid%em_znw(k)
+    grid%em_rdnw(k) = 1./grid%em_dnw(k)
+    grid%em_znu(k) = 0.5*(grid%em_znw(k+1)+grid%em_znw(k))
+   ENDDO
+   DO k=2, kde-1
+    grid%em_dn(k) = 0.5*(grid%em_dnw(k)+grid%em_dnw(k-1))
+    grid%em_rdn(k) = 1./grid%em_dn(k)
+    grid%em_fnp(k) = .5* grid%em_dnw(k  )/grid%em_dn(k)
+    grid%em_fnm(k) = .5* grid%em_dnw(k-1)/grid%em_dn(k)
+   ENDDO
+
+   cof1 = (2.*grid%em_dn(2)+grid%em_dn(3))/(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(2) 
+   cof2 =     grid%em_dn(2)        /(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(3) 
+   grid%cf1  = grid%em_fnp(2) + cof1
+   grid%cf2  = grid%em_fnm(2) - cof1 - cof2
+   grid%cf3  = cof2       
+
+   grid%cfn  = (.5*grid%em_dnw(kde-1)+grid%em_dn(kde-1))/grid%em_dn(kde-1)
+   grid%cfn1 = -.5*grid%em_dnw(kde-1)/grid%em_dn(kde-1)
+   grid%rdx = 1./config_flags%dx
+   grid%rdy = 1./config_flags%dy
+
+!  get the sounding from the ascii sounding file, first get dry sounding and 
+!  calculate base state
+
+     !!!!
+     !!!! user-modified wind speed
+     !!!!
+     mulu = 1.  !! default 
+     mulv = 1.  !! default
+     addu = 0.  !! default
+     addv = 0.  !! default
+     IF (config_flags%init_MU .ne. 0.) mulu = config_flags%init_MU
+     IF (config_flags%init_MV .ne. 0.) mulv = config_flags%init_MV
+     IF (config_flags%init_U  .ne. 0.) addu = config_flags%init_U
+     IF (config_flags%init_V  .ne. 0.) addv = config_flags%init_V
+     write(6,*) ' coeff for winds: ', mulu, mulv, addu, addv
+
+  dry_sounding = .true.
+  IF ( wrf_dm_on_monitor() ) THEN
+  write(6,*) ' getting dry sounding for base state '
+
+  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in, &
+                        mulu, mulv, addu, addv )
+  ENDIF
+  CALL wrf_dm_bcast_real( zk , nl_max )
+  CALL wrf_dm_bcast_real( p_in , nl_max )
+  CALL wrf_dm_bcast_real( pd_in , nl_max )
+  CALL wrf_dm_bcast_real( theta , nl_max )
+  CALL wrf_dm_bcast_real( rho , nl_max )
+  CALL wrf_dm_bcast_real( u , nl_max )
+  CALL wrf_dm_bcast_real( v , nl_max )
+  CALL wrf_dm_bcast_real( qv , nl_max )
+  CALL wrf_dm_bcast_integer ( nl_in , 1 ) 
+
+  write(6,*) ' returned from reading sounding, nl_in is ',nl_in
+
+!  find ptop for the desired ztop (ztop is input from the namelist),
+!  and find surface pressure
+
+  grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
+
+!!MARS
+!!MARS
+  open(unit=14,file='input_coord',form='formatted',status='old')
+  rewind(14)
+  read(14,*) lon_input
+  read(14,*) lat_input
+  close(14)
+  write(6,*) ' lon is ',lon_input
+  write(6,*) ' lat is ',lat_input
+!!MARS
+!!MARS
+
+!!MARS
+!!MARS
+  open(unit=18,file='input_more',form='formatted',status='old')
+  rewind(18)
+  read(18,*) alt_input, tsurf_input
+  close(18)
+  write(6,*) ' alt is ',alt_input
+  write(6,*) ' tsurf is ',tsurf_input
+!!MARS
+!!MARS
+
+  DO j=jts,jte
+  DO i=its,ite
+!!MARS
+    grid%ht(i,j) = alt_input
+      ! flat surface
+      !!    grid%ht(i,j) = 0.
+      !    grid%ht(i,j) = hm1*exp(-(( float(i-icm)/xa1)**2))   &
+      !               *( (cos(pii*float(i-icm)/xal1))**2 )
+      !****Mars
+      !!3D hill
+      !        grid%ht(i,j) = hm/(1.+(float(i-icm)/xa)**2+(float(j-jcm)/xa)**2)
+
+!!2D hill
+grid%ht(i,j) = alt_input + hm/(1.+(float(i-icm)/xa)**2)
+
+      !!!3D crater    
+      !!        grid%ht(i,j) = hm - hm/(1.+(float(i-icm)/xa)**2+(float(j-jcm)/xa)**2)
+      !!3D crater w/ rims
+      !  x_param = float(i-icm)
+      !  y_param = float(j-jcm)
+      !  dilat = xa/2
+      !  rho_param = sqrt(x_param**2 + y_param**2)
+      !  ! revolution surface ; seed is a fourth order polynom
+      !  grid%ht(i,j) = (rho_param+6*dilat)*(rho_param+10*dilat)
+      !  grid%ht(i,j) = (rho_param-6*dilat)*(rho_param-10*dilat)*grid%ht(i,j)
+      !  ! flat terrain elsewhere - smooth gradient (no abrupt fall)
+      !  grid%ht(i,j) = grid%ht(i,j)*(tanh(rho_param+7*dilat)/2 - tanh(rho_param-7*dilat)/2)
+      !  grid%ht(i,j) = hm - (hm*.4/1500)*grid%ht(i,j)/(dilat**4)
+      !           !NONONONONON  
+      !           !grid%ht(i,j) = grid%ht(i,j) + alt_input  
+      !   !if (rho_param .GE. dilat*10) ht(i,j) = hm
+
+    grid%tsk(i,j) = tsurf_input
+!!MARS
+    grid%xlat(i,j) = lat_input
+    grid%xlong(i,j) = lon_input!+float(i)*config_flags%dx/59000.
+    grid%mars_emiss(i,j)=0.95
+    grid%mars_cice(i,j)=0.
+    grid%slpx(i,j) = 0.
+    grid%slpy(i,j) = 0.
+!!MARS
+  ENDDO
+  ENDDO
+
+  xs=ide/2 -3
+  xs=ids   -3
+  xe=xs + 6
+  ys=jde/2 -3
+  ye=ys + 6
+  mtn_ht = 500
+#ifdef MTN
+  DO j=max(ys,jds),min(ye,jde-1)
+  DO i=max(xs,ids),min(xe,ide-1)
+     grid%ht(i,j) = mtn_ht * 0.25 * &
+               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * &
+               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
+  ENDDO
+  ENDDO
+#endif
+#ifdef EW_RIDGE
+  DO j=max(ys,jds),min(ye,jde-1)
+  DO i=ids,ide
+     grid%ht(i,j) = mtn_ht * 0.50 * &
+               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
+  ENDDO
+  ENDDO
+#endif
+#ifdef NS_RIDGE
+  DO j=jds,jde
+  DO i=max(xs,ids),min(xe,ide-1)
+     grid%ht(i,j) = mtn_ht * 0.50 * &
+               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) )
+  ENDDO
+  ENDDO
+#endif
+  DO j=jts,jte
+  DO i=its,ite
+    grid%em_phb(i,1,j) = g * grid%ht(i,j)
+    grid%em_ph0(i,1,j) = g * grid%ht(i,j)
+  ENDDO
+  ENDDO
+
+!!!dans hill_2d OK
+!    grid%em_phb(i,1,j) = g*grid%ht(i,j)
+!    grid%em_php(i,1,j) = 0.
+!    grid%em_ph0(i,1,j) = grid%em_phb(i,1,j)
+
+  DO J = jts, jte
+  DO I = its, ite
+
+    p_surf = interp_0( p_in, zk, grid%em_phb(i,1,j)/g, nl_in )
+    grid%em_mub(i,j) = p_surf-grid%p_top
+
+!  this is dry hydrostatic sounding (base state), so given grid%em_p (coordinate),
+!  interp theta (from interp) and compute 1/rho from eqn. of state
+
+    DO K = 1, kte-1
+      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
+      grid%em_pb(i,k,j) = p_level
+      grid%em_t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
+      grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm
+    ENDDO
+
+!  calc hydrostatic balance (alternatively we could interp the geopotential from the
+!  sounding, but this assures that the base state is in exact hydrostatic balance with
+!  respect to the model eqns.
+
+    DO k  = 2,kte
+      grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
+    ENDDO
+
+  ENDDO
+  ENDDO
+
+  IF ( wrf_dm_on_monitor() ) THEN
+    write(6,*) ' ptop is ',grid%p_top
+    write(6,*) ' base state grid%em_mub(1,1), p_surf is ',grid%em_mub(1,1),grid%em_mub(1,1)+grid%p_top
+  ENDIF
+
+!  calculate full state for each column - this includes moisture.
+
+!!!!!MARS MARS
+!  write(6,*) ' getting moist sounding for full state '
+!  dry_sounding = .false.
+  dry_sounding = .true.
+  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in, &
+                              mulu, mulv, addu, addv )
+
+  DO J = jts, min(jde-1,jte)
+  DO I = its, min(ide-1,ite)
+
+!  At this point grid%p_top is already set. find the DRY mass in the column 
+!  by interpolating the DRY pressure.  
+
+   pd_surf = interp_0( pd_in, zk, grid%em_phb(i,1,j)/g, nl_in )
+
+!  compute the perturbation mass and the full mass
+
+    grid%em_mu_1(i,j) = pd_surf-grid%p_top - grid%em_mub(i,j)
+    grid%em_mu_2(i,j) = grid%em_mu_1(i,j)
+    grid%em_mu0(i,j) = grid%em_mu_1(i,j) + grid%em_mub(i,j)
+
+! given the dry pressure and coordinate system, interp the potential
+! temperature and qv
+
+    do k=1,kde-1
+
+      p_level = grid%em_znu(k)*(pd_surf - grid%p_top) + grid%p_top
+
+      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
+      grid%em_t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
+      grid%em_t_2(i,k,j)          = grid%em_t_1(i,k,j)
+      
+
+    enddo
+
+!  integrate the hydrostatic equation (from the RHS of the bigstep
+!  vertical momentum equation) down from the top to get grid%em_p.
+!  first from the top of the model to the top pressure
+
+    k = kte-1  ! top level
+
+    qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
+    qvf2 = 1./(1.+qvf1)
+    qvf1 = qvf1*qvf2
+
+!    grid%em_p(i,k,j) = - 0.5*grid%em_mu_1(i,j)/grid%em_rdnw(k)
+    grid%em_p(i,k,j) = - 0.5*(grid%em_mu_1(i,j)+qvf1*grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2
+    qvf = 1. + rvovrd*moist(i,k,j,P_QV)
+    grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
+                (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
+    grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
+
+!  down the column
+
+    do k=kte-2,1,-1
+      qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
+      qvf2 = 1./(1.+qvf1)
+      qvf1 = qvf1*qvf2
+      grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_1(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1)
+      qvf = 1. + rvovrd*moist(i,k,j,P_QV)
+      grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
+                  (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
+      grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
+    enddo
+
+!  this is the hydrostatic equation used in the model after the
+!  small timesteps.  In the model, grid%em_al (inverse density)
+!  is computed from the geopotential.
+
+
+    grid%em_ph_1(i,1,j) = 0.
+    DO k  = 2,kte
+      grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
+                   (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ &
+                    grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j)  )
+                                                   
+      grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j) 
+      grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
+    ENDDO
+
+    IF ( wrf_dm_on_monitor() ) THEN
+    if((i==2) .and. (j==2)) then
+     write(6,*) ' grid%em_ph_1 calc ',grid%em_ph_1(2,1,2),grid%em_ph_1(2,2,2),&
+                              grid%em_mu_1(2,2)+grid%em_mub(2,2),grid%em_mu_1(2,2), &
+                              grid%em_alb(2,1,2),grid%em_al(1,2,1),grid%em_rdnw(1)
+    endif
+    ENDIF
+
+  ENDDO
+  ENDDO
+
+!#if 0
+
+!  thermal perturbation to kick off convection
+
+  write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
+  write(6,*) ' delt for perturbation ',delt
+
+  DO J = jts, min(jde-1,jte)
+    yrad = config_flags%dy*float(j-nyc)/10000.
+yrad = 0.
+    DO I = its, min(ide-1,ite)
+      xrad = config_flags%dx*float(i-nxc)/10000.
+xrad = 0.
+      DO K = 1, kte-1
+
+!  put in preturbation theta (bubble) and recalc density.  note,
+!  the mass in the column is not changing, so when theta changes,
+!  we recompute density and geopotential
+
+        zrad = 0.5*(grid%em_ph_1(i,k,j)+grid%em_ph_1(i,k+1,j)  &
+                   +grid%em_phb(i,k,j)+grid%em_phb(i,k+1,j))/g
+        zrad = (zrad-1500.)/1500.
+        RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
+        IF(RAD <= 1.) THEN
+           grid%em_t_1(i,k,j)=grid%em_t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
+           grid%em_t_2(i,k,j)=grid%em_t_1(i,k,j)
+           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
+           grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
+                        (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
+           grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
+        ENDIF
+      ENDDO
+
+!  rebalance hydrostatically
+
+      DO k  = 2,kte
+        grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
+                     (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ &
+                      grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j)  )
+                                                   
+        grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j) 
+        grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
+      ENDDO
+
+    ENDDO
+  ENDDO
+
+!#endif
+
+   IF ( wrf_dm_on_monitor() ) THEN
+   write(6,*) ' grid%em_mu_1 from comp ', grid%em_mu_1(1,1)
+   write(6,*) ' full state sounding from comp, ph, grid%em_p, grid%em_al, grid%em_t_1, qv '
+   do k=1,kde-1
+     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1)+grid%em_phb(1,k,1), &
+                                      grid%em_p(1,k,1)+grid%em_pb(1,k,1), grid%em_alt(1,k,1), &
+                                      grid%em_t_1(1,k,1)+t0, moist(1,k,1,P_QV)
+   enddo
+
+   write(6,*) ' pert state sounding from comp, grid%em_ph_1, pp, alp, grid%em_t_1, qv '
+   do k=1,kde-1
+     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1), &
+                                      grid%em_p(1,k,1), grid%em_al(1,k,1), &
+                                      grid%em_t_1(1,k,1), moist(1,k,1,P_QV)
+   enddo
+   ENDIF
+
+! interp v
+
+  DO J = jts, jte
+  DO I = its, min(ide-1,ite)
+
+    IF (j == jds) THEN
+      z_at_v = grid%em_phb(i,1,j)/g
+    ELSE IF (j == jde) THEN
+      z_at_v = grid%em_phb(i,1,j-1)/g
+    ELSE
+      z_at_v = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i,1,j-1))/g
+    END IF
+
+    p_surf = interp_0( p_in, zk, z_at_v, nl_in )
+
+    DO K = 1, kte-1
+      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
+      grid%em_v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
+      grid%em_v_2(i,k,j) = grid%em_v_1(i,k,j)
+    ENDDO
+
+  ENDDO
+  ENDDO
+
+! interp u
+
+  DO J = jts, min(jde-1,jte)
+  DO I = its, ite
+
+    IF (i == ids) THEN
+      z_at_u = grid%em_phb(i,1,j)/g
+    ELSE IF (i == ide) THEN
+      z_at_u = grid%em_phb(i-1,1,j)/g
+    ELSE
+      z_at_u = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i-1,1,j))/g
+    END IF
+
+    p_surf = interp_0( p_in, zk, z_at_u, nl_in )
+
+    DO K = 1, kte-1
+      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
+      grid%em_u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
+      grid%em_u_2(i,k,j) = grid%em_u_1(i,k,j)
+    ENDDO
+
+  ENDDO
+  ENDDO
+
+!  set w
+
+  DO J = jts, min(jde-1,jte)
+  DO K = kts, kte
+  DO I = its, min(ide-1,ite)
+    grid%em_w_1(i,k,j) = 0.
+    grid%em_w_2(i,k,j) = 0.
+  ENDDO
+  ENDDO
+  ENDDO
+
+!  set a few more things
+
+  DO J = jts, min(jde-1,jte)
+  DO K = kts, kte-1
+  DO I = its, min(ide-1,ite)
+    grid%h_diabatic(i,k,j) = 0.
+  ENDDO
+  ENDDO
+  ENDDO
+
+  IF ( wrf_dm_on_monitor() ) THEN
+  DO k=1,kte-1
+    grid%em_t_base(k) = grid%em_t_1(1,k,1)
+    grid%qv_base(k) = moist(1,k,1,P_QV)
+    grid%u_base(k) = grid%em_u_1(1,k,1)
+    grid%v_base(k) = grid%em_v_1(1,k,1)
+    grid%z_base(k) = 0.5*(grid%em_phb(1,k,1)+grid%em_phb(1,k+1,1)+grid%em_ph_1(1,k,1)+grid%em_ph_1(1,k+1,1))/g
+  ENDDO
+  ENDIF
+  CALL wrf_dm_bcast_real( grid%em_t_base , kte )
+  CALL wrf_dm_bcast_real( grid%qv_base , kte )
+  CALL wrf_dm_bcast_real( grid%u_base , kte )
+  CALL wrf_dm_bcast_real( grid%v_base , kte )
+  CALL wrf_dm_bcast_real( grid%z_base , kte )
+
+  DO J = jts, min(jde-1,jte)
+  DO I = its, min(ide-1,ite)
+     thtmp   = grid%em_t_2(i,1,j)+t0
+     ptmp    = grid%em_p(i,1,j)+grid%em_pb(i,1,j)
+     temp(1) = thtmp * (ptmp/p1000mb)**rcp
+     thtmp   = grid%em_t_2(i,2,j)+t0
+     ptmp    = grid%em_p(i,2,j)+grid%em_pb(i,2,j)
+     temp(2) = thtmp * (ptmp/p1000mb)**rcp
+     thtmp   = grid%em_t_2(i,3,j)+t0
+     ptmp    = grid%em_p(i,3,j)+grid%em_pb(i,3,j)
+     temp(3) = thtmp * (ptmp/p1000mb)**rcp
+
+!!MARS
+!     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
+     grid%tmn(I,J)=grid%tsk(I,J)-0.5
+!!!MARS
+!!TODO: passer la valeur a partir des donnees
+!grid%mars_tsoil(I,:,J)=grid%tsk(I,J)
+!!!MARS
+  ENDDO
+  ENDDO
+
+ END SUBROUTINE init_domain_rk
+
+   SUBROUTINE init_module_initialize
+   END SUBROUTINE init_module_initialize
+
+!---------------------------------------------------------------------
+
+!  test driver for get_sounding
+!
+!      implicit none
+!      integer n
+!      parameter(n = 1000)
+!      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
+!      logical dry
+!      integer nl,k
+!
+!      dry = .false.
+!      dry = .true.
+!      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
+!      write(6,*) ' input levels ',nl
+!      write(6,*) ' sounding '
+!      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
+!      do k=1,nl
+!        write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), pd(k), theta(k), rho(k), u(k), v(k), qv(k)
+!      enddo
+!      end
+!
+!---------------------------------------------------------------------------
+
+      subroutine get_sounding( zk, p, p_dry, theta, rho, &
+                               u, v, qv, dry, nl_max, nl_in, &
+                               mulu, mulv, addu, addv )
+      implicit none
+
+      integer nl_max, nl_in
+      real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
+           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
+      logical dry
+
+      integer n
+      parameter(n=1000)
+      logical debug
+
+!      parameter( debug = .false.)
+!****Mars
+      parameter( debug = .true.)
+      real mulu, mulv, addu, addv
+
+
+! input sounding data
+
+      real p_surf, th_surf, qv_surf
+      real pi_surf, pi(n)
+      real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
+
+      !! special MARS
+      real r_input(n)
+      real cp_input(n)
+      real cv_input(n)
+      real cvpm_input(n)
+      real pfile_input(n)
+      real t_input(n)
+      real rhofile_input(n)
+      !! special MARS
+
+! diagnostics
+
+      real rho_surf, p_input(n), rho_input(n)
+      real pm_input(n)  !  this are for full moist sounding
+
+! local data
+
+      real p1000mb,cv,cp,r,cvpm,g
+!****Mars
+!      parameter (p1000mb = 1.e+05, r = 287, cp = 1003., cv = cp-r, cvpm = -cv/cp, g=9.81 )
+!      parameter (p1000mb = 610., r = 192., cp = 844.6, cv = cp-r, cvpm = -cv/cp, g=3.72)
+      parameter (p1000mb = 610., r = 191., cp = 744.5, cv = cp-r, cvpm = -cv/cp, g=3.72)
+!****Mars
+      integer k, it, nl
+      real qvf, qvf1, dz
+
+!  first, read the sounding
+
+      call read_sounding( p_surf, th_surf, qv_surf, &
+                          h_input, th_input, qv_input, u_input, v_input, r_input, cp_input, pfile_input, t_input, rhofile_input, n, nl, debug )
+
+
+      !! special MARS
+      do k=1,nl
+        cv_input(k) = cp_input(k) - r_input(k)
+        cvpm_input(k) = - cv_input(k) / cp_input(k)
+      enddo
+      !! special MARS
+
+      if(dry) then
+       do k=1,nl
+         qv_input(k) = 0.
+       enddo
+      endif
+
+      if(debug) write(6,*) ' number of input levels = ',nl
+
+        nl_in = nl
+        if(nl_in .gt. nl_max ) then
+          write(6,*) ' too many levels for input arrays ',nl_in,nl_max
+          call wrf_error_fatal ( ' too many levels for input arrays ' )
+        end if
+
+!  compute diagnostics,
+!  first, convert qv(g/kg) to qv(g/g)
+
+      do k=1,nl
+        qv_input(k) = 0.001*qv_input(k)
+      enddo
+
+      p_surf = 100.*p_surf  ! convert to pascals
+      qvf = 1. + rvovrd*qv_input(1) 
+      rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
+      pi_surf = (p_surf/p1000mb)**(r/cp)
+          !!!!!! rcp variable
+          !rho_surf =  1./((r_input(1)/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm_input(1)))
+          !pi_surf = (p_surf/p1000mb)**(r_input(1)/cp_input(1))
+
+
+      if(debug) then
+        write(6,*) ' surface density is ',rho_surf
+        write(6,*) ' surface pi is      ',pi_surf
+      end if
+
+
+!  integrate moist sounding hydrostatically, starting from the
+!  specified surface pressure
+!  -> first, integrate from surface to lowest level
+
+          qvf = 1. + rvovrd*qv_input(1) 
+          qvf1 = 1. + qv_input(1)
+          rho_input(1) = rho_surf
+          dz = h_input(1)
+          do it=1,10
+!!MARS MARS
+            pm_input(1) = p_surf !&
+!                    - dz*(0.25*rho_surf+0.75*rho_input(1))*g*qvf1  !!! BEURK
+!                    - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1  !! parce que couche 1 tres proche
+            rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
+               !!!!!!! rcp variable
+               !rho_input(1) = 1./((r_input(1)/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm_input(1)))
+          enddo
+
+! integrate up the column
+
+          do k=2,nl
+            rho_input(k) = rho_input(k-1)
+            dz = h_input(k)-h_input(k-1)
+               !!!!!!! rcp variable
+               !dz = r_input(k) * t_input(k) * (- p_input(k) + p_input(k-1)) / p_input(k) / g
+               !!dz = - cp_input(k) * (- t_input(k) + t_input(k-1)) / g
+            qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
+            qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
+
+print *, 'input', pfile_input(k), rhofile_input(k)  
+
+            do it=1,10 !!ou moins??? non. !! trop de rho(k-1) donne une pression trop faible puis crash
+                                          !! mais coeff ci-dessous vont varier la pression calculÃe
+              pm_input(k) = pm_input(k-1) &
+                      - dz*(0.75*rho_input(k)+0.25*rho_input(k-1))*g*qvf1
+                      !- 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
+              rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
+!!
+!! marche pas 
+!!
+!pm_input(k) = pm_input(k-1) &
+!        - 0.5*dz*(1./rho_input(k)+1./rho_input(k-1))*g*qvf1
+!rho_input(k) = (r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm)
+              !!!!!!! rcp variable
+              !rho_input(k) = 1./((r_input(k)/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm_input(k)))
+              !print *, p_input(k), pm_input(k),((r_input(k)/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm_input(k))),k
+print *, it, pm_input(k), rho_input(k), dz
+            enddo
+          enddo
+
+
+!  we have the moist sounding
+
+!  next, compute the dry sounding using p at the highest level from the
+!  moist sounding and integrating down.
+
+        p_input(nl) = pm_input(nl)
+
+          do k=nl-1,1,-1
+            dz = h_input(k+1)-h_input(k)
+            p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
+          enddo
+
+
+        do k=1,nl
+
+          zk(k) = h_input(k)
+          p(k) = pm_input(k)
+          p_dry(k) = p_input(k)
+          theta(k) = th_input(k)
+          rho(k) = rho_input(k)
+          u(k) = mulu*u_input(k) + addu
+          v(k) = mulv*v_input(k) + addv
+          qv(k) = qv_input(k)
+
+         !!!! direct input from file
+         write(6,*) '*** DIRECT INPUT FROM FILE ***' 
+         p(k) = pfile_input(k)
+         p_dry(k) = pfile_input(k)
+         rho(k) = rhofile_input(k)
+
+        enddo
+
+     if(debug) then
+      write(6,*) ' sounding '
+      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
+      do k=1,nl
+        write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), p_dry(k), theta(k), rho(k), u(k), v(k), qv(k)
+      enddo
+
+     end if
+
+      end subroutine get_sounding
+
+!-------------------------------------------------------
+
+      subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,r,cp,p,t,rho,n,nl,debug )
+      implicit none
+      integer n,nl
+      real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n),r(n),cp(n),p(n),t(n),rho(n) 
+      logical end_of_file
+      logical debug
+
+      integer k
+
+      open(unit=10,file='input_sounding',form='formatted',status='old')
+      rewind(10)
+      read(10,*) ps, ts, qvs
+      if(debug) then
+        write(6,*) ' input sounding surface parameters '
+        write(6,*) ' surface pressure (mb) ',ps
+        write(6,*) ' surface pot. temp (K) ',ts
+        write(6,*) ' surface mixing ratio (g/kg) ',qvs
+      end if
+
+      end_of_file = .false.
+      k = 0
+
+      do while (.not. end_of_file)
+
+        read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
+        k = k+1
+        if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
+        go to 110
+ 100    end_of_file = .true.
+ 110    continue
+      enddo
+
+
+      !!! special MARS
+      open(unit=11,file='input_therm',form='formatted',status='old')
+      rewind(11)
+      end_of_file = .false.
+      k = 0
+      do while (.not. end_of_file)
+
+        read(11,*,end=101) r(k+1), cp(k+1), p(k+1), rho(k+1), t(k+1) 
+        write(*,*) k, r(k+1), cp(k+1), p(k+1), rho(k+1), t(k+1)
+        k = k+1
+        go to 112
+ 101    end_of_file = .true.
+ 112    continue
+      enddo
+      !!! special MARS
+
+
+
+      nl = k
+
+      close(unit=10,status = 'keep')
+
+      end subroutine read_sounding
+
+END MODULE module_initialize
Index: /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F_ref
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F_ref	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F_ref	(revision 31)
@@ -0,0 +1,904 @@
+!IDEAL:MODEL_LAYER:INITIALIZATION
+!
+
+!  This MODULE holds the routines which are used to perform various initializations
+!  for the individual domains.  
+
+!  This MODULE CONTAINS the following routines:
+
+!  initialize_field_test - 1. Set different fields to different constant
+!                             values.  This is only a test.  If the correct
+!                             domain is not found (based upon the "id")
+!                             then a fatal error is issued.               
+
+!-----------------------------------------------------------------------
+
+MODULE module_initialize
+
+   USE module_domain
+   USE module_io_domain
+   USE module_state_description
+   USE module_model_constants
+   USE module_bc
+   USE module_timing
+   USE module_configure
+   USE module_init_utilities
+#ifdef DM_PARALLEL
+   USE module_dm
+#endif
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------
+! this is a wrapper for the solver-specific init_domain routines.
+! Also dereferences the grid variables and passes them down as arguments.
+! This is crucial, since the lower level routines may do message passing
+! and this will get fouled up on machines that insist on passing down
+! copies of assumed-shape arrays (by passing down as arguments, the 
+! data are treated as assumed-size -- ie. f77 -- arrays and the copying
+! business is avoided).  Fie on the F90 designers.  Fie and a pox.
+
+   SUBROUTINE init_domain ( grid )
+
+   IMPLICIT NONE
+
+   !  Input data.
+   TYPE (domain), POINTER :: grid 
+   !  Local data.
+   INTEGER                :: dyn_opt 
+   INTEGER :: idum1, idum2
+
+   CALL nl_get_dyn_opt( 1, dyn_opt )
+   
+   CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
+
+   IF (      dyn_opt .eq. 1 &
+        .or. dyn_opt .eq. 2 &
+        .or. dyn_opt .eq. 3 &
+                                       ) THEN
+     CALL init_domain_rk( grid &
+!
+#include <em_actual_new_args.inc>
+!
+                        )
+
+   ELSE
+     WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt
+     CALL wrf_error_fatal ( ' init_domain: unknown or unimplemented dyn_opt ' )
+   ENDIF
+
+   END SUBROUTINE init_domain
+
+!-------------------------------------------------------------------
+
+   SUBROUTINE init_domain_rk ( grid &
+!
+# include <em_dummy_new_args.inc>
+!
+)
+   IMPLICIT NONE
+
+   !  Input data.
+   TYPE (domain), POINTER :: grid
+
+# include <em_dummy_new_decl.inc>
+
+   TYPE (grid_config_rec_type)              :: config_flags
+
+   !  Local data
+   INTEGER                             ::                       &
+                                  ids, ide, jds, jde, kds, kde, &
+                                  ims, ime, jms, jme, kms, kme, &
+                                  its, ite, jts, jte, kts, kte, &
+                                  i, j, k
+
+   ! Local data
+
+   INTEGER, PARAMETER :: nl_max = 1000
+   REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
+   INTEGER :: nl_in
+
+
+   INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc
+   REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
+   REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
+!   REAL, EXTERNAL :: interp_0
+   REAL    :: hm
+   REAL    :: pi
+
+!  stuff from original initialization that has been dropped from the Registry 
+   REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
+   REAL    :: qvf1, qvf2, pd_surf
+   INTEGER :: it
+   real :: thtmp, ptmp, temp(3)
+
+   LOGICAL :: moisture_init
+   LOGICAL :: stretch_grid, dry_sounding
+
+  INTEGER :: xs , xe , ys , ye
+  REAL :: mtn_ht
+   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
+
+!!MARS
+ REAL :: lon_input, lat_input, alt_input, tsurf_input
+!!MARS
+
+
+#ifdef DM_PARALLEL
+#    include <em_data_calls.inc>
+#endif
+
+
+   SELECT CASE ( model_data_order )
+         CASE ( DATA_ORDER_ZXY )
+   kds = grid%sd31 ; kde = grid%ed31 ;
+   ids = grid%sd32 ; ide = grid%ed32 ;
+   jds = grid%sd33 ; jde = grid%ed33 ;
+
+   kms = grid%sm31 ; kme = grid%em31 ;
+   ims = grid%sm32 ; ime = grid%em32 ;
+   jms = grid%sm33 ; jme = grid%em33 ;
+
+   kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
+   its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
+   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
+         CASE ( DATA_ORDER_XYZ )
+   ids = grid%sd31 ; ide = grid%ed31 ;
+   jds = grid%sd32 ; jde = grid%ed32 ;
+   kds = grid%sd33 ; kde = grid%ed33 ;
+
+   ims = grid%sm31 ; ime = grid%em31 ;
+   jms = grid%sm32 ; jme = grid%em32 ;
+   kms = grid%sm33 ; kme = grid%em33 ;
+
+   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
+   jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
+   kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
+         CASE ( DATA_ORDER_XZY )
+   ids = grid%sd31 ; ide = grid%ed31 ;
+   kds = grid%sd32 ; kde = grid%ed32 ;
+   jds = grid%sd33 ; jde = grid%ed33 ;
+
+   ims = grid%sm31 ; ime = grid%em31 ;
+   kms = grid%sm32 ; kme = grid%em32 ;
+   jms = grid%sm33 ; jme = grid%em33 ;
+
+   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
+   kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
+   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
+
+   END SELECT
+
+
+   stretch_grid = .true.
+!   stretch_grid = .false.
+
+   delt = 3.
+!   delt = 10.
+!   z_scale = .50
+   z_scale = .40
+
+!!!MARS
+!!!MARS
+!  open(unit=16,file='input_vert',form='formatted',status='old')
+!  rewind(16) 
+!  read(16,*) delt, z_scale
+!  write(6,*) 'delt, z_scale are  ', delt, z_scale
+!  close(16)
+!!!MARS
+!!!MARS
+
+   pi = 2.*asin(1.0)
+   write(6,*) ' pi is ',pi
+   nxc = (ide-ids)/2
+   nyc = (jde-jds)/2
+
+   CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
+
+! here we check to see if the boundary conditions are set properly
+
+   CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
+
+   moisture_init = .true.
+
+    grid%itimestep=0
+
+#ifdef DM_PARALLEL
+   CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
+   CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
+#endif
+
+    CALL nl_set_mminlu(1, '    ')
+    CALL nl_set_iswater(1,0)
+    CALL nl_set_cen_lat(1,40.)
+    CALL nl_set_cen_lon(1,-105.)
+    CALL nl_set_truelat1(1,0.)
+    CALL nl_set_truelat2(1,0.)
+    CALL nl_set_moad_cen_lat (1,0.)
+    CALL nl_set_stand_lon (1,0.)
+    CALL nl_set_map_proj(1,0)
+
+
+!  here we initialize data we currently is not initialized 
+!  in the input data
+
+    DO j = jts, jte
+      DO i = its, ite
+         grid%msft(i,j)     = 1.
+         grid%msfu(i,j)     = 1.
+         grid%msfv(i,j)     = 1.
+         grid%sina(i,j)     = 0.
+         grid%cosa(i,j)     = 1.
+         grid%e(i,j)        = 0.
+         grid%f(i,j)        = 0.         !! MARS: put coriolis here if needed
+
+      END DO
+   END DO
+
+    DO j = jts, jte
+    DO k = kts, kte
+      DO i = its, ite
+         grid%em_ww(i,k,j)     = 0.
+      END DO
+   END DO
+   END DO
+
+   grid%step_number = 0
+
+! set up the grid
+
+   IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
+     DO k=1, kde
+      grid%em_znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
+                                (1.-exp(-1./z_scale))
+     ENDDO
+   ELSE
+     DO k=1, kde
+      grid%em_znw(k) = 1. - float(k-1)/float(kde-1)
+     ENDDO
+   ENDIF
+
+!!!MARS
+!!!MARS
+!  open(unit=12,file='levels',form='formatted',status='old')
+!  rewind(12)
+!  DO k=1, kde
+!  read(12,*) grid%em_znw(k)
+!  write(6,*) 'read level ', k,grid%em_znw(k)
+!  ENDDO
+!  close(12)
+!!!MARS
+!!!MARS
+   
+
+   DO k=1, kde-1
+    grid%em_dnw(k) = grid%em_znw(k+1) - grid%em_znw(k)
+    grid%em_rdnw(k) = 1./grid%em_dnw(k)
+    grid%em_znu(k) = 0.5*(grid%em_znw(k+1)+grid%em_znw(k))
+   ENDDO
+   DO k=2, kde-1
+    grid%em_dn(k) = 0.5*(grid%em_dnw(k)+grid%em_dnw(k-1))
+    grid%em_rdn(k) = 1./grid%em_dn(k)
+    grid%em_fnp(k) = .5* grid%em_dnw(k  )/grid%em_dn(k)
+    grid%em_fnm(k) = .5* grid%em_dnw(k-1)/grid%em_dn(k)
+   ENDDO
+
+   cof1 = (2.*grid%em_dn(2)+grid%em_dn(3))/(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(2) 
+   cof2 =     grid%em_dn(2)        /(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(3) 
+   grid%cf1  = grid%em_fnp(2) + cof1
+   grid%cf2  = grid%em_fnm(2) - cof1 - cof2
+   grid%cf3  = cof2       
+
+   grid%cfn  = (.5*grid%em_dnw(kde-1)+grid%em_dn(kde-1))/grid%em_dn(kde-1)
+   grid%cfn1 = -.5*grid%em_dnw(kde-1)/grid%em_dn(kde-1)
+   grid%rdx = 1./config_flags%dx
+   grid%rdy = 1./config_flags%dy
+
+!  get the sounding from the ascii sounding file, first get dry sounding and 
+!  calculate base state
+
+  dry_sounding = .true.
+  IF ( wrf_dm_on_monitor() ) THEN
+  write(6,*) ' getting dry sounding for base state '
+
+  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
+  ENDIF
+  CALL wrf_dm_bcast_real( zk , nl_max )
+  CALL wrf_dm_bcast_real( p_in , nl_max )
+  CALL wrf_dm_bcast_real( pd_in , nl_max )
+  CALL wrf_dm_bcast_real( theta , nl_max )
+  CALL wrf_dm_bcast_real( rho , nl_max )
+  CALL wrf_dm_bcast_real( u , nl_max )
+  CALL wrf_dm_bcast_real( v , nl_max )
+  CALL wrf_dm_bcast_real( qv , nl_max )
+  CALL wrf_dm_bcast_integer ( nl_in , 1 ) 
+
+  write(6,*) ' returned from reading sounding, nl_in is ',nl_in
+
+!  find ptop for the desired ztop (ztop is input from the namelist),
+!  and find surface pressure
+
+  grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
+
+!!MARS
+!!MARS
+  open(unit=14,file='input_coord',form='formatted',status='old')
+  rewind(14)
+  read(14,*) lon_input
+  read(14,*) lat_input
+  close(14)
+  write(6,*) ' lon is ',lon_input
+  write(6,*) ' lat is ',lat_input
+!!MARS
+!!MARS
+
+!!MARS
+!!MARS
+  open(unit=18,file='input_more',form='formatted',status='old')
+  rewind(18)
+  read(18,*) alt_input, tsurf_input
+  close(18)
+  write(6,*) ' alt is ',alt_input
+  write(6,*) ' tsurf is ',tsurf_input
+!!MARS
+!!MARS
+
+  DO j=jts,jte
+  DO i=its,ite
+!!MARS
+    grid%ht(i,j) = alt_input
+    grid%tsk(i,j) = tsurf_input
+!!MARS
+    grid%xlat(i,j) = lat_input
+    grid%xlong(i,j) = lon_input!+float(i)*config_flags%dx/59000.
+    grid%mars_emiss(i,j)=0.95
+    grid%mars_cice(i,j)=0.
+    grid%slpx(i,j) = 0.
+    grid%slpy(i,j) = 0.
+!!MARS
+  ENDDO
+  ENDDO
+
+  xs=ide/2 -3
+  xs=ids   -3
+  xe=xs + 6
+  ys=jde/2 -3
+  ye=ys + 6
+  mtn_ht = 500
+#ifdef MTN
+  DO j=max(ys,jds),min(ye,jde-1)
+  DO i=max(xs,ids),min(xe,ide-1)
+     grid%ht(i,j) = mtn_ht * 0.25 * &
+               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * &
+               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
+  ENDDO
+  ENDDO
+#endif
+#ifdef EW_RIDGE
+  DO j=max(ys,jds),min(ye,jde-1)
+  DO i=ids,ide
+     grid%ht(i,j) = mtn_ht * 0.50 * &
+               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
+  ENDDO
+  ENDDO
+#endif
+#ifdef NS_RIDGE
+  DO j=jds,jde
+  DO i=max(xs,ids),min(xe,ide-1)
+     grid%ht(i,j) = mtn_ht * 0.50 * &
+               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) )
+  ENDDO
+  ENDDO
+#endif
+  DO j=jts,jte
+  DO i=its,ite
+    grid%em_phb(i,1,j) = g * grid%ht(i,j)
+    grid%em_ph0(i,1,j) = g * grid%ht(i,j)
+  ENDDO
+  ENDDO
+
+  DO J = jts, jte
+  DO I = its, ite
+
+    p_surf = interp_0( p_in, zk, grid%em_phb(i,1,j)/g, nl_in )
+    grid%em_mub(i,j) = p_surf-grid%p_top
+
+!  this is dry hydrostatic sounding (base state), so given grid%em_p (coordinate),
+!  interp theta (from interp) and compute 1/rho from eqn. of state
+
+    DO K = 1, kte-1
+      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
+      grid%em_pb(i,k,j) = p_level
+      grid%em_t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
+      grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm
+    ENDDO
+
+!  calc hydrostatic balance (alternatively we could interp the geopotential from the
+!  sounding, but this assures that the base state is in exact hydrostatic balance with
+!  respect to the model eqns.
+
+    DO k  = 2,kte
+      grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
+    ENDDO
+
+  ENDDO
+  ENDDO
+
+  IF ( wrf_dm_on_monitor() ) THEN
+    write(6,*) ' ptop is ',grid%p_top
+    write(6,*) ' base state grid%em_mub(1,1), p_surf is ',grid%em_mub(1,1),grid%em_mub(1,1)+grid%p_top
+  ENDIF
+
+!  calculate full state for each column - this includes moisture.
+
+!!!!!MARS MARS
+!  write(6,*) ' getting moist sounding for full state '
+!  dry_sounding = .false.
+  dry_sounding = .true.
+  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
+
+  DO J = jts, min(jde-1,jte)
+  DO I = its, min(ide-1,ite)
+
+!  At this point grid%p_top is already set. find the DRY mass in the column 
+!  by interpolating the DRY pressure.  
+
+   pd_surf = interp_0( pd_in, zk, grid%em_phb(i,1,j)/g, nl_in )
+
+!  compute the perturbation mass and the full mass
+
+    grid%em_mu_1(i,j) = pd_surf-grid%p_top - grid%em_mub(i,j)
+    grid%em_mu_2(i,j) = grid%em_mu_1(i,j)
+    grid%em_mu0(i,j) = grid%em_mu_1(i,j) + grid%em_mub(i,j)
+
+! given the dry pressure and coordinate system, interp the potential
+! temperature and qv
+
+    do k=1,kde-1
+
+      p_level = grid%em_znu(k)*(pd_surf - grid%p_top) + grid%p_top
+
+      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
+      grid%em_t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
+      grid%em_t_2(i,k,j)          = grid%em_t_1(i,k,j)
+      
+
+    enddo
+
+!  integrate the hydrostatic equation (from the RHS of the bigstep
+!  vertical momentum equation) down from the top to get grid%em_p.
+!  first from the top of the model to the top pressure
+
+    k = kte-1  ! top level
+
+    qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
+    qvf2 = 1./(1.+qvf1)
+    qvf1 = qvf1*qvf2
+
+!    grid%em_p(i,k,j) = - 0.5*grid%em_mu_1(i,j)/grid%em_rdnw(k)
+    grid%em_p(i,k,j) = - 0.5*(grid%em_mu_1(i,j)+qvf1*grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2
+    qvf = 1. + rvovrd*moist(i,k,j,P_QV)
+    grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
+                (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
+    grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
+
+!  down the column
+
+    do k=kte-2,1,-1
+      qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
+      qvf2 = 1./(1.+qvf1)
+      qvf1 = qvf1*qvf2
+      grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_1(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1)
+      qvf = 1. + rvovrd*moist(i,k,j,P_QV)
+      grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
+                  (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
+      grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
+    enddo
+
+!  this is the hydrostatic equation used in the model after the
+!  small timesteps.  In the model, grid%em_al (inverse density)
+!  is computed from the geopotential.
+
+
+    grid%em_ph_1(i,1,j) = 0.
+    DO k  = 2,kte
+      grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
+                   (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ &
+                    grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j)  )
+                                                   
+      grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j) 
+      grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
+    ENDDO
+
+    IF ( wrf_dm_on_monitor() ) THEN
+    if((i==2) .and. (j==2)) then
+     write(6,*) ' grid%em_ph_1 calc ',grid%em_ph_1(2,1,2),grid%em_ph_1(2,2,2),&
+                              grid%em_mu_1(2,2)+grid%em_mub(2,2),grid%em_mu_1(2,2), &
+                              grid%em_alb(2,1,2),grid%em_al(1,2,1),grid%em_rdnw(1)
+    endif
+    ENDIF
+
+  ENDDO
+  ENDDO
+
+!#if 0
+
+!  thermal perturbation to kick off convection
+
+  write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
+  write(6,*) ' delt for perturbation ',delt
+
+  DO J = jts, min(jde-1,jte)
+    yrad = config_flags%dy*float(j-nyc)/10000.
+!   yrad = 0.
+    DO I = its, min(ide-1,ite)
+      xrad = config_flags%dx*float(i-nxc)/10000.
+!     xrad = 0.
+      DO K = 1, kte-1
+
+!  put in preturbation theta (bubble) and recalc density.  note,
+!  the mass in the column is not changing, so when theta changes,
+!  we recompute density and geopotential
+
+        zrad = 0.5*(grid%em_ph_1(i,k,j)+grid%em_ph_1(i,k+1,j)  &
+                   +grid%em_phb(i,k,j)+grid%em_phb(i,k+1,j))/g
+        zrad = (zrad-1500.)/1500.
+        RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
+        IF(RAD <= 1.) THEN
+           grid%em_t_1(i,k,j)=grid%em_t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
+           grid%em_t_2(i,k,j)=grid%em_t_1(i,k,j)
+           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
+           grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
+                        (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
+           grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
+        ENDIF
+      ENDDO
+
+!  rebalance hydrostatically
+
+      DO k  = 2,kte
+        grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
+                     (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ &
+                      grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j)  )
+                                                   
+        grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j) 
+        grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
+      ENDDO
+
+    ENDDO
+  ENDDO
+
+!#endif
+
+   IF ( wrf_dm_on_monitor() ) THEN
+   write(6,*) ' grid%em_mu_1 from comp ', grid%em_mu_1(1,1)
+   write(6,*) ' full state sounding from comp, ph, grid%em_p, grid%em_al, grid%em_t_1, qv '
+   do k=1,kde-1
+     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1)+grid%em_phb(1,k,1), &
+                                      grid%em_p(1,k,1)+grid%em_pb(1,k,1), grid%em_alt(1,k,1), &
+                                      grid%em_t_1(1,k,1)+t0, moist(1,k,1,P_QV)
+   enddo
+
+   write(6,*) ' pert state sounding from comp, grid%em_ph_1, pp, alp, grid%em_t_1, qv '
+   do k=1,kde-1
+     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1), &
+                                      grid%em_p(1,k,1), grid%em_al(1,k,1), &
+                                      grid%em_t_1(1,k,1), moist(1,k,1,P_QV)
+   enddo
+   ENDIF
+
+! interp v
+
+  DO J = jts, jte
+  DO I = its, min(ide-1,ite)
+
+    IF (j == jds) THEN
+      z_at_v = grid%em_phb(i,1,j)/g
+    ELSE IF (j == jde) THEN
+      z_at_v = grid%em_phb(i,1,j-1)/g
+    ELSE
+      z_at_v = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i,1,j-1))/g
+    END IF
+    p_surf = interp_0( p_in, zk, z_at_v, nl_in )
+
+    DO K = 1, kte-1
+      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
+      grid%em_v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
+      grid%em_v_2(i,k,j) = grid%em_v_1(i,k,j)
+    ENDDO
+
+  ENDDO
+  ENDDO
+
+! interp u
+
+  DO J = jts, min(jde-1,jte)
+  DO I = its, ite
+
+    IF (i == ids) THEN
+      z_at_u = grid%em_phb(i,1,j)/g
+    ELSE IF (i == ide) THEN
+      z_at_u = grid%em_phb(i-1,1,j)/g
+    ELSE
+      z_at_u = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i-1,1,j))/g
+    END IF
+
+    p_surf = interp_0( p_in, zk, z_at_u, nl_in )
+
+    DO K = 1, kte-1
+      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
+      grid%em_u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
+      grid%em_u_2(i,k,j) = grid%em_u_1(i,k,j)
+    ENDDO
+
+  ENDDO
+  ENDDO
+
+!  set w
+
+  DO J = jts, min(jde-1,jte)
+  DO K = kts, kte
+  DO I = its, min(ide-1,ite)
+    grid%em_w_1(i,k,j) = 0.
+    grid%em_w_2(i,k,j) = 0.
+  ENDDO
+  ENDDO
+  ENDDO
+
+!  set a few more things
+
+  DO J = jts, min(jde-1,jte)
+  DO K = kts, kte-1
+  DO I = its, min(ide-1,ite)
+    grid%h_diabatic(i,k,j) = 0.
+  ENDDO
+  ENDDO
+  ENDDO
+
+  IF ( wrf_dm_on_monitor() ) THEN
+  DO k=1,kte-1
+    grid%em_t_base(k) = grid%em_t_1(1,k,1)
+    grid%qv_base(k) = moist(1,k,1,P_QV)
+    grid%u_base(k) = grid%em_u_1(1,k,1)
+    grid%v_base(k) = grid%em_v_1(1,k,1)
+    grid%z_base(k) = 0.5*(grid%em_phb(1,k,1)+grid%em_phb(1,k+1,1)+grid%em_ph_1(1,k,1)+grid%em_ph_1(1,k+1,1))/g
+  ENDDO
+  ENDIF
+  CALL wrf_dm_bcast_real( grid%em_t_base , kte )
+  CALL wrf_dm_bcast_real( grid%qv_base , kte )
+  CALL wrf_dm_bcast_real( grid%u_base , kte )
+  CALL wrf_dm_bcast_real( grid%v_base , kte )
+  CALL wrf_dm_bcast_real( grid%z_base , kte )
+
+  DO J = jts, min(jde-1,jte)
+  DO I = its, min(ide-1,ite)
+     thtmp   = grid%em_t_2(i,1,j)+t0
+     ptmp    = grid%em_p(i,1,j)+grid%em_pb(i,1,j)
+     temp(1) = thtmp * (ptmp/p1000mb)**rcp
+     thtmp   = grid%em_t_2(i,2,j)+t0
+     ptmp    = grid%em_p(i,2,j)+grid%em_pb(i,2,j)
+     temp(2) = thtmp * (ptmp/p1000mb)**rcp
+     thtmp   = grid%em_t_2(i,3,j)+t0
+     ptmp    = grid%em_p(i,3,j)+grid%em_pb(i,3,j)
+     temp(3) = thtmp * (ptmp/p1000mb)**rcp
+
+!!MARS
+!     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
+     grid%tmn(I,J)=grid%tsk(I,J)-0.5
+!!!MARS
+!!TODO: passer la valeur a partir des donnees
+!grid%mars_tsoil(I,:,J)=grid%tsk(I,J)
+!!!MARS
+  ENDDO
+  ENDDO
+
+ END SUBROUTINE init_domain_rk
+
+   SUBROUTINE init_module_initialize
+   END SUBROUTINE init_module_initialize
+
+!---------------------------------------------------------------------
+
+!  test driver for get_sounding
+!
+!      implicit none
+!      integer n
+!      parameter(n = 1000)
+!      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
+!      logical dry
+!      integer nl,k
+!
+!      dry = .false.
+!      dry = .true.
+!      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
+!      write(6,*) ' input levels ',nl
+!      write(6,*) ' sounding '
+!      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
+!      do k=1,nl
+!        write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), pd(k), theta(k), rho(k), u(k), v(k), qv(k)
+!      enddo
+!      end
+!
+!---------------------------------------------------------------------------
+
+      subroutine get_sounding( zk, p, p_dry, theta, rho, &
+                               u, v, qv, dry, nl_max, nl_in )
+      implicit none
+
+      integer nl_max, nl_in
+      real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
+           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
+      logical dry
+
+      integer n
+      parameter(n=1000)
+      logical debug
+      parameter( debug = .true.)
+
+! input sounding data
+
+      real p_surf, th_surf, qv_surf
+      real pi_surf, pi(n)
+      real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
+
+! diagnostics
+
+      real rho_surf, p_input(n), rho_input(n)
+      real pm_input(n)  !  this are for full moist sounding
+
+! local data
+
+      real p1000mb,cv,cp,r,cvpm,g
+!      parameter (p1000mb = 1.e+05, r = 287, cp = 1003., cv = cp-r, cvpm = -cv/cp, g=9.81 )
+      parameter (p1000mb = 610., r = 192., cp = 844.6, cv = cp-r, cvpm = -cv/cp, g=3.72)
+      integer k, it, nl
+      real qvf, qvf1, dz
+
+!  first, read the sounding
+
+      call read_sounding( p_surf, th_surf, qv_surf, &
+                          h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
+
+      if(dry) then
+       do k=1,nl
+         qv_input(k) = 0.
+       enddo
+      endif
+
+      if(debug) write(6,*) ' number of input levels = ',nl
+
+        nl_in = nl
+        if(nl_in .gt. nl_max ) then
+          write(6,*) ' too many levels for input arrays ',nl_in,nl_max
+          call wrf_error_fatal ( ' too many levels for input arrays ' )
+        end if
+
+!  compute diagnostics,
+!  first, convert qv(g/kg) to qv(g/g)
+
+      do k=1,nl
+        qv_input(k) = 0.001*qv_input(k)
+      enddo
+
+      p_surf = 100.*p_surf  ! convert to pascals
+      qvf = 1. + rvovrd*qv_input(1) 
+      rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
+      pi_surf = (p_surf/p1000mb)**(r/cp)
+
+      if(debug) then
+        write(6,*) ' surface density is ',rho_surf
+        write(6,*) ' surface pi is      ',pi_surf
+      end if
+
+
+!  integrate moist sounding hydrostatically, starting from the
+!  specified surface pressure
+!  -> first, integrate from surface to lowest level
+
+          qvf = 1. + rvovrd*qv_input(1) 
+          qvf1 = 1. + qv_input(1)
+          rho_input(1) = rho_surf
+          dz = h_input(1)
+          do it=1,10
+!            pm_input(1) = p_surf &
+!                    - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
+!!!MARS MARS MARS
+            pm_input(1) = p_surf 
+            rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
+          enddo
+
+! integrate up the column
+
+          do k=2,nl
+            rho_input(k) = rho_input(k-1)
+            dz = h_input(k)-h_input(k-1)
+            qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
+            qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
+ 
+            do it=1,10
+              pm_input(k) = pm_input(k-1) &
+                      - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
+              rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
+            enddo
+          enddo
+
+!  we have the moist sounding
+
+!  next, compute the dry sounding using p at the highest level from the
+!  moist sounding and integrating down.
+
+        p_input(nl) = pm_input(nl)
+
+          do k=nl-1,1,-1
+            dz = h_input(k+1)-h_input(k)
+            p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
+          enddo
+
+
+        do k=1,nl
+
+          zk(k) = h_input(k)
+          p(k) = pm_input(k)
+          p_dry(k) = p_input(k)
+          theta(k) = th_input(k)
+          rho(k) = rho_input(k)
+          u(k) = u_input(k)
+          v(k) = v_input(k)
+          qv(k) = qv_input(k)
+
+        enddo
+
+     if(debug) then
+      write(6,*) ' sounding '
+      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
+      do k=1,nl
+        write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), p_dry(k), theta(k), rho(k), u(k), v(k), qv(k)
+      enddo
+
+     end if
+
+      end subroutine get_sounding
+
+!-------------------------------------------------------
+
+      subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
+      implicit none
+      integer n,nl
+      real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
+      logical end_of_file
+      logical debug
+
+      integer k
+
+      open(unit=10,file='input_sounding',form='formatted',status='old')
+      rewind(10)
+      read(10,*) ps, ts, qvs
+      if(debug) then
+        write(6,*) ' input sounding surface parameters '
+        write(6,*) ' surface pressure (mb) ',ps
+        write(6,*) ' surface pot. temp (K) ',ts
+        write(6,*) ' surface mixing ratio (g/kg) ',qvs
+      end if
+
+      end_of_file = .false.
+      k = 0
+
+      do while (.not. end_of_file)
+
+        read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
+        k = k+1
+        if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
+        go to 110
+ 100    end_of_file = .true.
+ 110    continue
+      enddo
+
+      nl = k
+
+      close(unit=10,status = 'keep')
+
+      end subroutine read_sounding
+
+END MODULE module_initialize
Index: /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F_gcm
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F_gcm	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F_gcm	(revision 31)
@@ -0,0 +1,178 @@
+!WRF:MODEL_LAYER:CONSTANTS
+!
+!PLANET MARS
+
+ MODULE module_model_constants
+
+   !  2. Following are constants for use in defining real number bounds.
+
+   !  A really small number.
+
+   REAL    , PARAMETER :: epsilon         = 1.E-15
+
+   !  4. Following is information related to the physical constants.
+
+   !  These are the physical constants used within the model.
+
+! JM NOTE -- can we name this grav instead?
+   REAL    , PARAMETER :: g = 3.72  ! acceleration due to gravity (m {s}^-2)
+
+!#if ( NMM_CORE == 1 )
+!   REAL    , PARAMETER :: r_d          = 192.  ! gas constant m2 s-2 K-1
+!   REAL    , PARAMETER :: cp           = 844.6   ! r= 8.314511E+0 *1000.E+0/mugaz
+!#else
+!   REAL    , PARAMETER :: r_d          = 192.
+!   REAL    , PARAMETER :: cp           = 844.6
+!#endif
+
+   REAL    , PARAMETER :: r_d          = 191.
+   REAL    , PARAMETER :: cp           = 744.5
+
+   REAL    , PARAMETER :: r_v          = 461.6  ! gas constant for water vapor
+   REAL    , PARAMETER :: cv           = cp-r_d
+!   REAL    , PARAMETER :: cpv          = 4.4*r_v  ! cp / R = 4.4 ?? ! rcp = r/cp  ~0.256793
+   REAL    , PARAMETER :: cpv          = 3.9*r_v
+   REAL    , PARAMETER :: cvv          = cpv-r_v
+   REAL    , PARAMETER :: cvpm         = -cv/cp
+   REAL    , PARAMETER :: cliq         = 4190.
+   REAL    , PARAMETER :: cice         = 2106.
+   REAL    , PARAMETER :: psat         = 610.78
+   REAL    , PARAMETER :: rcv          = r_d/cv
+   REAL    , PARAMETER :: rcp          = r_d/cp
+   REAL    , PARAMETER :: rovg         = r_d/g
+   REAL    , PARAMETER :: c2           = cp * rcv
+
+   REAL    , PARAMETER :: p1000mb      = 610.
+   REAL    , PARAMETER :: t0           = 220. ! earth : 300 ... à remplacer par 220 ?
+   REAL    , PARAMETER :: p0           = 610.
+   REAL    , PARAMETER :: cpovcv       = cp/(cp-r_d)
+   REAL    , PARAMETER :: cvovcp       = 1./cpovcv
+   REAL    , PARAMETER :: rvovrd       = r_v/r_d
+
+   REAL    , PARAMETER :: reradius     = 1./3397200.	! Rayon de mars (m)  ~3397200 m
+
+!-------------------------------
+!****MARS
+! here some LMD stuff
+! >> TODO: possible modification with startfi.nc values ?
+   REAL    , PARAMETER :: wdaysec       = 88800.     ! duree du sol (s)  ~88775 s
+   REAL    , PARAMETER :: wmugaz        = 43.49      ! Masse molaire de l'atm (g.mol-1) ~43.49
+   REAL    , PARAMETER :: womeg         = 7.0721E-5  ! omega (rad.s-1)
+   REAL    , PARAMETER :: wyear_day     = 669        ! Duree de l'annee (sols) ~668.6
+   REAL    , PARAMETER :: wperiheli     = 206.66     ! Dist.min. soleil-mars (Mkm) ~206.66      
+   REAL    , PARAMETER :: waphelie      = 249.22     ! Dist.max. soleil-mars (Mkm) ~249.22
+   REAL    , PARAMETER :: wperi_day     = 485.       ! Date du perihelie (sols depuis printemps)
+   REAL    , PARAMETER :: wobliquit     = 25.2       ! Obliquite de la planete (deg) ~25.2
+   REAL    , PARAMETER :: wz0           = 1.e-2      ! surface roughness (m) ~0.01
+   REAL    , PARAMETER :: wlmixmin      = 30         ! longueur de melange ~100
+   REAL    , PARAMETER :: wemin_turb    = 1.e-6      ! energie minimale ~1.e-8 
+   REAL    , PARAMETER :: wemissiv      = 0.95       ! Emissivite du sol martien ~.95
+   REAL    , PARAMETER :: wemissiceN    = 0.95       ! Emissivite calotte nord
+   REAL    , PARAMETER :: wemissiceS    = 0.95       ! Emissivite calotte sud
+   REAL    , PARAMETER :: walbediceN    = 0.65        ! Albedo calotte nord !0.5
+   REAL    , PARAMETER :: walbediceS    = 0.65        ! Albedo calotte sud  !0.5
+   REAL    , PARAMETER :: wiceradiusN   = 100.e-6    ! mean scat radius of CO2 snow (north)
+   REAL    , PARAMETER :: wiceradiusS   = 100.e-6    ! mean scat radius of CO2 snow (south)
+   REAL    , PARAMETER :: wdtemisiceN   = 0.4  ! time scale for snow metamorphism (north) !2  
+   REAL    , PARAMETER :: wdtemisiceS   = 0.4  ! time scale for snow metamorphism (south) !2 
+#ifdef NEWPHYS
+   REAL    , PARAMETER :: wvolcapa      = 1.e6       ! volumetric capacity of soil (new soil model)
+#endif
+!
+!-------------------------------
+
+!   REAL    , PARAMETER :: asselin      = .025
+   REAL    , PARAMETER :: asselin      = .0
+   REAL    , PARAMETER :: cb           = 25.
+
+   REAL    , PARAMETER :: XLV0         = 3.15E6
+   REAL    , PARAMETER :: XLV1         = 2370.
+   REAL    , PARAMETER :: XLS0         = 2.905E6
+   REAL    , PARAMETER :: XLS1         = 259.532
+
+   REAL    , PARAMETER :: XLS          = 2.85E6
+   REAL    , PARAMETER :: XLV          = 2.5E6
+   REAL    , PARAMETER :: XLF          = 3.50E5
+
+   REAL    , PARAMETER :: rhowater     = 1000.
+   REAL    , PARAMETER :: rhosnow      = 100.
+   REAL    , PARAMETER :: rhoair0      = 0.02   !Earth Surface density: 1.217 kg/m3
+                                                !Mars Surface density: 0.020 kg/m3
+   REAL    , PARAMETER :: DEGRAD       = 3.1415926/180.
+   REAL    , PARAMETER :: DPD          = 360./669. ! longitude solaire d'un jour?
+
+   REAL    , PARAMETER ::  SVP1=0.6112
+   REAL    , PARAMETER ::  SVP2=17.67
+   REAL    , PARAMETER ::  SVP3=29.65
+   REAL    , PARAMETER ::  SVPT0=273.15
+   REAL    , PARAMETER ::  EP_1=R_v/R_d-1.
+   REAL    , PARAMETER ::  EP_2=R_d/R_v
+   REAL    , PARAMETER ::  KARMAN=0.4  ! von karman constant 
+   REAL    , PARAMETER ::  EOMEG=7.0721E-5 ! angular rotation rate (rad.s-1)
+   REAL    , PARAMETER ::  STBOLT=5.67051E-8 ! stefan-boltzmann constant
+
+                                      ! proportionality constants for eddy viscosity coefficient calc
+   REAL    , PARAMETER ::  c_s = .25  ! turbulence parameterization constant, for smagorinsky
+!   REAL    , PARAMETER ::  c_s = .1  ! turbulence parameterization constant, for smagorinsky
+   REAL    , PARAMETER ::  c_k = .15  ! turbulence parameterization constant, for TKE
+   REAL    , PARAMETER ::  prandtl = 1./3.0
+                                         ! constants for w-damping option
+!   REAL    , PARAMETER ::  w_alpha = 0.3 ! strength m/s/s
+   REAL    , PARAMETER ::  w_alpha = 2.0 ! strength m/s/s
+!   REAL    , PARAMETER ::  w_beta  = 1.0 ! activation cfl number
+   REAL    , PARAMETER ::  w_beta  = 0.8 ! activation cfl number
+
+       REAL , PARAMETER ::  pq0=379.90516
+       REAL , PARAMETER ::  epsq2=0.2
+!try both ?
+!       REAL , PARAMETER ::  epsq2=0.02
+       REAL , PARAMETER ::  a2=17.2693882
+       REAL , PARAMETER ::  a3=273.16
+       REAL , PARAMETER ::  a4=35.86
+       REAL , PARAMETER ::  epsq=1.e-12
+       REAL , PARAMETER ::  p608=rvovrd-1.
+!#if ( NMM_CORE == 1 )
+       REAL , PARAMETER ::  climit=1.e-20
+       REAL , PARAMETER ::  cm1=2937.4
+       REAL , PARAMETER ::  cm2=4.9283
+       REAL , PARAMETER ::  cm3=23.5518
+!       REAL , PARAMETER ::  defc=8.0
+!       REAL , PARAMETER ::  defm=32.0
+       REAL , PARAMETER ::  defc=0.0
+       REAL , PARAMETER ::  defm=99999.0
+       REAL , PARAMETER ::  epsfc=1./1.05
+       REAL , PARAMETER ::  epswet=0.0
+       REAL , PARAMETER ::  fcdif=1./3.
+!       REAL , PARAMETER ::  fcm=0.003
+       REAL , PARAMETER ::  fcm=0.00003
+!       REAL , PARAMETER ::  fcm=0.0
+       REAL , PARAMETER ::  gma=-r_d*(1.-rcp)*0.5
+       REAL , PARAMETER ::  p400=40000.0
+       REAL , PARAMETER ::  phitp=15000.0
+       REAL , PARAMETER ::  pi2=2.*3.1415926
+       REAL , PARAMETER ::  plbtm=105000.0
+       REAL , PARAMETER ::  plomd=64200.0
+       REAL , PARAMETER ::  pmdhi=35000.0
+       REAL , PARAMETER ::  q2ini=0.50
+       REAL , PARAMETER ::  rfcp=0.25/cp
+       REAL , PARAMETER ::  rhcrit_land=0.75
+       REAL , PARAMETER ::  rhcrit_sea=0.80
+       REAL , PARAMETER ::  rlag=14.8125
+       REAL , PARAMETER ::  rlx=0.90
+       REAL , PARAMETER ::  scq2=50.0
+       REAL , PARAMETER ::  slopht=0.001
+       REAL , PARAMETER ::  tlc=2.*0.703972477
+       REAL , PARAMETER ::  wa=0.15
+       REAL , PARAMETER ::  wght=0.35
+       REAL , PARAMETER ::  wpc=0.075
+       REAL , PARAMETER ::  z0land=0.10
+!       REAL , PARAMETER ::  z0max=0.01
+       REAL , PARAMETER ::  z0max=0.008
+       REAL , PARAMETER ::  z0sea=0.001
+!#endif
+
+
+ CONTAINS
+   SUBROUTINE init_module_model_constants
+   END SUBROUTINE init_module_model_constants
+ END MODULE module_model_constants
Index: /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F_ref
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F_ref	(revision 31)
+++ /trunk/mesoscale/LMD_MM_MARS/SRC/WRFV2/share/module_model_constants.F_ref	(revision 31)
@@ -0,0 +1,174 @@
+!WRF:MODEL_LAYER:CONSTANTS
+!
+!PLANET MARS
+
+ MODULE module_model_constants
+
+   !  2. Following are constants for use in defining real number bounds.
+
+   !  A really small number.
+
+   REAL    , PARAMETER :: epsilon         = 1.E-15
+
+   !  4. Following is information related to the physical constants.
+
+   !  These are the physical constants used within the model.
+
+! JM NOTE -- can we name this grav instead?
+   REAL    , PARAMETER :: g = 3.72  ! acceleration due to gravity (m {s}^-2)
+
+#if ( NMM_CORE == 1 )
+   REAL    , PARAMETER :: r_d          = 192.  ! gas constant m2 s-2 K-1
+   REAL    , PARAMETER :: cp           = 844.6   ! r= 8.314511E+0 *1000.E+0/mugaz
+#else
+   REAL    , PARAMETER :: r_d          = 192.
+   REAL    , PARAMETER :: cp           = 844.6
+#endif
+
+   REAL    , PARAMETER :: r_v          = 461.6  ! gas constant for water vapor
+   REAL    , PARAMETER :: cv           = cp-r_d
+   REAL    , PARAMETER :: cpv          = 4.4*r_v  ! cp / R = 4.4 ?? ! rcp = r/cp  ~0.256793
+   REAL    , PARAMETER :: cvv          = cpv-r_v
+   REAL    , PARAMETER :: cvpm         = -cv/cp
+   REAL    , PARAMETER :: cliq         = 4190.
+   REAL    , PARAMETER :: cice         = 2106.
+   REAL    , PARAMETER :: psat         = 610.78
+   REAL    , PARAMETER :: rcv          = r_d/cv
+   REAL    , PARAMETER :: rcp          = r_d/cp
+   REAL    , PARAMETER :: rovg         = r_d/g
+   REAL    , PARAMETER :: c2           = cp * rcv
+
+   REAL    , PARAMETER :: p1000mb      = 610.
+   REAL    , PARAMETER :: t0           = 220. ! earth : 300 ... à remplacer par 220 ?
+   REAL    , PARAMETER :: p0           = 610.
+   REAL    , PARAMETER :: cpovcv       = cp/(cp-r_d)
+   REAL    , PARAMETER :: cvovcp       = 1./cpovcv
+   REAL    , PARAMETER :: rvovrd       = r_v/r_d
+
+   REAL    , PARAMETER :: reradius     = 1./3397200.	! Rayon de mars (m)  ~3397200 m
+
+!-------------------------------
+!****MARS
+! here some LMD stuff
+! >> TODO: possible modification with startfi.nc values ?
+   REAL    , PARAMETER :: wdaysec       = 88800.     ! duree du sol (s)  ~88775 s
+   REAL    , PARAMETER :: wmugaz        = 43.49      ! Masse molaire de l'atm (g.mol-1) ~43.49
+   REAL    , PARAMETER :: womeg         = 7.0721E-5  ! omega (rad.s-1)
+   REAL    , PARAMETER :: wyear_day     = 669        ! Duree de l'annee (sols) ~668.6
+   REAL    , PARAMETER :: wperiheli     = 206.66     ! Dist.min. soleil-mars (Mkm) ~206.66      
+   REAL    , PARAMETER :: waphelie      = 249.22     ! Dist.max. soleil-mars (Mkm) ~249.22
+   REAL    , PARAMETER :: wperi_day     = 485.       ! Date du perihelie (sols depuis printemps)
+   REAL    , PARAMETER :: wobliquit     = 25.2       ! Obliquite de la planete (deg) ~25.2
+   REAL    , PARAMETER :: wz0           = 1.e-2      ! surface roughness (m) ~0.01
+   REAL    , PARAMETER :: wlmixmin      = 30         ! longueur de melange ~100
+   REAL    , PARAMETER :: wemin_turb    = 1.e-6      ! energie minimale ~1.e-8 
+   REAL    , PARAMETER :: wemissiv      = 0.95       ! Emissivite du sol martien ~.95
+   REAL    , PARAMETER :: wemissiceN    = 0.95       ! Emissivite calotte nord
+   REAL    , PARAMETER :: wemissiceS    = 0.95       ! Emissivite calotte sud
+   REAL    , PARAMETER :: walbediceN    = 0.65        ! Albedo calotte nord !0.5
+   REAL    , PARAMETER :: walbediceS    = 0.65        ! Albedo calotte sud  !0.5
+   REAL    , PARAMETER :: wiceradiusN   = 100.e-6    ! mean scat radius of CO2 snow (north)
+   REAL    , PARAMETER :: wiceradiusS   = 100.e-6    ! mean scat radius of CO2 snow (south)
+   REAL    , PARAMETER :: wdtemisiceN   = 0.4  ! time scale for snow metamorphism (north) !2  
+   REAL    , PARAMETER :: wdtemisiceS   = 0.4  ! time scale for snow metamorphism (south) !2 
+#ifdef NEWPHYS
+   REAL    , PARAMETER :: wvolcapa      = 1.e6       ! volumetric capacity of soil (new soil model)
+#endif
+!
+!-------------------------------
+
+!   REAL    , PARAMETER :: asselin      = .025
+   REAL    , PARAMETER :: asselin      = .0
+   REAL    , PARAMETER :: cb           = 25.
+
+   REAL    , PARAMETER :: XLV0         = 3.15E6
+   REAL    , PARAMETER :: XLV1         = 2370.
+   REAL    , PARAMETER :: XLS0         = 2.905E6
+   REAL    , PARAMETER :: XLS1         = 259.532
+
+   REAL    , PARAMETER :: XLS          = 2.85E6
+   REAL    , PARAMETER :: XLV          = 2.5E6
+   REAL    , PARAMETER :: XLF          = 3.50E5
+
+   REAL    , PARAMETER :: rhowater     = 1000.
+   REAL    , PARAMETER :: rhosnow      = 100.
+   REAL    , PARAMETER :: rhoair0      = 0.02   !Earth Surface density: 1.217 kg/m3
+                                                !Mars Surface density: 0.020 kg/m3
+   REAL    , PARAMETER :: DEGRAD       = 3.1415926/180.
+   REAL    , PARAMETER :: DPD          = 360./669. ! longitude solaire d'un jour?
+
+   REAL    , PARAMETER ::  SVP1=0.6112
+   REAL    , PARAMETER ::  SVP2=17.67
+   REAL    , PARAMETER ::  SVP3=29.65
+   REAL    , PARAMETER ::  SVPT0=273.15
+   REAL    , PARAMETER ::  EP_1=R_v/R_d-1.
+   REAL    , PARAMETER ::  EP_2=R_d/R_v
+   REAL    , PARAMETER ::  KARMAN=0.4  ! von karman constant 
+   REAL    , PARAMETER ::  EOMEG=7.0721E-5 ! angular rotation rate (rad.s-1)
+   REAL    , PARAMETER ::  STBOLT=5.67051E-8 ! stefan-boltzmann constant
+
+                                      ! proportionality constants for eddy viscosity coefficient calc
+   REAL    , PARAMETER ::  c_s = .25  ! turbulence parameterization constant, for smagorinsky
+!   REAL    , PARAMETER ::  c_s = .1  ! turbulence parameterization constant, for smagorinsky
+   REAL    , PARAMETER ::  c_k = .15  ! turbulence parameterization constant, for TKE
+   REAL    , PARAMETER ::  prandtl = 1./3.0
+                                         ! constants for w-damping option
+!   REAL    , PARAMETER ::  w_alpha = 0.3 ! strength m/s/s
+   REAL    , PARAMETER ::  w_alpha = 2.0 ! strength m/s/s
+!   REAL    , PARAMETER ::  w_beta  = 1.0 ! activation cfl number
+   REAL    , PARAMETER ::  w_beta  = 0.8 ! activation cfl number
+
+       REAL , PARAMETER ::  pq0=379.90516
+       REAL , PARAMETER ::  epsq2=0.2
+!try both ?
+!       REAL , PARAMETER ::  epsq2=0.02
+       REAL , PARAMETER ::  a2=17.2693882
+       REAL , PARAMETER ::  a3=273.16
+       REAL , PARAMETER ::  a4=35.86
+       REAL , PARAMETER ::  epsq=1.e-12
+       REAL , PARAMETER ::  p608=rvovrd-1.
+!#if ( NMM_CORE == 1 )
+       REAL , PARAMETER ::  climit=1.e-20
+       REAL , PARAMETER ::  cm1=2937.4
+       REAL , PARAMETER ::  cm2=4.9283
+       REAL , PARAMETER ::  cm3=23.5518
+!       REAL , PARAMETER ::  defc=8.0
+!       REAL , PARAMETER ::  defm=32.0
+       REAL , PARAMETER ::  defc=0.0
+       REAL , PARAMETER ::  defm=99999.0
+       REAL , PARAMETER ::  epsfc=1./1.05
+       REAL , PARAMETER ::  epswet=0.0
+       REAL , PARAMETER ::  fcdif=1./3.
+!       REAL , PARAMETER ::  fcm=0.003
+       REAL , PARAMETER ::  fcm=0.00003
+!       REAL , PARAMETER ::  fcm=0.0
+       REAL , PARAMETER ::  gma=-r_d*(1.-rcp)*0.5
+       REAL , PARAMETER ::  p400=40000.0
+       REAL , PARAMETER ::  phitp=15000.0
+       REAL , PARAMETER ::  pi2=2.*3.1415926
+       REAL , PARAMETER ::  plbtm=105000.0
+       REAL , PARAMETER ::  plomd=64200.0
+       REAL , PARAMETER ::  pmdhi=35000.0
+       REAL , PARAMETER ::  q2ini=0.50
+       REAL , PARAMETER ::  rfcp=0.25/cp
+       REAL , PARAMETER ::  rhcrit_land=0.75
+       REAL , PARAMETER ::  rhcrit_sea=0.80
+       REAL , PARAMETER ::  rlag=14.8125
+       REAL , PARAMETER ::  rlx=0.90
+       REAL , PARAMETER ::  scq2=50.0
+       REAL , PARAMETER ::  slopht=0.001
+       REAL , PARAMETER ::  tlc=2.*0.703972477
+       REAL , PARAMETER ::  wa=0.15
+       REAL , PARAMETER ::  wght=0.35
+       REAL , PARAMETER ::  wpc=0.075
+       REAL , PARAMETER ::  z0land=0.10
+!       REAL , PARAMETER ::  z0max=0.01
+       REAL , PARAMETER ::  z0max=0.008
+       REAL , PARAMETER ::  z0sea=0.001
+!#endif
+
+
+ CONTAINS
+   SUBROUTINE init_module_model_constants
+   END SUBROUTINE init_module_model_constants
+ END MODULE module_model_constants
Index: /trunk/mesoscale/LMD_MM_MARS/makemeso
===================================================================
--- /trunk/mesoscale/LMD_MM_MARS/makemeso	(revision 30)
+++ /trunk/mesoscale/LMD_MM_MARS/makemeso	(revision 31)
@@ -354,5 +354,5 @@
   if [[ "${phys}" == "newphys_" ]]
   then
-     sed s+"CFLAGS          =       "+"CFLAGS          =       -DNEWPHYS "+g configure.wrf > yeah 
+     sed s+"ARCHFLAGS       =       "+"ARCHFLAGS       =       -DNEWPHYS "+g configure.wrf > yeah 
      mv -f yeah configure.wrf
   fi
