Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 1973)
+++ trunk/LMDZ.MARS/README	(revision 1974)
@@ -2613,2 +2613,7 @@
 Cosmetic changes:
 - variables suffixes for sub-grid scale cloud scheme changed from --1 to --clf
+
+== 18/07/2018 == EM+MV
+Integration of the detached dust layer parametrizations (rocket dust storm, slope wind lifting, CW, and dust injection scheme, DB).
+Still experimental, default behaviour (rdstorm=.false., dustinjection=0) identical to previous revision.
+NB: Updated newstart requires an updated "surface.nc" containing the "hmons" field.
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/avg_horiz_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/avg_horiz_mod.F	(revision 1974)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/avg_horiz_mod.F	(revision 1974)
@@ -0,0 +1,94 @@
+      MODULE avg_horiz_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+      
+      subroutine avg_horiz(imo,jmo,imn,jmn,rlonuo,rlatvo,rlonv,
+     .             rlatu,zdata,outdata)
+
+c=======================================================================
+c Maximum/minimum value composite approach for interpolate the variables
+c  from surface.nc (summit, base... ) 
+c  Chao Wang
+c=======================================================================
+
+      implicit none
+      !inputs/outputs
+      integer, intent(in)  ::  imo,jmo
+      integer, intent(in)  :: imn,jmn
+      real, intent(in)  :: rlonv(imn+1)    
+      real, intent(in)  :: rlatu(jmn+1)   
+      real, intent(in)  :: zdata(imo*jmo)
+      real, intent(out)  :: outdata(imn+1,jmn+1)
+
+      !local variables:
+      real pi
+      integer i,j,l,m,n
+      real tmplon(imn+2)
+      real tmplat(jmn+2)
+      real rlonuo(imo)    
+      real rlatvo(jmo)   
+
+      integer iik(imo,jmo),jjk(imo,jmo)
+      integer tmp
+      integer nmesh
+      
+c     initialize 
+      iik(:,:)=0
+      jjk(:,:)=0
+      outdata(:,:)=0.
+      pi=2.*ASIN(1.)
+
+      tmplon(1)=rlonv(1)
+      do i=2,imn+1
+        tmplon(i)= (rlonv(i-1)+rlonv(i))/2.
+      enddo
+      tmplon(imn+2)=rlonv(imn+1)
+
+      tmplat(1)=rlatu(1)
+      do j=2,jmn+1
+        tmplat(j)=(rlatu(j-1)+rlatu(j))/2.
+      enddo
+      tmplat(jmn+2)=rlatu(jmn+1)
+
+
+      do i=1,imo 
+        do j=1,jmo
+          do m=1,imn+1
+            do n=1,jmn+1
+              
+              if ( ((rlonuo(i) .ge. tmplon(m)) .and.
+     &                       (rlonuo(i) .lt. tmplon(m+1))) .and.
+     &           ((rlatvo(j) .lt. tmplat(n) ) .and.
+     &            (rlatvo(j) .ge. tmplat(n+1) ))) then
+                iik(i,j)=m
+                jjk(i,j)=n
+                goto 1000
+              endif
+              
+            enddo 
+          enddo
+1000      continue
+        enddo
+      enddo
+
+      do m=1,imn+1
+        do n=1,jmn+1
+          nmesh=0
+          do i=1,imo 
+            do j=1,jmo
+               if ((iik(i,j) .eq. m) .and. (jjk(i,j) .eq. n)) then
+                  outdata(m,n)=outdata(m,n)+zdata(i+imo*(j-1)) 
+                  nmesh=nmesh+1
+               endif       
+
+            enddo 
+          enddo
+          outdata(m,n)=outdata(m,n)/nmesh
+        enddo
+      enddo
+      
+      END SUBROUTINE avg_horiz
+
+      END MODULE avg_horiz_mod
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F	(revision 1974)
@@ -1,5 +1,6 @@
 c=======================================================================
       SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0,
-     &                    zmea,zstd,zsig,zgam,zthe)
+     &                    zmea,zstd,zsig,zgam,zthe,
+     &                    hmons,summit,zavg)
 c=======================================================================
 c
@@ -42,7 +43,10 @@
 c=======================================================================
 
-      use ioipsl_getincom, only: getin
-      USE comconst_mod, ONLY: g,pi
+! to use  'getin'
+      use ioipsl_getincom, only: getin 
+      use comconst_mod, only: g,pi
       use datafile_mod, only: datadir
+      use avg_horiz_mod, only: avg_horiz
+      use mvc_horiz_mod, only: mvc_horiz
 
       implicit none
@@ -74,5 +78,8 @@
       REAL,intent(out) :: zgam(imdp1*jmdp1)
       REAL,intent(out) :: zthe(imdp1*jmdp1)
-
+      REAL,intent(out) :: hmons(imdp1*jmdp1) !CW17,361*180 hmons
+      REAL,intent(out) :: summit(imdp1*jmdp1) !CW17,361*180 hmons
+      REAL,intent(out) :: zavg(imdp1*jmdp1) !CW17,361*180 hmons
+      
 ! Local variables:
       REAL        zdata(imd*jmdp1)
@@ -95,5 +102,5 @@
 
       CHARACTER*20 string
-      DIMENSION string(0:4)
+      DIMENSION string(0:6)
 
 
@@ -233,4 +240,19 @@
           call grid_noro1(360, 180, longitude, latitude, zdata,
      .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
+     
+          !CW17
+          call avg_horiz(imd,jmdp1,iim,jjm,longitude,
+     .                latitude,rlonv,rlatu,zdata,pfield)
+
+          do j=1,jjp1 !CW17 49, iimp1=65, the last column = first column
+             pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
+          enddo 
+          do i=1,iimp1*jjp1
+c             if (pfield(i) .ne. -999999.) then
+                zavg(i) = pfield(i)
+c             else
+c               zavg(i)=-999999.
+c             endif
+          enddo
 
       endif
@@ -297,4 +319,95 @@
 c    FIN
 c-----------------------------------------------------------------------
-
+ 
+c=======================================================================
+c<<<<<<<<<<<<<<<<<<<<<<<add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! name of dataset
+      string(5) = 'hmons' !subgrid hmons
+c     the following data could be useful in future, but currently,
+c     we don't need to put them into restartfi.nc
+      string(6) = 'summit' !subgrid summit
+c     string(7) = 'base'   !base of summit (not subgrid)
+c     string(8) = 'bottom' !subgrid base
+      do k=5,6
+          write(*,*) 'string',k,string(k)  
+c-----------------------------------------------------------------------
+c    Lecture NetCDF  
+c-----------------------------------------------------------------------
+
+          ierr = NF_INQ_VARID (unit, string(k), nvarid)
+          if (ierr.ne.nf_noerr) then
+            write(*,*) 'datareadnc error, cannot find ',trim(string(k))
+            write(*,*) ' in file ',trim(datadir),'/surface.nc'
+            stop
+          endif
+
+c-----------------------------------------------------------------------
+c    initialisation
+c-----------------------------------------------------------------------
+          call initial0(iimp1*jjp1,pfield)
+          call initial0(imd*jmdp1,zdata)
+          call initial0(imdp1*jmdp1,zdataS)
+
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
+#else
+          ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
+#endif
+          if (ierr.ne.nf_noerr) then
+            write(*,*) 'datareadnc: error failed loading ',
+     .                  trim(string(k))
+            stop
+          endif
+
+c-----------------------------------------------------------------------
+c   Passage de zdata en grille (imdp1 jmdp1)
+c-----------------------------------------------------------------------
+          do j=1,jmdp1     !  180
+              do i=1,imd   !  360
+                  !copy zdata to zdataS, line by line
+                  !      i+ 361 *(j-1)          i+  360*(j-1)
+                  zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1))
+              enddo
+              ! the last column = the first column of zdata
+              zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1))
+          enddo
+
+          if (k .eq. 5 .or. k .eq. 6) then 
+              ! hmons, summit, keep the maximum of subgrids 
+              call mvc_horiz(imd,jmdp1,iim,jjm,longitude,
+     .                latitude,rlonv,rlatu,zdata,pfield)
+          endif
+
+
+c-----------------------------------------------------------------------
+c    Periodicite    
+c-----------------------------------------------------------------------
+
+          do j=1,jjp1 ! 49, iimp1=65, the last column = first column
+             pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
+          enddo 
+ 
+c-----------------------------------------------------------------------
+c    Sauvegarde des champs    
+c-----------------------------------------------------------------------
+
+          if (k.eq.5) then ! hmons
+             do i=1,iimp1*jjp1
+                  if (pfield(i) .ne. -999999.) then
+                    hmons(i) = pfield(i)
+                  else
+                    hmons(i)=0.
+                  endif
+             enddo
+          endif
+          
+          if (k.eq.6) then  ! summit
+             do i=1,iimp1*jjp1
+                    summit(i) = pfield(i)
+             enddo
+          endif
+
+      enddo
+c<<<<<<<<<<<<<<<<<<<<<<<<<done add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>>
+c=======================================================================
       END
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/mvc_horiz_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/mvc_horiz_mod.F	(revision 1974)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/mvc_horiz_mod.F	(revision 1974)
@@ -0,0 +1,93 @@
+      MODULE mvc_horiz_mod
+      
+      implicit none
+      
+      contains      
+
+      subroutine mvc_horiz(imo,jmo,imn,jmn,rlonuo,rlatvo,rlonv,
+     .             rlatu,zdata,outdata)
+
+c=======================================================================
+c Maximum/minimum value composite approach for interpolate the variables
+c  from surface.nc (summit, base... ) 
+c  Chao Wang
+c=======================================================================
+
+      implicit none
+      
+      !inputs/ouputs
+      integer, intent(in) :: imo,jmo
+      integer, intent(in) :: imn,jmn
+      real, intent(in)    ::  rlonv(imn+1)    
+      real, intent(in)    ::  rlatu(jmn+1)   
+      real, intent(in)    ::  zdata(imo*jmo)
+      real, intent(out)   ::   outdata(imn+1,jmn+1)
+      
+      !local variables:
+      real pi
+      integer i,j,l,m,n
+      real tmplon(imn+2)
+      real tmplat(jmn+2)
+      REAL rlonuo(imo)    
+      REAL rlatvo(jmo)   
+
+      integer iik(imo,jmo),jjk(imo,jmo)
+      integer tmp
+
+c     initialize 
+      iik(:,:)=0
+      jjk(:,:)=0
+      outdata(:,:)=-999999. !lowest point on Mars ~8200m
+
+      pi=2.*ASIN(1.)
+
+      tmplon(1)=rlonv(1)
+      do i=2,imn+1
+        tmplon(i)= (rlonv(i-1)+rlonv(i))/2.
+      enddo
+      tmplon(imn+2)=rlonv(imn+1)
+
+      tmplat(1)=rlatu(1)
+      do j=2,jmn+1
+        tmplat(j)=(rlatu(j-1)+rlatu(j))/2.
+      enddo
+      tmplat(jmn+2)=rlatu(jmn+1)
+
+
+      do i=1,imo 
+        do j=1,jmo
+          do m=1,imn+1
+            do n=1,jmn+1
+              
+              if ( ((rlonuo(i) .ge. tmplon(m)) .and.
+     &                       (rlonuo(i) .lt. tmplon(m+1))) .and.
+     &           ((rlatvo(j) .lt. tmplat(n) ) .and.
+     &            (rlatvo(j) .ge. tmplat(n+1) ))) then
+                iik(i,j)=m
+                jjk(i,j)=n
+                goto 1000
+              endif
+              
+            enddo 
+          enddo
+1000      continue
+        enddo
+      enddo
+
+      do m=1,imn+1
+        do n=1,jmn+1
+          do i=1,imo 
+            do j=1,jmo
+
+               if ((iik(i,j) .eq. m) .and. (jjk(i,j) .eq. n)) then
+                  outdata(m,n)=max(outdata(m,n),zdata(i+imo*(j-1))  )
+               endif       
+
+            enddo 
+          enddo
+        enddo
+      enddo
+
+      end subroutine mvc_horiz
+      
+      END MODULE mvc_horiz_mod
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 1974)
@@ -27,5 +27,5 @@
       use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
      &                     albedodat, z0_default, qsurf, tsurf,
-     &                     co2ice, emis
+     &                     co2ice, emis, hmons
       use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
       use control_mod, only: day_step, iphysiq, anneeref, planet_type
@@ -97,4 +97,7 @@
       real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
       real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
+      real hmonsS(iip1,jjp1)
+      real summitS(iip1,jjp1)
+      real zavgS(iip1,jjp1)
       real z0S(iip1,jjp1)
 
@@ -380,5 +383,6 @@
 
       CALL datareadnc(relief,phis,alb,surfith,z0S,
-     &          zmeaS,zstdS,zsigS,zgamS,ztheS)
+     &          zmeaS,zstdS,zsigS,zgamS,ztheS,
+     &          hmonsS,summitS,zavgS)
 
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
@@ -392,4 +396,7 @@
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
+      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,hmonsS,hmons)
+!      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,summitS,summit)
+!      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zavgS,zavg)
 
       endif ! of if (choix_1.ne.1)
@@ -1642,5 +1649,6 @@
       call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
      .              nqtot,dtphys,real(day_ini),0.0,
-     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
+     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
+     .              hmons)
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      &              dtphys,hour_ini,
Index: trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 1974)
@@ -7,12 +7,13 @@
       SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
      &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
-     &    nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
+     &    QREFvis3d,QREFir3d,omegaREFir3d,
+     &    totstormfract,clearatm,dsords,
      &    clearsky,totcloudfrac)
-                                                   
+                                                         
 ! to use  'getin'
       USE ioipsl_getincom, only: getin
       use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
      &                      igcm_dust_submicron, rho_dust, rho_ice,
-     &                      nqdust
+     &                      nqdust, igcm_stormdust_mass
       use geometry_mod, only: latitude ! grid point latitudes (rad)
       use comgeomfi_h, only: sinlat ! sines of grid point latitudes
@@ -26,5 +27,7 @@
      &            iaerdust,tauvis,
      &            iaer_dust_conrath,iaer_dust_doubleq,
-     &            iaer_dust_submicron,iaer_h2o_ice
+     &            iaer_dust_submicron,iaer_h2o_ice,
+     &            iaer_stormdust_doubleq
+      USE calcstormfract_mod
        IMPLICIT NONE
 c=======================================================================
@@ -56,5 +59,4 @@
 c   QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
 c   QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
-c   omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
 c   omegaREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
 c
@@ -65,8 +67,8 @@
 c   aerosol      aerosol(ig,l,1) is the dust optical
 c                depth in layer l, grid point ig
-
+c   taualldust   CW17 total opacity for all dust scatterer stormdust included
 c
 c=======================================================================
-#include "callkeys.h"
+      include "callkeys.h"
 
 c-----------------------------------------------------------------------
@@ -77,25 +79,27 @@
 c    Input/Output
 c    ------------
-      INTEGER ngrid,nlayer,nq
-
-      REAL ls,zday,expfactor    
-      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
-      REAL pq(ngrid,nlayer,nq)
-      REAL tauref(ngrid), tau(ngrid,naerkind)
-      REAL aerosol(ngrid,nlayer,naerkind)
-      REAL dsodust(ngrid,nlayer)
-      REAL reffrad(ngrid,nlayer,naerkind)
-      REAL nueffrad(ngrid,nlayer,naerkind)
-      REAL QREFvis3d(ngrid,nlayer,naerkind)
-      REAL QREFir3d(ngrid,nlayer,naerkind)
-      REAL omegaREFvis3d(ngrid,nlayer,naerkind)
-      REAL omegaREFir3d(ngrid,nlayer,naerkind)
-      real,intent(in) :: totcloudfrac(ngrid) ! total cloud fraction
-      logical,intent(in) :: clearsky ! true for part without clouds, 
-                                     ! false for part with clouds (total or sub-grid clouds)
-      real :: CLFtot ! total cloud fraction
+      INTEGER, INTENT(IN) ::  ngrid,nlayer,nq
+      REAL, INTENT(IN) :: ls,zday
+      REAL, INTENT(IN) ::  pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
+      REAL, INTENT(IN) ::  pq(ngrid,nlayer,nq)
+      REAL, INTENT(OUT) :: tauref(ngrid)
+      REAL, INTENT(OUT) :: tau(ngrid,naerkind)
+      REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
+      REAL, INTENT(OUT) ::  dsodust(ngrid,nlayer)
+      REAL, INTENT(OUT) ::  dsords(ngrid,nlayer) !dso of stormdust  
+      REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
+      REAL, INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind)
+      REAL, INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
+      REAL, INTENT(IN) ::  omegaREFir3d(ngrid,nlayer,naerkind)
+      LOGICAL, INTENT(IN) :: clearatm
+      REAL, INTENT(IN) :: totstormfract(ngrid)
+      REAL, INTENT(OUT) ::  tauscaling(ngrid) ! Scaling factor for qdust and Ndust
+      REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total cloud fraction
+      LOGICAL,INTENT(IN) :: clearsky ! true for part without clouds,false for part with clouds (total or sub-grid clouds)
 c
 c    Local variables :
 c    -----------------
+      REAL CLFtot ! total cloud fraction
+      real expfactor 
       INTEGER l,ig,iq,i,j
       INTEGER iaer           ! Aerosol index
@@ -110,7 +114,10 @@
 c       for dust or water ice particles in the radiative transfer
 c       (see callradite.F for more information).
-      REAL taudusttmp(ngrid)! Temporary dust opacity
-                               !   used before scaling
-      REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
+      REAL taudusttmp(ngrid)! Temporary dust opacity used before scaling
+      REAL taubackdusttmp(ngrid)! Temporary background dust opacity used before scaling
+      REAL taualldust(ngrid)! dust opacity all dust
+      REAL taudust(ngrid)! dust opacity dust doubleq
+      REAL taustormdust(ngrid)! dust opacity stormdust doubleq
+      REAL taustormdusttmp(ngrid)! dust opacity stormdust doubleq before tauscaling
       REAL taudustvis(ngrid) ! Dust opacity after scaling
       REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
@@ -157,4 +164,6 @@
 
       tau(1:ngrid,1:naerkind)=0
+      dsords(:,:)=0. !CW17: initialize dsords
+      dsodust(:,:)=0.
 
 ! identify tracers
@@ -166,5 +175,6 @@
         DO iaer=1,naerkind
           txt=name_iaer(iaer)
-          IF (txt(1:4).eq."dust") THEN
+        ! CW17: choice tauscaling for stormdust or not
+          IF ((txt(1:4).eq."dust").OR.(txt(1:5).eq."storm")) THEN
             naerdust=naerdust+1
             iaerdust(naerdust)=iaer
