Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90	(revision 613)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PREP_MARS/readmeteo.F90	(revision 665)
@@ -662,5 +662,5 @@
     if (ierr.ne.NF_NOERR) then
       write(*,*) "...No ccn - CCN mass set to 0"
-      dustfile(:,:,:,:)=0.
+      ccnfile(:,:,:,:)=0.
     else
       ierr=NF_GET_VAR_REAL(nid,nvarid,ccnfile)
@@ -671,5 +671,5 @@
     if (ierr.ne.NF_NOERR) then
       write(*,*) "...No ccnN - CCN number set to 0"
-      dustnfile(:,:,:,:)=0.
+      ccnnfile(:,:,:,:)=0.
     else
       ierr=NF_GET_VAR_REAL(nid,nvarid,ccnnfile)
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM	(revision 613)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/Registry/Registry.EM	(revision 665)
@@ -87,27 +87,27 @@
 state  real  THETA     ij   misc  1  -  rd   "THETA"     "SLOPE INCLINATION"               "deg"     #SAVEMARS2 theta_sl
 state  real  PSI       ij   misc  1  -  rd   "PSI"       "SLOPE ORIENTATION"               "deg"     #SAVEMARS2 psi_sl
-state  real  TAU_DUST  ij   misc  1  -  rd   "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""        #SAVEMARS2 tauref
+state  real  TAU_DUST  ij   misc  1  -  rhd   "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""        #SAVEMARS2 tauref
 state  real  SWDOWNZ   ij   misc  1  -  rd   "SWDOWNZ"   "DOWNWARD SW FLUX AT SURFACE"     "W m-2"   #SAVEMARS2 fluxsurf_sw_tot      
 state  real  LWDOWNZ   ij   misc  1  -  rd   "LWDOWNZ"   "DOWNWARD LW FLUX AT SURFACE"     "W m-2"   #SAVEMARS2 fluxsurf_lw
 state  real  SWUP      ij   misc  1  -  rd   "SWUP"      "UPWARD SW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_sw_tot    
 state  real  LWUP      ij   misc  1  -  rd   "LWUP"      "UPWARD LW FLUX AT TOP"           "W m-2"   #SAVEMARS2 fluxtop_lw
-state  real  MTOT      ij   misc  1  -  rd   "MTOT"      "TOTAL MASS WATER VAPOR in pmic"  "pmic"    #SAVEMARS2 mtot
-state  real  ICETOT    ij   misc  1  -  rd   "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"  #SAVEMARS2 icetot
-state  real  RAVE      ij   misc  1  -  rd   "RAVE"      "MEAN ICE RADIUS"                 "m"       #SAVEMARS2 rave
+state  real  MTOT      ij   misc  1  -  rhd   "MTOT"      "TOTAL MASS WATER VAPOR in pmic"  "pmic"    #SAVEMARS2 mtot
+state  real  ICETOT    ij   misc  1  -  rhd   "ICETOT"    "TOTAL MASS WATER ICE"            "kg m-2"  #SAVEMARS2 icetot
+state  real  RAVE      ij   misc  1  -  rhd   "RAVE"      "MEAN ICE RADIUS"                 "m"       #SAVEMARS2 rave
 state  real  RICE      ikj  misc  1  -  rd   "RICE"      "ICE RADIUS"                      "m"       #SAVEMARS3 rice
 state  real  HR_SW     ikj  misc  1  -  rd   "HR_SW"     "HEATING RATE SW"                 "K/s"     #SAVEMARS3 zdtsw
 state  real  HR_LW     ikj  misc  1  -  rd   "HR_LW"     "HEATING RATE LW"                 "K/s"     #SAVEMARS3 zdtlw
 state  real  HR_SH     ikj  misc  1  -  rd   "HR_SH"     "HEATING RATE sens. heat"         "K/s"     #SAVEMARS3 zdtdif
-state  real  QSURFICE  ij   misc  1  -  rd   "QSURFICE"  "WATER ICE AT SURFACE"            "kg m-2"  #SAVEMARS2 qsurfice
-state  real  RDUST     ikj  misc  1  -  rd   "RDUST"     "DUST RADIUS"                     "m"       #SAVEMARS3 rdust
+state  real  QSURFICE  ij   misc  1  -  rhd   "QSURFICE"  "WATER ICE AT SURFACE"            "kg m-2"  #SAVEMARS2 qsurfice
+state  real  RDUST     ikj  misc  1  -  rhd   "RDUST"     "DUST RADIUS"                     "m"       #SAVEMARS3 rdust
 state  real  HR_NIR    ikj  misc  1  -  rd   "HR_NIR"    "HEATING RATE nirco2"             "K/s"     #SAVEMARS3 zdtnirco2
 state  real  HR_NLTE   ikj  misc  1  -  rd   "HR_NLTE"   "HEATING RATE nlte"               "K/s"     #SAVEMARS3 zdtnlte
-state  real  ALBBARE   ij   misc  1  -  rd   "ALBBARE"   "SOIL ALBEDO"                     ""        #SAVEMARS2 albedodat
-state  real  VMR_ICE   ikj  misc  1  -  rd   "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"     #SAVEMARS3 vmr
-state  real  TAU_ICE   ij   misc  1  -  rd   "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""        #SAVEMARS2 tauTES
+state  real  ALBBARE   ij   misc  1  -  rhd   "ALBBARE"   "SOIL ALBEDO"                     ""        #SAVEMARS2 albedodat
+state  real  VMR_ICE   ikj  misc  1  -  rhd   "VMR_ICE"   "VOL. MIXING RATIO ICE"           "ppm"     #SAVEMARS3 vmr
+state  real  TAU_ICE   ij   misc  1  -  rhd   "TAU_ICE"   "CLOUD OD at 825 cm-1 TES"        ""        #SAVEMARS2 tauTES
 state  real  PDTZ      ikj  misc  1  -  rd   "PDT"       "TEMP TENDENCY"                   "K s-1"   #SAVEMARS3 pdt
