Index: /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 3921)
+++ /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 3922)
@@ -5,5 +5,5 @@
 contains
 
-      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, &
+      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt,pq,zls, &
          aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col, &
          cloudfrac,totcloudfrac,clearsky)
@@ -20,5 +20,5 @@
        use geometry_mod, only: latitude
        use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl,	&
-		CLFvarying,CLFfixval,dusttau,			 	&
+		CLFvarying,CLFfixval,dusttau,timedepdust,		 	&
 		pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,	&
 		pres_bottom_strato,pres_top_strato,obs_tau_col_strato,  &
@@ -69,4 +69,5 @@
       REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
+      REAL,INTENT(IN) :: zls ! Stellar longitude (rad)
       REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K)
       REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
@@ -100,8 +101,12 @@
       CHARACTER(LEN=20) :: tracername ! to temporarily store text
 
-      ! for fixed dust profiles
+      ! for dust profiles
       real topdust, expfactor, zp
       REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
       REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
+
+      ! time-dependent dust (MM)
+      real zlsconst, odpref, taueq, tauS, tauN
+      real tau_pref_MGS(ngrid), tauscaling(ngrid)
 
       real CLFtot
@@ -311,4 +316,9 @@
 !==================================================================
 !             Dust 
+!             Either constant/homogeneous or
+!             following MGS scenario for
+!             present-day Mars as per:
+!             Montmessin et al., 2004
+!             (DOI: 10.1029/2004JE002284)
 !==================================================================
       if (iaero_dust.ne.0) then
@@ -316,12 +326,63 @@
 !         1. Initialization 
           aerosol(1:ngrid,1:nlayer,iaer)=0.0
-          
-          topdust=30.0 ! km  (used to be 10.0 km) LK
-
-!       2. Opacity calculation
+
+!       2. Opacity calculation
+
+          IF (timedepdust) THEN
+!           Time-dependent dust (MGS scenarion for present-day Mars)
+
+            zlsconst = sin(zls-2.76)
+            taudusttmp(:) = 0
+            odpref = 610. ! Reference pressure (Pa) of
+                          ! DOD (Dust optical Depth) tau_pref_*
+
+            DO l=1,nlayer-1
+              DO ig=1,ngrid
+
+                  topdust = 60.+18.*zlsconst                     & ! From
+                    - (32.+18.*zlsconst)*(sin(latitude(ig)))**4  & ! Montmessin
+                    -  8.*zlsconst*(sin(latitude(ig)))**5          ! et al. 2004
+                  if (pplay(ig,l).ge.odpref/(988.**(topdust/70.))) then ! What is the use of this line?
+                    zp = (odpref/pplay(ig,l))**(70./topdust)
+                    expfactor = max(exp(0.007*(1.-max(zp,1.))),1.e-3)
+                  else
+                    expfactor = 1.e-3
+                  endif
+
+!                 Vertical scaling function
+                  aerosol(ig,l,iaer) = (pplev(ig,l)-pplev(ig,l+1)) &
+                                     *  expfactor
+
+!                 Horizontal scaling of the dust opacity
+                  if (l==1) then
+
+                    taueq = 0.2 + (0.5-0.2) * (cos(0.5*(zls-4.363)))**14
+                    tauS  = 0.1 + (0.5-0.1) * (cos(0.5*(zls-4.363)))**14
+                    tauN  = 0.1
+
+                    if (latitude(ig).ge.0) then
+                    ! Northern hemisphere
+                      tau_pref_MGS(ig) = tauN + (taueq-tauN)*0.5 &
+                             *(1+tanh((45-latitude(ig)*180./pi)*6/60))
+                    else
+                    ! Southern hemisphere
+                      tau_pref_MGS(ig) = tauS + (taueq-tauS)*0.5 &
+                             *(1+tanh((45+latitude(ig)*180./pi)*6/60))
+                    endif
+                  endif
+
+              ENDDO
+            ENDDO
+
+          ELSE
+!           Fixed dust
 
 !           expfactor=0.
-           DO l=1,nlayer-1
-             DO ig=1,ngrid
+            topdust=30.0 ! km  (used to be 10.0 km) LK
+
+            DO l=1,nlayer-1
+              DO ig=1,ngrid
+
+            
 !             Typical mixing ratio profile
 