@@ -213,8 +223,9 @@
         call getin("tauvis",tauvis)
 
-        IF (freedust) THEN
+        IF (freedust.or.rdstorm) THEN ! if rdstorm no need to held opacity constant at the first levels
           cstdustlevel = 1
         ELSE
-          cstdustlevel = cstdustlevel0
+          cstdustlevel = cstdustlevel0 !Opacity in the first levels is held constant to 
+                                       !avoid unrealistic values due to constant lifting
         ENDIF
 
@@ -374,8 +385,8 @@
      &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
                 ! DENSITY SCALED OPACITY IN INFRARED:
-                dsodust(ig,l) =
-     &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
-     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
-     &          pq(ig,cstdustlevel,igcm_dust_mass)
+!                dsodust(ig,l) =
+!     &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
+!     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+!     &          pq(ig,cstdustlevel,igcm_dust_mass)
               ENDDO
             ELSE
@@ -387,8 +398,8 @@
      &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
                 ! DENSITY SCALED OPACITY IN INFRARED:
-                dsodust(ig,l) =
-     &          (  0.75 * QREFir3d(ig,l,iaer) /
-     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
-     &          pq(ig,l,igcm_dust_mass)
+!                dsodust(ig,l) =
+!     &          (  0.75 * QREFir3d(ig,l,iaer) /
+!     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+!     &          pq(ig,l,igcm_dust_mass)
               ENDDO
             ENDIF
@@ -473,19 +484,36 @@
         ENDIF ! end (clearsky)
 
-c       3. Outputs -- Now done in physiq.F
-!        IF (ngrid.NE.1) THEN
-!          CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
-!     &      ' ',2,taucloudvis)
-!          CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
-!     &      ' ',2,taucloudtes)
-!          IF (callstats) THEN
-!            CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
-!     &        ' ',2,taucloudvis)
-!            CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
-!     &        ' ',2,taucloudtes)
-!          ENDIF
-!        ELSE
-!         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
-!        ENDIF
+c==================================================================
+        CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for 
+c       stormdust  (transport of mass and number mixing ratio) 
+c==================================================================
+c       aerosol is calculated twice : once within the dust storm (clearatm=false)
+c       and once in the part of the mesh without dust storm (clearatm=true)
+        aerosol(1:ngrid,1:nlayer,iaer) = 0.
+        IF (clearatm) THEN  ! considering part of the mesh without storm
+          aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
+        ELSE  ! part of the mesh with concentred dust storm
+          DO l=1,nlayer
+             IF (l.LE.cstdustlevel) THEN
+c          Opacity in the first levels is held constant to 
+c           avoid unrealistic values due to constant lifting:
+               DO ig=1,ngrid
+                  aerosol(ig,l,iaer) = 
+     &           (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
+     &           ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &           pq(ig,cstdustlevel,igcm_stormdust_mass) *
+     &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
+               ENDDO
+             ELSE
+              DO ig=1,ngrid
+                 aerosol(ig,l,iaer) =
+     &           (  0.75 * QREFvis3d(ig,l,iaer) /
+     &           ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &           pq(ig,l,igcm_stormdust_mass) *
+     &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
+              ENDDO
+             ENDIF
+          ENDDO
+        ENDIF
 c==================================================================
         END SELECT aerkind
@@ -498,4 +526,5 @@
 c   is originally scaled to an equivalent odpref Pa pressure surface.
 c -----------------------------------------------------------------
+
 
 #ifdef DUSTSTORM
@@ -640,5 +669,22 @@
       IF (freedust) THEN
           tauscaling(:) = 1.
-
+c        opacity obtained with stormdust
+        IF (rdstorm) THEN
+           taustormdusttmp(1:ngrid)=0.
+           DO l=1,nlayer
+             DO ig=1,ngrid
+                taustormdusttmp(ig) = taustormdusttmp(ig)+ 
+     &            aerosol(ig,l,iaerdust(2))
+             ENDDO
+           ENDDO
+           !opacity obtained with background dust only
+           taubackdusttmp(1:ngrid)=0.  
+           DO l=1,nlayer
+             DO ig=1,ngrid
+                taubackdusttmp(ig) = taubackdusttmp(ig)+
+     &           aerosol(ig,l,iaerdust(1))
+             ENDDO
+           ENDDO
+        ENDIF !rdsstorm
       ELSE
 c       Temporary scaling factor
@@ -669,10 +715,4 @@
      &                aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
           ENDDO
-        ENDDO
-      ENDDO
-
-      DO l=1,nlayer
-        DO ig=1,ngrid
-          dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
         ENDDO
       ENDDO
@@ -700,17 +740,5 @@
         tauref(:) = tauref(:) * odpref / pplev(:,1)
       ENDIF
-      
-c output for debug
-c        IF (ngrid.NE.1) THEN
-c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
-c     &      '#',2,taudusttmp)
-c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
-c     &      '#',2,tauscaling)
-c        ELSE
-c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
-c     &      '#',0,taudusttmp)
-c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
-c     &      '#',0,tauscaling)
-c        ENDIF
+
 c -----------------------------------------------------------------
 c Column integrated visible optical depth in each point
@@ -724,4 +752,38 @@
       ENDDO
 
+c     for diagnostics: opacity for all dust scatterers stormdust included 
+      taualldust(1:ngrid)=0.
+      DO iaer=1,naerdust
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            taualldust(ig) = taualldust(ig) + 
+     &                         aerosol(ig,l,iaerdust(iaer))
+          ENDDO
+        ENDDO
+      ENDDO
+     
+      IF (rdstorm) THEN
+ 
+c     for diagnostics: opacity for dust in background only  
+       taudust(1:ngrid)=0.
+        DO l=1,nlayer
+         DO ig=1,ngrid
+           taudust(ig) = taudust(ig) + 
+     &                       aerosol(ig,l,iaer_dust_doubleq)
+         ENDDO
+        ENDDO
+
+c     for diagnostics: opacity for dust in storm only  
+       taustormdust(1:ngrid)=0.
+        DO l=1,nlayer
+         DO ig=1,ngrid
+           taustormdust(ig) = taustormdust(ig) + 
+     &                       aerosol(ig,l,iaer_stormdust_doubleq)
+         ENDDO
+        ENDDO
+ 
+      ENDIF
+      
+
 #ifdef DUSTSTORM
       IF (firstcall) THEN
@@ -733,39 +795,42 @@
 c Density scaled opacity and column opacity output
 c -----------------------------------------------------------------
-c      dsodust(1:ngrid,1:nlayer) = 0.
-c      DO iaer=1,naerdust
-c        DO l=1,nlayer
-c          DO ig=1,ngrid
-c            dsodust(ig,l) = dsodust(ig,l) +
-c     &                      aerosol(ig,l,iaerdust(iaer)) * g /
-c     &                      (pplev(ig,l) - pplev(ig,l+1))
-c          ENDDO
-c        ENDDO
-c        IF (ngrid.NE.1) THEN
-c          write(txt2,'(i1.1)') iaer
-c          call WRITEDIAGFI(ngrid,'taudust'//txt2,
-c     &                    'Dust col opacity',
-c     &                    ' ',2,tau(1,iaerdust(iaer)))
-c          IF (callstats) THEN
-c            CALL wstats(ngrid,'taudust'//txt2,
-c     &                 'Dust col opacity',
-c     &                 ' ',2,tau(1,iaerdust(iaer)))
-c          ENDIF
-c        ENDIF
-c      ENDDO
-
-c      IF (ngrid.NE.1) THEN
-c       CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
-c    &                    'm2.kg-1',3,dsodust)
-c        IF (callstats) THEN
-c          CALL wstats(ngrid,'dsodust',
-c     &               'tau*g/dp',
-c     &               'm2.kg-1',3,dsodust)
-c        ENDIF
-c      ELSE
-c        CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
-c     &                   dsodust)
-c      ENDIF ! of IF (ngrid.NE.1)
-c -----------------------------------------------------------------
+          IF (rdstorm) then
+            DO l=1,nlayer
+              IF (l.LE.cstdustlevel) THEN
+                DO ig=1,ngrid 
+                  dsodust(ig,l)=dsodust(ig,l) +
+     &                      aerosol(ig,l,iaer_dust_doubleq) * g /
+     &                      (pplev(ig,l) - pplev(ig,l+1))
+
+                  dsords(ig,l) = dsords(ig,l) +
+     &              aerosol(ig,l,iaer_stormdust_doubleq)* g/
+     &              (pplev(ig,l) - pplev(ig,l+1))
+                ENDDO
+              ELSE
+                DO ig=1,ngrid
+                  dsodust(ig,l) =dsodust(ig,l) +
+     &                      aerosol(ig,l,iaer_dust_doubleq) * g /
+     &                      (pplev(ig,l) - pplev(ig,l+1))
+                   dsords(ig,l) = dsords(ig,l) +
+     &                        aerosol(ig,l,iaer_stormdust_doubleq)* g/
+     &                        (pplev(ig,l) - pplev(ig,l+1))
+                ENDDO
+              ENDIF
+            ENDDO
+          ENDIF
+
+c -----------------------------------------------------------------
+c -----------------------------------------------------------------
+c aerosol/X for stormdust to prepare calculation of radiative transfer
+c -----------------------------------------------------------------
+      if (rdstorm) then
+         DO l=1,nlayer
+           DO ig=1,ngrid
+               aerosol(ig,l,iaer_stormdust_doubleq) = 
+     &           aerosol(ig,l,iaer_stormdust_doubleq)/totstormfract(ig)
+           ENDDO
+         ENDDO
+      endif 
+
 
       END SUBROUTINE aeropacity
Index: trunk/LMDZ.MARS/libf/phymars/calcstormfract_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/calcstormfract_mod.F90	(revision 1974)
+++ trunk/LMDZ.MARS/libf/phymars/calcstormfract_mod.F90	(revision 1974)
@@ -0,0 +1,64 @@
+      MODULE calcstormfract_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+!=======================================================================
+! ROCKET DUST STORM - Defining dust storm fraction of the mesh
+!=======================================================================
+! before calling radiative transfer
+! SUBGRID PARAMETRIZATION 
+! -----------------------------------------------------------------------
+! Author : T. Bertrand 2014-05 
+! Institution : Laboratoire de Meteorologie Dynamique (LMD) Paris, France
+! -----------------------------------------------------------------------
+! defining the fraction of the mesh containing the rocket dust storm.
+! file using the fraction: aeropacity.F, before calling radiative transfer
+! the radiative transfer is performed in the mesh without the storm and in the mesh with the storm
+! the fraction is used to calculate the true opacity of dust within the storm
+
+      SUBROUTINE calcstormfract(ngrid,nlayer,nq,pq,totstormfract)
+                                 
+      use tracer_mod, only: igcm_stormdust_mass
+
+      implicit none
+
+!--------------------------------------------------------
+! Input Variables
+!--------------------------------------------------------
+
+      INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
+      INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
+      INTEGER, INTENT(IN) :: nq ! number of tracer species
+      REAL, INTENT(IN) ::  pq(ngrid,nlayer,nq) ! advected field nq
+
+!--------------------------------------------------------
+! Output Variables
+!--------------------------------------------------------
+
+      REAL, INTENT(OUT) :: totstormfract(ngrid) ! fraction of the mesh containing the dust storm
+
+!--------------------------------------------------------
+! Local variables 
+!--------------------------------------------------------
+
+      REAL, PARAMETER :: mmr_ref=5.e-4 ! mass reference mixing ratio (corresponding
+                                       ! to the tau=10 reference OMEGA local dust storm)
+      REAL, PARAMETER :: fracmax=0.6    ! maximal fraction of the storm
+      REAL, PARAMETER :: fracmin=1.e-2 ! minimal fraction of the storm
+      INTEGER ig     
+ 
+! **********************************************************************
+! Definition of stormfraction
+! **********************************************************************
+
+      DO ig=1, ngrid
+       totstormfract(ig)=abs(maxval(pq(ig,2:nlayer,  & !why from 2nd layer?
+     &                     igcm_stormdust_mass)))/mmr_ref
+       totstormfract(ig)=max(totstormfract(ig),fracmin)  
+       totstormfract(ig)=min(totstormfract(ig),fracmax) 
+      ENDDO
+
+      END SUBROUTINE calcstormfract
+
+      END MODULE calcstormfract_mod
Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 1974)
@@ -13,13 +13,14 @@
      &   ,lifting,freedust,callddevil,scavenging,sedimentation          &
      &   ,activice,water,tifeedback,microphys,supersat,caps,photochem   &
-     &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds,&
-     &   co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying,        &
-     &   satindexco2
+     &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
+     &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
+     &   ,satindexco2,rdstorm
      
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
-     &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear
+     &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear,dustinjection
      
       COMMON/callkeys_r/topdustref,semi,alphan,euveff,                  &
-     &   tke_heat_flux,dustrefir,fixed_euv_value,CLFfixval
+     &   tke_heat_flux,dustrefir,fixed_euv_value,CLFfixval,             &
+     &   coeff_injection
      
       LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
@@ -42,4 +43,5 @@
       real euveff
       real tke_heat_flux
+      real coeff_injection ! dust injection scheme coefficient
       real CLFfixval
 
@@ -53,10 +55,12 @@
       integer solvarmod   ! model for solar EUV variation
       integer solvaryear  ! mars year for realisticly varying solar EUV 
+      integer dustinjection ! dust injection scheme number 
 
       logical rayleigh
       logical tracer
       integer dustbin
-      logical freedust 
+      logical freedust  
       logical active,doubleq,submicron,lifting,callddevil,scavenging
+      logical rdstorm ! rocket dust storm parametrization
       logical sedimentation
       logical activice,tifeedback,supersat,caps
Index: trunk/LMDZ.MARS/libf/phymars/callradite_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 1974)
@@ -7,7 +7,9 @@
       SUBROUTINE callradite(icount,ngrid,nlayer,nq,zday,ls,pq,albedo,
      $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
-     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
-     &     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
-     &     nuice,co2ice,clearsky,totcloudfrac)
+     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
+     $     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
+     $     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,
+     $     totstormfract,clearatm,dsords,
+     $     clearsky,totcloudfrac)
 
       use aeropacity_mod, only: aeropacity
@@ -16,5 +18,6 @@
       use dimradmars_mod, only: naerkind, name_iaer,
      &            iaer_dust_conrath,iaer_dust_doubleq,
-     &            iaer_dust_submicron,iaer_h2o_ice
+     &            iaer_dust_submicron,iaer_h2o_ice,
+     &            iaer_stormdust_doubleq 
       use yomlw_h, only: gcp, nlaylte
       use comcstfi_h, only: g,cpp
@@ -154,5 +157,4 @@
 c                         the "naerkind" kind of aerosol optical
 c                         properties.
-
 c=======================================================================
 c
@@ -160,5 +162,5 @@
 c    -------------
 c
-#include "callkeys.h"
+      include "callkeys.h"
 
 c-----------------------------------------------------------------------
@@ -170,5 +172,5 @@
 
       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
-      REAL,INTENT(IN) :: tauscaling(ngrid) ! Conversion factor for 
+      REAL,INTENT(INOUT) :: tauscaling(ngrid) ! Conversion factor for 
                                ! qdust and Ndust
       REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
@@ -194,7 +196,14 @@
       REAL,INTENT(OUT) :: nuice(ngrid,nlayer)  ! Estimated effective variance
       REAL,INTENT(IN) :: co2ice(ngrid)           ! co2 ice surface layer (kg.m-2)
+
+c     rocket dust storm
+      LOGICAL,INTENT(IN) :: clearatm ! true for background dust
+      REAL,INTENT(IN) :: totstormfract(ngrid) ! dust storm mesh fraction
+      REAL,INTENT(OUT) :: rstormdust(ngrid,nlayer)  ! Storm dust geometric mean radius (m)
+      REAL dsords(ngrid,nlayer) ! density scaled opacity for rocket dust storm dust
+
 c     sub-grid scale water ice clouds
-      real,intent(out) :: totcloudfrac(ngrid)
-      logical,intent(in) :: clearsky
+      LOGICAL,INTENT(IN) :: clearsky
+      REAL,INTENT(IN) :: totcloudfrac(ngrid)
 
 c
@@ -292,4 +301,5 @@
          iaer_dust_submicron=0
          iaer_h2o_ice=0
+         iaer_stormdust_doubleq=0
 
          aer_count=0
@@ -326,4 +336,12 @@
            enddo
          endif
+         if (rdstorm.AND.active) then
+           do iaer=1,naerkind
+             if (name_iaer(iaer).eq."stormdust_doubleq") then
+               iaer_stormdust_doubleq = iaer
+               aer_count = aer_count + 1
+             endif
+           enddo
+         end if
 
 c        Check that we identified all tracers:
@@ -372,5 +390,5 @@
 c     Updating aerosol size distributions:
       CALL updatereffrad(ngrid,nlayer,
-     &                rdust,rice,nuice,
+     &                rdust,rstormdust,rice,nuice,
      &                reffrad,nueffrad,
      &                pq,tauscaling,tau,pplay)
@@ -385,7 +403,8 @@
 c     Computing aerosol optical depth in each layer:
       CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &     pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
-     &     nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
-     &     clearsky,totcloudfrac)
+     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
+     &    QREFvis3d,QREFir3d,omegaREFir3d,
+     &    totstormfract,clearatm,dsords,
+     &    clearsky,totcloudfrac)
 
 c     Starting loop on sub-domain
Index: trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F	(revision 1974)
@@ -4,11 +4,11 @@
 
       CONTAINS
-      
-      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
-     &                pplev,zlev, zlay, pt, pdt, rdust, rice,
+
+      SUBROUTINE callsedim(ngrid,nlay,ptimestep,
+     &                pplev,zlev,zlay,pt,pdt,rdust,rstormdust,rice,
      &                rsedcloud,rhocloud,
-     &                pq, pdqfi, pdqsed,pdqs_sed,nq, 
+     &                pq,pdqfi,pdqsed,pdqs_sed,nq, 
      &                tau,tauscaling)
-! to use  'getin'
+
       USE ioipsl_getincom, only: getin
       USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
@@ -17,6 +17,7 @@
      &                      igcm_ccn_mass, igcm_ccn_number,
      &                      igcm_h2o_ice, nuice_sed, nuice_ref,
-     &                      igcm_ccnco2_mass, igcm_ccnco2_number,
-     &                      igcm_co2_ice
+     &                      igcm_ccnco2_mass,igcm_ccnco2_number,
+     &                      igcm_co2_ice, igcm_stormdust_mass, 
+     &                      igcm_stormdust_number
       USE newsedim_mod, ONLY: newsedim
       USE comcstfi_h, ONLY: g