-state  real  ZMAX_TH   ij   misc  1  -  rd   "ZMAX_TH"   "MAXIMUM LEVEL REACHED IN TH"     "m"       #SAVEMARS2 zmax_th
+state  real  ZMAX_TH   ij   misc  1  -  rhd   "ZMAX_TH"   "MAXIMUM LEVEL REACHED IN TH"     "m"       #SAVEMARS2 zmax_th
 state  real  HFMAX_TH  ij   misc  1  -  rd   "HFMAX_TH"  "MAXIMUM TH HEAT FLUX"            "m.K/s"   #SAVEMARS2 hfmax_th
-state  real  WSTAR     ij   misc  1  -  rd   "WSTAR"     "FREE CONVECTION VELOCITY FROM TH" "m/s"    #SAVEMARS2 wstar
+state  real  WSTAR     ij   misc  1  -  rhd   "WSTAR"     "FREE CONVECTION VELOCITY FROM TH" "m/s"    #SAVEMARS2 wstar
 state  real  Z0SET     ij   misc  1  -  rd   "Z0SET"     "SET SURFACE ROUGHNESS"           "m"       #SAVEMARS2 z0
 
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/aeropacity_use.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/aeropacity_use.F	(revision 613)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/aeropacity_use.F	(revision 665)
@@ -189,5 +189,6 @@
         msolsir(1:nlayer,1:naerkind)=0
         mqextsqabs(1:nlayer,1:naerkind)=0
-        WRITE(*,*) "Typical profiles of solsir and Qext/Qabs(IR):"
+        WRITE(*,*) "Typical profiles of Qext(vis)/Qext(IR)"
+        WRITE(*,*) "  and Qext(IR)/Qabs(IR):"
         DO iaer = 1, naerkind ! Loop on aerosol kind
           WRITE(*,*) "Aerosol # ",iaer
@@ -210,4 +211,9 @@
 !       otherwise default value read from starfi.nc file will be used)
         call getin("tauvis",tauvis)
+
+!       Some information about the water cycle
+        IF (water) THEN
+          write(*,*) "water_param CCN reduc. fac. ", ccn_factor
+        ENDIF
 
         firstcall=.false.
@@ -445,5 +451,4 @@
             taudusttmp(ig) = taudusttmp(ig) + 
      &                       aerosol(ig,l,iaerdust(iaer))
-
           ENDDO
         ENDDO
@@ -590,5 +595,4 @@
 !         ccn(ig,l) = ccn(ig,l) / 8.
 c -----------------------------------------------------------------
-       write(*,*) "water_param CCN reduc. fac. ", ccn_factor
        DO l=1,nlayer
          DO ig=1,ngrid
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/callradite_use.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/callradite_use.F	(revision 613)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/callradite_use.F	(revision 665)
@@ -264,8 +264,13 @@
 c        PLEASE MAKE SURE that you set up the right number of
 c          scatterers in dimradmars.h (naerkind);
-c          name_iaer(1) = "dust_conrath"   !! poussiere classique
-          name_iaer(1) = "dust_doubleq"
-cc        name_iaer(2) = "dust_submicron" !! JB: experimental
-c          name_iaer(2) = "h2o_ice"
+          name_iaer(1) = "dust_conrath"   !! default choice is good old Conrath profile
+          IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
+          if (nqmx.gt.1) then
+           ! trick to avoid problems compiling with 1 tracer
+           ! and picky compilers who know name_iaer(2) is out of bounds
+           j=2
+          IF (water.AND.activice) name_iaer(j) = "h2o_ice"      !! radiatively-active clouds
+          IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
+          endif ! of if (nqmx.gt.1)
 c        ----------------------------------------------------------
 
@@ -389,6 +394,7 @@
 c     Computing aerosol optical depth in each layer:
       CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &            pq,ccn,tauref,tau,aerosol,reffrad,nueffrad, QREFvis3d,
-     &            QREFir3d,omegaREFvis3d,omegaREFir3d,zdqnorm)
+     &      pq,ccn,tauref,tau,aerosol,reffrad,nueffrad, QREFvis3d,
+     &      QREFir3d,omegaREFvis3d,omegaREFir3d,zdqnorm)
+
 c     Starting loop on sub-domain
 c     ----------------------------
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/initracer_use.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/initracer_use.F	(revision 613)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/initracer_use.F	(revision 665)
@@ -43,5 +43,4 @@
 #include "advtrac.h"
 #include "comgeomfi.h"
-#include "watercap.h"
 #include "chimiedata.h"
 
@@ -111,6 +110,11 @@
               write(*,*) 'initracer: error expected dustbin=2'
             else
-              noms(1)='dust_mass'   ! dust mass mixing ratio
-              noms(2)='dust_number' ! dust number mixing ratio
+!              noms(1)='dust_mass'   ! dust mass mixing ratio
+!              noms(2)='dust_number' ! dust number mixing ratio
+! same as above, but avoid explict possible out-of-bounds on array
+               noms(1)='dust_mass'   ! dust mass mixing ratio
+               do iq=2,2
+               noms(iq)='dust_number' ! dust number mixing ratio
+               enddo
             endif
           endif
@@ -120,6 +124,11 @@
           noms(nqmx)='h2o_vap'
           mmol(nqmx)=18.
-          noms(nqmx-1)='h2o_ice'
-          mmol(nqmx-1)=18.
+!            noms(nqmx-1)='h2o_ice'
+!            mmol(nqmx-1)=18.
+          !"loop" to avoid potential out-of-bounds in array
+          do iq=nqmx-1,nqmx-1
+            noms(iq)='h2o_ice'
+            mmol(iq)=18.
+          enddo
         endif
         ! 3. Chemistry
@@ -157,5 +166,8 @@
       if (oldnames.and.water) then
         write(*,*)'initracer: moving surface water ice to index ',nqmx-1
-        qsurf(1:ngridmx,nqmx-1)=qsurf(1:ngridmx,nqmx)
+        ! "loop" to avoid potential out-of-bounds on arrays
+        do iq=nqmx-1,nqmx-1
+          qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1)
+        enddo
         qsurf(1:ngridmx,nqmx)=0
       endif 
@@ -218,4 +230,19 @@
       igcm_ar=0
       igcm_ar_n2=0
+      igcm_n=0
+      igcm_no=0
+      igcm_no2=0
+      igcm_n2d=0
+      igcm_co2plus=0
+      igcm_oplus=0
+      igcm_o2plus=0
+      igcm_coplus=0
+      igcm_cplus=0
+      igcm_nplus=0
+      igcm_noplus=0
+      igcm_n2plus=0
+      igcm_hplus=0
+      igcm_elec=0
+
 
       ! 1. find dust tracers
