Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 226)
+++ /trunk/LMDZ.MARS/README	(revision 227)
@@ -814,2 +814,5 @@
 - Slight cosmetic changes to physiq.f and inifis.F
   --- some of them e.g. to prepare convergence between meso_physiq and physiq
+- Dropped the meso_ prefix from the slope routines
+  (currently those are only used in the mesoscale model,
+   but one could imagine using those in hi-res GCM runs)
Index: /trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_inifisini.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_inifisini.F	(revision 226)
+++ /trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_inifisini.F	(revision 227)
@@ -68,5 +68,5 @@
 
         !
-        ! in meso_slope.h
+        ! in slope.h
         !
         theta_sl(:)=wtheta(:)
Index: /trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_slope.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_slope.F	(revision 226)
+++ /trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_slope.F	(revision 227)
@@ -26,5 +26,5 @@
            if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
      !!!!!!!!!!!!!!!!!!!!!!!!!!
-        CALL meso_param_slope( mu0(ig), declin, sl_ra, sl_lat, 
+        CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 
      &            sl_tau, sl_alb, 
      &            sl_the, sl_psi, sl_di0, sl_fl0, sl_flu)
Index: /trunk/LMDZ.MARS/libf/phymars/meso_inifis.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_inifis.F	(revision 226)
+++ /trunk/LMDZ.MARS/libf/phymars/meso_inifis.F	(revision 227)
@@ -65,5 +65,5 @@
 #include "datafile.h"
 #ifdef MESOSCALE
-#include "meso_slope.h"
+#include "slope.h"
 #include "comsoil.h"     !!MESOSCALE -- needed to fill volcapa
 #include "meso_inc/meso_inc_inifisvar.F"