@@ -61,4 +62,5 @@
 c    Aerosol radius provided by the water ice microphysical scheme:
       real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
+      real,intent(out) :: rstormdust(ngrid,nlay) ! Stormdust geometric mean radius (m) 
       real,intent(out) :: rice(ngrid,nlay)  ! H2O Ice geometric mean radius (m)
 c     Sedimentation radius of water ice
@@ -89,8 +91,7 @@
       real r0dust(ngrid,nlay) ! geometric mean radius used for
                                     !   dust (m)
-!      real r0ccn(ngrid,nlay)  ! geometric mean radius used for
+      real r0stormdust(ngrid,nlay) ! Geometric mean radius used for stormdust (m)
 !                                    !   CCNs (m)
       real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
-
 c     for ice radius computation
       REAL Mo,No
@@ -133,4 +134,8 @@
       INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
                                   !   mix. ratio
+      INTEGER,SAVE :: istormdust_mass  !  index of tracer containing
+                                       !stormdust mass mix. ratio
+      INTEGER,SAVE :: istormdust_number !  index of tracer containing
+                                        !stormdust number mix. ratio                      
       INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
       INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
@@ -140,7 +145,8 @@
       LOGICAL,SAVE :: firstcall=.true.
 
+
+
 c    ** un petit test de coherence
 c       --------------------------
-
       ! AS: firstcall OK absolute
       IF (firstcall) THEN
@@ -240,6 +246,5 @@
        ENDIF                    !of if (co2clouds)
 
-
-        IF (water) THEN
+       IF (water) THEN
          write(*,*) "correction for the shape of the ice particles ?"
          beta=0.75 ! default value
@@ -251,5 +256,30 @@
             write(*,*) "water_param nueff Radiative:", nuice_ref
           ENDIF
-        ENDIF
+       ENDIF
+
+       IF (rdstorm) THEN ! identifying stormdust tracers for sedimentation
+           istormdust_mass=0      ! dummy initialization
+           istormdust_number=0    ! dummy initialization
+
+           do iq=1,nq
+             if (noms(iq).eq."stormdust_mass") then
+               istormdust_mass=iq
+               write(*,*)"callsedim: istormdust_mass=",istormdust_mass
+             endif
+             if (noms(iq).eq."stormdust_number") then
+               istormdust_number=iq
+               write(*,*)"callsedim: istormdust_number=",
+     &                                           istormdust_number
+             endif
+           enddo
+
+           ! check that we did find the tracers
+           if ((istormdust_mass.eq.0).or.(istormdust_number.eq.0)) then
+             write(*,*) 'callsedim: error! could not identify'
+             write(*,*) ' tracers for stormdust mass and number mixing'
+             write(*,*) ' ratio and rdstorm is activated!'
+             stop
+           endif
+       ENDIF !of if (rdstorm)
       
         firstcall=.false.
@@ -268,5 +298,4 @@
       zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
      &                         +pdt(1:ngrid,1:nlay)*ptimestep
-
 
 c    Computing the different layer properties
@@ -295,6 +324,16 @@
         end do
       endif
-
-
+      ! rocket dust storm
+      if (rdstorm) then
+        do l=1,nlay
+          do ig=1, ngrid
+     
+         call updaterdust(zqi(ig,l,igcm_stormdust_mass),
+     &               zqi(ig,l,igcm_stormdust_number),r0stormdust(ig,l),
+     &               tauscaling(ig))
+          
+          end do
+        end do
+      endif
 c =================================================================
       do iq=1,nq
@@ -307,11 +346,11 @@
 c -----------------------------------------------------------------
           if ((doubleq.and.
-     &        ((iq.eq.idust_mass).or.
-     &         (iq.eq.idust_number)))) then
+     &     ((iq.eq.idust_mass).or.(iq.eq.idust_number).or.
+     &     (iq.eq.istormdust_mass).or.(iq.eq.istormdust_number)))) then
      
 c           Computing size distribution:
 c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-c            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
+            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
               do  l=1,nlay
                 do ig=1, ngrid
@@ -319,9 +358,17 @@
                 end do
               end do
-              sigma0 = varian
+            else if ((iq.eq.istormdust_mass).or.
+     &                                (iq.eq.istormdust_number)) then
+              do  l=1,nlay
+                do ig=1, ngrid
+                  r0(ig,l)=r0stormdust(ig,l)
+                end do
+              end do
+            endif
+            sigma0 = varian
 
 c        Computing mass mixing ratio for each particle size
 c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
+          IF ((iq.EQ.idust_mass).or.(iq.EQ.istormdust_mass)) then
             radpower = 2
           ELSE  ! number
@@ -376,5 +423,4 @@
 
           do ir=1,nr
-          
                call newsedim(ngrid,nlay,1,1,ptimestep,
      &         pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
@@ -391,5 +437,5 @@
                  zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
                ENDDO
-             ENDDO
+             ENDDO            
           enddo ! of do ir=1,nr
 c -----------------------------------------------------------------
@@ -457,5 +503,15 @@
        ENDDO
       endif ! of if (doubleq)
-      
+
+      if (rdstorm) then
+       DO l = 1, nlay
+        DO ig=1,ngrid
+         call updaterdust(zqi(ig,l,igcm_stormdust_mass),
+     &                zqi(ig,l,igcm_stormdust_number),rstormdust(ig,l),
+     &                tauscaling(ig))    
+        ENDDO
+       ENDDO
+      endif ! of if (rdstorm)
+ 
 c     Update the ice particle size "rice"
 c     -------------------------------------
@@ -490,2 +546,3 @@
       
       END MODULE callsedim_mod
+
Index: trunk/LMDZ.MARS/libf/phymars/compute_dtau_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/compute_dtau_mod.F90	(revision 1974)
+++ trunk/LMDZ.MARS/libf/phymars/compute_dtau_mod.F90	(revision 1974)
@@ -0,0 +1,83 @@
+	MODULE compute_dtau_mod
+
+	IMPLICIT NONE
+
+	REAL,PARAMETER :: ti_injection=10/24. ! time of beginning injection
+	REAL,PARAMETER :: tf_injection=18/24. ! time of end injection
+        REAL,PARAMETER :: t_scenario=14/24.   ! time of tauref_scenario
+
+	CONTAINS
+
+	SUBROUTINE compute_dtau(ngrid,nlayer,            		&
+       & 			zday,pplev,tauref,			&
+       & 			ptimestep,dustliftday,local_time)
+
+	USE geometry_mod, only: longitude_deg
+        USE time_phylmdz_mod, only: dtphys, daysec
+        USE comcstfi_h, only: g
+        USE tracer_mod, only: alpha_lift,igcm_dust_mass,igcm_dust_number
+	
+	IMPLICIT NONE
+        
+        include "callkeys.h"
+	
+        INTEGER, INTENT(in) :: ngrid
+	INTEGER, INTENT(in) :: nlayer
+        REAL, INTENT(in) :: zday ! date at lon=0, in fraction of sols
+	REAL, INTENT(in) :: pplev(ngrid,nlayer+1) ! pressure (Pa)
+	REAL, INTENT(in) :: tauref(ngrid) ! Computed dust opacity at 610Pa
+	REAL, INTENT(in) :: ptimestep 
+        REAL, INTENT(in) :: local_time(ngrid)
+        REAL, INTENT(out) :: dustliftday(ngrid) ! Dust injection rate (s-1)
+        
+        INTEGER :: ig, l
+	INTEGER, SAVE :: nb_daystep ! nomber of step a day
+        REAL :: tauref_scenario(ngrid) ! from dust scenario
+        REAL :: zday_scenario
+        REAL,ALLOCATABLE,SAVE :: local_time_prev(:)
+        REAL,ALLOCATABLE,SAVE :: dtau(:) ! Dust opacity increment (at 610Pa)
+        REAL,PARAMETER :: odpref=610. !DOD reference pressure (Pa)
+        
+        LOGICAL, SAVE :: firstcall=.TRUE. ! signals first call to physics
+        
+        
+        IF(firstcall)THEN
+                ALLOCATE(local_time_prev(ngrid))
+                ALLOCATE(dtau(ngrid))
+                DO ig=1,ngrid
+                   local_time_prev(ig)=modulo(1.+(zday-ptimestep/daysec)&
+                                      -INT(zday-ptimestep/daysec) 	&
+                                      +(longitude_deg(ig)/15)/24,1.)
+                   dtau(ig)=0.
+                ENDDO
+                nb_daystep=(daysec/dtphys)
+                firstcall=.FALSE.
+        ENDIF
+        
+	! 1. Obtain tauref_scenario from dust scenario
+	zday_scenario=zday-modulo(zday,1.)
+        call read_dust_scenario(ngrid,nlayer,zday_scenario,pplev,	&
+       & 			tauref_scenario)
+
+	! 2. Compute dtau
+	DO ig=1,ngrid
+         IF ((local_time(ig).ge.t_scenario).and. 			&
+       &          (local_time_prev(ig).lt.(t_scenario)))THEN
+	         dtau(ig)=tauref_scenario(ig)-tauref(ig)
+         ENDIF
+         IF (dtau(ig).LT.0) THEN
+                  dtau(ig)=0.
+         ENDIF
+         
+        ! 3. Use dtau to compute dustliftday
+         dustliftday(ig)=coeff_injection*				&
+       &		(dtau(ig)*pplev(ig,1)/odpref)			&
+       &		/(daysec*(tf_injection-ti_injection))
+        ENDDO 
+         
+        ! 4. Save local time 
+        local_time_prev(1:ngrid)=local_time(1:ngrid)
+        
+        end subroutine compute_dtau
+
+	end module compute_dtau_mod
Index: trunk/LMDZ.MARS/libf/phymars/comsaison_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/comsaison_h.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/comsaison_h.F90	(revision 1974)
@@ -9,4 +9,5 @@
   real,save,allocatable :: mu0(:)
   real,save,allocatable :: fract(:)
+  real,save,allocatable :: local_time(:) ! local solar time as fraction of day (0,1)
 
 contains
@@ -19,4 +20,5 @@
     allocate(mu0(ngrid))
     allocate(fract(ngrid))
+    allocate(local_time(ngrid))
   
   end subroutine ini_comsaison_h
@@ -29,4 +31,5 @@
     if (allocated(mu0)) deallocate(mu0)
     if (allocated(fract)) deallocate(fract)
+    if (allocated(local_time)) deallocate(local_time)
 
   end subroutine end_comsaison_h
Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1974)
@@ -280,4 +280,10 @@
          write(*,*)" calllott = ",calllott
 
+! rocket dust storm injection scheme
+         write(*,*)"call rocket dust storm and slope lifting", 
+     &              " parametrization"
+         rdstorm=.false. ! default value
+         call getin("rdstorm",rdstorm)
+         write(*,*)" rdstorm = ",rdstorm
 
          write(*,*)"rad.transfer is computed every iradia",
@@ -359,4 +365,13 @@
          endif
 
+! dust injection scheme
+        dustinjection=0 ! default: no injection scheme
+        call getin("dustinjection",dustinjection)
+        write(*,*)" dustinjection = ",dustinjection
+! dust injection scheme coefficient        
+        coeff_injection=1. ! default value
+        call getin("coeff_injection",coeff_injection)
+        write(*,*)" coeff_in,jection = ",coeff_injection
+
 ! free evolving dust
 ! freedust=true just says that there is no lifting and no dust opacity scaling.
@@ -372,11 +387,46 @@
          ! this test is valid in GCM case
          ! ... not in mesoscale case, for we want to activate mesoscale lifting
-         if (freedust.and.lifting) then
-           print*,'if freedust is used, then lifting should not be used'
-           print*,'lifting forced to false !!'
-           lifting=.false.
+         if (freedust.and.dustinjection.eq.0)then
+           if(lifting) then
+             print*,'if freedust is used and dustinjection = 0, 
+     &      then lifting should not be used'
+             stop
+           endif
          endif
 #endif
-
+         if (dustinjection.eq.1)then
+           if(.not.lifting) then
+             print*,"if dustinjection=1, then lifting should be true"
+             stop
+           endif
+           if(.not.freedust) then
+             print*,"if dustinjection=1, then freedust should be true"
+             stop
+           endif
+         endif
+! rocket dust storm
+! Test of incompatibility:
+! if rdstorm is used, then doubleq should be true
+         if (rdstorm.and..not.doubleq) then
+           print*,'if rdstorm is used, then doubleq should be used !'
+           stop
+         endif
+         if (rdstorm.and..not.active) then
+           print*,'if rdstorm is used, then active should be used !'
+           stop
+         endif
+         if (rdstorm.and..not.lifting) then
+           print*,'if rdstorm is used, then lifting should be used !'
+           stop
+         endif
+         if (rdstorm.and..not.freedust) then
+           print*,'if rdstorm is used, then freedust should be used !'
+           stop
+         endif
+         if (rdstorm.and.(dustinjection.eq.0)) then
+           print*,'if rdstorm is used, then dustinjection should 
+     &             be used !'
+           stop
+         endif
 ! Dust IR opacity
          write(*,*)" Wavelength for infrared opacity of dust ?"
@@ -642,9 +692,4 @@
            WRITE(*,*) 'equal to 2.'
            CALL ABORT
-         ELSE IF ( (.NOT.activice).AND.(naerkind.GT.1) ) THEN
-           WRITE(*,*) 'naerkind is greater than unity, but'
-           WRITE(*,*) 'activice has not been set to .true.'
-           WRITE(*,*) 'in callphys.def; this is not logical!'
-           CALL ABORT
          ENDIF
 
@@ -666,7 +711,12 @@
           ! and picky compilers who know name_iaer(2) is out of bounds
           j=2
+         IF (rdstorm.AND..NOT.activice) name_iaer(2) = 
+     &                       "stormdust_doubleq" !! storm dust two-moment scheme
+         IF (rdstorm.AND.water.AND.activice) name_iaer(3) = 
+     &                                             "stormdust_doubleq" 
          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 (nq.gt.1)
+
 c        ----------------------------------------------------------
 
Index: trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90	(revision 1974)
@@ -32,4 +32,7 @@
                               ! submicron population of dust
                               ! particles
+  integer iaer_stormdust_doubleq ! Storm dust profile is given by the
+                              ! mass mixing ratio of the two moment scheme 
+                              ! method (doubleq)
   integer iaer_h2o_ice ! Water ice particles
 
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1974)
@@ -11,5 +11,5 @@
      &                     albedice, iceradius, dtemisice, z0,
      &                     zmea, zstd, zsig, zgam, zthe, phisfi,
-     &                     watercaptag
+     &                     watercaptag, hmons
       use slope_mod, only: theta_sl, psi_sl
       use phyredem, only: physdem0,physdem1
@@ -557,4 +557,6 @@
       zgam(1)=0.E+0
       zthe(1)=0.E+0
+c  for rocket dust storm scheme      
+      hmons(1)=0.E+0
 
 
@@ -716,5 +718,5 @@
       call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
      &              nq,dtphys,float(day0),time,cell_area,
-     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,hmons)
       call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
      &              dtphys,time,
Index: trunk/LMDZ.MARS/libf/phymars/initracer.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/initracer.F	(revision 1974)
@@ -24,5 +24,5 @@
 
 
-#include "callkeys.h"
+      include "callkeys.h"
 
       integer,intent(in) :: ngrid ! number of atmospheric columns
@@ -32,9 +32,10 @@
       integer iq,ig,count
       real r0_lift , reff_lift, nueff_lift
+      real r0_storm,reff_storm
 c     Ratio of small over large dust particles (used when both 
 c       doubleq and the submicron mode are active); In Montmessin
 c       et al. (2002), a value of 25 has been deduced;
       real, parameter :: popratio = 25.
-      character(len=20) :: txt ! to store some text
+      character(len=30) :: txt ! to store some text
 
 c-----------------------------------------------------------------------
@@ -69,4 +70,6 @@
       igcm_h2o_vap=0
       igcm_h2o_ice=0
+      igcm_stormdust_mass=0
+      igcm_stormdust_number=0
       igcm_co2=0
       igcm_co=0
@@ -147,4 +150,16 @@
         enddo
       endif ! of if (submicron)
+       if (rdstorm) then
+        do iq=1,nq
+          if (noms(iq).eq."stormdust_mass") then
+            igcm_stormdust_mass=iq
+            count=count+1
+          endif
+          if (noms(iq).eq."stormdust_number") then
+            igcm_stormdust_number=iq
+            count=count+1
+          endif
+        enddo
+      endif ! of if (rdstorm)     
       ! 2. find chemistry and water tracers
       do iq=1,nq
@@ -422,4 +437,9 @@
 #else
         alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
+        IF (dustinjection.ge.1) THEN
+        	reff_lift = 3.0e-6 ! Effective radius of lifted dust (m)
+        	alpha_lift(igcm_dust_mass)=(4/3.)*reff_lift*rho_dust
+     &           				/2.4
+     	ENDIF
 #endif
 
@@ -437,4 +457,29 @@
         write(*,*) "initracer: doubleq_param alpha_lift:",
      &    alpha_lift(igcm_dust_mass)
+!c ----------------------------------------------------------------------
+!c rocket dust storm scheme
+!c lifting tracer stormdust using same distribution than
+!c normal dust 
+        if (rdstorm) then
+          reff_storm=3.e-6 ! reff_lift !3.e-6
+          r0_storm=reff_storm/ref_r0
+          rho_q(igcm_stormdust_mass)=rho_dust
+          rho_q(igcm_stormdust_number)=rho_dust
+
+          alpha_devil(igcm_stormdust_mass)=9.e-9   ! dust devil lift mass coeff
+          alpha_lift(igcm_stormdust_mass)=4./3./2.4*reff_storm*rho_dust
+
+          write(*,*) 'alpha_lift(rds):',alpha_lift(igcm_stormdust_mass)
+  
+          alpha_devil(igcm_stormdust_number)=r3n_q*
+     &                      alpha_devil(igcm_stormdust_mass)/r0_storm**3
+          alpha_lift(igcm_stormdust_number)=r3n_q*
+     &                       alpha_lift(igcm_stormdust_mass)/r0_storm**3
+  
+          radius(igcm_stormdust_mass) = reff_storm
+          radius(igcm_stormdust_number) = reff_storm
+        end if !(rdstorm)
+!c ----------------------------------------------------------------------
+      
       else
 
@@ -618,5 +663,21 @@
           endif 
        endif