@@ -320,4 +347,74 @@
           count=count+1
         endif
+ 	if (noms(iq).eq."n") then
+          igcm_n=iq
+          mmol(igcm_n)=14.
+          count=count+1
+        endif
+        if (noms(iq).eq."no") then
+          igcm_no=iq
+          mmol(igcm_no)=30.
+          count=count+1
+        endif
+        if (noms(iq).eq."no2") then
+          igcm_no2=iq
+          mmol(igcm_no2)=46.
+          count=count+1
+        endif
+        if (noms(iq).eq."n2d") then
+          igcm_n2d=iq
+          mmol(igcm_n2d)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."co2plus") then
+          igcm_co2plus=iq
+          mmol(igcm_co2plus)=44.
+          count=count+1
+        endif
+        if (noms(iq).eq."oplus") then
+          igcm_oplus=iq
+          mmol(igcm_oplus)=16.
+          count=count+1
+        endif
+        if (noms(iq).eq."o2plus") then
+          igcm_o2plus=iq
+          mmol(igcm_o2plus)=32.
+          count=count+1
+        endif
+        if (noms(iq).eq."coplus") then
+          igcm_coplus=iq
+          mmol(igcm_coplus)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."cplus") then
+          igcm_cplus=iq
+          mmol(igcm_cplus)=12.
+          count=count+1
+        endif
+        if (noms(iq).eq."nplus") then
+          igcm_nplus=iq
+          mmol(igcm_nplus)=14.
+          count=count+1
+        endif
+        if (noms(iq).eq."noplus") then
+          igcm_noplus=iq
+          mmol(igcm_noplus)=30.
+          count=count+1
+        endif
+ 	if (noms(iq).eq."n2plus") then
+          igcm_n2plus=iq
+          mmol(igcm_n2plus)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."hplus") then
+          igcm_hplus=iq
+          mmol(igcm_hplus)=1.
+          count=count+1
+        endif
+        if (noms(iq).eq."elec") then
+          igcm_elec=iq
+          mmol(igcm_elec)=1./1822.89
+          count=count+1
+        endif
         if (noms(iq).eq."h2o_vap") then
           igcm_h2o_vap=iq
@@ -336,4 +433,6 @@
           count=count+1
         endif
+ 
+
       enddo ! of do iq=1,nqmx
 !      count=count+nbqchem
@@ -486,74 +585,4 @@
          alpha_lift(igcm_h2o_vap) =0.
          alpha_devil(igcm_h2o_vap)=0.
-
-c       "Dryness coefficient" controlling the evaporation and
-c        sublimation from the ground water ice (close to 1)
-c        HERE, the goal is to correct for the fact
-c        that the simulated permanent water ice polar caps
-c        is larger than the actual cap and the atmospheric
-c        opacity not always realistic.
-
-         do ig=1,ngridmx
-           if (ngridmx.ne.1) watercaptag(ig)=.false.
-           dryness(ig) = 1.
-         enddo
-
-         IF (caps) THEN
-c Perennial H20 north cap defined by watercaptag=true (allows surface to be
-c hollowed by sublimation in vdifc).
-         yeyey = 0
-         do ig=1,ngridmx
-!          !!! TESTS TESTS outliers
-!          !!! TESTS TESTS outliers
-!          if ( ( lati(ig)*180./pi      .ge.  75 ) .and.
-!     .         ( lati(ig)*180./pi      .le.  77 ) .and.
-!     .         ( ( ( long(ig)*180./pi .ge. 000. ) .and.
-!     .              ( long(ig)*180./pi .le. 120. ) ) 
-!     .             .or.
-!     .             ( ( long(ig)*180./pi .ge. -130. ) .and.
-!     .             ( long(ig)*180./pi .le. -115. ) ) ) ) then
-!             if (yeyey .eq. 0) then  !!! 1/2 en 64x48 sinon trop large en lat
-!              write(*,*) "outliers ", lati(ig)*180./pi, long(ig)*180./pi
-!              if (ngridmx.ne.1) watercaptag(ig)=.true.
-!              dryness(ig) = 1.
-!              albedodat(ig) = 0.45 !! comme alb_surfice
-!              yeyey = 1
-!             else
-!              yeyey = 0
-!             endif
-!          endif
-!          !!! TESTS TESTS outliers
-!          !!! TESTS TESTS outliers
-!
-!          !!! TESTS TESTS addcap
-!          !!! TESTS TESTS addcap     
-!          if ( ( lati(ig)*180./pi      .ge.  82 ) .and.
-!     .         ( lati(ig)*180./pi      .le.  84 ) .and.
-!     .         ( ( long(ig)*180./pi .gt. -030. ) .and.
-!     .              ( long(ig)*180./pi .lt. 090. ) ) ) then
-!              write(*,*) "capadd ", lati(ig)*180./pi, long(ig)*180./pi
-!              if (ngridmx.ne.1) watercaptag(ig)=.true.
-!              albedodat(ig) = 0.45 !! comme alb_surfice
-!              dryness(ig) = 1.
-!          endif
-!          !!! TESTS TESTS addcap
-!          !!! TESTS TESTS addcap 
-    
-           if (lati(ig)*180./pi.gt.84) then
-             if (ngridmx.ne.1) watercaptag(ig)=.true.
-             dryness(ig) = 1.
-c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
-c	      if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
-c	        if (ngridmx.ne.1) watercaptag(ig)=.true.
-c               dryness(ig) = 1.
-c	      endif
-c             if (lati(ig)*180./pi.ge.85) then
-c               if (ngridmx.ne.1) watercaptag(ig)=.true.
-c               dryness(ig) = 1.
-c	      endif
-           endif  ! (lati>80 deg)
-         end do ! (ngridmx)
-        ENDIF ! (caps)
-
          if(water.and.(nqmx.ge.2)) then
            radius(igcm_h2o_ice)=3.e-6
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/physiq_modif_use.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/physiq_modif_use.F	(revision 613)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/GCM_modif/physiq_modif_use.F	(revision 665)
@@ -1,10 +1,14 @@
-      SUBROUTINE physiq(ngrid,nlayer,nq,
-     $            firstcall,lastcall,
-     $            pday,ptime,ptimestep,
-     $            pplev,pplay,pphi,
-     $            pu,pv,pt,pq,
-     $            pw,
-     $            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
-
+      SUBROUTINE physiq(
+     $            ngrid,nlayer,nq
+     $            ,firstcall,lastcall
+     $            ,pday,ptime,ptimestep
+     $            ,pplev,pplay,pphi
+     $            ,pu,pv,pt,pq
+     $            ,pw
+     $            ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_invar.F"
+#endif
+     $            )
 
       IMPLICIT NONE
@@ -42,11 +46,10 @@
 c      8. Contribution to tendencies due to thermosphere
 c      9. Surface and sub-surface temperature calculations
-c     10. Renormalisation
-c     11. Write outputs :
+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     12. Diagnostic: mass conservation of tracers
+c     11. Diagnostic: mass conservation of tracers
 c 
 c   author: 
@@ -61,4 +64,5 @@
 c           Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009)
 c             Nb: See callradite.F for more information.