Index: unk/LMDZ.MARS/libf/phymars/meso_param_slope.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_param_slope.F90	(revision 226)
+++ 	(revision )
@@ -1,220 +1,0 @@
-SUBROUTINE meso_param_slope( &
-  !
-  ! INPUTS
-  !
-               &  csza, declin, rho, latitude &
-               &  ,taudust, albedo  &
-               &  ,theta_s, psi_s &
-               &  ,fdir_0, ftot_0  &
-  !
-  ! OUTPUTS
-  !
-               &  ,ftot &
- )
-
-
-!!*****************************************************************************************
-!
-! SUBROUTINE:
-! param_slope
-! 
-!
-! PURPOSE: 
-! computes total solar irradiance on a given Martian slope 
-!
-!
-! INPUTS:
-! csza		cosine solar zenith angle
-! declin	sun declination (rad)
-! rho		sun right ascension (rad)
-! latitude	latitude (deg)
-! taudust	dust optical depth at reference wavelength 0.67 mic. 
-! albedo	spectrally integrated surface Lambertian reflection albedo 
-! theta_s	slope inclination angle (deg)
-!			0 is horizontal, 90 is vertical
-! phi_s		slope azimuth (deg) 
-!       		0 >> Northward
-!       		90 >> Eastward
-!       		180 >> Southward
-!       		270 >> Westward
-! ftot_0	spectrally integrated total irradiance on an horizontal surface (W/m2)
-! fdir_0	spectrally integrated direct irradiance on an horizontal surface (W/m2)
-!
-!
-! OUTPUTS:
-! ftot		spectrally integrated total irradiance on the slope (W/m2)
-!
-! REFERENCE: 
-! "Fast and accurate estimation of irradiance on Martian slopes"
-! A. Spiga & F. Forget
-! .....
-!
-! AUTHOR: 
-! A. Spiga (spiga@lmd.jussieu.fr)
-! March 2008
-!
-!!*****************************************************************************************
-
-IMPLICIT NONE
-
-!!
-!! INPUT
-!!
-REAL, INTENT(IN) :: csza, declin, rho, latitude, taudust, theta_s, psi_s, albedo, ftot_0 , fdir_0 
-
-!!
-!! LOCAL
-!!
-REAL :: pi, deg2rad
-REAL :: a
-REAL :: mu_s, sigma_s
-REAL :: fdir, fscat, fscat_0, fref
-REAL, DIMENSION(4,2) :: mat_M, mat_N, mat_T
-REAL, DIMENSION(2) :: g_vector
-REAL, DIMENSION(4) :: s_vector
-REAL :: ratio
-
-!!
-!! OUTPUT
-!!
-REAL, INTENT(OUT) :: ftot
-
-!!*****************************************************************************************
-
-!
-! Prerequisite
-!
-  pi = 2.*asin(1.)
-  deg2rad = pi/180.
-  if ((theta_s > 90.) .or. (theta_s < 0.)) then
-        print *, 'please set theta_s between 0 and 90', theta_s
-        stop 
-  endif
-
-!
-! Solar Zenith angle (radian) 
-!
-  if (csza .lt. 0.01) then
-      !print *, 'sun below horizon' 
-      !fdir_0=0.
-      fdir=0.
-      fscat_0=0.
-      fscat=0.   
-      fref=0.
-  else
-
-!!
-!! Low incidence fix
-!! 
-!  if (csza .lt. 0.15) csza = 0.15
-
-!
-! 'Slope vs Sun' azimuth (radian)
-!
-  if ( ( (cos(declin)*sin(rho)) .eq. 0.0 ) &
-               .and. &
-     ( ( sin(deg2rad*latitude)*cos(declin)*cos(rho) &
-              -cos(deg2rad*latitude)*sin(declin) ) .eq. 0.0 ) ) then
-    a = deg2rad*psi_s ! some compilator need specfying value for atan2(0,0)  
-  else
-    a = deg2rad*psi_s + atan2(cos(declin)*sin(rho), &
-               sin(deg2rad*latitude)*cos(declin)*cos(rho)-cos(deg2rad*latitude)*sin(declin))
-  end if
-  
-!
-! Cosine of slope-sun phase angle 
-!
-  mu_s = csza*cos(deg2rad*theta_s) - cos(a)*sin(deg2rad*theta_s)*sqrt(1-csza**2)
-  if (mu_s .le. 0.) mu_s=0.
-
-!
-! Sky-view factor
-!
-  sigma_s=0.5*(1.+cos(deg2rad*theta_s))
-
-!
-! Direct flux on the slope
-!
-  fdir = fdir_0 * mu_s/csza
-
-!
-! Reflected flux on the slope
-!
-  fref = albedo * (1-sigma_s) * ftot_0
-
-!
-! Scattered flux on a flat surface
-!
-  fscat_0 = ftot_0 - fdir_0
-
-!
-! Scattering vector (slope vs sky)
-!
-  s_vector=(/ 1., exp(-taudust) , sin(deg2rad*theta_s), sin(deg2rad*theta_s)*exp(-taudust) /)
-
-!
-! Geometry vector (slope vs sun)
-!
-  g_vector=(/ mu_s/csza, 1. /)
-
-!
-! Coupling matrix
-!
-  if (csza .ge. 0.5) then
-        mat_M(:,1) = (/ -0.264,  1.309,  0.208, -0.828 /)
-        mat_M(:,2) = (/  1.291*sigma_s, -1.371*sigma_s, -0.581,  1.641 /)
-        mat_N(:,1) = (/  0.911, -0.777, -0.223,  0.623 /)
-        mat_N(:,2) = (/ -0.933*sigma_s,  0.822*sigma_s,  0.514, -1.195 /)
-
-  else
-        mat_M(:,1) = (/ -0.373,  0.792, -0.095,  0.398 /)
-        mat_M(:,2) = (/  1.389*sigma_s, -0.794*sigma_s, -0.325,  0.183 /)
-        mat_N(:,1) = (/  1.079,  0.275,  0.419, -1.855 /)
-        mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075,  1.844 /)
-  endif
-  !
-  mat_T = mat_M + csza*mat_N
-
-
-!
-! Scattered flux slope ratio
-!
-  if (deg2rad*theta_s <= 0.0872664626) then
-  !
-  ! low angles
-  !
-        s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /)
-        ratio = DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector )
-        ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626
-  else
-  !
-  ! general case
-  !
-        ratio= DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector )  
-                  !
-                  ! NB: ratio= DOT_PRODUCT ( s_vector, MATMUL( mat_T, g_vector ) ) is equivalent
-  endif
-
-!
-! Scattered flux on the slope
-!
-  fscat = ratio * fscat_0
-
-
-endif !! if (csza < 0.01)
-
-!
-! Total flux on the slope
-!
-  ftot = fdir + fref + fscat
-
-!!
-!! Display results
-!!
-!  print *, 'sca component 0 ', ftot_0-fdir_0
-!  print *, 'dir component 0 ', fdir_0
-!  print *, 'scattered component ', fscat
-!  print *, 'direct component ', fdir
-!  print *, 'reflected component ', fref
-
-END SUBROUTINE meso_param_slope
Index: /trunk/LMDZ.MARS/libf/phymars/meso_physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_physiq.F	(revision 226)
+++ /trunk/LMDZ.MARS/libf/phymars/meso_physiq.F	(revision 227)
@@ -136,5 +136,5 @@
 
 #ifdef MESOSCALE
-#include "meso_slope.h"
+#include "slope.h"
 #include "wrf_output_2d.h"
 #include "wrf_output_3d.h"