-       
+ 
+       if (rdstorm) then 
+       ! verify that we indeed have stormdust_mass and stormdust_number tracers 
+         if (igcm_stormdust_mass.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use rdstorm option without ",
+     &                "a stormdust_mass tracer !"
+           stop
+         endif
+         if (igcm_stormdust_number.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use rdstorm option without ",
+     &                "a stormdust_number tracer !"
+           stop
+         endif
+       endif
+      
        if (callnlte) then ! NLTE requirements
          if (nltemodel.ge.1) then
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 1974)
@@ -12,5 +12,5 @@
   use tracer_mod, only: noms ! tracer names
   use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
-                       zmea, zstd, zsig, zgam, zthe
+                       zmea, zstd, zsig, zgam, zthe, hmons
   use iostart, only: nid_start, open_startphy, close_startphy, &
                      get_field, get_var, inquire_field, &
@@ -19,4 +19,6 @@
 
   implicit none
+  
+  include "callkeys.h"
 !======================================================================
 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
@@ -179,4 +181,22 @@
 endif
 
+! hmons 
+call get_field("hmons",hmons,found)
+if (.not.found) then
+  write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
+  if (rdstorm) then
+    call abort
+  else
+    write(*,*) "will continue anyway..."
+    write(*,*) "because you may not need it."
+    hmons(:)=0.
+  end if
+else
+  do ig=1,ngrid
+    if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
+  enddo
+  write(*,*) "phyetat0: <hmons> range:", &
+             minval(hmons), maxval(hmons)
+endif
       
 ! Time axis
@@ -277,5 +297,5 @@
 if (.not.found) then
   write(*,*) "phyetat0: <tauscaling> not in file"
-  tauscaling(:) = -1
+  tauscaling(:) = 1
 else
   write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
@@ -296,5 +316,5 @@
 call get_field("wstar",wstar,found,indextime)
 if (.not.found) then
-  write(*,*) "phyetat0: <totcloudfrac> not in file! Set to zero"
+  write(*,*) "phyetat0: <wstar> not in file! Set to zero"
   wstar(:)=0
 else
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1974)
@@ -1,5 +1,2 @@
-!
-! $Id: $
-!
 module phyredem
 
@@ -10,5 +7,6 @@
 subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
                          phystep,day_ini,time,airefi, &
-                         alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
+                         alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe, &
+                         phmons)
 ! create physics restart file and write time-independent variables
   use comsoil_h, only: inertiedat, volcapa, mlayer
@@ -16,5 +14,6 @@
   use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, &
                        z0_default, albedice, emisice, emissiv, &
-                       iceradius, dtemisice, phisfi, z0
+                       iceradius, dtemisice, phisfi, z0,   &
+                       hmons
   use dimradmars_mod, only: tauvis
   use iostart, only : open_restartphy, close_restartphy, & 
@@ -46,5 +45,6 @@
   real,intent(in) :: pzgam(ngrid)
   real,intent(in) :: pzthe(ngrid)
- 
+  real,intent(in) :: phmons(ngrid)
+
   real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
   
@@ -131,5 +131,9 @@
   call put_field("ZGAM","Relief: gamma parameter",zgam)
   call put_field("ZTHE","Relief: theta parameter",zthe)
-  
+  ! Write hmons if it was read in startfi
+  if (maxval(hmons).ne.minval(hmons)) then
+    call put_field("hmons","Relief: hmons parameter",hmons)
+  end if
+     
   ! Soil thermal inertia
   call put_field("inertiedat","Soil thermal inertia",ith)
Index: trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 1974)
@@ -47,4 +47,6 @@
       use time_phylmdz_mod, only: init_time
       use co2cloud_mod, only: ini_co2cloud,end_co2cloud
+      use rocketduststorm_mod, only: ini_rocketduststorm_mod, &
+                                     end_rocketduststorm_mod
 
       IMPLICIT NONE
@@ -112,4 +114,8 @@
       call ini_co2cloud(ngrid,nlayer)
       
+      ! allocate arrays in "rocketduststorm_mod":
+      call end_rocketduststorm_mod
+      call ini_rocketduststorm_mod(ngrid)
+      
       END SUBROUTINE phys_state_var_init
 
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1974)
@@ -17,7 +17,8 @@
       use co2cloud_mod, only: co2cloud, mem_Mccn_co2, mem_Mh2o_co2,
      &                        mem_Nccn_co2
-      use aeropacity_mod
-      use callradite_mod
-      use callsedim_mod
+      use callradite_mod, only: callradite
+      use callsedim_mod, only: callsedim
+      use rocketduststorm_mod, only: rocketduststorm, dustliftday
+      use calcstormfract_mod, only: calcstormfract
       use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
      &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
@@ -27,8 +28,10 @@
      &                      igcm_dust_mass, igcm_dust_number, igcm_h2o2,
      &                      nuice_ref, rho_ice, rho_dust, ref_r0,
-     &                      igcm_he
+     &                      igcm_he, igcm_stormdust_mass,
+     &                      igcm_stormdust_number 
       use comsoil_h, only: inertiedat, ! soil thermal inertia
      &                     tsoil, nsoilmx ! number of subsurface layers
-      use geometry_mod, only: longitude, latitude, cell_area
+      use geometry_mod, only: longitude, latitude, cell_area,
+     &                        longitude_deg 
       use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
       use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
@@ -36,6 +39,7 @@
      &                     frost_albedo_threshold,
      &                     tsurf, co2ice, emis,
-     &                     capcal, fluxgrd, qsurf
-      use comsaison_h, only: dist_sol, declin, mu0, fract
+     &                     capcal, fluxgrd, qsurf,
+     &                     hmons
+      use comsaison_h, only: dist_sol, declin, mu0, fract, local_time
       use slope_mod, only: theta_sl, psi_sl
       use conc_mod, only: rnew, cpnew, mmean
@@ -48,5 +52,5 @@
       use planete_h, only: aphelie, periheli, year_day, peri_day,
      &                     obliquit
-      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
+      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad 
       USE calldrag_noro_mod, ONLY: calldrag_noro
       USE vdifc_mod, ONLY: vdifc
@@ -54,4 +58,5 @@
      &                      fill_data_thermos, allocate_param_thermos
       use iono_h, only: allocate_param_iono
+      use compute_dtau_mod, only: compute_dtau
 #ifdef MESOSCALE
       use comsoil_h, only: mlayer,layer
@@ -64,4 +69,5 @@
       use eofdump_mod, only: eofdump
       USE vertical_layers_mod, ONLY: ap,bp,aps,bps
+     
 #endif
 
@@ -207,4 +213,5 @@
       REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
 
+      
 
 c Local saved variables:
@@ -212,7 +219,9 @@
       INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
       INTEGER,SAVE :: icount     ! counter of calls to physiq during the run.
+
 #ifdef DUSTSTORM
       REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
 #endif
+            
 c     Variables used by the water ice microphysical scheme:
       REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
@@ -257,8 +266,14 @@
       REAL fluxtop_sw(ngrid,2)     !Outgoing SW (solar) flux to space (W.m-2)
       REAL tauref(ngrid)           ! Reference column optical depth at odpref
+c     rocket dust storm
+      REAL totstormfract(ngrid)     ! fraction of the mesh where the dust storm is contained
+      logical clearatm              ! clearatm used to calculate twice the radiative 
+                                    ! transfer when rdstorm is active : 
+                                    !            - in a mesh with stormdust and background dust (false)
+                                    !            - in a mesh with background dust only (true)
+      
       real,parameter :: odpref=610. ! DOD reference pressure (Pa)
       REAL tau(ngrid,naerkind)     ! Column dust optical depth at each point
                                    ! AS: TBD: this one should be in a module !
-      REAL dsodust(ngrid,nlayer)   ! density-scaled opacity (in infrared) 
       REAL zls                       !  solar longitude (rad)
       REAL zday                      ! date (time since Ls=0, in martian days)
@@ -271,6 +286,6 @@
       REAL zdtlw(ngrid,nlayer)     ! (K/s)
       REAL zdtsw(ngrid,nlayer)     ! (K/s)
-!      REAL cldtlw(ngrid,nlayer)     ! (K/s) LW heating rate for clear area
-!      REAL cldtsw(ngrid,nlayer)     ! (K/s) SW heating rate for clear area
+      REAL pdqrds(ngrid,nlayer,nq) ! tendency for dust after rocketduststorm
+
       REAL zdtnirco2(ngrid,nlayer) ! (K/s)
       REAL zdtnlte(ngrid,nlayer)   ! (K/s)
@@ -315,4 +330,9 @@
       PARAMETER (length=100)
 
+c     Variables for the total dust for diagnostics
+      REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust
+
+      INTEGER iaer
+
 c local variables only used for diagnostic (output in file "diagfi" or "stats")
 c -----------------------------------------------------------------------------
@@ -320,4 +340,5 @@
       REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
       REAL zq(ngrid,nlayer,nq)
+
       REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
       character*2 str2
@@ -325,4 +346,5 @@
       real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
       real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
+      real rstormdust(ngrid,nlayer) ! stormdust geometric mean radius (m)
       integer igmin, lmin
       logical tdiag
@@ -340,4 +362,6 @@
       real mass(nq)             ! global mass of tracers (g)
       REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
+      REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
+      REAL mdusttot(ngrid)    ! Total mass of dust tracer (kg/m2)
       REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
       REAL mtotco2(ngrid)       ! Total mass of co2 vapor (kg/m2)
@@ -360,4 +384,13 @@
       REAL nccn(ngrid,nlayer)  ! true n ccn (kg/kg)
       REAL qccn(ngrid,nlayer)  ! true q ccn (kg/kg)
+c     definition tendancies of stormdust tracers
+      REAL rdsdqdustsurf(ngrid) ! surface q stormdust flux (kg/m2/s)
+      REAL rdsdndustsurf(ngrid) ! surface n stormdust flux (number/m2/s)
+      REAL rdsndust(ngrid,nlayer) ! true n stormdust (kg/kg)
+      REAL rdsqdust(ngrid,nlayer) ! true q stormdust (kg/kg)
+      REAL wspeed(ngrid,nlayer+1) ! vertical speed tracer stormdust
+      REAL dsodust(ngrid,nlayer) 
+      REAL dsords(ngrid,nlayer)
+
       REAL nccnco2(ngrid,nlayer)   ! true n ccnco2 (kg/kg)
       REAL qccnco2(ngrid,nlayer)  ! true q ccnco2 (kg/kg)
@@ -396,5 +429,7 @@
       REAL tstar(ngrid)  ! friction velocity and friction potential temp
       REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
-!      REAL zu2(ngrid)
+      real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer)
+      real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer)   
+      real qdusttotal0(ngrid),qdusttotal1(ngrid)
 
 c sub-grid scale water ice clouds (A. Pottier 2013)
@@ -543,5 +578,6 @@
      &                 nsoilmx,ngrid,nlayer,nq,
      &                 ptimestep,pday,time_phys,cell_area,
-     &                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+     &                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
+     &                 hmons)
          endif
 #endif
@@ -565,4 +601,5 @@
       zdtsurf(:)=0
       dqsurf(:,:)=0
+      
 #ifdef DUSTSTORM
       pq_tmp(:,:,:)=0
@@ -572,5 +609,9 @@
 
       zday=pday+ptime ! compute time, in sols (and fraction thereof)
-    
+      ! Compute local time at each grid point
+      DO ig=1,ngrid
+       local_time(ig)=modulo(1.+(zday-INT(zday))
+     &                +(longitude_deg(ig)/15)/24,1.)
+      ENDDO
 c     Compute Solar Longitude (Ls) :
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -702,11 +743,17 @@
      &                             zdtnlte(1:ngrid,1:nlayer)/86400.
               endif
-           else
+           ELSE
               zdtnlte(:,:)=0.
-           endif
+           ENDIF !end callnlte
 
 c          Find number of layers for LTE radiation calculations
            IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
      &          CALL nlthermeq(ngrid,nlayer,zplev,zplay)
+
+c          rocketstorm : compute dust storm mesh fraction
+           IF (rdstorm) THEN
+              CALL calcstormfract(ngrid,nlayer,nq,pq,
+     &                 totstormfract)
+           ENDIF
 
 c          Note: Dustopacity.F has been transferred to callradite.F
@@ -724,5 +771,6 @@
 c          Transfer through CO2 (except NIR CO2 absorption)
 c            and aerosols (dust and water ice)
-
+           ! callradite for background dust
+           clearatm=.true.
 c          Radiative transfer
 c          ------------------
@@ -730,8 +778,11 @@
            clearsky=.false. ! part with clouds for both cases CLFvarying true/false 
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
-     $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
-     $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
-     $     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
-     $     nuice,co2ice,clearsky,totcloudfrac)
+     &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
+     &     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
+     &     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
+     &     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,
+     &     totstormfract,clearatm,dsords,
+     &     clearsky,totcloudfrac)
+
            ! case of sub-grid water ice clouds: callradite for the clear case
             IF (CLFvarying) THEN
@@ -746,5 +797,6 @@
      &              fluxsurf_swclf,fluxtop_lwclf,fluxtop_swclf,tauref,
      &              tau,aerosol,dsodust,tauscaling,taucloudtesclf,rdust,
-     &              rice,nuice,co2ice, clearsky, totcloudfrac)
+     &              rice,nuice,co2ice,rstormdust,totstormfract,
+     &              clearatm,dsords,clearsky,totcloudfrac)
                clearsky = .false.  ! just in case.
                ! Sum the fluxes and heating rates from cloudy/clear
@@ -774,10 +826,18 @@
 
             ENDIF ! (CLFvarying)
-
+            
+           ! Dustinjection
+           if (dustinjection.gt.0) then
+             CALL compute_dtau(ngrid,nlayer,
+     &        			zday,pplev,tauref,
+     &        			ptimestep,dustliftday,local_time)
+           endif
+c============================================================================
+           
 #ifdef DUSTSTORM
 !! specific case: compute the added quantity of dust for perturbation
        if (firstcall) then
          pdq(1:ngrid,1:nlayer,igcm_dust_mass)= 
-     &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)
+     &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)     
      &      - pq_tmp(1:ngrid,1:nlayer,1)
      &      + pq(1:ngrid,1:nlayer,igcm_dust_mass)
@@ -897,6 +957,86 @@
       ENDIF ! of IF (callrad)
 
+c     3. Rocket dust storm
+c    -------------------------------------------
+      IF (rdstorm) THEN
+         clearatm=.false.
+         pdqrds(:,:,:)=0.
+         qdusttotal0(:)=0.
+         qdusttotal1(:)=0.
+         do ig=1,ngrid
+           do l=1,nlayer
+            zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential
+                                             ! temperature tendency
+            ! for diagnostics
+            qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+
+     &           pdq(ig,l,igcm_dust_mass)*ptimestep
+            qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+
+     &           pdq(ig,l,igcm_stormdust_mass)*ptimestep
+            qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+
+     &          qstormrds0(ig,l))*(zplev(ig,l)-
+     &           zplev(ig,l+1))/g
+           enddo
+         enddo
+         call writediagfi(ngrid,'qdustrds0','qdust before rds',
+     &           'kg/kg ',3,qdustrds0)
+         call writediagfi(ngrid,'qstormrds0','qstorm before rds',
+     &           'kg/kg ',3,qstormrds0)
+
+         CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,
+     &                      pq,pdq,pt,pdt,zplev,zplay,zzlev,
+     &                      zzlay,zdtsw,zdtlw,
+c               for radiative transfer
+     &                      clearatm,icount,zday,zls,
+     &                      tsurf,igout,totstormfract,
+c               input sub-grid scale cloud
+     &                      clearsky,totcloudfrac,
+c               output
+     &                      pdqrds,wspeed,dsodust,dsords)
+                   
+c      update the tendencies of both dust after vertical transport
+         DO l=1,nlayer
+          DO ig=1,ngrid
+           pdq(ig,l,igcm_stormdust_mass)=
+     &              pdq(ig,l,igcm_stormdust_mass)+ 
+     &                      pdqrds(ig,l,igcm_stormdust_mass)
+           pdq(ig,l,igcm_stormdust_number)=
+     &              pdq(ig,l,igcm_stormdust_number)+
+     &                      pdqrds(ig,l,igcm_stormdust_number)
+
+           pdq(ig,l,igcm_dust_mass)=
+     &       pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
+           pdq(ig,l,igcm_dust_number)=
+     &           pdq(ig,l,igcm_dust_number)+ 
+     &                  pdqrds(ig,l,igcm_dust_number)
+
+          ENDDO
+         ENDDO
+         do l=1,nlayer
+           do ig=1,ngrid
+             qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+
+     &           pdq(ig,l,igcm_dust_mass)*ptimestep
+             qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+
+     &           pdq(ig,l,igcm_stormdust_mass)*ptimestep
+             qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+
+     &          qstormrds1(ig,l))*(zplev(ig,l)-
+     &           zplev(ig,l+1))/g
+           enddo
+         enddo
+
+c        for diagnostics
+         call writediagfi(ngrid,'qdustrds1','qdust after rds',
+     &           'kg/kg ',3,qdustrds1)
+         call writediagfi(ngrid,'qstormrds1','qstorm after rds',
+     &           'kg/kg ',3,qstormrds1)
+         
+         call writediagfi(ngrid,'qdusttotal0','q sum before rds',
+     &           'kg/kg ',2,qdusttotal0)
+         call writediagfi(ngrid,'qdusttotal1','q sum after rds',
+     &           'kg/kg ',2,qdusttotal1)
+
+      ENDIF ! end of if(rdstorm)
+
 c-----------------------------------------------------------------------
-c    3. Gravity wave and subgrid scale topography drag :
+c    4. Gravity wave and subgrid scale topography drag :
 c    -------------------------------------------------
 
@@ -917,5 +1057,5 @@
 
 c-----------------------------------------------------------------------
-c    4. Vertical diffusion (turbulent mixing):
+c    5. Vertical diffusion (turbulent mixing):
 c    -----------------------------------------
 
@@ -968,6 +1108,6 @@
      $        zdum1,zdum2,zdh,pdq,zflubid,
      $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
-     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux)
-
+     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux,
+     &	      dustliftday,local_time)
 
           DO ig=1,ngrid
@@ -1018,8 +1158,8 @@
              zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
              DO iq=1, nq
-              IF ((iq .ne. igcm_dust_mass)
-     &        .and. (iq .ne. igcm_dust_number)) THEN
-                zdqsdif(:,iq)=0.
-              ENDIF
+               IF ((iq .ne. igcm_dust_mass)
+     &          .and. (iq .ne. igcm_dust_number)) THEN
+                 zdqsdif(:,iq)=0.
+               ENDIF
              ENDDO
            ELSE
@@ -1041,5 +1181,5 @@
 
 c-----------------------------------------------------------------------
-c   5. Thermals :
+c   6. Thermals :
 c   -----------------------------
 
@@ -1085,7 +1225,6 @@
         lmax_th_out(:)=0.
       end if