+c           Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
 c           
 c   arguments:
@@ -124,5 +128,4 @@
 
 #include "chimiedata.h"
-#include "watercap.h"
 #include "param.h"
 #include "param_v3.h"
@@ -131,5 +134,12 @@
 #include "netcdf.inc"
 
-
+#include "slope.h"
+
+#ifdef MESOSCALE
+#include "wrf_output_2d.h"
+#include "wrf_output_3d.h"
+#include "advtrac.h"   !!! this is necessary for tracers (in dyn3d)
+#include "meso_inc/meso_inc_var.F"
+#endif
 
 c Arguments :
@@ -180,6 +190,7 @@
       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) 
+      REAL q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
+      
+      REAL watercapflag(ngridmx)     ! water cap flag
 
 c     Variables used by the water ice microphysical scheme:
@@ -187,7 +198,11 @@
       REAL nuice(ngridmx,nlayermx)   ! Estimated effective variance
                                      !   of the size distribution
-c     Albedo of deposited surface ice
-      !!REAL, PARAMETER :: alb_surfice = 0.4 ! 0.45
-      REAL, PARAMETER :: alb_surfice = 0.45 !!TESTS_JB 
+
+c     Variables used by the slope model
+      REAL sl_ls, sl_lct, sl_lat
+      REAL sl_tau, sl_alb, sl_the, sl_psi
+      REAL sl_fl0, sl_flu
+      REAL sl_ra, sl_di0
+      REAL sky
 
       SAVE day_ini, icount
@@ -195,5 +210,4 @@
       SAVE co2ice,albedo,emis, q2
       SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
-      SAVE ig_vl1
 
       REAL stephan   
@@ -251,5 +265,5 @@
       REAL zdqchim(ngridmx,nlayermx,nqmx)
       REAL zdqschim(ngridmx,nqmx)
-      REAL zdqnorm(ngridmx,nlayermx,2)                     !quantity of dust which have to be added by the dynamical core
+      REAL zdqnorm(ngridmx,nlayermx,2) !quantity of dust which have to be added
 
 
@@ -304,4 +318,18 @@
       REAL time_phys
 
+c Variables for PBL
+
+      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
+      REAL, SAVE :: wmax_th(ngridmx)
+      REAL hfmax_th(ngridmx)
+      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
+      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
+      INTEGER lmax_th(ngridmx)
+      REAL dtke_th(ngridmx,nlayermx+1)
+      REAL zcdv(ngridmx), zcdh(ngridmx)
+      REAL Teta_out(ngridmx),u_out(ngridmx)  ! Interpolated teta and u at z_out
+      REAL z_out                          ! height of interpolation between z0 and z1
+      REAL ustar(ngridmx),tstar(ngridmx)  ! friction velocity and friction potential temp
+      REAL zu2(ngridmx)
 c=======================================================================
 
@@ -315,11 +343,13 @@
 c        variables set to 0
 c        ~~~~~~~~~~~~~~~~~~
-         call zerophys(ngrid*nlayer*naerkind,aerosol)
-         call zerophys(ngrid*nlayer,dtrad)
-         call zerophys(ngrid,fluxrad)
+         aerosol(:,:,:)=0
+         dtrad(:,:)=0
+         fluxrad(:)=0
+
+         wmax_th(:)=0.
 
 c        read startfi 
 c        ~~~~~~~~~~~~
-
+#ifndef MESOSCALE
 ! Read netcdf initial physical parameters.
          CALL phyetat0 ("startfi.nc",0,0,
@@ -327,4 +357,7 @@
      &         day_ini,time_phys,
      &         tsurf,tsoil,emis,q2,qsurf,co2ice)
+#else
+#include "meso_inc/meso_inc_ini.F"
+#endif
 
          if (pday.ne.day_ini) then
@@ -337,4 +370,11 @@
 
          write (*,*) 'In physiq day_ini =', day_ini
+
+c        initialize tracers
+c        ~~~~~~~~~~~~~~~~~~
+         tracerdyn=tracer
+         IF (tracer) THEN
+            CALL initracer(qsurf,co2ice)
+         ENDIF  ! end tracer
 
 c        Initialize albedo and orbital calculation
@@ -358,29 +398,10 @@
          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 
-
+#ifndef MESOSCALE
 c        Initialize thermospheric parameters
 c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
          if (callthermos) call param_read
-
+#endif
 c        Initialize R and Cp as constant
 
@@ -396,5 +417,6 @@
 
         IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
-          write(*,*)"physiq: water_param Surface ice alb:",alb_surfice
+          write(*,*)"physiq: water_param Surface water ice albedo:", 
+     .                  albedo_h2o_ice
         ENDIF
                    
@@ -415,12 +437,12 @@
 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)
+      pdv(:,:)=0
+      pdu(:,:)=0
+      pdt(:,:)=0
+      pdq(:,:,:)=0
+      pdpsrf(:)=0
+      zflubid(:)=0
+      zdtsurf(:)=0
+      dqsurf(:,:)=0
       igout=ngrid/2+1 
 
@@ -469,4 +491,5 @@
       ENDDO
 
