Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/LMD_MMM_d1_10km_PSFC_UV_Ls244.8.sh
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/LMD_MMM_d1_10km_PSFC_UV_Ls244.8.sh	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/LMD_MMM_d1_10km_PSFC_UV_Ls244.8.sh	(revision 559)
@@ -0,0 +1,1 @@
+/donnees/aslmd/MODELES//MESOSCALE//LMD_MM_MARS//SRC/PYTHON/pp.py -f wrfout_d01_9999-09-09_09:36:00 -v PSFC -W -s 2 -i 4 -l 1. --div 25 -z 33 -m 300. -M 500. -S png 
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/callphys.def
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/callphys.def	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/callphys.def	(revision 559)
@@ -0,0 +1,94 @@
+General options
+~~~~~~~~~~~~~~~
+tracer    (Run with or without tracer transport ?)
+F
+diurnal   (Diurnal cycle ?  if diurnal=F, diurnal averaged solar heating)
+T
+season    (Seasonal cycle ? if season=F, Ls stays constant like in "start")
+T
+lwrite    (want some more output on the screen ?) 
+F
+stats     (Saving statistics in file "cumul" ?) -- noWRF
+F
+calleofdump (Saving EOF profiles in file "profiles" for Climate Database ?) -- noWRF
+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 
+1        (=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 ?)
+F T
+callsoil  (call thermal conduction in the soil ?)
+T
+calllott  (call Lott's gravity wave/subgrid topography scheme ?) -- noWRF
+F
+Radiative transfer options :
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+iradia    (the rad.transfer is computed every "iradia" physical timestep)
+1 37 75 
+callg2d   (Output of the exchange coefficient mattrix ? for diagnostic only)
+F
+rayleigh  (Rayleigh scattering : should be =F for now)
+F
+Tracer (dust water, ice and/or chemical species) options (use if tracer=T) :
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dustbin    (DUST: Transported dust ? (if >0, uses q(1) to q(dustbin))
+0
+active     (DUST: Radiatively active dust ? (uses q(1) to q(dustbin))
+F
+doubleq    (DUST: needs dustbin=1, use mass (q(1)) and number (q(2)) mixing ratio to predict dust size ?)
+F
+lifting    (DUST: lifted by GCM surface winds ?)
+F
+dustdevil  (DUST: lifted by dust devils ?)
+F
+scavenging (DUST: Scavenging by CO2 snowfall ?)
+F
+sedimentation (DUST/WATERICE: Gravitationnal sedimentation ?)
+F
+iceparty   (WATERICE: Water cycle includes water ice mixing ratio q(nqmx-1))
+F
+activice   (WATERICE: Radiatively active transported atmospheric water ice ?)
+F
+water      (WATER: Compute water cycle using q(nqmx) )
+F
+caps       (WATER: put the current permanent caps at both poles)
+F
+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/DEF/TESTCRATER/dustopacity.def
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/dustopacity.def	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/dustopacity.def	(revision 559)
@@ -0,0 +1,3 @@
+0.3
+
+0.6
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/levels
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/levels	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/levels	(revision 559)
@@ -0,0 +1,1 @@
+link vert_level/levels
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/namelist.input
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/namelist.input	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/namelist.input	(revision 559)
@@ -0,0 +1,65 @@
+ &time_control
+ run_seconds         = 14800   !! Duration of simulation                       [LES: !=0]  
+ history_interval_s  = 600     !! Frequency of outputs in seconds              [LES: !=0]
+ frames_per_outfile  = 6       !! Size of time dimension in files
+ start_year          = 9999    !! Idealized simulation                         [LES: 9999]
+ end_year            = 9999    !! Idealized simulation                         [LES: 9999]
+ restart             = .false. !! (*) Output restart files ?
+ restart_interval    = 8880    !! (*) Frequency of output restart files ?
+ io_form_history     = 2       !! (*) Choice of NETCDF for ouputs
+ io_form_restart     = 2       !! (*) Choice of NETCDF for ouputs
+ io_form_input       = 2       !! (*) Choice of NETCDF for ouputs
+ io_form_boundary    = 2       !! (*) Choice of NETCDF for ouputs
+ debug_level         = 0       !! (*) Verbose level
+ input_from_file     = F,      !! (n)(p2) Initialize a given domain with an input file
+ /
+
+ &domains
+ ztop                = 25000.  !! Altitude above ground of model top           [LES: !=0]       
+ time_step           = 10      !! Dynamical timestep
+ time_step_fract_num = 0       !! Additional fraction to time_step: numerator
+ time_step_fract_den = 1       !! Additional fraction to time_step: denominator
+ dx                  = 10000   !! (p2) Horizontal resolution
+ dy                  = 10000   !! (p2) Horizontal resolution (should be equal to dx)
+ e_we                = 61      !! (r)(p2) Number of longitude grid points
+ e_sn                = 61      !! (r)(p2) Number of latitude grid points
+ e_vert              = 61      !! (r)(p2) Number of vertical levels
+ /
+
+ &physics
+ init_TI             = 250.    !! Define constant thermal inertia value        [LES: !=0]
+ init_AL             = 0.25    !! Define constant albedo value                 [LES: !=0]
+ !!
+ !! OPTIONAL
+ !!
+ radt                = 1       !! Ratio between physical and dynamical time step
+ mars                = 0       !! (r)(p2) Configuration of tracers:
+                               !!         0: no tracers, 1: water vapor + ice, 2: dust
+ init_MU             = 0.001   !! (p3) Multiply ini & bdy zonal wind by init_U
+ init_MV             = 0.001   !! (p3) Multiply ini & bdy meridional wind by init_V
+ init_U = 15.
+ init_V = 0.
+ /
+
+ &dynamics
+ pert_coriolis         = T        !! Coriolis only acts on wind perturbation (idealized)
+ /
+
+ &bdy_control
+ specified = F,              !! (n)(p3) Boundary conditions specified by GCM
+ open_xs = T,                !! (p3) Open boundary conditions @ western boundary
+ open_xe = T,                !! (p3) Open boundary conditions @ eastern boundary
+ open_ys = T,                !! (p3) Open boundary conditions @ southern boundary
+ open_ye = T,                !! (p3) Open boundary conditions @ northern boundary
+ /
+
+ &grib2
+ /
+
+ &fdda
+ /
+
+ &namelist_quilt              !! (*)
+ nio_tasks_per_group = 0,     !! (*)
+ nio_groups = 1,              !! (*)
+ /                            !! (*)
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/levels
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/levels	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/levels	(revision 559)
@@ -0,0 +1,61 @@
+      1.00000
+     0.999244
+     0.996981
+     0.993219
+     0.987974
+     0.981268
+     0.973127
+     0.963586
+     0.952684
+     0.940463
+     0.926975
+     0.912271
+     0.896411
+     0.879456
+     0.861471
+     0.842525
+     0.822690
+     0.802037
+     0.780644
+     0.758585
+     0.735938
+     0.712781
+     0.689192
+     0.665247
+     0.641024
+     0.616596
+     0.592038
+     0.567421
+     0.542814
+     0.518284
+     0.493895
+     0.469707
+     0.445778
+     0.422161
+     0.398907
+     0.376063
+     0.353671
+     0.331770
+     0.310396
+     0.289579
+     0.269347
+     0.249722
+     0.230726
+     0.212374
+     0.194679
+     0.177649
+     0.161291
+     0.145608
+     0.130600
+     0.116264
+     0.102595
+    0.0895848
+    0.0772241
+    0.0655012
+    0.0544028
+    0.0439138
+    0.0340178
+    0.0246973
+    0.0159339
+   0.00770813
+      0.00000
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/levspe.pro
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/levspe.pro	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/levspe.pro	(revision 559)
@@ -0,0 +1,95 @@
+pro levspe
+
+
+;
+; TWEAK PARAM
+;
+@param.idl
+
+
+;;
+;; HARD PARAM
+;;
+;!p.multi=[0,2,2]
+set_plot, 'ps'
+device, file='plot.ps';, /landscape
+tinv=5      ;; intervalle plot
+hache=10.   ;; scale height
+psurf=610.  ;; pression surface -- ne change pas les eta levels
+
+;
+; PRELIM
+;
+param = 'e_vert='+string(nlev,'(I0)')+' ztop='+string(altmax,'(I0)')+' e='+string(epsilon,'(I0)')+' c='+string(elong_cos,'(F4.2)')
+;
+epsilon = epsilon / (nlev-1) / 100.
+;
+exposant = !pi/nlev/elong_cos
+;
+; alpha (plus grand ecart en km) est determine pour que max(altitudes)=altmax
+;
+alpha =  altmax / ( (nlev-1)/2. - sin(2.*exposant*(nlev-1))/4./exposant + epsilon*(nlev-1)^2/2. )
+print, 'alpha (km)', alpha
+;
+x=[findgen(nlev)]
+
+
+;
+; CALC
+;
+ecart = alpha - alpha*cos(exposant*x)^2 + epsilon*alpha*x 
+altitudes = alpha*x/2. - alpha*sin(2.*exposant*x)/4./exposant + epsilon*alpha*x^2/2. 
+logpressions = alog(psurf) - altitudes/hache
+pressions=exp(logpressions)
+ptop=psurf*exp(-altmax/hache) & print, ptop  
+etas=(pressions-ptop)/(psurf-ptop)
+etas(nlev-1)=0.
+
+
+;
+; PRINT
+;
+print, 'eta levels'
+print, etas
+;
+print, 'pressure'
+pi=etas*(psurf-ptop)+ptop
+print, pi
+print, 'ptop', ptop
+;
+print, 'pseudo-altitude (meters)'
+pseudo=10.*alog(610./pi)
+print, 1000.*10.*alog(610./pi)
+
+
+;
+; PLOT
+;
+plot, x, etas, title='ETA '+param, xtickinterval=tinv 
+!psym=7 & oplot, x, etas
+!psym=0
+plot, x, pseudo, title='ALTITUDE (km)', xtickinterval=tinv
+!psym=7 & oplot, x, pseudo 
+!psym=0
+plot, x, pi, title='PRESSURE (Pa)', ylog=1, xtickinterval=tinv
+!psym=7 & oplot, x, pi 
+!psym=0
+;
+diff=pseudo - shift(pseudo,-1) & diff=-diff(0:nlev-2)
+;
+plot, x, ecart, title='ECART (km, L=theor x=calc)', xtickinterval=tinv
+!psym=7 & oplot, x, diff
+!psym=0
+
+
+;
+; PRINT
+;
+;openw, 1, 'levels'
+for i=0, nlev-1 do begin
+;        printf, 1, etas(i);,',' 
+print, etas(i);,','
+endfor
+close, 1
+
+end
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/param.idl
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/param.idl	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/param.idl	(revision 559)
@@ -0,0 +1,89 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=61
+altmax=120.    ;; altitude maximale en km
+;;
+epsilon=30.    ;; augmentation totale en pourcentage de l'ecart max (modifie inflexion)
+               ;; -- epsilon = 0   : point d'inflexion parfaitement plat
+               ;; -- epsilon = 100 : pas de point d'inflexion
+elong_cos=2./3.  ;; plus petit rapproche l'inflexion du sol ;; 2/3 parfait pour un nlev divisible par 3
+elong_cos=0.64
+;; plus petit ajoute des inflexions...
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;nlev=61
+;altmax=150.
+;epsilon=30.
+;elong_cos=2./3.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=61
+altmax=150.
+epsilon=100.
+elong_cos=2.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=201
+altmax=150.
+epsilon=0.
+elong_cos=2.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=101
+altmax=80. ;;100.
+epsilon=100.
+elong_cos=60.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=201
+altmax=150.
+epsilon=0.
+elong_cos=2. ;; pas inflexion  
+elong_cos=2.5
+elong_cos=3.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;nlev=201
+;altmax=220.
+;epsilon=0.
+;elong_cos=2. ;; pas inflexion
+;elong_cos=2.5
+;elong_cos=3.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=61
+altmax=60.
+epsilon=0.
+elong_cos=2. ;; pas inflexion
+elong_cos=2.5
+elong_cos=3.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+epsilon=100.
+elong_cos=60.
+
+;altmax=100.
+altmax=130.
+altmax=150.
+
+;;; ici OK avec simu
+
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=201
+altmax=150.
+epsilon=0.
+elong_cos=3.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+nlev=151
+altmax=150.
+epsilon=100.
+elong_cos=60.
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;nlev=61
+nlev = 181
+altmax = 180.
+
+
+
+nlev=61
+altmax=25.
+epsilon=100
+elong_cos=100.
+
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/simple.pro
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/simple.pro	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/vert_level/simple.pro	(revision 559)
@@ -0,0 +1,94 @@
+pro simple, z_scale
+
+
+nlev=101 
+;
+; z-scale doit etre suffisamment bas pour que l avant-dernier niveau soit assez haut
+; mais pas trop sinon instable (il faut avoir des niveaux espaces en haut semble-t-il)
+;
+z_scale = .25
+;z_scale = .15 marche pas (niveaux pas assez espaces en haut ?)
+ptop = 0.2
+npbl = 5
+	;z_scale = .15
+	;ptop = 0.05 ;marche pas
+
+z_scale = .20 & ptop = 0.1 & npbl = 6
+
+
+psurf=610.
+
+pbl = [1.0000,0.9995,0.9980,0.9950,0.9850,0.9700,0.9400,0.9000] 
+x=[findgen(nlev)]
+nlev = nlev - npbl + 1 
+
+etas=[findgen(nlev)]
+FOR k=0, nlev-1 DO BEGIN
+     etas(k) = (exp(-float(k)/float(nlev-1)/z_scale) - exp(-1./z_scale)) / (1.-exp(-1./z_scale)) 
+ENDFOR
+etas(nlev-1)=0.
+
+etas = [pbl(0:npbl-1),etas(1:nlev-1)]
+nlev = nlev + npbl - 1 
+print, nlev
+
+;pressions = ptop + etas * (psurf-ptop)
+;zzz = 10. * alog ( psurf / pressions )
+;print, zzz
+
+;
+; HARD PARAM
+;
+!p.multi=[0,2,2]
+;set_plot, 'ps'
+;device, file='plot.ps', /landscape
+
+
+
+;
+; PRINT
+;
+print, 'eta levels'
+print, etas
+;
+print, 'pressure'
+pi=etas*(psurf-ptop)+ptop
+print, pi
+print, 'ptop', ptop
+;
+print, 'pseudo-altitude'
+pseudo=10.*alog(610./pi)
+print, 1000.*10.*alog(610./pi)
+
+
+;
+; PLOT
+;
+plot, x, etas, title='ETA', xtickinterval=tinv 
+!psym=7 & oplot, x, etas
+!psym=0
+plot, x, pseudo, title='ALTITUDE (km)', xtickinterval=tinv
+!psym=7 & oplot, x, pseudo 
+!psym=0
+plot, x, pi, title='PRESSURE (Pa)', ylog=1, xtickinterval=tinv
+!psym=7 & oplot, x, pi 
+!psym=0
+;
+diff=pseudo - shift(pseudo,-1) & diff=-diff(0:nlev-2)
+;
+plot, x, diff, title='ECART (km, L=theor x=calc)', xtickinterval=tinv
+!psym=7 & oplot, x, diff
+!psym=0
+
+;stop
+
+;
+; PRINT
+;
+openw, 1, 'levels'
+for i=0, nlev-1 do begin
+        printf, 1, etas(i);,',' 
+endfor
+close, 1
+
+end
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/ze_hill
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/ze_hill	(revision 559)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/TESTCRATER/ze_hill	(revision 559)
@@ -0,0 +1,38 @@
+-4000.
+2.0
+
+1000.
+4.0
+
+4000.
+2.0
+
+1000.
+4.0
+
+6000.
+4.0
+
+4000.
+8.0
+
+50.
+8.0
+
+4000.
+8.0
+
+500.
+5.0
+
+4000.
+8.0
+
+100.
+4.0
+
+1000.
+2.0
+
+2000.
+6.0
Index: /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/newphys_THARSIS_WATER/callphys.def.mars11_radact_therm
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/newphys_THARSIS_WATER/callphys.def.mars11_radact_therm	(revision 558)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SIMU/DEF/newphys_THARSIS_WATER/callphys.def.mars11_radact_therm	(revision 559)
@@ -54,5 +54,5 @@
 callnlte = .false.
 # call CO2 NIR absorption ?   matters only if callrad=T
-callnirco2 = .true.
+callnirco2 = .false. 
 # call turbulent vertical diffusion ?
 calldifv  = .true.
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F	(revision 558)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/dyn_em/module_initialize_quarter_ss.F	(revision 559)
@@ -264,5 +264,5 @@
          grid%cosa(i,j)     = 1.
          grid%e(i,j)        = 0.
-         grid%f(i,j)        = 0.         !! MARS: put coriolis here if needed
+         grid%f(i,j)        = 0.
 
       END DO
@@ -401,11 +401,30 @@
       !    grid%ht(i,j) = hm1*exp(-(( float(i-icm)/xa1)**2))   &
       !               *( (cos(pii*float(i-icm)/xal1))**2 )
-    write(6,*) '3D hill or crater '
-    write(6,*) 'height, width ',hm,xa
-    grid%ht(i,j) = alt_input + hm/(1.+(float(i-icm)/xa)**2+(float(j-jcm)/xa)**2)
+    IF (hm .gt. 0.) THEN
+      write(6,*) '3D hill. height, width: ',hm,xa
+      write(6,*) 'input sounding is out of the mountain'
+      grid%ht(i,j) = alt_input + hm/(1.+(float(i-icm)/xa)**2+(float(j-jcm)/xa)**2)
+    ELSE IF (hm .lt. 0.) THEN
+      write(6,*) '3D crater. height, width: ',hm,xa
+      write(6,*) 'input sounding is at the bottom of crater'
+      grid%ht(i,j) = (alt_input - hm) + hm/(1.+(float(i-icm)/xa)**2+(float(j-jcm)/xa)**2)
+      !! AS: cannot use same formula as hill because it would force ideal.exe to extrapolate
+      !!       which is not possible given how interp_0 is written
+    ELSE
+      write(6,*) 'Nothing. Height is 0. Flat topography'
+      grid%ht(i,j) = alt_input
+    ENDIF
   ELSE
     write(6,*) 'linear slope '
     write(6,*) 'height ',hm
-    grid%ht(i,j) = alt_input + hm * float(i)  
+    IF (hm .gt. 0.) THEN  
+      grid%ht(i,j) = alt_input + hm * float(i)  
+    ELSE IF (hm .lt. 0.) THEN
+      !! see above, crater case
+      grid%ht(i,j) = (alt_input - hm) + hm * float(i)
+    ELSE
+      write(6,*) 'Nothing. Height is 0. Flat topography'
+      grid%ht(i,j) = alt_input
+    ENDIF
   ENDIF
       !!!3D crater    
@@ -437,6 +456,12 @@
     grid%mars_cice(i,j)=0.
     grid%mars_wice(i,j)=0.
+    write(6,*) 'NOTE TO SELF. slpx and slpy set to 0 which means no slope insolation.'
     grid%slpx(i,j) = 0.
     grid%slpy(i,j) = 0.
+    !!! COMMENT THE LINES BELOW IF YOU DON'T WANT CORIOLIS TERMS
+    !!! cf. doc WRF2008 page 11 for e and f expressions
+    grid%e(i,j) = 2. * EOMEG * COS(pi*lat_input/180.)
+    grid%f(i,j) = 2. * EOMEG * SIN(pi*lat_input/180.)
+    write(6,*) 'CALCULATE CORIOLIS TERM', grid%f(i,j),grid%e(i,j)
 !!MARS
   ENDDO