-
 c-----------------------------------------------------------------------
-c   5. Dry convective adjustment:
+c   7. Dry convective adjustment:
 c   -----------------------------
 
@@ -1107,5 +1246,4 @@
      $                zduadj,zdvadj,zdhadj,
      $                zdqadj)
-
 
          DO l=1,nlayer
@@ -1133,10 +1271,10 @@
 
 c-----------------------------------------------------------------------
-c   6. Specific parameterizations for tracers 
+c   8. Specific parameterizations for tracers 
 c:   -----------------------------------------
 
       if (tracer) then 
 
-c   6a. Water and ice
+c   8a. Water and ice
 c     ---------------
 
@@ -1219,5 +1357,5 @@
          END IF  ! of IF (water)
 
-c   6a bis. CO2 clouds (CL & JA)
+c   8a bis. CO2 clouds (CL & JA)
 c        ---------------------------------------
 c        CO2 ice cloud condensation in the atmosphere
@@ -1342,7 +1480,6 @@
 
 
-c   6b. Aerosol particles
+c   8b. Aerosol particles
 c     -------------------
-
 c        ----------
 c        Dust devil :
@@ -1371,5 +1508,5 @@
 c        ------------- 
 c        Sedimentation :   acts also on water ice
-c        ------------- 
+c        -------------
          IF (sedimentation) THEN 
            zdqsed(1:ngrid,1:nlayer,1:nq)=0
@@ -1378,8 +1515,8 @@
 c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
 c Zdqssed isn't
-           call callsedim(ngrid,nlayer, ptimestep,
-     &                zplev,zzlev, zzlay, pt, pdt, rdust, rice,
-     &                rsedcloud,rhocloud,
-     &                pq, pdq, zdqsed, zdqssed,nq, 
+           call callsedim(ngrid,nlayer,ptimestep,
+     &                zplev,zzlev,zzlay,pt,pdt,rdust,rstormdust,
+     &                rice,rsedcloud,rhocloud,
+     &                pq,pdq,zdqsed,zdqssed,nq, 
      &                tau,tauscaling)
 c Flux at the surface of co2 ice computed in co2cloud microtimestep
@@ -1387,4 +1524,20 @@
               zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
            ENDIF
+
+           IF (rdstorm) THEN
+c             Storm dust cannot sediment to the surface
+              DO ig=1,ngrid  
+                 zdqsed(ig,1,igcm_stormdust_mass)=
+     &             zdqsed(ig,1,igcm_stormdust_mass)+
+     &             zdqssed(ig,igcm_stormdust_mass) /
+     &             ((pplev(ig,1)-pplev(ig,2))/g)
+                 zdqsed(ig,1,igcm_stormdust_number)=
+     &             zdqsed(ig,1,igcm_stormdust_number)+
+     &             zdqssed(ig,igcm_stormdust_number) /
+     &               ((pplev(ig,1)-pplev(ig,2))/g) 
+                 zdqssed(ig,igcm_stormdust_mass)=0.
+                 zdqssed(ig,igcm_stormdust_number)=0.
+              ENDDO
+           ENDIF !rdstorm
 
            DO iq=1, nq
@@ -1400,6 +1553,7 @@
              ENDDO
            ENDDO
+
          END IF   ! of IF (sedimentation)
-        
+
 c Add lifted dust to tendancies after sedimentation in the LES (AC)
       IF (turb_resolved) THEN
@@ -1417,7 +1571,6 @@
            ENDDO
       ENDIF
-
 c
-c   6c. Chemical species
+c   8c. Chemical species
 c     ------------------
 
@@ -1493,5 +1646,5 @@
 #endif
 
-c   6d. Updates
+c   8d. Updates
 c     ---------
 
@@ -1511,5 +1664,5 @@
 #ifndef MESOSCALE
 c-----------------------------------------------------------------------
-c   7. THERMOSPHERE CALCULATION
+c   9. THERMOSPHERE CALCULATION
 c-----------------------------------------------------------------------
 
@@ -1535,7 +1688,6 @@
       endif ! of if (callthermos)
 #endif
-
 c-----------------------------------------------------------------------
-c   8. Carbon dioxide condensation-sublimation:
+c   10. Carbon dioxide condensation-sublimation:
 c     (should be the last atmospherical physical process to be computed)
 c   -------------------------------------------
@@ -1607,13 +1759,11 @@
 #endif
       
-      ENDIF  ! of IF (callcond)
-      
-
+      ENDIF  ! of IF (callcond)      
 c-----------------------------------------------------------------------
-c   9. Surface  and sub-surface soil temperature
+c   11. Surface  and sub-surface soil temperature
 c-----------------------------------------------------------------------
 c
 c
-c   9.1 Increment Surface temperature:
+c   11.1 Increment Surface temperature:
 c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1656,5 +1806,5 @@
 
 c
-c   9.2 Compute soil temperatures and subsurface heat flux:
+c   11.2 Compute soil temperatures and subsurface heat flux:
 c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       IF (callsoil) THEN
@@ -1669,8 +1819,40 @@
         ENDIF
       ENDIF
-      
+
+c     To avoid negative values
+      IF (rdstorm) THEN
+           where (pq(:,:,igcm_stormdust_mass) + 
+     &      ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.)
+             pdq(:,:,igcm_stormdust_mass) = 
+     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
+             pdq(:,:,igcm_stormdust_number) = 
+     &       - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30
+           end where
+           where (pq(:,:,igcm_stormdust_number) + 
+     &      ptimestep*pdq(:,:,igcm_stormdust_number) < 0.)
+             pdq(:,:,igcm_stormdust_mass) = 
+     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
+             pdq(:,:,igcm_stormdust_number) = 
+     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+           end where
+
+            where (pq(:,:,igcm_dust_mass) + 
+     &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
+             pdq(:,:,igcm_dust_mass) = 
+     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+             pdq(:,:,igcm_dust_number) = 
+     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+           end where
+           where (pq(:,:,igcm_dust_number) + 
+     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
+             pdq(:,:,igcm_dust_mass) = 
+     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
+             pdq(:,:,igcm_dust_number) = 
+     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
+           end where
+      ENDIF !(rdstorm)   
 
 c-----------------------------------------------------------------------
-c  10. Write output files
+c  12. Write output files
 c  ----------------------
 
@@ -1885,4 +2067,30 @@
               endif
            endif ! of (doubleq)
+
+           if (rdstorm) then   ! diagnostics of stormdust tendancies for 1D and 3D
+              mstormdtot(:)=0
+              mdusttot(:)=0
+              qdusttotal(:,:)=0
+              do ig=1,ngrid 
+                rdsdqdustsurf(ig) = 
+     &                zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
+                rdsdndustsurf(ig) = 
+     &                zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
+                rdsndust(ig,:) =
+     &                pq(ig,:,igcm_stormdust_number)*tauscaling(ig)
+                rdsqdust(ig,:) =
+     &                pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
+                do l=1,nlayer
+                    mstormdtot(ig) = mstormdtot(ig) + 
+     &                      zq(ig,l,igcm_stormdust_mass) * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
+                    mdusttot(ig) = mdusttot(ig) + 
+     &                      zq(ig,l,igcm_dust_mass) * 
+     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
+                    qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
+                enddo
+              enddo
+           endif !(rdstorm)
+
            if (co2clouds) then
               do ig=1,ngrid 
@@ -1906,6 +2114,5 @@
               enddo
            endif ! of if (co2clouds)
-           
-                    
+                              
            if (water) then
              mtot(:)=0
@@ -2138,4 +2345,12 @@
                call wstats(ngrid,'dustN','Dust number',
      &                        'part/kg',3,ndust)
+               if (rdstorm) then
+                 call wstats(ngrid,'reffstormdust','reffdust',
+     &                          'm',3,rstormdust*ref_r0)
+                 call wstats(ngrid,'rdsdustq','Dust mass mr',
+     &                          'kg/kg',3,rdsqdust)
+                 call wstats(ngrid,'rdsdustN','Dust number',
+     &                        'part/kg',3,rdsndust)
+               end if
              else
                do iq=1,dustbin
@@ -2345,5 +2560,5 @@
 c        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
 c    &                   'w.m-2',3,zdtlw)
-  
+ 
             if (.not.activice) then
                CALL WRITEDIAGFI(ngrid,'tauTESap',
@@ -2434,5 +2649,4 @@
      &                       qsurf(1:ngrid,igcm_h2o_ice))
 #endif
-
             CALL WRITEDIAGFI(ngrid,'mtot',
      &                       'total mass of water vapor',
@@ -2443,9 +2657,9 @@
             vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
      &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
-            call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
+            CALL WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
      &                       'mol/mol',3,vmr)
             vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
      &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
-            call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
+            CALL WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
      &                       'mol/mol',3,vmr)
             CALL WRITEDIAGFI(ngrid,'reffice',
@@ -2555,8 +2769,9 @@
      &                        'part/kg',3,ndust)
              call WRITEDIAGFI(ngrid,'dsodust',
-     &                        'dust density scaled opacity',
-     &                        'm2.kg-1',3,dsodust)
-c             call WRITEDIAGFI(ngrid,"tauscaling",
-c     &                    "dust conversion factor"," ",2,tauscaling)
+     &                       'density scaled optcial depth',
+     &                       'm2.kg-1',3,dsodust)
+             call WRITEDIAGFI(ngrid,'dso',
+     &                       'density scaled optcial depth',
+     &                       'm2.kg-1',3,dsodust+dsords)
            else
              do iq=1,dustbin
@@ -2569,4 +2784,55 @@
            endif ! (doubleq)
 
+           if (rdstorm) then  ! writediagfi tendencies stormdust tracers
+             call WRITEDIAGFI(ngrid,'reffdust','reffdust',
+     &                        'm',3,rdust*ref_r0)
+             call WRITEDIAGFI(ngrid,'reffstormdust','reffstormdust',
+     &                        'm',3,rstormdust*ref_r0)
+             call WRITEDIAGFI(ngrid,'mstormdtot',
+     &                        'total mass of stormdust only',
+     &                        'kg.m-2',2,mstormdtot)
+             call WRITEDIAGFI(ngrid,'mdusttot',
+     &                        'total mass of dust only',
+     &                        'kg.m-2',2,mdusttot)
+             call WRITEDIAGFI(ngrid,'rdsdqsdust',
+     &                        'deposited surface stormdust mass',
+     &                        'kg.m-2.s-1',2,rdsdqdustsurf)
+             call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
+     &                        'kg/kg',3,rdsqdust)
+             call WRITEDIAGFI(ngrid,'rdsdustqmodel','storm Dust massmr',
+     &                        'kg/kg',3,pq(:,:,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'rdsdustN','storm Dust number',
+     &                        'part/kg',3,rdsndust)
+             call WRITEDIAGFI(ngrid,"stormfract",
+     &                "fraction of the mesh, with stormdust","none",
+     &                                               2,totstormfract)
+             call WRITEDIAGFI(ngrid,'qsurf',
+     &                 'stormdust injection',
+     &                 'kg.m-2',2,qsurf(:,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'pdqsurf',
+     &                  'tendancy stormdust mass at surface',
+     &                  'kg.m-2',2,dqsurf(:,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
+     &                        'm/s',3,wspeed(:,1:nlayer))
+             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
+     &          ,'sedimentation q','kg.m-2.s-1',3,
+     &          zdqsed(:,:,igcm_dust_mass))
+             call WRITEDIAGFI(ngrid,'zdqssed_dustq'
+     &          ,'sedimentation q','kg.m-2.s-1',2,
+     &          zdqssed(:,igcm_dust_mass))
+             call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
+     &          ,'sedimentation q','kg.m-2.s-1',3,
+     &          zdqsed(:,:,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'rdust','rdust',
+     &                        'm',3,rdust)
+             call WRITEDIAGFI(ngrid,'rstormdust','rstormdust',
+     &                        'm',3,rstormdust)
+             call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
+     &                        'kg/kg',3,qdusttotal)
+             call WRITEDIAGFI(ngrid,'dsords',
+     &                       'density scaled opacity of stormdust',
+     &                       'm2.kg-1',3,dsords)
+           endif ! (rdstorm)
+
            if (scavenging) then
              call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
@@ -2574,4 +2840,8 @@
              call WRITEDIAGFI(ngrid,'ccnN','CCN number',
      &                        'part/kg',3,nccn)
+             call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr',
+     &                        'kg.m-2',2,qsurf(1,igcm_ccn_mass))
+             call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number',
+     &                        'kg.m-2',2,qsurf(1,igcm_ccn_number))
            endif ! (scavenging)
 
@@ -2804,5 +3074,52 @@
              call WRITEDIAGFI(ngrid,'rdust','rdust',
      &                        'm',1,rdust)