+#ifndef MESOSCALE
 c-----------------------------------------------------------------------
 c    1.2.5 Compute mean mass, cp, and R
@@ -476,4 +499,5 @@
          call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
       endif
+#endif
 c-----------------------------------------------------------------------
 c    2. Compute radiative tendencies :
@@ -518,4 +542,5 @@
 c          Radiative transfer
 c          ------------------
+
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
      $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
@@ -523,7 +548,55 @@
      &     tauref,tau,aerosol,ccn,rdust,rice,nuice,zdqnorm)
 
+c          Outputs for basic check (middle of domain)
+c          ------------------------------------------
+           print*, 'Ls =',zls*180./pi,
+     &             'check lat lon', lati(igout)*180/pi,
+     &                              long(igout)*180/pi
+           print*, 'tauref(700 Pa) =',tauref(igout),
+     &             ' tau(700 Pa) =',tau(igout,1)*700./pplev(igout,1)
+
+c          ---------------------------------------------------------
+c          Call slope parameterization for direct and scattered flux
+c          ---------------------------------------------------------
+           IF(callslope) THEN
+            print *, 'Slope scheme is on and computing...'
+            DO ig=1,ngrid  
+              sl_the = theta_sl(ig)
+              IF (sl_the .ne. 0.) THEN
+                ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2)
+                DO l=1,2
+                 sl_lct = ptime*24. + 180.*long(ig)/pi/15.
+                 sl_ra  = pi*(1.0-sl_lct/12.)
+                 sl_lat = 180.*lati(ig)/pi
+                 sl_tau = tau(ig,1)
+                 sl_alb = albedo(ig,l)
+                 sl_psi = psi_sl(ig)
+                 sl_fl0 = fluxsurf_sw(ig,l)
+                 sl_di0 = 0.
+                 if (mu0(ig) .gt. 0.) then
+                  sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
+                  sl_di0 = sl_di0*1370./dist_sol/dist_sol
+                  sl_di0 = sl_di0/ztim1
+                  sl_di0 = fluxsurf_sw(ig,l)*sl_di0
+                 endif
+                 ! you never know (roundup concern...)
+                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
+                 !!!!!!!!!!!!!!!!!!!!!!!!!!
+                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 
+     &                             sl_tau, sl_alb, sl_the, sl_psi,
+     &                             sl_di0, sl_fl0, sl_flu )
+                 !!!!!!!!!!!!!!!!!!!!!!!!!!
+                 fluxsurf_sw(ig,l) = sl_flu
+                ENDDO
+              !!! compute correction on IR flux as well
+              sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
+              fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
+              ENDIF
+            ENDDO
+           ENDIF
+
 c          CO2 near infrared absorption
 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-           call zerophys(ngrid*nlayer,zdtnirco2)
+           zdtnirco2(:,:)=0
            if (callnirco2) then
               call nirco2abs (ngrid,nlayer,pplay,dist_sol,
@@ -554,6 +627,4 @@
            ENDIF
 
-
-
         ENDIF ! of if(mod(icount-1,iradia).eq.0)
 
@@ -569,6 +640,9 @@
      $         stephan*zplanck(ig)*zplanck(ig)
                fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
+               IF(callslope) THEN
+                 sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
+                 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
+               ENDIF
            ENDDO
-
 
          DO l=1,nlayer
@@ -602,7 +676,6 @@
 c    4. Vertical diffusion (turbulent mixing):
 c    -----------------------------------------
-c
+
       IF (calldifv) THEN
-
 
          DO ig=1,ngrid
@@ -610,6 +683,6 @@
          ENDDO
 
-         CALL zerophys(ngrid*nlayer,zdum1)
-         CALL zerophys(ngrid*nlayer,zdum2)
+         zdum1(:,:)=0
+         zdum2(:,:)=0
          do l=1,nlayer
             do ig=1,ngrid
@@ -617,5 +690,28 @@
             enddo
          enddo
-         
+
+
+#ifdef MESOSCALE
+      IF (.not.flag_LES) THEN
+#endif
+c ----------------------
+c Treatment of a special case : using new surface layer (Richardson based)
+c without using the thermals in gcm and mesoscale can yield problems in
+c weakly unstable situations when winds are near to 0. For those cases, we add
+c a unit subgrid gustiness. Remember that thermals should be used we using the 
+c Richardson based surface layer model.
+        IF ( .not.calltherm .and. callrichsl ) THEN
+          DO ig=1, ngridmx
+             IF (zh(ig,1) .lt. tsurf(ig)) THEN
+               wmax_th(ig)=1.
+             ENDIF        
+          ENDDO
+        ENDIF
+c ----------------------
+#ifdef MESOSCALE
+      ENDIF
+#endif
+
+
 c        Calling vdif (Martian version WITH CO2 condensation)
          CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
@@ -625,6 +721,9 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif)
-
+     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
+
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_les.F"
+#endif
          DO l=1,nlayer
             DO ig=1,ngrid
@@ -638,7 +737,7 @@
          ENDDO
 
-         DO ig=1,ngrid
-            zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
-         ENDDO
+          DO ig=1,ngrid
+             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
+          ENDDO
 
          if (tracer) then 