Index: unk/LMDZ.MARS/libf/phymars/meso_readtesassim.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_readtesassim.F90	(revision 226)
+++ 	(revision )
@@ -1,297 +1,0 @@
-subroutine readtesassim(ngrid,nlayer,zday,pplev,tauref)
-
-! Reading of the dust assimilation file
-
-implicit none
-
-#include "dimensions.h"
-#include "dimphys.h"
-#include "comgeomfi.h"
-#include "netcdf.inc"
-#include "datafile.h"
-#include "callkeys.h"
-
-integer, intent(in) :: ngrid,nlayer
-real, intent(in) :: zday ! date in martian days
-real, dimension(ngrid,nlayer+1), intent(in) :: pplev
-real, dimension(ngrid), intent(out) :: tauref
-
-! Local variables
-
-real :: realday
-integer nid,nvarid,ierr
-integer ltloop,lsloop,iloop,jloop,varloop,ig
-real, dimension(2) :: taubuf
-real tau1(4),tau
-real alt(4)
-integer latp(4),lonp(4)
-integer yinf,ysup,xinf,xsup,tinf,tsup
-real latinf,latsup,loninf,lonsup
-real latintmp,lonintmp,latdeg,londeg
-real colat,dlat,dlon,colattmp
-logical, save :: firstcall=.true.
-logical :: timeflag
-real,save :: radeg,pi
-integer :: timedim,londim,latdim
-real, dimension(:), allocatable, save :: lat,lon,time
-real, dimension(:,:,:), allocatable, save :: tautes
-integer, save :: timelen,lonlen,latlen
-character(len=33) :: filename
-
-realday=mod(zday,669.)
-
-if (firstcall) then
-   firstcall=.false.
-
-   pi=acos(-1.)
-   radeg=180/pi
-   
-   ! assimilated dust file: (NB: iaervar is a common in "callkeys.h")
-   ! iaervar=4 means read dust_tes.nc file
-   ! iaervar=24 means read dust_tes_MY24.nc file
-   ! iaervar=25 means read dust_tes_MY25.nc file
-   ! iaervar=26 means read dust_tes_MY26.nc file
-   if (iaervar.eq.4) then
-     filename="dust_tes.nc"
-   else if (iaervar.eq.24) then
-     filename="dust_tes_MY24.nc" ! which currently is identical to dust_tes.nc
-   else if (iaervar.eq.25) then
-     filename="dust_tes_MY25.nc" 
-   else if (iaervar.eq.26) then
-     filename="dust_tes_MY26.nc"
-   else if ((iaervar.eq.6).or.(iaervar.eq.7)) then ! modified TES scenario
-     filename="modifieddust.nc"
-   endif
-   
-   ! Note: datafile() is defined in "datafile.h"
-   ierr=NF_OPEN(trim(datafile)//"/"//trim(filename),NF_NOWRITE,nid)
-   IF (ierr.NE.NF_NOERR) THEN
-     write(*,*)'Problem opening ',trim(filename),' (in phymars/readtesassim.F90)'
-     write(*,*)'It should be in :',trim(datafile),'/'
-     write(*,*)'1) You can change this directory address in '
-     write(*,*)'   file phymars/datafile.h'
-     write(*,*)'2) If necessary, dust.nc (and other datafiles)'
-     write(*,*)'   can be obtained online on:'
-     write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
-     CALL ABORT
-   ENDIF
-
-   ierr=NF_INQ_DIMID(nid,"Time",timedim)
-   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
-   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
-   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
-   ierr=NF_INQ_DIMID(nid,"longitude",londim)
-   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
-
-   allocate(tautes(lonlen,latlen,timelen))
-   allocate(lat(latlen), lon(lonlen), time(timelen))
-
-   ierr = NF_INQ_VARID (nid, "dustop",nvarid)
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tautes)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, tautes)
-#endif
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT*, "Error: Readtesassim <dustop> not found"
-      stop
-   ENDIF
-   ierr = NF_INQ_VARID (nid, "Time",nvarid)
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
-#endif
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT*, "Error: Readtesassim <Time> not found"
-      stop
-   ENDIF
-
-   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
-#endif
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT*, "Error: Readtesassim <latitude> not found"
-      stop
-   ENDIF
-
-   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
-#ifdef NC_DOUBLE
-   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
-#else
-   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
-#endif
-   IF (ierr .NE. NF_NOERR) THEN
-      PRINT*, "Error: Readtesassim <longitude> not found"
-      stop
-   ENDIF
-
-endif ! of if (firstcall)
-
-do ig=1,ngrid
-
-! Find the four nearest points, arranged as follows:
-!                               1 2
-!                               3 4
-
-   latdeg=lati(ig)*radeg  ! latitude, in degrees
-   londeg=long(ig)*radeg  ! longitude, in degrees east
-   colat=90-latdeg        ! colatitude, in degrees
-
- ! Find encompassing latitudes
-   if (colat<(90-lat(1))) then ! between north pole and lat(1)
-      ysup=1
-      yinf=1
-   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
-      ysup=latlen
-      yinf=latlen
-   else ! general case
-      do iloop=2,latlen
-         if(colat<(90-lat(iloop))) then
-           ysup=iloop-1
-           yinf=iloop
-           exit
-         endif
-      enddo
-   endif
-
-   latinf=lat(yinf)
-   latsup=lat(ysup)
-
-
- ! Find encompassing longitudes
-   ! Note: in input file, lon(1)=-180.
-   if (londeg>lon(lonlen)) then
-      xsup=1
-      xinf=lonlen
-      loninf=lon(xsup)
-      lonsup=180.0 ! since lon(1)=-180.0
-   else
-      do iloop=2,lonlen
-         if(londeg<lon(iloop)) then
-              xsup=iloop
-              xinf=iloop-1
-              exit
-         endif
-      enddo
-      loninf=lon(xinf)
-      lonsup=lon(xsup)
-   endif
-
-   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
-     .OR.(ysup.lt.1)) then
-      write (*,*) "Readtesassim: SYSTEM ERROR on x or y in ts_gcm"
-      write (*,*) "xinf: ",xinf
-      write (*,*) "xsup: ",xsup
-      write (*,*) "yinf: ",yinf
-      write (*,*) "ysup: ",ysup
-      stop
-   endif
-
-!   loninf=lon(xinf)
-!   lonsup=lon(xsup)
-!   latinf=lat(yinf)
-!   latsup=lat(ysup)
-
-! The four neighbouring points are arranged as follows:
-!                               1 2
-!                               3 4
-
-   latp(1)=ysup
-   latp(2)=ysup
-   latp(3)=yinf
-   latp(4)=yinf
-
-   lonp(1)=xinf
-   lonp(2)=xsup
-   lonp(3)=xinf
-   lonp(4)=xsup
-
-! Linear interpolation on time, for all four neighbouring points
-
-   if ((realday<time(1)).or.(realday>time(timelen))) then
-      tinf=timelen
-      tsup=1
-      timeflag=.true.
-   else 
-      do iloop=2,timelen
-         if (realday<time(iloop)) then
-            tinf=iloop-1
-            tsup=iloop
-            exit
-         endif
-      enddo
-   endif
-
-! Bilinear interpolation on the four nearest points
-
-   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
-      dlat=0
-   else
-      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
-   endif
-
-   if (lonsup==loninf) then
-      dlon=0
-   else
-      dlon=(londeg-loninf)/(lonsup-loninf)
-   endif
-
-   do iloop=1,4
-      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
-      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
-      if (timeflag) then
-         if (realday>time(timelen)) then
-            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf))) 
-         else
-            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
-         endif
-      else
-         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
-      endif
-      if (tau1(iloop)<0) then
-          write (*,*) "Readtesassim: SYSTEM ERROR on tau"
-          write (*,*) "utime ",realday
-          write (*,*) "time(tinf) ",time(tinf)
-          write (*,*) "time(tsup) ",time(tsup)
-          write (*,*) "tau1 ",taubuf(1)
-          write (*,*) "tau2 ",taubuf(2)
-          write (*,*) "tau ",tau1(iloop)
-          stop
-      endif
-   enddo
-   timeflag=.false.
-
-   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
-      write (*,*) "Readtesassim: SYSTEM ERROR on dlat or dlon in ts_gcm"
-      write (*,*) "dlat: ",dlat
-      write (*,*) "lat: ",latdeg
-      write (*,*) "dlon: ",dlon
-      write (*,*) "lon: ",londeg
-      stop
-   endif
-
-   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
-
-! Multiply by correction coefficient (value depends on SW radiative transfer)
-! NB: swrtype is set in callkeys.h
-   if (swrtype.eq.1) then ! Fouquart 
-     tauref(ig)=tau*0.825
-   else
-     if (swrtype.eq.2) then ! Toon
-       tauref(ig)=tau*1.3
-        ! special case for enhanced flushing storm scenario
-        if ((iaervar.eq.7).and.(tauref(ig).gt.1.1*1.3)) then
-          tauref(ig)=2.0*tauref(ig)
-        endif
-     else
-       write(*,*) "readtesassim: wrong value for flag swrtype !"
-       stop
-     endif
-   endif ! of if (swrtype.eq.1)
-enddo
-
-end
Index: unk/LMDZ.MARS/libf/phymars/meso_slope.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/meso_slope.h	(revision 226)
+++ 	(revision )
@@ -1,5 +1,0 @@
-      COMMON/comslope/theta_sl(ngridmx),psi_sl(ngridmx)
-
-      REAL theta_sl   ! slope angle versus horizontale (deg)
-      REAL psi_sl     ! orientation of the slope (deg)
-
Index: /trunk/LMDZ.MARS/libf/phymars/param_slope.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/param_slope.F90	(revision 227)
+++ /trunk/LMDZ.MARS/libf/phymars/param_slope.F90	(revision 227)
@@ -0,0 +1,220 @@
+SUBROUTINE param_slope( &
+  !
+  ! INPUTS
+  !
+               &  csza, declin, rho, latitude &
+               &  ,taudust, albedo  &
+               &  ,theta_s, psi_s &
+               &  ,fdir_0, ftot_0  &
+  !
+  ! OUTPUTS
+  !
+               &  ,ftot &
+ )
+
+
+!!*****************************************************************************************
+!
+! SUBROUTINE:
+! param_slope
+! 
+!
+! PURPOSE: 
+! computes total solar irradiance on a given Martian slope 
+!
+!
+! INPUTS:
+! csza		cosine solar zenith angle
+! declin	sun declination (rad)
+! rho		sun right ascension (rad)
+! latitude	latitude (deg)
+! taudust	dust optical depth at reference wavelength 0.67 mic. 
+! albedo	spectrally integrated surface Lambertian reflection albedo 
+! theta_s	slope inclination angle (deg)
+!			0 is horizontal, 90 is vertical
+! phi_s		slope azimuth (deg) 
+!       		0 >> Northward
+!       		90 >> Eastward
+!       		180 >> Southward
+!       		270 >> Westward
+! ftot_0	spectrally integrated total irradiance on an horizontal surface (W/m2)
+! fdir_0	spectrally integrated direct irradiance on an horizontal surface (W/m2)
+!
+!
+! OUTPUTS:
+! ftot		spectrally integrated total irradiance on the slope (W/m2)
+!
+! REFERENCE: 
+! "Fast and accurate estimation of irradiance on Martian slopes"
+! A. Spiga & F. Forget
+! .....
+!
+! AUTHOR: 
+! A. Spiga (spiga@lmd.jussieu.fr)
+! March 2008
+!
+!!*****************************************************************************************
+
+IMPLICIT NONE
+
+!!
+!! INPUT
+!!
+REAL, INTENT(IN) :: csza, declin, rho, latitude, taudust, theta_s, psi_s, albedo, ftot_0 , fdir_0 
+
+!!
+!! LOCAL
+!!
+REAL :: pi, deg2rad
+REAL :: a
+REAL :: mu_s, sigma_s
+REAL :: fdir, fscat, fscat_0, fref
+REAL, DIMENSION(4,2) :: mat_M, mat_N, mat_T
+REAL, DIMENSION(2) :: g_vector
+REAL, DIMENSION(4) :: s_vector
+REAL :: ratio
+
+!!
+!! OUTPUT
+!!
+REAL, INTENT(OUT) :: ftot
+
+!!*****************************************************************************************
+
+!
+! Prerequisite
+!
+  pi = 2.*asin(1.)
+  deg2rad = pi/180.
+  if ((theta_s > 90.) .or. (theta_s < 0.)) then
+        print *, 'please set theta_s between 0 and 90', theta_s
+        stop 
+  endif
+
+!
+! Solar Zenith angle (radian) 
+!
+  if (csza .lt. 0.01) then
+      !print *, 'sun below horizon' 
+      !fdir_0=0.
+      fdir=0.
+      fscat_0=0.
+      fscat=0.   
+      fref=0.
+  else
+
+!!
+!! Low incidence fix
+!! 
+!  if (csza .lt. 0.15) csza = 0.15
+
+!
+! 'Slope vs Sun' azimuth (radian)
+!
+  if ( ( (cos(declin)*sin(rho)) .eq. 0.0 ) &
+               .and. &
+     ( ( sin(deg2rad*latitude)*cos(declin)*cos(rho) &
+              -cos(deg2rad*latitude)*sin(declin) ) .eq. 0.0 ) ) then
+    a = deg2rad*psi_s ! some compilator need specfying value for atan2(0,0)  
+  else
+    a = deg2rad*psi_s + atan2(cos(declin)*sin(rho), &
+               sin(deg2rad*latitude)*cos(declin)*cos(rho)-cos(deg2rad*latitude)*sin(declin))
+  end if
+  
+!
+! Cosine of slope-sun phase angle 
+!
+  mu_s = csza*cos(deg2rad*theta_s) - cos(a)*sin(deg2rad*theta_s)*sqrt(1-csza**2)
+  if (mu_s .le. 0.) mu_s=0.
+
+!
+! Sky-view factor
+!
+  sigma_s=0.5*(1.+cos(deg2rad*theta_s))
+
+!
+! Direct flux on the slope
+!
+  fdir = fdir_0 * mu_s/csza
+
+!
+! Reflected flux on the slope
+!
+  fref = albedo * (1-sigma_s) * ftot_0
+
+!
+! Scattered flux on a flat surface
+!
+  fscat_0 = ftot_0 - fdir_0
+
+!
+! Scattering vector (slope vs sky)
+!
+  s_vector=(/ 1., exp(-taudust) , sin(deg2rad*theta_s), sin(deg2rad*theta_s)*exp(-taudust) /)
+
+!
+! Geometry vector (slope vs sun)
+!
+  g_vector=(/ mu_s/csza, 1. /)
+
+!
+! Coupling matrix
+!
+  if (csza .ge. 0.5) then
+        mat_M(:,1) = (/ -0.264,  1.309,  0.208, -0.828 /)
+        mat_M(:,2) = (/  1.291*sigma_s, -1.371*sigma_s, -0.581,  1.641 /)
+        mat_N(:,1) = (/  0.911, -0.777, -0.223,  0.623 /)
+        mat_N(:,2) = (/ -0.933*sigma_s,  0.822*sigma_s,  0.514, -1.195 /)
+
+  else
+        mat_M(:,1) = (/ -0.373,  0.792, -0.095,  0.398 /)
+        mat_M(:,2) = (/  1.389*sigma_s, -0.794*sigma_s, -0.325,  0.183 /)
+        mat_N(:,1) = (/  1.079,  0.275,  0.419, -1.855 /)
+        mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075,  1.844 /)
+  endif
+  !
+  mat_T = mat_M + csza*mat_N
+
+
+!
+! Scattered flux slope ratio
+!
+  if (deg2rad*theta_s <= 0.0872664626) then
+  !
+  ! low angles
+  !
+        s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /)
+        ratio = DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector )
+        ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626
+  else
+  !
+  ! general case
+  !
+        ratio= DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector )  
+                  !
+                  ! NB: ratio= DOT_PRODUCT ( s_vector, MATMUL( mat_T, g_vector ) ) is equivalent
+  endif
+
+!
+! Scattered flux on the slope
+!
+  fscat = ratio * fscat_0
+
+
+endif !! if (csza < 0.01)
+
+!
+! Total flux on the slope
+!
+  ftot = fdir + fref + fscat
+
+!!
+!! Display results
+!!
+!  print *, 'sca component 0 ', ftot_0-fdir_0
+!  print *, 'dir component 0 ', fdir_0
+!  print *, 'scattered component ', fscat
+!  print *, 'direct component ', fdir
+!  print *, 'reflected component ', fref
+
+END SUBROUTINE param_slope
Index: /trunk/LMDZ.MARS/libf/phymars/slope.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/slope.h	(revision 227)
+++ /trunk/LMDZ.MARS/libf/phymars/slope.h	(revision 227)
@@ -0,0 +1,5 @@
+      COMMON/comslope/theta_sl(ngridmx),psi_sl(ngridmx)
+
+      REAL theta_sl   ! slope angle versus horizontale (deg)
+      REAL psi_sl     ! orientation of the slope (deg)
+
Index: /trunk/MESOSCALE/LMD_MM_MARS/makemeso
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/makemeso	(revision 226)
+++ /trunk/MESOSCALE/LMD_MM_MARS/makemeso	(revision 227)
@@ -436,4 +436,6 @@
                        echo 3 | configure > log_compile 2> log_error 
                        sed -f mars.sed configure.wrf > yeah ; mv -f yeah configure.wrf 