-           endif
+           endif ! doubleq 1D
+           if (rdstorm) then
+             call writediagfi(1,'aerosol_dust','opacity of env. dust',''
+     &                        ,1,aerosol(:,:,igcm_dust_mass))
+             call writediagfi(1,'aerosol_stormdust',
+     &                         'opacity of storm dust',''
+     &                        ,1,aerosol(:,:,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
+     &                       'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass))
+             call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
+     &                  'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'mstormdtot',
+     &                        'total mass of stormdust only',
+     &                        'kg.m-2',0,mstormdtot)
+             call WRITEDIAGFI(ngrid,'mdusttot',
+     &                        'total mass of dust only',
+     &                        'kg.m-2',0,mdusttot)
+             call WRITEDIAGFI(ngrid,'tauref',
+     &                    'Dust ref opt depth','NU',0,tauref)
+             call WRITEDIAGFI(ngrid,'rdsdqsdust',
+     &                        'deposited surface stormdust mass',
+     &                        'kg.m-2.s-1',0,rdsdqdustsurf)
+             call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
+     &                        'kg/kg',1,rdsqdust)
+             call WRITEDIAGFI(ngrid,"stormfract", 
+     &                         "fraction of the mesh,with stormdust",
+     &                         "none",0,totstormfract)
+             call WRITEDIAGFI(ngrid,'rdsqsurf',
+     &                 'stormdust at surface',
+     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'qsurf',
+     &                  'dust mass at surface',
+     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass))
+             call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
+     &                        'm/s',1,wspeed)
+             call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
+     &                        'kg/kg',1,qdusttotal)
+             call WRITEDIAGFI(ngrid,'dsords',
+     &                       'density scaled opacity of stormdust',
+     &                       'm2.kg-1',1,dsords)
+             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
+     &          ,'sedimentation q','kg.m-2.s-1',1,
+     &          zdqsed(1,:,igcm_dust_mass))
+             call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
+     &          ,'sedimentation q','kg.m-2.s-1',1,
+     &          zdqsed(1,:,igcm_stormdust_mass))
+           endif !(rdstorm 1D)
+
            if (water.AND.tifeedback) then
              call WRITEDIAGFI(ngrid,"soiltemp",
@@ -2975,4 +3292,5 @@
 
       icount=icount+1
+	
 
       END SUBROUTINE physiq
Index: trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 1974)
+++ trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 1974)
@@ -0,0 +1,1075 @@
+      MODULE rocketduststorm_mod
+
+      IMPLICIT NONE
+
+      REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1)
+      REAL, SAVE, ALLOCATABLE :: alpha_hmons(:) ! slope winds lifting mesh fraction
+      
+      CONTAINS
+
+!=======================================================================
+! ROCKET DUST STORM - vertical transport and detrainment 
+!=======================================================================
+! calculating the vertical flux
+! calling vl_storm : transport scheme of stormdust tracers
+! detrainement of stormdust into nomal background dust
+! -----------------------------------------------------------------------
+! Authors: C. Wang; F. Forget; T. Bertrand 
+! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
+! -----------------------------------------------------------------------
+
+      SUBROUTINE rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,      &
+                                 pq,pdqfi,pt,pdtfi,pplev,pplay,pzlev,  &
+                                 pzlay,pdtsw,pdtlw,                    &
+!             input for radiative transfer
+                                 clearatm,icount,zday,zls,             &
+                                 tsurf,igout,totstormfract,            &
+!             input sub-grid scale cloud
+                                 clearsky,totcloudfrac,                &
+!                   output
+                                 pdqrds,wspeed,dsodust,dsords)
+
+      use tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
+                            ,igcm_dust_mass,igcm_dust_number          &
+                            ,rho_dust 
+      USE comcstfi_h, only: r,g,cpp,rcp
+      use dimradmars_mod, only: albedo,naerkind
+      use comsaison_h, only: dist_sol,mu0,fract
+      use surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons
+      use planetwide_mod, only: planetwide_maxval,planetwide_minval
+!      use rocketstorm_h, only: rdsinject
+      use callradite_mod
+      implicit none
+
+!--------------------------------------------------------
+! Input Variables
+!--------------------------------------------------------
+
+      INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
+      INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
+      INTEGER, INTENT(IN) :: nq ! number of tracer species
+      REAL, INTENT(IN) :: ptime
+      REAL, INTENT(IN) :: ptimestep
+
+      REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
+      REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq
+      REAL, INTENT(IN) :: pt(ngrid,nlayer)      ! temperature at mid-layer (K)
+      REAL, INTENT(IN) :: pdtfi(ngrid,nlayer)   ! tendancy temperature at mid-layer (K)
+
+      REAL, INTENT(IN) :: pplay(ngrid,nlayer)     ! pressure at middle of the layers
+      REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels
+      REAL, INTENT(IN) :: pzlay(ngrid,nlayer)     ! altitude at the middle of the layers
+      REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1)   ! altitude at layer boundaries
+
+      REAL, INTENT(IN) :: pdtsw(ngrid,nlayer)     ! (K/s) env
+      REAL, INTENT(IN) :: pdtlw(ngrid,nlayer)     ! (K/s) env
+
+!     input for second radiative transfer 
+      logical, INTENT(IN) :: clearatm
+      INTEGER, INTENT(INOUT) :: icount
+      real, intent(in) :: zday
+      real, intent(in) :: zls
+      real, intent(in) :: tsurf(ngrid)
+      integer, intent(in) :: igout
+      real, intent(in) :: totstormfract(ngrid)
+!     sbgrid scale water ice clouds
+      logical, intent(in) :: clearsky
+      real, intent(in) :: totcloudfrac(ngrid)      
+ 
+!--------------------------------------------------------
+! Output Variables
+!--------------------------------------------------------
+
+      REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining
+      REAL, INTENT(OUT) :: wspeed(ngrid,nlayer+1)   ! vertical speed within the rocket dust storm
+      REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
+      REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
+
+!--------------------------------------------------------
+! Local variables 
+!--------------------------------------------------------
+      INTEGER l,ig,tsub,iq,ll
+! chao local variables from callradite.F       
+      REAL zdtlw1(ngrid,nlayer)    ! (K/s) storm
+      REAL zdtsw1(ngrid,nlayer)    ! (K/s) storm
+      REAL zt(ngrid,nlayer)        ! actual temperature at mid-layer (K)
+      REAL zdtvert(nlayer)   ! dT/dz , lapse rate
+      REAL ztlev(nlayer)     ! temperature at intermediate levels l+1/2 without last level
+
+      REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2
+      REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer)   ! rad. heating rate at intermediate levels l+1/2
+
+      REAL zqstorm_mass(ngrid,nlayer) ! tracer pq mass intermediate 
+      REAL zqstorm_mass_col(nlayer) 
+      REAL zqstorm_number(ngrid,nlayer) ! tracer field pq number intermediate
+      REAL zqstorm_number_col(nlayer) 
+
+      REAL zqi_mass(ngrid,nlayer)           ! tracer pq mass intermediate 
+      REAL zqi_number(ngrid,nlayer)         ! tracer pq mass intermediate 
+      REAL zdqvlstorm_mass(ngrid,nlayer)    ! tendancy pdq mass after vertical transport
+      REAL zdqvlstorm_number(ngrid,nlayer)  ! tendancy pdq number after vertical transport
+      REAL zdqdetstorm_mass(ngrid,nlayer)   ! tendancy field pq mass after detrainment only
+      REAL zdqdetstorm_number(ngrid,nlayer) ! tendancy field pq number after detrainment only
+
+      REAL zdqenv_mass(ngrid,nlayer)    ! tendancy pdq mass from dust->
+                                        !  stormdust in slp
+      REAL zdqenv_number(ngrid,nlayer)  ! tendancy pdq number from dust->
+                                        !  stormdust in slp
+
+      REAL masse(nlayer)
+      REAL zq(ngrid,nlayer,nq)   ! updated tracers
+      
+      REAL wrds(nlayer)          ! vertical flux within the rocket dust storm
+      REAL wqrdsmass(nlayer+1)   ! flux mass from vl_storm 
+      REAL wqrdsnumber(nlayer+1) ! flux number from vl_storm
+
+      INTEGER nsubtimestep    !number of subtimestep when calling vl_storm 
+      REAL subtimestep        !ptimestep/nsubtimestep
+      REAL dtmax              !considered time needed for dust to cross one layer 
+                              !minimal value over a column
+      logical storm(ngrid)    !logical : true if you have some storm dust in the column
+!     real slope(ngrid)    !logical : true if you don't have storm and have
+                              !a slope
+!     real wslplev(ngrid,nlayer)
+!     real wslp(ngrid,nlayer)
+      REAL coefdetrain           !coefficient for detrainment : % of stormdust detrained
+
+      REAL,PARAMETER:: coefmin =0.025 !C 0<c<1 Minimum fraction of stormdust detrained  
+      REAL,PARAMETER:: detrainlim =0.1!0.25 !L  stormdust detrained if wspeed < detrainlim  
+      REAL,PARAMETER:: wlim =10.   ! maximum vertical speed of rocket storms (m/s)
+      REAL,PARAMETER:: secu=3.      !coefficient on wspeed to avoid dust crossing many layers during subtimestep
+ 
+      ! terms for buoyancy and W^2 in equation:
+      !   w*dw/dz = k1*g*(T'-T)/T - k2*w^2
+      real,parameter:: k1=0.00001
+      real,parameter:: k2=0.25
+      real,parameter:: mu0lim=0.1
+      
+      ! diagnose 
+      REAL detrainment(ngrid,nlayer)
+      real lapserate(ngrid,nlayer) 
+      real deltahr(ngrid,nlayer+1)
+      real stormdust_m0(ngrid,nlayer)
+      real stormdust_m1(ngrid,nlayer)
+      real stormdust_m2(ngrid,nlayer)
+   
+      real hmax,hmin
+!     real zh(ngrid,nlayer)
+    
+      logical,save :: firstcall=.true.
+      real alpha(ngrid) ! scale of the vertical motion (applicable for
+                        ! rds and slp
+      ! variables for radiative transfer
+      REAL  fluxsurf_lw1(ngrid)
+      REAL  fluxsurf_sw1(ngrid,2)
+      REAL  fluxtop_lw1(ngrid)
+      REAL  fluxtop_sw1(ngrid,2)
+      REAL  tauref(ngrid)
+      REAL  tau(ngrid,naerkind)
+      REAL  aerosol(ngrid,nlayer,naerkind)
+      REAL  tauscaling(ngrid)
+      REAL  taucloudtes(ngrid)
+      REAL  rdust(ngrid,nlayer)
+      REAL  rstormdust(ngrid,nlayer)
+      REAL  rice(ngrid,nlayer)
+      REAL  nuice(ngrid,nlayer)
+
+      !variables related to slope,reference layer...
+!     integer lref(ngrid),llref ! the reference layer of slopewind
+!     real buoyt(nlayer) ! buoyancy term when there is a subgrid mountain
+      real slpbg(ngrid)  ! temperature difference at half height of a mountain
+
+      real zqdustslp(ngrid,nlayer),zndustslp(ngrid,nlayer)
+      real zqstormdustslp(ngrid,nlayer),znstormdustslp(ngrid,nlayer)
+
+      real rdsdustqvl0(ngrid,nlayer)
+      real rdsdustqvl1(ngrid,nlayer)
+!     real q2rds(ngrid,nlayer)
+
+      !for second formule of wslp
+      real wtemp(ngrid,nlayer) ! a intermediate variable for wspeed
+ 
+      !merge rds and slp
+      real newzt(ngrid,nlayer) !temperature with perturbation (integrated from
+                               !  vetical motion)
+      real w0(ngrid) !prescribed slope winds at the first layer of atmosphere
+!     real w1(ngrid) !prescribed slope winds at the first layer of atmosphere
+!     real ztb1(ngrid)
+      real wadiabatic(ngrid,nlayer) !for diagnosis
+!     real czt(nlayer),czlay(nlayer),czlev(nlayer+1) !temporary variables
+      real tnew(nlayer) !interpolated temperature profile above the top of Mons
+      real envtemp(nlayer) !interpolated background temperature profile 
+                           ! as the same height as tnew
+      real envt(ngrid,nlayer) ! output,for diagnosing
+      integer scheme(ngrid)   ! diagnose
+
+      ! Chao: for checking conservation of dust
+!     real totdust0(ngrid)
+!     real totdust1(ngrid)
+
+
+      ! **********************************************************************
+      ! **********************************************************************
+      ! Detached dust layers parametrization: two processes are included
+      ! 1) rocket dust storm
+      !     The radiative warming due to the presence of storm dust is 
+      !     balanced by the adiabatic cooling. The tracer "storm dust"  
+      !     is transported by the upward/downward flow. 
+      ! 2) daytime slope winds 
+      !     The daytime thermally driven upslope wind blows dust from the
+      !     bottom to the top of the mountain, upward flow keeps rising untill
+      !     the velocity becomes zero. Both the storm dust and environment dust
+      !     will be transported.
+      ! **********************************************************************
+      ! **********************************************************************      
+      !!    1. Radiative transfer in storm dust
+      !!    2. Compute vertical velocity for storm dust
+      !!      case 1 storm = false and nighttime: nothing to do
+      !!      case 2 daytime slope wind scheme: (mu0(ig) > mu0lim and with storm=false)
+      !!      case 3 rocket dust storm (storm=true)
+      !!    3. Vertical transport
+      !!    4. Detrainment
+      ! **********************************************************************
+      ! **********************************************************************
+      
+
+      !! **********************************************************************
+      !! Firstcall: Evaluate slope wind mesh fraction  
+      IF (firstcall) then   
+        call planetwide_maxval(hmons,hmax )
+        call planetwide_minval(hmons,hmin )
+        do ig=1,ngrid
+          ! It's hard to know the exact the scale of upward flow,
+          !  we assume that the maximun is 10% of the mesh area.
+          alpha_hmons(ig)= 0.1*(hmons(ig)-hmin)/(hmax-hmin)
+        enddo
+        firstcall = .false.
+      ENDIF !firstcall
+      
+      ! **********************************************************************
+      ! 0. Initializations
+      ! **********************************************************************
+      storm(:)=.false.
+      pdqrds(:,:,:) = 0.
+      zdqdetstorm_mass(:,:)=0.
+      zdqdetstorm_number(:,:)=0.
+      wspeed(:,:)=0.
+      detrainment(:,:)=0.
+      zqstorm_mass_col(:)=0.
+      zqstorm_number_col(:)=0.
+      lapserate(:,:)=0.
+      deltahr(:,:)=0.
+      rdsdustqvl0(:,:)=0.
+      rdsdustqvl1(:,:)=0.
+      zqstormdustslp(:,:)=0.
+      znstormdustslp(:,:)=0.
+      zqdustslp(:,:)=0.
+      zndustslp(:,:)=0.
+      zq(:,:,:) = 0.
+      
+      w0(:)=0.
+!     w1(:)=0.
+!     ztb1(:)=0.
+      newzt(:,:)=0
+      wtemp(:,:)=0.
+      wadiabatic(:,:)=0.
+      slpbg(:)=0.
+!     buoyt(:)=0.
+      tnew(:)=0.
+      envtemp(:)=0.
+      envt(:,:)=0.
+      scheme(:)=0
+      alpha(:)=0.
+      stormdust_m0(:,:)=0.
+      stormdust_m1(:,:)=0.
+      stormdust_m2(:,:)=0.
+!     totdust0(:)=0.
+!     totdust1(:)=0.
+      
+      !! no update for the stormdust tracer and temperature fields 
+      !! because previous callradite was for background dust
+      zq(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq)
+      zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)
+
+      !! get actual q for stormdust and dust tracers
+      do l=1,nlayer
+           do ig=1, ngrid     
+             zqi_mass(ig,l)=zq(ig,l,igcm_dust_mass)
+             zqi_number(ig,l)=zq(ig,l,igcm_dust_number)
+             zqstorm_mass(ig,l)=zq(ig,l,igcm_stormdust_mass)
+             zqstorm_number(ig,l)=zq(ig,l,igcm_stormdust_number)      
+             !for diagnostics:
+             stormdust_m0(ig,l)=zqstorm_mass(ig,l)
+           enddo
+      enddo ! of do l=1,nlayer
+
+      !! Check if there is a rocket dust storm 
+      do ig=1,ngrid
+        storm(ig)=.false.
+        do l=1,nlayer
+          if (zqstorm_mass(ig,l)/zqi_mass(ig,l) .gt. 1.E-4) then 
+            storm(ig)=.true.
+            exit
+          endif
+        enddo       
+      enddo
+      
+      ! *********************************************************************
+      ! 1. Call the second radiative transfer for stormdust, obtain the extra heating
+      ! *********************************************************************
+      CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,       &
+                 emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,   &
+                 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1,  &
+                 fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling,    &
+                 taucloudtes,rdust,rice,nuice,co2ice,rstormdust,       &
+                 totstormfract,clearatm,dsords,                 &
+                 clearsky,totcloudfrac)
+
+      ! **********************************************************************
+      ! 2. Compute vertical velocity for storm dust
+      ! **********************************************************************     
+      DO ig=1,ngrid
+        !! **********************************************************************
+        !! 2.1 case 1: Nothing to do when no storm and no slope
+!        IF ((mu0(ig) .LE. mu0lim) .AND. .NOT.(storm(ig)) ) then
+!          scheme(ig)=1
+!          cycle
+!        endif
+        IF ((alpha_hmons(ig) .EQ. 0.) .AND. .NOT.(storm(ig))) then
+          scheme(ig)=1
+          cycle !!no slope 
+        endif
+        
+        
+        ! whatever the situation is, we need the vertical velocity computed by
+        !  the rds scheme
+        zdtvert(1)=0. !This is the env. lapse rate 
+        DO l=1,nlayer-1
+          zdtvert(l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
+          lapserate(ig,l+1)=zdtvert(l+1) !for diagnosing
+        ENDDO
+      
+        ! computing heating rates gradient at boundraies of each layer
+        ! start from surface
+        zdtlw1_lev(1)=0.
+        zdtsw1_lev(1)=0.
+        zdtlw_lev(1)=0.
+        zdtsw_lev(1)=0.
+        ztlev(1)=zt(ig,1)
+
+        DO l=1,nlayer-1
+          ! Calculation for the dust storm fraction
+          zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
+                       zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
+                          (pzlay(ig,l+1)-pzlay(ig,l))
+         
+          zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
+                       zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
+                          (pzlay(ig,l+1)-pzlay(ig,l))
+          !MV18: calculation for the background dust fraction
+          zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
+                       pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
+                          (pzlay(ig,l+1)-pzlay(ig,l))
+         
+          zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
+                       pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
+                          (pzlay(ig,l+1)-pzlay(ig,l))
+         
+          ztlev(l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+          &
+                       zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /      &
+                          (pzlay(ig,l+1)-pzlay(ig,l))
+        ENDDO
+      
+        DO l=1,nlayer
+          deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l))  &
+                                      -(zdtlw_lev(l)+zdtsw_lev(l))
+          wadiabatic(ig,l)=-deltahr(ig,l)/(g/cpp+   &
+                                     max(zdtvert(l),-0.99*g/cpp))
+      
+          !limit vertical wind in case of lapse rate close to adiabat
+          wadiabatic(ig,l)=max(wadiabatic(ig,l),-wlim) 
+          wadiabatic(ig,l)=min(wadiabatic(ig,l),wlim)
+        ENDDO
+        
+          !! **********************************************************************
+          !! 2.2 case 2: daytime slope wind scheme      
+          IF ((mu0(ig) .gt. mu0lim) .and. .not. storm(ig)) then
+            scheme(ig)=2
+            alpha(ig) = alpha_hmons(ig)
+            ! interpolate the env. temperature above the mountain top 
+            call intep_vtemp(nlayer,hmons(ig),zt(ig,:),pzlay(ig,:), &
+                               envtemp,slpbg(ig))
+            envt(ig,:)=envtemp(:) !for diagnosis
+ 
+            ! second: estimate the vertical velocity at boundraies of each layer
+            !wspeed(ig,1)=0. ! at surface, already initialized
+ 
+            !!!!!!!the first layer of atmosphere!!!!!!!!!!
+            IF (slpbg(ig) .gt. 0.) THEN !only postive buoyancy
+               ! if slpbg(ig) lt 0, means the slope flow is colder than env. (night or
+               ! early morning ?)  !!!!!!!!!!!
+               ! ideal method
+              !w1(ig)=-sqrt(g*slpbg(ig)/zt(ig,1)*hmons(ig))*sin(zsig(ig))
+               ! new scheme, simply proportional to temperature and
+               ! mountain height
+               w0(ig)=-1.5e-4*g*slpbg(ig)/zt(ig,1)*hmons(ig)
+               ! otherwise, w0(ig) =0.
+               wspeed(ig,2)=w0(ig)
+            ELSE
+               wspeed(ig,2)=wadiabatic(ig,2) !! for early morning ?
+            ENDIF
+       
+            ! prepare the integration, NOTE: if w is too small, may have artificials
+            IF (abs(wspeed(ig,2)) .lt. 0.01 ) & 
+                                 wspeed(ig,2)=sign(0.01,wspeed(ig,2))
+ 
+            newzt(ig,1)= zt(ig,1) !temperature of the first layer atmosphere
+                                  !  above the mountain top (radiative
+                                  !  equilibrium on Mars)
+
+            ! estimate the vertical velocities
+            ! if w0(ig) >= 0, means downward motion, no upslope winds, we switch to
+            ! assume that the extra heating integrally convert to
+            ! vertical motion.
+            if ( w0(ig) .ge. 0 ) then  !! normal, it is impossible,
+                                       !!  because mu(ig)>0.1 here
+               do l=3,nlayer
+                  wspeed(ig,l)=wadiabatic(ig,l)
+               enddo
+            else
+            ! estimate the velocities by taking into account the heating due
+            ! to storm dust, the cooling due to vertical motion ...
+            !!!!!!!!!!!the simple scheme!!!!!!!!!
+               do l=2,nlayer-1
+                 !if superadiabatic layer
+                 if ( zdtvert(l) .lt. -g/cpp) then
+                    !case 1
+                    ! test, also decrease adiabatically ?
+                   !newzt(ig,l)= &
+                   !        zt(ig,l-1)-g/cpp*(pzlay(ig,l)-pzlay(ig,l-1))
+                    newzt(ig,l)=zt(ig,l)
+                   !wspeed(ig,l+1)=wspeed(ig,l)
+                 else
+                    !not superadiabatic
+                    newzt(ig,l)=newzt(ig,l-1)+(deltahr(ig,l)/ &
+                                  (-wspeed(ig,l))-g/cpp)*        &
+                                  (pzlay(ig,l)-pzlay(ig,l-1))
+                 endif ! end of if superadiabatic or not
+ 
+                !wtemp(ig,l+1)=wspeed(ig,l)**2+2.*g*(pzlev(ig,l+1) & 
+                !              -pzlev(ig,l))*(k1*(newzt(ig,l)       &
+                !              -envtemp(l))/envtemp(l))
+                 wtemp(ig,l+1)=(1.-2.*k1*(pzlev(ig,l+1)-pzlev(ig,l)))*&
+                                wspeed(ig,l)**2+2.*k2*g*(pzlev(ig,l+1) & 
+                               -pzlev(ig,l))*((newzt(ig,l)       &
+                               -envtemp(l))/envtemp(l))
+ 
+                 if (wtemp(ig,l+1) .gt. 0.) then
+                   !case 2
+                   wspeed(ig,l+1)=-sqrt(wtemp(ig,l+1))
+
+                   ! if |wspeed| < |wadiabatic| then go to wadiabatic
+                   if (wspeed(ig,l+1) .gt. wadiabatic(ig,l+1)) then
+                     do ll=l,nlayer-1
+                       newzt(ig,ll)=envtemp(ll)
+                       wspeed(ig,ll+1)=wadiabatic(ig,ll+1)
+                     enddo 
+                     exit
+                   endif
+
+                   ! avoid artificials
+                   if (abs(wspeed(ig,l+1)) .lt. 0.01 ) &
+                              wspeed(ig,l+1)=sign(0.01,wspeed(ig,l+1))
+       
+                 else if (l .lt. nlayer) then
+                   !case 3
+                   do ll=l,nlayer-1
+                     newzt(ig,ll)=envtemp(ll)
+                     wspeed(ig,ll+1)=wadiabatic(ig,ll+1)
+                   enddo !overshot
+                   exit
+          
+                 else
+                   wspeed(ig,l+1)=0.
+                 endif
+       
+               enddo
+       
+            endif !w0
+ 
+         ELSE
+          !! **********************************************************************
+          !! 2.3 case 3: storm=true
+          if (storm(ig)) then
+              scheme(ig)=3
+              alpha(ig) = totstormfract(ig)
+              do l=1,nlayer
+                  wspeed(ig,l)=wadiabatic(ig,l)
+              enddo
+           endif ! storm=1
+ 
+         ENDIF ! rds or slp
+
+
+        !!!!!!!! estimate the amount of dust for diagnostics
+        DO l=1,nlayer
+            ! transfer background dust + storm dust (concentrated)
+            zqstormdustslp(ig,l) =zqi_mass(ig,l)+ &
+                                      zqstorm_mass(ig,l)/alpha(ig)
+            znstormdustslp(ig,l) =zqi_number(ig,l)+ &
+                                    zqstorm_number(ig,l)/alpha(ig)
+            zqdustslp(ig,l) = zqi_mass(ig,l)
+            zndustslp(ig,l) = zqi_number(ig,l)
+      
+            rdsdustqvl0(ig,l)=zqstormdustslp(ig,l) !for diagnosis
+        ENDDO
+
+      ! **********************************************************************
+      ! 3. Vertical transport
+      ! **********************************************************************
+        do l=1,nlayer
+             masse(l)=(pplev(ig,l)-pplev(ig,l+1))/g
+        enddo
+        !Estimation of "dtmax" (s) to be used for vertical transport    
+        dtmax=ptimestep
+        !secu is a margin on subtimstep to avoid dust crossing many layers
+        do l=2,nlayer
+          if (wspeed(ig,l).lt.0.) then       ! case up
+             dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/  &
+                               (secu*abs(wspeed(ig,l))))
+          else if (wspeed(ig,l).gt.0.) then
+             dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/  &
+                               (secu*abs(wspeed(ig,l))))
+          endif
+        enddo
+
+        nsubtimestep= int(ptimestep/dtmax) 
+        subtimestep=ptimestep/float(nsubtimestep)
+        
+        do l=1,nlayer
+           wrds(l)=wspeed(ig,l)*pplev(ig,l)/(r*ztlev(l)) &
+                    *subtimestep
+        enddo
+      
+        do l=1,nlayer
+           zqstorm_mass_col(l)= zqstormdustslp(ig,l) !zqstorm_mass(ig,l)
+           zqstorm_number_col(l)=znstormdustslp(ig,l) ! zqstorm_number(ig,l)
+        enddo
+        
+        do tsub=1,nsubtimestep
+           wqrdsmass(:)=0.
+           wqrdsnumber(:)=0.
+           CALL vl_storm(nlayer,zqstorm_mass_col,2.,   &
+                 masse,wrds ,wqrdsmass)
+           CALL vl_storm(nlayer,zqstorm_number_col,2.,  &
+                 masse,wrds ,wqrdsnumber)
+        enddo
+      
+        !!!!!generate the "extra" dust
+        do l=1,nlayer
+           rdsdustqvl1(ig,l)=zqstorm_mass_col(l) ! for diagnosis
+
+          ! extra dust = storm dust
+          !zqdustslp(ig,l)=zqi_mass(ig,l) !(1.-alpha(ig))*zqi_mass(ig,l)
+          !zndustslp(ig,l)=zqi_number(ig,l) !(1.-alpha(ig))*zqi_number(ig,l)
+          !zqstorm_mass_col(l)=alpha(ig)*zqstorm_mass_col(l)
+          !zqstorm_number_col(l)=alpha(ig)*zqstorm_number_col(l)
+
+           !with compensation
+           if (zqstorm_mass_col(l) .lt. zqi_mass(ig,l)  ) then
+              ! the following two equations are easier to understand
+              zqdustslp(ig,l)=(1.-alpha(ig))*zqi_mass(ig,l)+alpha(ig)* &
+                         zqstorm_mass_col(l)
+              zndustslp(ig,l)=(1.-alpha(ig))*zqi_number(ig,l)+alpha(ig)&
+                         *zqstorm_number_col(l)
+              !with a bug, should be  zqi+alpha****
+              !zqdustslp(ig,l)=zqi_mass(ig,l)-alpha(ig)* &
+              !           (zqstorm_mass_col(l)-zqi_mass(ig,l))
+              !zndustslp(ig,l)=zqi_number(ig,l)-alpha(ig)* &
+              !           (zqstorm_number_col(l)-zqi_number(ig,l))
+              zqstorm_mass_col(l)=0.  
+              zqstorm_number_col(l)=0.
+           else
+              zqstorm_mass_col(l)=alpha(ig)* &
+                         (zqstorm_mass_col(l)-zqi_mass(ig,l))
+              zqstorm_number_col(l)=alpha(ig)* &
+                         (zqstorm_number_col(l)-zqi_number(ig,l))
+              ! the mass mixing ratio of environmental dust doesn't change.
+           endif
+           !diagnose
+           stormdust_m1(ig,l)=zqstorm_mass_col(l)
+        enddo
+
+        !=======================================================================
+        ! calculate the tendencies due to vertical transport
+        do l=1,nlayer
+           ! tendencies due to vertical transport 
+           zdqvlstorm_mass(ig,l)= (zqstorm_mass_col(l)- &
+                                   zqstorm_mass(ig,l)) /ptimestep
+           zdqvlstorm_number(ig,l)= (zqstorm_number_col(l)- &
+                                 zqstorm_number(ig,l)) /ptimestep
+
+           zdqenv_mass(ig,l)=(zqdustslp(ig,l)-zqi_mass(ig,l))/ptimestep
+           zdqenv_number(ig,l)=(zndustslp(ig,l)-zqi_number(ig,l))  &
+                                /ptimestep
+
+           ! chao for output only
+          !qstormdustvl1(ig,l)=zqstorm_mass_col(l)
+          !nstormdustvl1(ig,l)=zqstorm_number_col(l)
+          !stormdust_m_col1(ig)=stormdust_m_col1(ig)+zqstorm_mass_col(l)  &
+          !                      *(pplev(ig,l)-pplev(ig,l+1))/g
+          !rdsdustqvl1(ig,l)=zqstorm_mass_col(l)
+        enddo
+
+      ! **********************************************************************
+      ! 4. Detrainment: convert dust storm to background dust
+      ! **********************************************************************
+        do l=1,nlayer
+             ! compute the coefficient of detrainment
+          if ((max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) .lt.  & 
+                        detrainlim) .or. (zqdustslp(ig,l) .gt.  &
+                                   10000.*zqstorm_mass_col(l))) then
+             coefdetrain=1.
+          else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))   &
+                                                     .le. wlim) then
+                  ! case where detrainment depends on vertical wind 
+           ! coefdetrain=0.5*(((1-coefmin)/(detrainlim-3.)**2)*     &
+           !     (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-3.)**2 &
+           !                                               +coefmin)
+             coefdetrain=1.*(((1-coefmin)/(detrainlim-wlim)**2)*     &
+                 (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-wlim)**2 &
+                                                             +coefmin)
+           !coefdetrain=0.5
+          else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))).gt. 10. )&
+                  then
+             coefdetrain=0.025
+          else
+             coefdetrain=coefmin
+          endif
+
+          detrainment(ig,l)=coefdetrain !for diagnose
+
+          ! Calculate tendancies corresponding to pdq after detrainement
+          ! pdqdet = tendancy corresponding to detrainment only 
+          zdqdetstorm_mass(ig,l)=-coefdetrain*zqstorm_mass_col(l) &
+                                                        /ptimestep
+          zdqdetstorm_number(ig,l)=-coefdetrain*zqstorm_number_col(l) &
+                                                        /ptimestep
+      
+          ! pdqrds ( tendancy corresponding to vertical transport and
+          ! detrainment) = zdqvlstorm + pdqdet
+          pdqrds(ig,l,igcm_stormdust_mass)=zdqdetstorm_mass(ig,l) &
+                                                 +zdqvlstorm_mass(ig,l)
+          pdqrds(ig,l,igcm_stormdust_number)=zdqdetstorm_number(ig,l) &
+                                               +zdqvlstorm_number(ig,l)
+          pdqrds(ig,l,igcm_dust_mass)= zdqenv_mass(ig,l) &
+                             -zdqdetstorm_mass(ig,l) 
+          pdqrds(ig,l,igcm_dust_number)= zdqenv_number(ig,l) &
+                             -zdqdetstorm_number(ig,l)
+
+          !diagnose 
+          stormdust_m2(ig,l)=zqstorm_mass_col(l)-coefdetrain*zqstorm_mass_col(l) 
+        enddo ! nlayer
+!       endif
+      !=======================================================================
+      enddo ! end do ig=1,ngrid
+ 
+!     !chao check conservation here
+!     do l=1,nlayer
+!       do ig=1,ngrid
+!         totdust0(ig)=totdust0(ig)+                &
+!            zq(ig,l,igcm_stormdust_mass)*        &
+!            ((pplev(ig,l) - pplev(ig,l+1)) / g)  &
+!            +  zq(ig,l,igcm_dust_mass)*          &
+!            ((pplev(ig,l) - pplev(ig,l+1)) / g)
+
+!         totdust1(ig)=totdust1(ig)+                &
+!            (zq(ig,l,igcm_stormdust_mass) + & 
+!             pdqrds(ig,l,igcm_stormdust_mass)*ptimestep)*        &
+!            ((pplev(ig,l) - pplev(ig,l+1)) / g)  &
+!            +  ( zq(ig,l,igcm_dust_mass)+          &
+!             pdqrds(ig,l,igcm_dust_mass)*ptimestep)*        &
+!            ((pplev(ig,l) - pplev(ig,l+1)) / g)
+!       enddo
+!     enddo
+     
+!       call writediagfi(ngrid,'totdust0','total dust before rds', &
+!              ' ',2,totdust0)
+!       call writediagfi(ngrid,'totdust1','total dust after rds', &
+!              ' ',2,totdust1)
+        !output for diagnosis
+        call WRITEDIAGFI(ngrid,'detrainment',         &
+                        'coefficient of detrainment',' ',3,detrainment)
+       !call WRITEDIAGFI(ngrid,'qstormvl1','mmr of stormdust after rds_vl', &
+       !   &                        'kg/kg',3,qstormdustvl1)
+        call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', &
+           &                        'k/m',3,lapserate)
+        call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', &
+           &                        'K/s',3,deltahr)
+        call WRITEDIAGFI(ngrid,'wold', &
+                             'wind generated from adiabatic colling ', &
+           &                        'm/s',3,wadiabatic)
+        call WRITEDIAGFI(ngrid,'newzt','perturbated temperature', &
+           &                        'K/s',3,newzt)
+        call WRITEDIAGFI(ngrid,'zt','unperturbated temperature', &
+           &                        'K/s',3,zt)
+        call WRITEDIAGFI(ngrid,'wtemp','under square root', &
+           &                        'K/s',3,wtemp)
+       !call WRITEDIAGFI(ngrid,'stormdust_m_col1','mass of stormdust after rds_vl', &
+       !   &                        'kg/kg',2,stormdust_m_col1)
+       !call WRITEDIAGFI(ngrid,'temprds','temp for calculating zdtvert', &
+       !   &                        'k',3,temprds)
+        call WRITEDIAGFI(ngrid,'stormdust_m0','mass col  of stormdust before rds_vl', &
+           &                        'kg/kg',3,stormdust_m0)
+        call WRITEDIAGFI(ngrid,'stormdust_m1','mass col  of stormdust after rds_vl', &
+           &                        'kg/kg',3,stormdust_m1)
+        call WRITEDIAGFI(ngrid,'stormdust_m2','mass col  of stormdust after rds_vl', &
+           &                        'kg/kg',3,stormdust_m2)
+      
+      ! call writediagfi(ngrid,'wslp','estimated slope winds','m/s',3,wslp)
+      ! call writediagfi(ngrid,'wslp2','estimated slope winds2','m/s',3,wslp2)
+      ! call writediagfi(ngrid,'zhb','estimated slope winds2','m/s',3,zhb)
+      ! call writediagfi(ngrid,'bruntf','bouyancy frequency',' ',3,bruntf)
+      ! call writediagfi(ngrid,'slpdepth','slope depth','m',2,slpdepth)
+      ! call writediagfi(ngrid,'slpu','perbulation u','m/s',3,slpu)
+      ! call writediagfi(ngrid,'slpzh','perbulation zh',' ',3,slpzh)
+      ! call writediagfi(ngrid,'zqslp','zq in rocketduststorm','ikg/kg',3, &
+      !                     zq(:,:,igcm_dust_mass))
+      ! call writediagfi(ngrid,'zrdsqslp','zq in rocketduststorm','ikg/kg',3, &
+      !                     zq(:,:,igcm_stormdust_mass))
+      ! call writediagfi(ngrid,'wslplev','estimated slope winds','m/s',3,wslplev)
+      ! call writediagfi(ngrid,'slope','identified slope wind effect',' ',2,slope)
+        call writediagfi(ngrid,'w0','max of slope wind',' ',2,w0)
+!       call writediagfi(ngrid,'w1','max of slope wind',' ',2,w1)
+        call writediagfi(ngrid,'mu0','cosine of solar incidence angle',&
+                                                           ' ',2,mu0)
+!       call writediagfi(ngrid,'storm','identified rocket dust storm',&
+!                                                  ' ',2,real(storm))
+        call writediagfi(ngrid,'scheme','which scheme',&
+                                                   ' ',2,real(scheme))
+        call writediagfi(ngrid,'alpha','coefficient alpha',' ',2,alpha)
+      ! call writediagfi(ngrid,'q2rds','alpha zq',' ', &
+      !                         3,q2rds)
+        call writediagfi(ngrid,'rdsdustqvl0','not vl storm slp',       &
+                                              'kg/kg',3,zqstormdustslp)
+        call writediagfi(ngrid,'rdsdustqvl1','vled storm slp',         &
+                                                 'kg/kg',3,rdsdustqvl1)
+        call writediagfi(ngrid,'dustqvl0','not vl  slp',               &
+                                                    'kg/kg',3,zqi_mass)
+        call writediagfi(ngrid,'dustqvl1','vled  slp',                 &
+                                                   'kg/kg',3,zqdustslp)
+      ! call WRITEDIAGFI(ngrid,'lmax_th2', &
+      !                  'hauteur du thermique','point', &
+      !                             2,real(lmax_th(:)))
+      ! call WRITEDIAGFI(ngrid,'zmax_th2', &
+      !                  'hauteur du thermique','m', &
+      !                             2,zmax_th)
+      ! call writediagfi(ngrid,'lslope','lenght of slope',' ',2,lslope)
+      ! call writediagfi(ngrid,'hmon','identified slope wind effect',' ',2,hmon)
+        call writediagfi(ngrid,'envt','interpolated env. temp.',       &
+                                                           'K',3,envt)
+        call writediagfi(ngrid,'hmons','identified slope wind effect', &
+                                                           ' ',2,hmons)
+      ! call writediagfi(ngrid,'slpbg','temp. diff along slope',       &
+      !                                                    ' ',2,slpbg)
+      ! call writediagfi(ngrid,'zmea','identified slope wind effect',  &
+      !                                                    ' ',2,zmea)
+      ! call writediagfi(ngrid,'zsig','identified slope wind effect',  &
+      !                                                    ' ',2,zsig)
+      ! call writediagfi(ngrid,'zhslpenv','difference of zh above mons',' ',2,zhslpenv)
+      ! call writediagfi(ngrid,'lref','identified slope wind effect',' ',2,real(lref))
+        END SUBROUTINE rocketduststorm
+
+!********************************************************************************
+!********************************************************************************
+      SUBROUTINE vl_storm(nlay,q,pente_max,masse,w,wq)
+!
+!     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
+!
+!    ********************************************************************
+!     Shema  d'advection " pseudo amont " dans la verticale
+!    pour appel dans la physique (sedimentation)
+!    ********************************************************************
+!    q rapport de melange (kg/kg)...
+!    masse : masse de la couche Dp/g
+!    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
+!    pente_max = 2 conseillee
+!
+!
+!   --------------------------------------------------------------------
+      IMPLICIT NONE
+!
+
+!   Arguments:
+!   ----------
+      integer,intent(in) :: nlay ! number of atmospheric layers
+      real masse(nlay),pente_max
+      REAL q(nlay)
+      REAL w(nlay)
+      REAL wq(nlay+1)
+!
+!      Local 
+!   ---------
+!
+      INTEGER i,l,j,ii
+!
+
+      real dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
+      real newmasse
+      real sigw, Mtot, MQtot
+      integer m
+
+      REAL      SSUM,CVMGP,CVMGT
+      integer ismax,ismin
+
+
+!    On oriente tout dans le sens de la pression c'est a dire dans le
+!    sens de W
+
+      do l=2,nlay
+            dzqw(l)=q(l-1)-q(l)
+            adzqw(l)=abs(dzqw(l))
+      enddo
+
+      do l=2,nlay-1
+#ifdef CRAY
+            dzq(l)=0.5*
+     ,      cvmgp(dzqw(l)+dzqw(l+1),0.,dzqw(l)*dzqw(l+1))
+#else
+            if(dzqw(l)*dzqw(l+1).gt.0.) then
+                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
+            else
+                dzq(l)=0.
+            endif
+#endif
+            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
+            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
+      enddo
+
+      dzq(1)=0.
+      dzq(nlay)=0.
+       
+!      write(*,*),'TB14 wq before up',wq(1,:)
+!      write(*,*),'TB14 q before up',q(1,:)
+! ---------------------------------------------------------------
+!   .... calcul des termes d'advection verticale  .......
+! ---------------------------------------------------------------
+
+! calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
+!
+!      No flux at the model top:
+       wq(nlay+1)=0.
+
+!      1) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION)     
+!      ===============================
+
+!      Surface flux up:
+         if(w(1).lt.0.) wq(1)=0. ! warning : not always valid
+
+       do l = 1,nlay-1  ! loop different than when w>0
+         if(w(l+1).le.0)then
+
+!         Regular scheme (transfered mass < 1 layer)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          if(-w(l+1).le.masse(l))then
+            sigw=w(l+1)/masse(l)
+            wq(l+1)=w(l+1)*(q(l)-0.5*(1.+sigw)*dzq(l))
+!         Extended scheme (transfered mass > 1 layer)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          else 
+             m = l-1
+             Mtot = masse(m+1)
+             MQtot = masse(m+1)*q(m+1)
+             if (m.le.0)goto 77
+             do while(-w(l+1).gt.(Mtot+masse(m)))
+!             do while(-w(l+1).gt.Mtot)
+                m=m-1
+                Mtot = Mtot + masse(m+1)
+                MQtot = MQtot + masse(m+1)*q(m+1)
+                if (m.le.0)goto 77
+             end do
+ 77          continue
+
+             if (m.gt.0) then
+                sigw=(w(l+1)+Mtot)/masse(m)
+                wq(l+1)= (MQtot + (-w(l+1)-Mtot)*         &
+                (q(m)-0.5*(1.+sigw)*dzq(m))  )
+             else
+! new
+                w(l+1) = -Mtot
+                wq(l+1) = -MQtot 
+!                write(*,*) 'TB14 MQtot = ',MQtot,l
+
+! old
+!                wq(l+1)= (MQtot + (-w(l+1)-Mtot)*qm(1))
+!                write(*,*) 'a rather weird situation in vlz_fi !'
+!                stop
+             end if
+          endif
+         endif  ! w<0 (up)
+       enddo
+
+       do l = 1,nlay-1  ! loop different than when w>0
+
+          q(l)=q(l) +  (wq(l+1)-wq(l))/masse(l)
+
+       enddo
+
+!       write(*,*),'TB14 masse',masse(1,:)
+!       write(*,*),'TB14 wq before down after up',wq(1,:)
+!       write(*,*),'TB14 q before down',q(1,:)
+     
+!      2) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)     
+!      ===============================
+
+!      Initialisation wq = 0 to consider now only downward flux
+         wq(:)=0. ! 
+
+       do l = 1,nlay          ! loop different than when w<0
+
+         if(w(l).gt.0.)then
+
+!         Regular scheme (transfered mass < 1 layer)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          if(w(l).le.masse(l))then
+            sigw=w(l)/masse(l)
+            wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l))
+!            write(*,*),'TB14 wq after up',wq(1,:)
+            
+
+!         Extended scheme (transfered mass > 1 layer)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          else 
+            m=l
+            Mtot = masse(m)
+            MQtot = masse(m)*q(m)
+            if(m.ge.nlay)goto 88
+            do while(w(l).gt.(Mtot+masse(m+1)))
+                m=m+1
+                Mtot = Mtot + masse(m)
+                MQtot = MQtot + masse(m)*q(m)
+                if(m.ge.nlay)goto 88
+            end do
+ 88         continue
+            if (m.lt.nlay) then
+                sigw=(w(l)-Mtot)/masse(m+1)
+                wq(l)=(MQtot + (w(l)-Mtot)* &
+                          (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
+            else
+                w(l) = Mtot
+                wq(l) = MQtot 
+            end if
+          end if
+         end if ! w>0 (down)
+       enddo
+       
+       do l = 1,nlay          ! loop different than when w<0
+
+          q(l)=q(l) +  (wq(l+1)-wq(l))/masse(l)
+
+       enddo
+      end subroutine vl_storm
+!********************************************************************************
+      SUBROUTINE intep_vtemp(nlayer,hm,temp,zlay,envtemp,slpb)
+
+       USE comcstfi_h, only: g,cpp
+
+       implicit none
+       ! this subroutine is using for obtaining the environmental
+       ! temperature profile when a subgrid mountain exists.
+
+
+!      input:
+       integer,intent(in) :: nlayer
+       real,intent(in) :: hm  ! the height of mountain
+       real,intent(in) :: temp(nlayer)   !large scale temp. profile of one mesh
+       real,intent(in) :: zlay(nlayer)   ! height at the middle of each layer
+
+!      output:
+       real,intent(out) :: envtemp(nlayer) ! environment temperature 
+       real,intent(out) :: slpb !the temperature difference between slope and
+                                ! env. at the half height of mountain
+
+!      local variables
+       integer l,il,tmpl
+       integer lnew  !the layer of atmosphere above the mountain
+                             ! corresponding to the env. (for buoyancy
+                             ! calc. ) 
+       real newh(nlayer)     !height at the middle of each layer
+                             ! account for the exist of mountain
+!      real g,cpp 
+       real halfh !  half the height of a mountain
+
+       !initilize the array
+       lnew=0
+       newh(:)=0.
+       envtemp(:)=0.
+       tmpl=1  
+
+       do l=1,nlayer
+         newh(l)=hm+zlay(l)
+         
+         do il=tmpl,nlayer-1 !MV18: added the -1
+            if (newh(l) .ge. zlay(il) .and. newh(l) .lt. zlay(il+1))then
+              ! find the corresponding layer
+              lnew=il
+
+              ! interpolate
+              envtemp(l) = temp(il)+(newh(l)-zlay(lnew))*&
+                          (temp(il+1)-temp(il))/(zlay(il+1)-zlay(il))
+
+              exit !go to the next layer
+            else if (newh(l) .ge. zlay(nlayer)) then
+              ! higher than the highest layer
+              lnew=nlayer
+              envtemp(l)=temp(nlayer)   
+
+            endif
+         enddo
+         ! this can accelerate the loop
+         tmpl=lnew
+
+       enddo
+
+       halfh=0.5*hm
+       if (halfh .le. zlay(1) ) then 
+         slpb=0.
+       else if (halfh .gt. zlay(nlayer)) then 
+         !normally, impossible for halfh gt zlay(l), anyway...
+         tmpl=nlayer
+         !difference between surface and atmosphere (env.)
+         slpb=temp(1)-(temp(nlayer-1)+((halfh-zlay(nlayer-1))* &
+              (temp(tmpl)-temp(tmpl-1))/(zlay(tmpl)-zlay(tmpl-1)))) 
+       else
+         do l=1,nlayer-1
+           if ((halfh .gt. zlay(l)) .and. (halfh .le. zlay(l+1)))then
+             tmpl= l
+             slpb=temp(1)-(temp(tmpl)+(halfh-zlay(tmpl))*  &
+                 (temp(tmpl+1)-temp(tmpl))/(zlay(tmpl+1)-zlay(tmpl)))
+           endif
+         enddo
+       endif          
+       end subroutine intep_vtemp
+
+       ! initialization module variables
+       subroutine ini_rocketduststorm_mod(ngrid)
+       
+       implicit none
+       
+       integer, intent(in) :: ngrid
+       
+       allocate(dustliftday(ngrid))
+       allocate(alpha_hmons(ngrid))
+       
+       end subroutine ini_rocketduststorm_mod
+       
+       subroutine end_rocketduststorm_mod
+       
+       implicit none
+       
+       if (allocated(dustliftday)) deallocate(dustliftday)
+       if (allocated(alpha_hmons)) deallocate(alpha_hmons)
+
+       end subroutine end_rocketduststorm_mod       
+      
+      END MODULE rocketduststorm_mod
Index: trunk/LMDZ.MARS/libf/phymars/suaer.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/suaer.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/suaer.F90	(revision 1974)
@@ -6,4 +6,5 @@
                     iaer_dust_conrath,iaer_dust_doubleq,&
                     iaer_dust_submicron,iaer_h2o_ice,&
+                    iaer_stormdust_doubleq,&
                     file_id,radiustab, gvis, omegavis, &
                     QVISsQREF, gIR, omegaIR, &
@@ -153,4 +154,16 @@
 !       Reference wavelength in the infrared:
         longrefir(iaer)=12.1E-6  ! 825cm-1 TES/MGS
+!==================================================================
+        CASE("stormdust_doubleq") aerkind   ! Two-moment scheme for stormdust - radiative properties
+!==================================================================
+!       Visible domain:
+        file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix
+!       Infrared domain:
+        file_id(iaer,2) = 'optprop_dustir_n50.dat'     !Mie
+!       Reference wavelength in the visible:
+        longrefvis(iaer)=0.67E-6
+!       If not equal to 0.67e-6 -> change readtesassim accordingly;
+!       Reference wavelength in the infrared:
+        longrefir(iaer)=dustrefir
 !==================================================================
         END SELECT aerkind
Index: trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 1974)
@@ -20,4 +20,5 @@
   real,save :: iceradius(2) , dtemisice(2)
   real,save,allocatable :: zmea(:),zstd(:),zsig(:),zgam(:),zthe(:)
+  real,save,allocatable :: hmons(:)
   real,save,allocatable :: z0(:) ! surface roughness length (m)
   real,save :: z0_default ! default (constant over planet) surface roughness (m)
@@ -55,5 +56,6 @@
     allocate(capcal(ngrid))
     allocate(fluxgrd(ngrid))
- 
+    allocate(hmons(ngrid))
+    
   end subroutine ini_surfdat_h
 
@@ -79,4 +81,5 @@
     if (allocated(capcal))      deallocate(capcal)
     if (allocated(fluxgrd))     deallocate(fluxgrd)
+    if (allocated(hmons))       deallocate(hmons)
 
   end subroutine end_surfdat_h
Index: trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/tracer_mod.F90	(revision 1974)
@@ -6,5 +6,5 @@
       integer,save :: nqmx ! initialized in conf_phys
    
-      character*20,allocatable,save ::  noms(:)  ! name of the tracer
+      character*30,allocatable,save ::  noms(:)  ! name of the tracer
       real,allocatable,save :: mmol(:)           ! mole mass of tracer (g/mol-1) 
       real,allocatable,save :: radius(:)   ! dust and ice particle radius (m)
@@ -23,5 +23,5 @@
       real,save :: nuiceco2_sed   ! Sedimentation effective variance of the co2 ice dist.
       real,save :: nuiceco2_ref   ! Effective variance of the co2 ice dist.
-
+      
       real,save :: ccn_factor  ! ratio of nuclei for water ice particles
 
@@ -42,5 +42,7 @@
       integer,save :: igcm_ccn_number ! CCN number mixing ratio
       integer,save :: igcm_dust_submicron ! submicron dust mixing ratio
-    
+      integer,save :: igcm_stormdust_mass   !  storm dust mass mixing ratio
+      integer,save :: igcm_stormdust_number !  storm dust number mixing ratio
+   
       integer,save :: igcm_ccnco2_mass   ! CCN (dust and/or water ice) for CO2 mass mixing ratio
       integer,save :: igcm_ccnco2_number ! CCN (dust and/or water ice) for CO2 number mixing ratio
Index: trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F	(revision 1974)
@@ -6,5 +6,5 @@
       
       SUBROUTINE updatereffrad(ngrid,nlayer,
-     &                rdust,rice,nuice,
+     &                rdust,rstormdust,rice,nuice,
      &                reffrad,nueffrad,
      &                pq,tauscaling,tau,pplay)
@@ -14,9 +14,11 @@
      &                       igcm_h2o_ice, igcm_ccn_mass, radius,
      &                       igcm_ccn_number, nuice_ref, varian,
-     &                       ref_r0, igcm_dust_submicron
+     &                       ref_r0, igcm_dust_submicron,
+     &                       igcm_stormdust_mass,igcm_stormdust_number
        USE dimradmars_mod, only: nueffdust,naerkind,
      &            name_iaer,
      &            iaer_dust_conrath,iaer_dust_doubleq,
-     &            iaer_dust_submicron,iaer_h2o_ice
+     &            iaer_dust_submicron,iaer_h2o_ice,
+     &            iaer_stormdust_doubleq 
 
        IMPLICIT NONE
@@ -45,28 +47,24 @@
 
 c-----------------------------------------------------------------------
-c     Inputs:
+c     Inputs/outputs:
 c     ------
 
-      INTEGER ngrid,nlayer
+      INTEGER, INTENT(in) :: ngrid,nlayer
 c     Ice geometric mean radius (m)
-      REAL :: rice(ngrid,nlayer)
+      REAL, INTENT(out) :: rice(ngrid,nlayer)
 c     Estimated effective variance of the size distribution (n.u.)
-      REAL :: nuice(ngrid,nlayer)
+      REAL, INTENT(out) :: nuice(ngrid,nlayer)
 c     Tracer mass mixing ratio (kg/kg)
-      REAL pq(ngrid,nlayer,nqmx)
-      REAL rdust(ngrid,nlayer) ! Dust geometric mean radius (m)
-      
-      REAL pplay(ngrid,nlayer) ! altitude at the middle of the layers
-      REAL tau(ngrid,naerkind)
-
-
-c     Outputs:
-c     -------
-
+      REAL, INTENT(in) :: pq(ngrid,nlayer,nqmx)
+      REAL, INTENT(out) :: rdust(ngrid,nlayer) ! Dust geometric mean radius (m)
+      REAL, INTENT(out) :: rstormdust(ngrid,nlayer) ! Dust geometric mean radius (m) 
+      REAL, INTENT(in) :: pplay(ngrid,nlayer) ! altitude at the middle of the layers
+      REAL, INTENT(in) :: tau(ngrid,naerkind)
 c     Aerosol effective radius used for radiative transfer (meter)
-      REAL :: reffrad(ngrid,nlayer,naerkind)
+      REAL, INTENT(out) :: reffrad(ngrid,nlayer,naerkind)
 c     Aerosol effective variance used for radiative transfer (n.u.)
-      REAL :: nueffrad(ngrid,nlayer,naerkind)
-
+      REAL, INTENT(out) :: nueffrad(ngrid,nlayer,naerkind)
+      REAL, INTENT(in) :: tauscaling(ngrid)         ! Convertion factor for qccn and Nccn
+      
 c     Local variables:
 c     ---------------
@@ -85,5 +83,4 @@
       REAL Mo,No                       ! Mass and number of ccn
       REAL rhocloud(ngrid,nlayer)  ! Cloud density (kg.m-3)
-      REAL tauscaling(ngrid)         ! Convertion factor for qccn and Nccn
 
       LOGICAL,SAVE :: firstcall=.true.
@@ -114,4 +111,15 @@
           ENDDO
         ENDIF
+
+        ! updating radius of stormdust particles
+        IF (rdstorm.AND.active) THEN
+          DO l=1,nlayer
+            DO ig=1, ngrid
+              call updaterdust(pq(ig,l,igcm_stormdust_mass),
+     &                 pq(ig,l,igcm_stormdust_number),rstormdust(ig,l))
+              nueffdust(ig,l) = exp(varian**2.)-1.
+             ENDDO
+           ENDDO
+        ENDIF
         
 c       1.2 Water-ice particles
@@ -126,6 +134,6 @@
 
           IF (firstcall) THEN
-            IF (minval(tauscaling).lt.0) tauscaling(:) = 1.e-3 ! default value when non-read in startfi is -1
-            IF (freedust)                tauscaling(:) = 1.    ! if freedust, enforce no rescaling at all
+            !IF (minval(tauscaling).lt.0) tauscaling(:) = 1.e-3 ! default value when non-read in startfi is -1
+            !IF (freedust)                tauscaling(:) = 1.    ! if freedust, enforce no rescaling at all
             firstcall = .false.
           ENDIF
@@ -206,4 +214,14 @@
           ENDDO
 c==================================================================
+        CASE("stormdust_doubleq") aerkind! Two-moment scheme for
+c       stormdust; same distribution than normal dust
+c==================================================================
+          DO l=1,nlayer
+            DO ig=1,ngrid
+              reffrad(ig,l,iaer) = rstormdust(ig,l) * ref_r0
+              nueffrad(ig,l,iaer) = nueffdust(ig,l)
+            ENDDO
+          ENDDO
+c==================================================================
         END SELECT aerkind
       ENDDO ! iaer (loop on aerosol kind)
Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 1973)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F	(revision 1974)
@@ -12,12 +12,14 @@
      $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
      $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