@@ -662,6 +761,57 @@
      s      (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
          ENDDO
+#ifdef MESOSCALE
+         IF (flag_LES) THEN
+            write(*,*) 'LES mode !' 
+            write(*,*) 'Please set calldifv to T in callphys.def'
+            STOP
+         ENDIF
+#endif
       ENDIF ! of IF (calldifv)
 
+c-----------------------------------------------------------------------
+c   TEST. Thermals :
+c HIGHLY EXPERIMENTAL, BEWARE !!
+c   -----------------------------
+ 
+      if(calltherm) then
+ 
+        call calltherm_interface(firstcall,
+     $ long,lati,zzlev,zzlay,
+     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
+     $ pplay,pplev,pphi,zpopsk,
+     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
+     $ dtke_th,hfmax_th,wmax_th)
+ 
+         DO l=1,nlayer
+           DO ig=1,ngrid
+              pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
+              pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
+              pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
+              q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
+           ENDDO
+        ENDDO
+ 
+        DO ig=1,ngrid
+          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
+        ENDDO      
+    
+        if (tracer) then
+        DO iq=1,nq
+         DO l=1,nlayer
+           DO ig=1,ngrid
+             pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
+           ENDDO
+         ENDDO
+        ENDDO
+        endif
+
+        lmax_th_out(:)=real(lmax_th(:))
+
+        else   !of if calltherm
+        lmax_th(:)=0
+        wmax_th(:)=0.
+        lmax_th_out(:)=0.
+        end if
 
 c-----------------------------------------------------------------------
@@ -676,14 +826,15 @@
             ENDDO
          ENDDO
-         CALL zerophys(ngrid*nlayer,zduadj)
-         CALL zerophys(ngrid*nlayer,zdvadj)
-         CALL zerophys(ngrid*nlayer,zdhadj)
+         zduadj(:,:)=0
+         zdvadj(:,:)=0
+         zdhadj(:,:)=0
 
          CALL convadj(ngrid,nlayer,nq,ptimestep,
-     $                pplay,pplev,zpopsk,
+     $                pplay,pplev,zpopsk,lmax_th,
      $                pu,pv,zh,pq,
      $                pdu,pdv,zdh,pdq,
      $                zduadj,zdvadj,zdhadj,
      $                zdqadj)
+
 
          DO l=1,nlayer
@@ -718,5 +869,5 @@
      $              co2ice,albedo,emis,
      $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
-     $	            fluxsurf_sw,zls) 
+     $              fluxsurf_sw,zls) 
 
          DO l=1,nlayer
@@ -793,4 +944,5 @@
 c     ------------------
 
+#ifndef MESOSCALE
 c        --------------
 c        photochemistry :
@@ -837,4 +989,5 @@
 
          END IF  ! of IF (photochem.or.thermochem)
+#endif
 
 c   7c. Aerosol particles
@@ -908,5 +1061,5 @@
       endif !  of if (tracer) 
 
-
+#ifndef MESOSCALE
 c-----------------------------------------------------------------------
 c   8. THERMOSPHERE CALCULATION
@@ -933,4 +1086,5 @@
 
       endif ! of if (callthermos)
+#endif
 
 c-----------------------------------------------------------------------
@@ -951,4 +1105,5 @@
 
       IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
+#ifndef MESOSCALE
          if (caps.and.(obliquit.lt.27.)) then
            ! NB: Updated surface pressure, at grid point 'ngrid', is
@@ -957,6 +1112,7 @@
      &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
          endif
+#endif
 c       -------------------------------------------------------------
-c       Change of surface albedo (set to 0.4) in case of ground frost
+c       Change of surface albedo in case of ground frost
 c       everywhere except on the north permanent cap and in regions
 c       covered by dry ice. 
@@ -965,7 +1121,10 @@
          do ig=1,ngrid
            if ((co2ice(ig).eq.0).and.
-     &        (qsurf(ig,igcm_h2o_ice).gt.0.005)) then
-              albedo(ig,1) = alb_surfice
-              albedo(ig,2) = alb_surfice
+     &        (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
+              albedo(ig,1) = albedo_h2o_ice
+              albedo(ig,2) = albedo_h2o_ice
+c              write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice)
+c              write(*,*) "physiq.F frost :"
+c     &        ,lati(ig)*180./pi, long(ig)*180./pi
            endif
          enddo  ! of do ig=1,ngrid
@@ -994,5 +1153,5 @@
 
 c-----------------------------------------------------------------------
-c  11. Write output files
+c  10. Write output files
 c  ----------------------
 
@@ -1039,7 +1198,16 @@
       ENDDO
 
-c    Compute surface stress : (NB: z0 is a common in planete.h)
+      ! Potential Temperature
+
+       DO ig=1,ngridmx
+          DO l=1,nlayermx
+              zh(ig,l) = zt(ig,l)*(zplay(ig,l)/zplev(ig,1))**rcp
+          ENDDO
+       ENDDO
+
+
+c    Compute surface stress : (NB: z0 is a common in surfdat.h)
 c     DO ig=1,ngrid
-c        cd = (0.4/log(zzlay(ig,1)/z0))**2
+c        cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
 c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
 c     ENDDO
@@ -1047,6 +1215,6 @@
 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)
+             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 ******************************************************
@@ -1061,5 +1229,5 @@
         ENDDO
       ENDDO
-      if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then	
+      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),
@@ -1088,9 +1256,6 @@
 
       IF (ngrid.NE.1) THEN
-         print*,'Ls =',zls*180./pi,
-     &   ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),
-     &   ' tau(Viking1) =',tau(ig_vl1,1)
-
-
+
+#ifndef MESOSCALE
 c        -------------------------------------------------------------------
 c        Writing NetCDF file  "RESTARTFI" at the end of the run
@@ -1111,4 +1276,5 @@
      .              zgam,zthe)
          ENDIF
+#endif
 
 c        -------------------------------------------------------------------
@@ -1120,8 +1286,8 @@
            if (water) then
 
-             call zerophys(ngrid,mtot)
-             call zerophys(ngrid,icetot)
-             call zerophys(ngrid,rave)
-             call zerophys(ngrid,tauTES)
+             mtot(:)=0
+             icetot(:)=0
+             rave(:)=0
+             tauTES(:)=0
              do ig=1,ngrid 
                do l=1,nlayermx
@@ -1160,5 +1326,4 @@
 c        which can later be used to make the statistic files of the run:
 c        "stats")          only possible in 3D runs !
-
          
          IF (callstats) THEN
@@ -1268,4 +1433,35 @@
          ENDIF
 
+
+#ifdef MESOSCALE
+      !!!
+      !!! OUTPUT FIELDS
+      !!!
+      wtsurf(1:ngrid) = tsurf(1:ngrid)    !! surface temperature
+      wco2ice(1:ngrid) = co2ice(1:ngrid)  !! co2 ice
+      mtot(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice
+      !! JF 
+      TAU_lay(:)=tau(:,1)!!true opacity (not a reference like tauref)
+      IF (igcm_dust_mass .ne. 0) THEN
+        qsurfice_dust(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass)
+      ENDIF
+      IF (igcm_h2o_ice .ne. 0) THEN      
+        qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
+        vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
+     .           * mugaz / mmol(igcm_h2o_ice)
+      ENDIF
+      !! Dust quantity integration along the vertical axe
+      dustot(:)=0
+      do ig=1,ngrid
+       do l=1,nlayermx
+        dustot(ig) = dustot(ig) +
+     &               zq(ig,l,igcm_dust_mass)
+     &               *  (pplev(ig,l) - pplev(ig,l+1)) / g
+       enddo
+      enddo
+c AUTOMATICALLY GENERATED FROM REGISTRY
+#include "fill_save.inc"
+#else
+
 c        ==========================================================
 c        WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing
@@ -1275,6 +1471,8 @@
 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)