+                       ### the following line is necessary because readtesassim has a line "use netcdf" in it!
+                       sed s+"-L../mars_lmd/libo -llmd"+"-L../mars_lmd/libo -llmd $NETCDF/lib/libnetcdf.a"+g configure.wrf > yeah ; mv -f yeah configure.wrf
                        sed s+"mpif90"+"${WHERE_MPI}/mpif90"+g configure.wrf | sed s+"mpicc"+"${WHERE_MPI}/mpicc"+g | sed s+"mpicc"+"mpicc -DMPI2_SUPPORT"+g > yeah
                                 #| sed s+"fastsse"+"fast"+g > yeah
@@ -582,17 +584,9 @@
   sed s/--xsize--/${physx}/g meso_inc/meso_dimphys.h_ref | sed s/--ysize--/${physy}/g | sed s/--physize--/${physize}/g | sed s/--zsize--/${physz}/g | sed s/--soilsize--/${soilsize}/g > dimphys.h
   head -15 dimphys.h
-     ###
      ### TEST new new phys
      if [[ "${phys}" == "newphys_" ]]
      then
-       #touch callkeys.h.meso  
-       #  ln -sf meso_callkeys.h callkeys.h
-       #touch dustlift.F.meso 
-       #  ln -sf meso_dustlift.F dustlift.F  ## attention avec "ancienne nouvelle physique"
-       touch readtesassim.F90.meso
-         ln -sf meso_readtesassim.F90 readtesassim.F90  ## attention idem
        touch gr_fi_dyn.F.lien