-     $                hfmax,sensibFlux)
+     $                hfmax,sensibFlux,dustliftday,local_time)
 
       use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
      &                      igcm_dust_submicron, igcm_h2o_vap,
-     &                      igcm_h2o_ice, alpha_lift
+     &                      igcm_h2o_ice, alpha_lift,
+     &                      igcm_stormdust_mass, igcm_stormdust_number
       use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
       USE comcstfi_h, ONLY: cpp, r, rcp, g
       use turb_mod, only: turb_resolved, ustar, tstar
+      use compute_dtau_mod, only: ti_injection,tf_injection
       IMPLICIT NONE
 
@@ -79,5 +81,7 @@
       real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
       real,intent(out) :: pdqdif(ngrid,nlay,nq) 
-      real,intent(out) :: pdqsdif(ngrid,nq) 
+      real,intent(out) :: pdqsdif(ngrid,nq)
+      REAL,INTENT(in) :: dustliftday(ngrid)
+      REAL,INTENT(in) :: local_time(ngrid)
       
 c   local:
@@ -691,5 +695,7 @@
              end do
            else if (doubleq) then
-             do ig=1,ngrid
+             if (dustinjection.eq.0) then !injection scheme 0 (old) 
+             				  !or 2 (injection in CL)
+              do ig=1,ngrid
                if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2
                  pdqsdif(ig,igcm_dust_mass) =