+c        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+c    &                  emis)
+!         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
+!         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
      &                  tsurf)
@@ -1282,4 +1480,7 @@
         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
      &                  co2ice)
+
+c         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
+c     &                  "K",2,zt(1,7))
 c         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
 c     &                  fluxsurf_lw)
@@ -1291,11 +1492,10 @@
 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)
+         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
 c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
-c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
+!        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
 c        call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
 c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
@@ -1305,11 +1505,18 @@
 c        call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
 c    &                   'w.m-2',3,zdtlw)
-
-!!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL
+c        CALL WRITEDIAGFI(ngridmx,'tauTESap',
+c     &                         'tau abs 825 cm-1',
+c     &                         '',2,tauTES)
+
+#ifdef MESOINI
+        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
+     &                  emis)
+        call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
         call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
      &                       "K",3,tsoil)
         call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
      &                       "K",3,inertiedat)
-!!!!!!!!!!!!!!!!!!!!!!!!SOIL SOIL SOIL
+#endif
+
 
 c        ----------------------------------------------------------
@@ -1324,5 +1531,5 @@
         
          ! Compute co2 column
-         call zerophys(ngrid,co2col)
+         co2col(:)=0
          do l=1,nlayermx
            do ig=1,ngrid
@@ -1331,6 +1538,6 @@
            enddo
          enddo
-c         call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
-c     &                  co2col)
+         call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2,
+     &                  co2col)
          endif ! of if (tracer.and.(igcm_co2.ne.0))
 
@@ -1341,4 +1548,5 @@
            if (water) then
 
+#ifdef MESOINI
             !!!! waterice = q01, voir readmeteo.F90
             call WRITEDIAGFI(ngridmx,'q01',noms(igcm_h2o_ice),
@@ -1353,4 +1561,5 @@
      &                      'kg.m-2',2,
      &                       qsurf(1:ngridmx,igcm_h2o_ice))
+#endif
 
              CALL WRITEDIAGFI(ngridmx,'mtot',
@@ -1360,22 +1569,34 @@
      &                       'total mass of water ice',
      &                       'kg/m2',2,icetot)
-cc            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
-cc    &                *mugaz/mmol(igcm_h2o_ice)
-cc            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
-cc    &                       'mol/mol',3,vmr)
-c             CALL WRITEDIAGFI(ngridmx,'reffice',
-c     &                       'Mean reff',
-c     &                       'm',2,rave)
-cc            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
-cc    &                       'm',3,rice)
-cc            If activice is true, tauTES is computed in aeropacity.F;
-c             if (.not.activice) then
-c               CALL WRITEDIAGFI(ngridmx,'tauTESap',
-c     &                         'tau abs 825 cm-1',
-c     &                         '',2,tauTES)
-c             endif
-c             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
-c     &                       'surface h2o_ice',
-c     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+c            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+c    &                *mugaz/mmol(igcm_h2o_ice)
+c            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
+c    &                       'mol/mol',3,vmr)
+             CALL WRITEDIAGFI(ngridmx,'reffice',
+     &                       'Mean reff',
+     &                       'm',2,rave)
+c            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
+c    &                       'm',3,rice)
+c            If activice is true, tauTES is computed in aeropacity.F;
+             if (.not.activice) then
+               CALL WRITEDIAGFI(ngridmx,'tauTESap',
+     &                         'tau abs 825 cm-1',
+     &                         '',2,tauTES)
+             endif
+             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
+     &                       'surface h2o_ice',
+     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+
+            if (caps) then
+             do ig=1,ngridmx
+                if (watercaptag(ig)) watercapflag(ig) = 1
+             enddo
+             CALL WRITEDIAGFI(ngridmx,'watercaptag',
+     &                         'Ice water caps',
+     &                         '',2,watercapflag)
+            endif
+            CALL WRITEDIAGFI(ngridmx,'albedo',
+     &                         'albedo',
+     &                         '',2,albedo(1:ngridmx,1))
            endif !(water)
 
@@ -1407,18 +1628,22 @@
 c          call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
            if (doubleq) then
-cc            call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
-cc    &                       'kg.m-2',2,qsurf(1,1))
-cc            call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
-cc    &                       'N.m-2',2,qsurf(1,2))
-cc            call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
-cc    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
-cc            call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
-cc    &                       'kg.m-2.s-1',2,zdqssed(1,1))
-c             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
-c     &                        'm',3,rdust*ref_r0)
+c            call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
+c    &                       'kg.m-2',2,qsurf(1,1))
+c            call WRITEDIAGFI(ngridmx,'Nsurf','N particles',
+c    &                       'N.m-2',2,qsurf(1,2))
+c            call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift',
+c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
+c            call WRITEDIAGFI(ngridmx,'dqssed','sedimentation',
+c    &                       'kg.m-2.s-1',2,zdqssed(1,1))
+             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
+     &                        'm',3,rdust*ref_r0)
              call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
      &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
-            call WRITEDIAGFI(ngridmx,'dustN','Dust number',
+c            call WRITEDIAGFI(ngridmx,'dustN','Dust number',
+c    &                        'part/kg',3,pq(1,1,igcm_dust_number))
+#ifdef MESOINI
+             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
      &                        'part/kg',3,pq(1,1,igcm_dust_number))
+#endif
            else
              do iq=1,dustbin
@@ -1437,4 +1662,87 @@
 
 c        ----------------------------------------------------------