-         ln -sf ../dyn3d/gr_fi_dyn.F .  ## dommage, a corriger ## attention idem
-     ###############################
+         ln -sf ../dyn3d/gr_fi_dyn.F .  ## dommage, a corriger 
      ###############################
 \rm scatterers.h
Index: /trunk/MESOSCALE_DEV/WORK/meso_readtesassim.F90
===================================================================
--- /trunk/MESOSCALE_DEV/WORK/meso_readtesassim.F90	(revision 227)
+++ /trunk/MESOSCALE_DEV/WORK/meso_readtesassim.F90	(revision 227)
@@ -0,0 +1,297 @@
+subroutine readtesassim(ngrid,nlayer,zday,pplev,tauref)
+
+! Reading of the dust assimilation file
+
+implicit none
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comgeomfi.h"
+#include "netcdf.inc"
+#include "datafile.h"
+#include "callkeys.h"
+
+integer, intent(in) :: ngrid,nlayer
+real, intent(in) :: zday ! date in martian days
+real, dimension(ngrid,nlayer+1), intent(in) :: pplev
+real, dimension(ngrid), intent(out) :: tauref
+
+! Local variables
+
+real :: realday
+integer nid,nvarid,ierr
+integer ltloop,lsloop,iloop,jloop,varloop,ig
+real, dimension(2) :: taubuf
+real tau1(4),tau
+real alt(4)
+integer latp(4),lonp(4)
+integer yinf,ysup,xinf,xsup,tinf,tsup
+real latinf,latsup,loninf,lonsup
+real latintmp,lonintmp,latdeg,londeg
+real colat,dlat,dlon,colattmp
+logical, save :: firstcall=.true.
+logical :: timeflag
+real,save :: radeg,pi
+integer :: timedim,londim,latdim
+real, dimension(:), allocatable, save :: lat,lon,time
+real, dimension(:,:,:), allocatable, save :: tautes
+integer, save :: timelen,lonlen,latlen
+character(len=33) :: filename
+
+realday=mod(zday,669.)
+
+if (firstcall) then
+   firstcall=.false.
+
+   pi=acos(-1.)
+   radeg=180/pi
+   
+   ! assimilated dust file: (NB: iaervar is a common in "callkeys.h")
+   ! iaervar=4 means read dust_tes.nc file
+   ! iaervar=24 means read dust_tes_MY24.nc file
+   ! iaervar=25 means read dust_tes_MY25.nc file
+   ! iaervar=26 means read dust_tes_MY26.nc file
+   if (iaervar.eq.4) then
+     filename="dust_tes.nc"
+   else if (iaervar.eq.24) then
+     filename="dust_tes_MY24.nc" ! which currently is identical to dust_tes.nc
+   else if (iaervar.eq.25) then
+     filename="dust_tes_MY25.nc" 
+   else if (iaervar.eq.26) then
+     filename="dust_tes_MY26.nc"
+   else if ((iaervar.eq.6).or.(iaervar.eq.7)) then ! modified TES scenario
+     filename="modifieddust.nc"
+   endif
+   
+   ! Note: datafile() is defined in "datafile.h"
+   ierr=NF_OPEN(trim(datafile)//"/"//trim(filename),NF_NOWRITE,nid)
+   IF (ierr.NE.NF_NOERR) THEN
+     write(*,*)'Problem opening ',trim(filename),' (in phymars/readtesassim.F90)'
+     write(*,*)'It should be in :',trim(datafile),'/'
+     write(*,*)'1) You can change this directory address in '
+     write(*,*)'   file phymars/datafile.h'
+     write(*,*)'2) If necessary, dust.nc (and other datafiles)'
+     write(*,*)'   can be obtained online on:'
+     write(*,*)'   http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
+     CALL ABORT
+   ENDIF
+
+   ierr=NF_INQ_DIMID(nid,"Time",timedim)
+   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
+   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
+   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
+   ierr=NF_INQ_DIMID(nid,"longitude",londim)
+   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
+
+   allocate(tautes(lonlen,latlen,timelen))
+   allocate(lat(latlen), lon(lonlen), time(timelen))
+
+   ierr = NF_INQ_VARID (nid, "dustop",nvarid)
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tautes)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, tautes)
+#endif
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT*, "Error: Readtesassim <dustop> not found"
+      stop
+   ENDIF
+   ierr = NF_INQ_VARID (nid, "Time",nvarid)
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, time)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, time)
+#endif
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT*, "Error: Readtesassim <Time> not found"
+      stop
+   ENDIF
+
+   ierr = NF_INQ_VARID (nid, "latitude",nvarid)
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lat)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lat)
+#endif
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT*, "Error: Readtesassim <latitude> not found"
+      stop
+   ENDIF
+
+   ierr = NF_INQ_VARID (nid, "longitude",nvarid)
+#ifdef NC_DOUBLE
+   ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lon)
+#else
+   ierr = NF_GET_VAR_REAL(nid, nvarid, lon)
+#endif
+   IF (ierr .NE. NF_NOERR) THEN
+      PRINT*, "Error: Readtesassim <longitude> not found"
+      stop
+   ENDIF
+
+endif ! of if (firstcall)
+
+do ig=1,ngrid
+
+! Find the four nearest points, arranged as follows:
+!                               1 2
+!                               3 4
+
+   latdeg=lati(ig)*radeg  ! latitude, in degrees
+   londeg=long(ig)*radeg  ! longitude, in degrees east
+   colat=90-latdeg        ! colatitude, in degrees
+
+ ! Find encompassing latitudes
+   if (colat<(90-lat(1))) then ! between north pole and lat(1)
+      ysup=1
+      yinf=1
+   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
+      ysup=latlen
+      yinf=latlen
+   else ! general case
+      do iloop=2,latlen
+         if(colat<(90-lat(iloop))) then
+           ysup=iloop-1
+           yinf=iloop
+           exit
+         endif
+      enddo
+   endif
+
+   latinf=lat(yinf)
+   latsup=lat(ysup)
+
+
+ ! Find encompassing longitudes
+   ! Note: in input file, lon(1)=-180.
+   if (londeg>lon(lonlen)) then
+      xsup=1
+      xinf=lonlen
+      loninf=lon(xsup)
+      lonsup=180.0 ! since lon(1)=-180.0
+   else
+      do iloop=2,lonlen
+         if(londeg<lon(iloop)) then
+              xsup=iloop
+              xinf=iloop-1
+              exit
+         endif
+      enddo
+      loninf=lon(xinf)
+      lonsup=lon(xsup)
+   endif
+
+   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
+     .OR.(ysup.lt.1)) then
+      write (*,*) "Readtesassim: SYSTEM ERROR on x or y in ts_gcm"
+      write (*,*) "xinf: ",xinf
+      write (*,*) "xsup: ",xsup
+      write (*,*) "yinf: ",yinf
+      write (*,*) "ysup: ",ysup
+      stop
+   endif
+
+!   loninf=lon(xinf)
+!   lonsup=lon(xsup)
+!   latinf=lat(yinf)
+!   latsup=lat(ysup)
+
+! The four neighbouring points are arranged as follows:
+!                               1 2
+!                               3 4
+
+   latp(1)=ysup
+   latp(2)=ysup
+   latp(3)=yinf
+   latp(4)=yinf
+
+   lonp(1)=xinf
+   lonp(2)=xsup
+   lonp(3)=xinf
+   lonp(4)=xsup
+
+! Linear interpolation on time, for all four neighbouring points
+
+   if ((realday<time(1)).or.(realday>time(timelen))) then
+      tinf=timelen
+      tsup=1
+      timeflag=.true.
+   else 
+      do iloop=2,timelen
+         if (realday<time(iloop)) then
+            tinf=iloop-1
+            tsup=iloop
+            exit
+         endif
+      enddo
+   endif
+
+! Bilinear interpolation on the four nearest points
+
+   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
+      dlat=0
+   else
+      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
+   endif
+
+   if (lonsup==loninf) then
+      dlon=0
+   else
+      dlon=(londeg-loninf)/(lonsup-loninf)
+   endif
+
+   do iloop=1,4
+      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
+      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
+      if (timeflag) then
+         if (realday>time(timelen)) then
+            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf))) 
+         else
+            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
+         endif
+      else
+         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
+      endif
+      if (tau1(iloop)<0) then
+          write (*,*) "Readtesassim: SYSTEM ERROR on tau"
+          write (*,*) "utime ",realday
+          write (*,*) "time(tinf) ",time(tinf)
+          write (*,*) "time(tsup) ",time(tsup)
+          write (*,*) "tau1 ",taubuf(1)
+          write (*,*) "tau2 ",taubuf(2)
+          write (*,*) "tau ",tau1(iloop)
+          stop
+      endif
+   enddo
+   timeflag=.false.
+
+   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
+      write (*,*) "Readtesassim: SYSTEM ERROR on dlat or dlon in ts_gcm"
+      write (*,*) "dlat: ",dlat
+      write (*,*) "lat: ",latdeg
+      write (*,*) "dlon: ",dlon
+      write (*,*) "lon: ",londeg
+      stop
+   endif
+
+   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
+
+! Multiply by correction coefficient (value depends on SW radiative transfer)
+! NB: swrtype is set in callkeys.h
+   if (swrtype.eq.1) then ! Fouquart 
+     tauref(ig)=tau*0.825
+   else
+     if (swrtype.eq.2) then ! Toon
+       tauref(ig)=tau*1.3
+        ! special case for enhanced flushing storm scenario
+        if ((iaervar.eq.7).and.(tauref(ig).gt.1.1*1.3)) then
+          tauref(ig)=2.0*tauref(ig)
+        endif
+     else
+       write(*,*) "readtesassim: wrong value for flag swrtype !"
+       stop
+     endif
+   endif ! of if (swrtype.eq.1)
+enddo
+
+end