@@ -698,5 +704,40 @@
      &             -alpha_lift(igcm_dust_number) 
                end if 
-             end do
+              end do
+             elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface
+              do ig=1,ngrid
+               if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2
+                IF((ti_injection.LE.local_time(ig)).and.
+     &            	(local_time(ig).LE.tf_injection)) THEN
+     		  if (rdstorm) then !Rocket dust storm scheme
+             	    	pdqsdif(ig,igcm_stormdust_mass) = 
+     &             		-alpha_lift(igcm_stormdust_mass)
+     &				*dustliftday(ig)
+               		pdqsdif(ig,igcm_stormdust_number) = 
+     &             		-alpha_lift(igcm_stormdust_number)
+     &				*dustliftday(ig)
+                 	pdqsdif(ig,igcm_dust_mass)= 0.
+                        pdqsdif(ig,igcm_dust_number)= 0.
+                  else
+                 	pdqsdif(ig,igcm_dust_mass)=
+     &            		-dustliftday(ig)*
+     &            		alpha_lift(igcm_dust_mass)               
+                 	pdqsdif(ig,igcm_dust_number)= 
+     &             		-dustliftday(ig)*
+     &            		alpha_lift(igcm_dust_number)
+		  endif
+             	  if (submicron) then
+                   	pdqsdif(ig,igcm_dust_submicron) = 0.
+                  endif ! if (submicron)
+     		ELSE ! outside dust injection time frame
+                 	pdqsdif(ig,igcm_dust_mass)= 0.
+                        pdqsdif(ig,igcm_dust_number)= 0.
+                 	pdqsdif(ig,igcm_stormdust_mass)= 0.
+                        pdqsdif(ig,igcm_stormdust_number)= 0.
+     		ENDIF
+               
+                end if ! of if(co2ice(ig).lt.1)
+              end do
+             endif ! end if dustinjection
            else if (submicron) then
              do ig=1,ngrid
@@ -723,5 +764,5 @@
 c      Inversion pour l'implicite sur q 
 c       --------------------------------
-        do iq=1,nq
+        do iq=1,nq  !for all tracers including stormdust
           CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
 