+c        ----------------------------------------------------------
+c        PBL OUTPUS
+c        ----------------------------------------------------------
+c        ----------------------------------------------------------
+
+
+c        ----------------------------------------------------------
+c        Outputs of surface layer
+c        ----------------------------------------------------------
+
+
+         z_out=0.
+         if (calltherm .and. (z_out .gt. 0.)) then
+         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
+     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar)
+
+         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
+         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
+     &              'horizontal velocity norm','m/s',
+     &                         2,zu2)
+
+         call WRITEDIAGFI(ngridmx,'Teta_out',
+     &              'potential temperature at z_out','K',
+     &                         2,Teta_out)
+         call WRITEDIAGFI(ngridmx,'u_out',
+     &              'horizontal velocity norm at z_out','m/s',
+     &                         2,u_out)
+         call WRITEDIAGFI(ngridmx,'u*',
+     &              'friction velocity','m/s',
+     &                         2,ustar)
+         call WRITEDIAGFI(ngridmx,'teta*',
+     &              'friction potential temperature','K',
+     &                         2,tstar)
+         else
+           if((.not. calltherm).and.(z_out .gt. 0.)) then
+            print*, 'WARNING : no interpolation in surface-layer :'
+            print*, 'Outputing surface-layer quantities without thermals
+     & does not make sense'
+           endif
+         endif
+
+c        ----------------------------------------------------------
+c        Outputs of thermals
+c        ----------------------------------------------------------
+         if (calltherm) then
+
+!        call WRITEDIAGFI(ngrid,'dtke',
+!     &              'tendance tke thermiques','m**2/s**2',
+!     &                         3,dtke_th)
+!        call WRITEDIAGFI(ngrid,'d_u_ajs',
+!     &              'tendance u thermiques','m/s',
+!     &                         3,pdu_th*ptimestep)
+!        call WRITEDIAGFI(ngrid,'d_v_ajs',
+!     &              'tendance v thermiques','m/s',
+!     &                         3,pdv_th*ptimestep)
+!        if (tracer) then
+!        if (nq .eq. 2) then
+!        call WRITEDIAGFI(ngrid,'deltaq_th',
+!     &              'delta q thermiques','kg/kg',
+!     &                         3,ptimestep*pdq_th(:,:,2))
+!        endif
+!        endif
+
+        call WRITEDIAGFI(ngridmx,'lmax_th',
+     &              'hauteur du thermique','K',
+     &                         2,lmax_th_out)
+        call WRITEDIAGFI(ngridmx,'hfmax_th',
+     &              'maximum TH heat flux','K.m/s',
+     &                         2,hfmax_th)
+        call WRITEDIAGFI(ngridmx,'wmax_th',
+     &              'maximum TH vertical velocity','m/s',
+     &                         2,wmax_th)
+
+         endif
+
+c        ----------------------------------------------------------
+c        ----------------------------------------------------------
+c        END OF PBL OUTPUS
+c        ----------------------------------------------------------
+c        ----------------------------------------------------------
+
+
+c        ----------------------------------------------------------
 c        Output in netcdf file "diagsoil.nc" for subterranean
 c          variables (output every "ecritphy", as for writediagfi)
@@ -1451,4 +1759,5 @@
 c        END OF WRITEDIAGFI
 c        ==========================================================
+#endif
 
       ELSE     ! if(ngrid.eq.1)
@@ -1470,4 +1779,71 @@
 c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
 
+! THERMALS STUFF 1D
+
+         z_out=0.
+         if (calltherm .and. (z_out .gt. 0.)) then
+         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
+     &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar)
+
+         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
+         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
+     &              'horizontal velocity norm','m/s',
+     &                         0,zu2)
+
+         call WRITEDIAGFI(ngridmx,'Teta_out',
+     &              'potential temperature at z_out','K',
+     &                         0,Teta_out)
+         call WRITEDIAGFI(ngridmx,'u_out',
+     &              'horizontal velocity norm at z_out','m/s',
+     &                         0,u_out)
+         call WRITEDIAGFI(ngridmx,'u*',
+     &              'friction velocity','m/s',
+     &                         0,ustar)
+         call WRITEDIAGFI(ngridmx,'teta*',
+     &              'friction potential temperature','K',
+     &                         0,tstar)
+         else
+           if((.not. calltherm).and.(z_out .gt. 0.)) then
+            print*, 'WARNING : no interpolation in surface-layer :'
+            print*, 'Outputing surface-layer quantities without thermals
+     & does not make sense'
+           endif
+         endif
+
+         if(calltherm) then
+
+        call WRITEDIAGFI(ngridmx,'lmax_th',
+     &              'hauteur du thermique','point',
+     &                         0,lmax_th_out)
+        call WRITEDIAGFI(ngridmx,'hfmax_th',
+     &              'maximum TH heat flux','K.m/s',
+     &                         0,hfmax_th)
+        call WRITEDIAGFI(ngridmx,'wmax_th',
+     &              'maximum TH vertical velocity','m/s',
+     &                         0,wmax_th)
+
+         co2col(:)=0.
+         if (tracer) then
+         do l=1,nlayermx
+           do ig=1,ngrid
+             co2col(ig)=co2col(ig)+
+     &                  zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
+         enddo
+         enddo
+
+         end if
+         call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
+     &                                      ,'kg/m-2',0,co2col)
+         endif
+         call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &
+     &                              ,'m/s',1,pw)
+         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
+         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
+     &                  tsurf)
+         call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
+         call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
+
+         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
+         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
 ! or output in diagfi.nc (for testphys1d)
          call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F	(revision 613)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F	(revision 665)
@@ -472,5 +472,9 @@
 c--------------------------------------------------------------------
 
-      justbackground=.false.
+      justbackground=.false. !! default value
+      write(*,*) "justbackground"
+      call getin("justbackground",justbackground)
+      write(*,*) "justbackground = ",justbackground
+
 
       IF (justbackground .eq. .true.)  THEN
@@ -775,7 +779,11 @@
         DO l=1,nlayer
           DO ig=1,ngrid
-            dsodust(ig,l) = dsodust(ig,l) +
-     &                      aerosol(ig,l,iaerdust(iaer)) * g /
-     &                      (pplev(ig,l) - pplev(ig,l+1))
+c            dsodust(ig,l) = dsodust(ig,l) +
+c     &                      aerosol(ig,l,iaerdust(iaer)) * g /
+c     &                      (pplev(ig,l) - pplev(ig,l+1))
+            dsodust(ig,l) = 
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &          pq(ig,l,igcm_dust_mass)
           ENDDO
         ENDDO