@@ -329,11 +390,13 @@
                  expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
 
-!             Vertical scaling function
-              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
-               *expfactor
+!               Vertical scaling function
+                aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
+                 *expfactor
 
 
              ENDDO
            ENDDO
+          ENDIF ! of if timedepdust
+
 
 !          Rescaling each layer to reproduce the choosen (or assimilated)
@@ -348,14 +411,26 @@
                 ENDDO
               ENDDO
+    
+            if (timedepdust) then
+!             Dust opacity scaling
+              tauscaling(:) = tau_pref_MGS(:) * pplev(:,1) / odpref
+            else
+              tauscaling(:) = 1
+            endif
+
             DO l=1,nlayer-1
                DO ig=1,ngrid
-                  aerosol(ig,l,iaer) = max(1E-20, &
-                          dusttau &
-                       *  pplev(ig,1) / pplev(ig,1) &
+                aerosol(ig,l,iaer) = max(1E-20, &
+                          dusttau * tauscaling(ig) &
+                       *  pplev(ig,1) / pplev(ig,1) & ! what is the use of this line ? (MM)
                        *  aerosol(ig,l,iaer) &
                        /  taudusttmp(ig))
+                
 
               ENDDO
             ENDDO
+            
+            call writediagfi(ngrid,"taudust","Optical depth at pref","-",2, dusttau * tauscaling)
+
       end if ! If dust aerosol   
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 3921)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 3922)
@@ -5,5 +5,5 @@
 CONTAINS
 
-      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
+      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zls,       &
           albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    & 
           tsurf,fract,dist_star,aerosol,muvar,                 &
@@ -80,4 +80,5 @@
       INTEGER,INTENT(IN) :: nq                     ! Number of tracers.
       REAL,INTENT(IN) :: qsurf(ngrid,nq)           ! Tracers on surface (kg.m-2).
+      REAL,INTENT(IN) :: zls                       ! Stellar longitude (rad).
       REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)   ! Spectral Short Wavelengths Albedo. By MT2015
       REAL,INTENT(IN) :: emis(ngrid)               ! Long Wave emissivity.
@@ -554,5 +555,5 @@
 
       ! Get aerosol optical depths.
-      call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
+      call aeropacity(ngrid,nlayer,nq,pplay,pplev,pt,pq,zls,aerosol,      &
            reffrad,nueffrad,QREFvis3d,QREFir3d,                             & 
            tau_col,cloudfrac,totcloudfrac,clearsky)                
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3921)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3922)
@@ -121,7 +121,8 @@
       real,save :: topdustref
       real,save :: dusttau
+      logical,save :: timedepdust
       real,save :: Fat1AU
       real,save :: stelTbb
-!$OMP THREADPRIVATE(topdustref,dusttau,Fat1AU,stelTbb)
+!$OMP THREADPRIVATE(topdustref,dusttau,timedepdust,Fat1AU,stelTbb)
       real,save :: Tstrat
       real,save :: tplanet
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3921)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3922)
@@ -723,4 +723,9 @@
      call getin_p("dusttau",dusttau)
      if (is_master) write(*,*)trim(rname)//": dusttau = ",dusttau
+     
+     if (is_master) write(*,*)trim(rname)//": Use time-dependent dust ?:"
+     timedepdust=.false. ! default value
+     call getin_p("timedepdust",timedepdust)
+     if (is_master) write(*,*)trim(rname)//": timedepdust = ",timedepdust
 
      if (is_master) write(*,*)trim(rname)//": Radiatively active CO2 aerosols?"
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3921)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3922)
@@ -1058,5 +1058,5 @@
                ! standard callcorrk
                clearsky=.false.
-               call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
+               call callcorrk(ngrid,nlayer,pq,nq,qsurf,zls,                        &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
                               tsurf,fract,dist_star,aerosol,muvar,                &
@@ -1089,5 +1089,5 @@
                   ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ...
                   clearsky=.true.
-                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
+                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,zls,                       &
                                  albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt,  &
                                  tsurf,fract,dist_star,aerosol,muvar,                &