Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile	(revision 226)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/compile	(revision 234)
@@ -34,6 +34,6 @@
 nz=25
 ##################
-
-
+scatterers=2
+scatterers=1
 
 
@@ -62,5 +62,5 @@
 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
+makegcm -t ${tracers} -p mars -d ${nx}x${ny}x${nz} -s ${scatterers} gcm | tee libo/log_compile_gcm | grep 'warnings' >> libo/log_compile_gcm
 
 echo well, OK.
Index: trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F
===================================================================
--- trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F	(revision 226)
+++ trunk/MESOSCALE/LMDZ.MARS.new/in_lmdz_mars_newphys/physiq.F	(revision 234)
@@ -7,4 +7,7 @@
      $            ,pw
      $            ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_invar.F"
+#endif
      $            )
 
@@ -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:
@@ -131,5 +135,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 :
@@ -181,5 +192,4 @@
       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) 
 
 c     Variables used by the water ice microphysical scheme:
@@ -191,9 +201,15 @@
       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
       SAVE aerosol, tsurf,tsoil
       SAVE co2ice,albedo,emis, q2
       SAVE capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
-      SAVE ig_vl1
 
       REAL stephan   
@@ -306,5 +322,5 @@
       REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
       REAL wmax_th(ngridmx)
-      REAL ,SAVE :: hfmax_th(ngridmx)
+      REAL, SAVE :: hfmax_th(ngridmx)
       REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
       REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
@@ -330,5 +346,5 @@
 c        read startfi 
 c        ~~~~~~~~~~~~
-
+#ifndef MESOSCALE
 ! Read netcdf initial physical parameters.
          CALL phyetat0 ("startfi.nc",0,0,
@@ -336,4 +352,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
@@ -375,21 +394,14 @@
          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 
-
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_caps.F"
+#endif
+
+#ifndef MESOSCALE
 c        Initialize thermospheric parameters
 c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
          if (callthermos) call param_read
-
+#endif
 c        Initialize R and Cp as constant
 
@@ -478,4 +490,5 @@
       ENDDO
 
+#ifndef MESOSCALE
 c-----------------------------------------------------------------------
 c    1.2.5 Compute mean mass, cp, and R
@@ -485,4 +498,5 @@
          call concentrations(pplay,pt,pdt,pq,pdq,ptimestep)
       endif
+#endif
 c-----------------------------------------------------------------------
 c    2. Compute radiative tendencies :
@@ -532,4 +546,52 @@
      $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
      &     tauref,tau,aerosol,ccn,rdust,rice,nuice)
+
+c          Outputs for basic check (middle of domain)
+c          ------------------------------------------
+           print*, 'check lat lon', lati(igout)*180/pi,
+     .                              long(igout)*180/pi
+           print*, 'Ls =',zls*180./pi
+           print*, 'tauref(700 Pa) =',tauref(igout)
+           print*, '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
@@ -577,4 +639,8 @@
      $         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
 
@@ -633,4 +699,7 @@
      &        zdqdif,zdqsdif)
 
+#ifdef MESOSCALE
+#include "meso_inc/meso_inc_les.F"
+#endif
          DO l=1,nlayer
             DO ig=1,ngrid
@@ -668,4 +737,11 @@
      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)
 
@@ -840,4 +916,5 @@
 c     ------------------
 
+#ifndef MESOSCALE
 c        --------------
 c        photochemistry :
@@ -884,4 +961,5 @@
 
          END IF  ! of IF (photochem.or.thermochem)
+#endif
 
 c   7c. Aerosol particles
@@ -955,5 +1033,5 @@
       endif !  of if (tracer) 
 
-
+#ifndef MESOSCALE
 c-----------------------------------------------------------------------
 c   8. THERMOSPHERE CALCULATION
@@ -980,4 +1058,5 @@
 
       endif ! of if (callthermos)
+#endif
 
 c-----------------------------------------------------------------------
@@ -998,4 +1077,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
@@ -1004,4 +1084,5 @@
      &                     (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
          endif
+#endif
 c       -------------------------------------------------------------
 c       Change of surface albedo (set to 0.4) in case of ground frost
@@ -1095,5 +1176,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),
@@ -1122,9 +1203,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
@@ -1145,4 +1223,5 @@
      .              zgam,zthe)
          ENDIF
+#endif
 
 c        -------------------------------------------------------------------
@@ -1300,4 +1379,21 @@
             CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps)
          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
+      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
+c AUTOMATICALLY GENERATED FROM REGISTRY
+#include "fill_save.inc"
+#else
 
 c        ==========================================================
@@ -1526,4 +1622,5 @@
 c        END OF WRITEDIAGFI
 c        ==========================================================
+#endif
 
       ELSE     ! if(ngrid.eq.1)
