Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F	(revision 2077)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F	(revision 2079)
@@ -2,5 +2,5 @@
       SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0,
      &                    zmea,zstd,zsig,zgam,zthe,
-     &                    hmons,summit,zavg)
+     &                    hmons,summit,base,zavg)
 c=======================================================================
 c
@@ -79,6 +79,7 @@
       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
+      REAL,intent(out) :: summit(imdp1*jmdp1)
+      REAL,intent(out) :: base(imdp1*jmdp1)
+      REAL,intent(out) :: zavg(imdp1*jmdp1)
       
 ! Local variables:
@@ -327,7 +328,7 @@
 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)
+      string(7) = 'base'   !base of summit (not subgrid)
 c     string(8) = 'bottom' !subgrid base
-      do k=5,6
+      do k=5,7
           write(*,*) 'string',k,string(k)  
 c-----------------------------------------------------------------------
@@ -373,5 +374,5 @@
           enddo
 
-          if (k .eq. 5 .or. k .eq. 6) then 
+          if (k .eq. 5 .or. k .eq. 6 .or. k .eq. 7) then 
               ! hmons, summit, keep the maximum of subgrids 
               call mvc_horiz(imd,jmdp1,iim,jjm,longitude,
@@ -404,5 +405,19 @@
           if (k.eq.6) then  ! summit
              do i=1,iimp1*jjp1
+                  if (pfield(i) .ne. -999999.) then             
                     summit(i) = pfield(i)
+                  else
+                    summit(i)=0.
+                  endif
+             enddo
+          endif
+
+          if (k.eq.7) then  ! base
+             do i=1,iimp1*jjp1
+                  if (pfield(i) .ne. -999999.) then             
+                    base(i) = pfield(i)
+                  else
+                    base(i)=0.
+                  endif
              enddo
           endif
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 2077)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 2079)
@@ -27,5 +27,5 @@
       use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
      &                     albedodat, z0_default, qsurf, tsurf,
-     &                     co2ice, emis, hmons
+     &                     co2ice, emis, hmons, summit, base
       use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
       use control_mod, only: day_step, iphysiq, anneeref, planet_type
@@ -100,4 +100,5 @@
       real hmonsS(iip1,jjp1)
       real summitS(iip1,jjp1)
+      real baseS(iip1,jjp1)
       real zavgS(iip1,jjp1)
       real z0S(iip1,jjp1)
@@ -385,5 +386,5 @@
       CALL datareadnc(relief,phis,alb,surfith,z0S,
      &          zmeaS,zstdS,zsigS,zgamS,ztheS,
-     &          hmonsS,summitS,zavgS)
+     &          hmonsS,summitS,baseS,zavgS)
 
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
@@ -398,5 +399,6 @@
       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,summitS,summit)
+      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,baseS,base)
 !      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zavgS,zavg)
 
@@ -1651,5 +1653,5 @@
      &              nqtot,dtphys,real(day_ini),0.0,cell_area,
      &              albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
-     &              hmons)
+     &              hmons,summit,base)
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      &              dtphys,hour_ini,
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2077)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2079)
@@ -11,5 +11,5 @@
      &                     albedice, iceradius, dtemisice, z0,
      &                     zmea, zstd, zsig, zgam, zthe, phisfi,
-     &                     watercaptag, hmons
+     &                     watercaptag, hmons, summit, base
       use slope_mod, only: theta_sl, psi_sl
       use phyredem, only: physdem0,physdem1
@@ -556,9 +556,9 @@
       zsig(1)=0.E+0
       zgam(1)=0.E+0
-      zthe(1)=0.E+0
-c  for rocket dust storm scheme      
+      zthe(1)=0.E+0      
       hmons(1)=0.E+0
-
-
+      summit(1)=0.E+0
+      base(1)=0.E+0
+      
 c   Specific initializations for "physiq"
 c   -------------------------------------
@@ -718,5 +718,6 @@
       call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
      &              nq,dtphys,float(day0),time,cell_area,
-     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,hmons)
+     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
+     &              hmons,summit,base)
       call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
      &              dtphys,time,
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2077)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2079)
@@ -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, hmons
+                       zmea, zstd, zsig, zgam, zthe, hmons, summit, base
   use iostart, only: nid_start, open_startphy, close_startphy, &
                      get_field, get_var, inquire_field, &
@@ -199,5 +199,43 @@
              minval(hmons), maxval(hmons)
 endif
-      
+
+! summit 
+call get_field("summit",summit,found)
+if (.not.found) then
+  write(*,*) "WARNING: phyetat0: Failed loading <summit>"
+  if (rdstorm) then
+    call abort
+  else
+    write(*,*) "will continue anyway..."
+    write(*,*) "because you may not need it."
+    summit(:)=0.
+  end if
+else
+  do ig=1,ngrid
+    if (summit(ig) .eq. -999999.)  summit(ig)=0.
+  enddo
+  write(*,*) "phyetat0: <summit> range:", &
+             minval(summit), maxval(summit)
+endif
+
+! summit 
+call get_field("base",base,found)
+if (.not.found) then
+  write(*,*) "WARNING: phyetat0: Failed loading <base>"
+  if (rdstorm) then
+    call abort
+  else
+    write(*,*) "will continue anyway..."
+    write(*,*) "because you may not need it."
+    base(:)=0.
+  end if
+else
+  do ig=1,ngrid
+    if (base(ig) .eq. -999999.)  base(ig)=0.
+  enddo
+  write(*,*) "phyetat0: <base> range:", &
+             minval(base), maxval(base)
+endif
+ 
 ! Time axis
 ! obtain timestart from run.def
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2077)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2079)
@@ -8,5 +8,5 @@
                          phystep,day_ini,time,airefi, &
                          alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe, &
-                         phmons)
+                         phmons,psummit,pbase)
 ! create physics restart file and write time-independent variables
   use comsoil_h, only: inertiedat, volcapa, mlayer
@@ -15,5 +15,5 @@
                        z0_default, albedice, emisice, emissiv, &
                        iceradius, dtemisice, phisfi, z0,   &
-                       hmons
+                       hmons,summit,base
   use dimradmars_mod, only: tauvis
   use iostart, only : open_restartphy, close_restartphy, & 
@@ -46,4 +46,6 @@
   real,intent(in) :: pzthe(ngrid)
   real,intent(in) :: phmons(ngrid)
+  real,intent(in) :: psummit(ngrid)
+  real,intent(in) :: pbase(ngrid)
 
   real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
@@ -131,8 +133,7 @@
   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
+  call put_field("hmons","Relief: hmons parameter (summit - base)",hmons)
+  call put_field("summit","Relief: altitude from the aeroid of the summit of the highest subgrid topography",summit)
+  call put_field("base","Relief: altitude from the aeroid of the base of the highest subgrid topography",base)
      
   ! Soil thermal inertia
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2077)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2079)
@@ -42,5 +42,5 @@
      &                     tsurf, co2ice, emis,
      &                     capcal, fluxgrd, qsurf,
-     &                     hmons
+     &                     hmons,summit,base
       use comsaison_h, only: dist_sol, declin, mu0, fract, local_time
       use slope_mod, only: theta_sl, psi_sl
@@ -579,5 +579,5 @@
      &                 ptimestep,pday,time_phys,cell_area,
      &                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
-     &                 hmons)
+     &                 hmons,summit,base)
          endif
 #endif
Index: trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2077)
+++ trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2079)
@@ -4,5 +4,4 @@
 
       REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1)
-      REAL, SAVE, ALLOCATABLE :: alpha_hmons(:) ! slope winds lifting mesh fraction
       
       CONTAINS
@@ -11,9 +10,9 @@
 ! 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
+! calculation of the vertical flux
+! call of vl_storm : Van Leer transport scheme of the dust tracers
+! detrainement of stormdust in the background dust
 ! -----------------------------------------------------------------------
-! Authors: C. Wang; F. Forget; T. Bertrand 
+! Authors: M. Vals; C. Wang; F. Forget; T. Bertrand 
 ! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
 ! -----------------------------------------------------------------------
@@ -27,18 +26,16 @@
 !             input sub-grid scale cloud
                                  clearsky,totcloudfrac,                &
-!                   output
-                                 pdqrds,wspeed,dsodust,dsords)
-
-      use tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
+!             output
+                                 pdqrds,wrad,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
+      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 callradite_mod
+      IMPLICIT NONE
 
 !--------------------------------------------------------
@@ -66,11 +63,12 @@
 
 !     input for second radiative transfer 
-      logical, INTENT(IN) :: clearatm
+      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)
+      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
@@ -82,5 +80,5 @@
 
       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) :: wrad(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
@@ -90,5 +88,5 @@
 !--------------------------------------------------------
       INTEGER l,ig,tsub,iq,ll
-! chao local variables from callradite.F       
+!     local variables from callradite.F       
       REAL zdtlw1(ngrid,nlayer)    ! (K/s) storm
       REAL zdtsw1(ngrid,nlayer)    ! (K/s) storm
@@ -97,68 +95,49 @@
       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 zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for stormdust
+      REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer)   ! rad. heating rate at intermediate levels l+1/2 for background dust
+
+      REAL zq_stormdust_mass(ngrid,nlayer) ! intermediate tracer stormdust mass 
+      REAL zq_stormdust_number(ngrid,nlayer) ! intermediate tracer stormdust number
+      REAL zq_dust_mass(ngrid,nlayer)           ! intermediate tracer dust mass 
+      REAL zq_dust_number(ngrid,nlayer)         ! intermediate tracer dust number
+
+      REAL mr_stormdust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust mass)
+      REAL mr_stormdust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust number)
+      REAL mr_dust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (dust mass)
+      REAL mr_dust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (sdust number)
+
+      REAL zq_vl_col(nlayer) ! column intermediate tracer used by Van Leer (number)
+      REAL zn_vl_col(nlayer) ! column intermediate tracer used by Van Leer (mass)
+                   
+      REAL dqvl_stormdust_mass(ngrid,nlayer)    ! tendancy of vertical transport (stormdust mass)
+      REAL dqvl_stormdust_number(ngrid,nlayer)  ! tendancy of vertical transport (stormdust number)
+      REAL dqvl_dust_mass(ngrid,nlayer)    ! tendancy of vertical transport (dust mass)
+      REAL dqvl_dust_number(ngrid,nlayer)  ! tendancy of vertical transport (dust number)
+      REAL dqdet_stormdust_mass(ngrid,nlayer)   ! tendancy of detrainement (stormdust mass)
+      REAL dqdet_stormdust_number(ngrid,nlayer) ! tendancy of detrainement (stormdust number)
+
+      REAL masse(nlayer)     ! mass of atmosphere (kg/m2)
       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)
+      REAL w(nlayer)          ! air mass flux (calculated with the vertical wind velocity profile) used as input in Van Leer (kgair/m2)
+      REAL wqmass(nlayer+1)   ! tracer (dust_mass) mass flux in Van Leer (kg/m2)
+      REAL wqnumber(nlayer+1) ! tracer (dust_number) mass flux in Van Leer (kg/m2)
+
+      LOGICAL storm(ngrid)    ! true when there is a dust storm (if the opacity is high): trigger the rocket dust storm scheme
+      REAL coefdetrain        ! coefficient for detrainment : % of stormdust detrained
+      INTEGER scheme(ngrid)   ! triggered scheme
+            
+      REAL,PARAMETER:: coefmin =0.025 ! 0<coefmin<1 Minimum fraction of stormdust detrained  
+      REAL,PARAMETER:: wmin =0.1!0.25 ! stormdust detrainment if wrad < wmin  
+      REAL,PARAMETER:: wmax =10.   ! maximum vertical velocity of the rocket dust storms (m/s)
+      
+!     diagnostics 
+      REAL lapserate(ngrid,nlayer) 
+      REAL deltahr(ngrid,nlayer+1)
     
-      logical,save :: firstcall=.true.
-      real alpha(ngrid) ! scale of the vertical motion (applicable for
-                        ! rds and slp
-      ! variables for radiative transfer
+      LOGICAL,SAVE :: firstcall=.true.
+
+!     variables for the radiative transfer
       REAL  fluxsurf_lw1(ngrid)
       REAL  fluxsurf_sw1(ngrid,2)
@@ -175,115 +154,43 @@
       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
+
+      ! **********************************************************************
+      ! **********************************************************************
+      ! Rocket dust storm parametrization to reproduce the detached dust layers
+      ! during the dust storm season:
       !     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
+      !!      case 1 storm = false: nothing to do
+      !!      case 2 rocket dust storm (storm=true)
+      !!    3. Vertical transport (Van Leer)
       !!    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
+
+      
+      ! **********************************************************************
+      ! Initializations
       ! **********************************************************************
       storm(:)=.false.
       pdqrds(:,:,:) = 0.
-      zdqdetstorm_mass(:,:)=0.
-      zdqdetstorm_number(:,:)=0.
-      wspeed(:,:)=0.
-      detrainment(:,:)=0.
-      zqstorm_mass_col(:)=0.
-      zqstorm_number_col(:)=0.
+      mr_dust_mass(:,:)=0.
+      mr_dust_number(:,:)=0.
+      mr_stormdust_mass(:,:)=0.
+      mr_stormdust_number(:,:)=0.
+      dqvl_dust_mass(:,:)=0.
+      dqvl_dust_number(:,:)=0.
+      dqvl_stormdust_mass(:,:)=0.
+      dqvl_stormdust_number(:,:)=0.
+      dqdet_stormdust_mass(:,:)=0.
+      dqdet_stormdust_number(:,:)=0.
+      wrad(:,:)=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 
@@ -292,26 +199,23 @@
       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
+
+      zq_dust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_mass)
+      zq_dust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_number)
+      zq_stormdust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_mass)
+      zq_stormdust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_number)
+
+      ! *********************************************************************
+      ! 0. Check if there is a 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 
+        DO l=1,nlayer
+          IF (zq(ig,l,igcm_stormdust_mass) &
+          .gt. zq(ig,l,igcm_dust_mass)*(1.E-4)) THEN
             storm(ig)=.true.
-            exit
-          endif
-        enddo       
-      enddo
+            EXIT
+          ENDIF
+        ENDDO     
+      ENDDO
       
       ! *********************************************************************
@@ -328,457 +232,260 @@
       ! **********************************************************************
       ! 2. Compute vertical velocity for storm dust
-      ! **********************************************************************     
-      DO ig=1,ngrid
-        !! **********************************************************************
-        !! 2.1 case 1: Nothing to do when no storm and no slope, or
-        !!             no storm and not daytime
-        IF (((alpha_hmons(ig) .EQ. 0.) .AND. .NOT.(storm(ig))) &
-             .OR. ((mu0(ig) .LE. mu0lim) .AND. .NOT.(storm(ig))) ) then
-          scheme(ig)=1
-          cycle
-        endif
+      ! **********************************************************************
+        !! **********************************************************************
+        !! 2.1 Nothing to do when no storm
+        !!             no storm
+        DO ig=1,ngrid        
+          IF (.NOT.(storm(ig))) then
+            scheme(ig)=1
+            cycle
+          ENDIF ! IF (storm(ig))
+        ENDDO ! DO ig=1,ngrid                 
         
+        !! **********************************************************************
+        !! 2.2 Calculation of the extra heating : computing heating rates
+        !! gradient at boundaries of each layer, start from surface
+        DO ig=1,ngrid
+          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))
+            !! 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-1
+        ENDDO ! DO ig=1,ngrid     
+
+        !! **********************************************************************
+        !! 2.3 Calculation of the vertical velocity : extra heating
+        !! balanced by adiabatic cooling 
+        DO ig=1,ngrid
+          IF (storm(ig)) THEN       
+            scheme(ig)=2         
+            DO l=1,nlayer
+              deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l))  &
+                                          -(zdtlw_lev(l)+zdtsw_lev(l))
+              wrad(ig,l)=-deltahr(ig,l)/(g/cpp+   &
+                                         max(zdtvert(l),-0.99*g/cpp))       
+              !! Limit vertical wind in case of lapse rate close to adiabatic
+              wrad(ig,l)=max(wrad(ig,l),-wmax) 
+              wrad(ig,l)=min(wrad(ig,l),wmax)
+            ENDDO
+          ENDIF ! IF (storm(ig))
+        ENDDO ! DO ig=1,ngrid
+
+      ! **********************************************************************
+      ! 3. Vertical transport
+      ! **********************************************************************
+        !! **********************************************************************
+        !! 3.1 Transport of background dust + storm dust (concentrated)
+        DO ig=1,ngrid
+          IF (storm(ig)) THEN 
+            DO l=1,nlayer
+              mr_dust_mass(ig,l) = zq_dust_mass(ig,l)
+              mr_dust_number(ig,l) = zq_dust_number(ig,l)
+              mr_stormdust_mass(ig,l) = zq_dust_mass(ig,l)+ &
+                                      zq_stormdust_mass(ig,l)/totstormfract(ig)
+              mr_stormdust_number(ig,l) = zq_dust_number(ig,l)+ &
+                                      zq_stormdust_number(ig,l)/totstormfract(ig)
+            ENDDO ! DO l=1,nlayer
+          ENDIF ! IF (storm(ig))
+        ENDDO ! DO ig=1,ngrid
+
+        !! **********************************************************************
+        !! 3.2 Vertical transport by a Van Leer scheme
+        DO ig=1,ngrid
+          IF (storm(ig)) THEN
+           
+            !! Mass of atmosphere in the layer
+            DO l=1,nlayer
+               masse(l)=(pplev(ig,l)-pplev(ig,l+1))/g
+            ENDDO
+            
+            !! Mass flux in kg/m2
+            DO l=1,nlayer
+               w(l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(l)))*ptimestep
+            ENDDO
+            
+            !! Vector column
+            DO l=1,nlayer
+               zq_vl_col(l)= mr_stormdust_mass(ig,l)
+               zn_vl_col(l)= mr_stormdust_number(ig,l)
+            ENDDO
+            
+            !! Van Leer scheme
+            wqmass(:)=0.
+            wqnumber(:)=0.
+            CALL vl_storm(nlayer,zq_vl_col,2.,   &
+                  masse,w,wqmass)
+            CALL vl_storm(nlayer,zn_vl_col,2.,  &
+                  masse,w,wqnumber)
+            !! Mass mixing ratio after transport      
+            mr_stormdust_mass(ig,:) = zq_vl_col(:)
+            mr_stormdust_number(ig,:) = zn_vl_col(:)
+                        
+          ENDIF ! IF storm(ig)
+        ENDDO ! DO ig=1,ngrid  
+
+        !! **********************************************************************
+        !! 3.3 Re-calculation of the mixing ratios after vertical transport
+        DO ig=1,ngrid
+         IF (storm(ig)) THEN
+           DO l=1,nlayer
+           
+             !! General and "healthy" case
+             IF (mr_stormdust_mass(ig,l).ge.mr_dust_mass(ig,l)) THEN
+               zq_dust_mass(ig,l) = mr_dust_mass(ig,l)
+               zq_dust_number(ig,l) = mr_dust_number(ig,l)
+               zq_stormdust_mass(ig,l) = totstormfract(ig)*(mr_stormdust_mass(ig,l)-mr_dust_mass(ig,l))
+               zq_stormdust_number(ig,l) = totstormfract(ig)*(mr_stormdust_number(ig,l)-mr_dust_number(ig,l))
+             !! Particular case: there is less than initial dust mixing ratio after the vertical transport
+             ELSE
+               zq_dust_mass(ig,l) = (1.-totstormfract(ig))*mr_dust_mass(ig,l)+totstormfract(ig)*mr_stormdust_mass(ig,l)
+               zq_dust_number(ig,l) = (1.-totstormfract(ig))*mr_dust_number(ig,l)+totstormfract(ig)*mr_stormdust_number(ig,l)
+               zq_stormdust_mass(ig,l) = 0.
+               zq_stormdust_number(ig,l) = 0.
+             ENDIF
+             
+           ENDDO ! DO l=1,nlayer            
+         ENDIF ! IF storm(ig)
+        ENDDO ! DO ig=1,ngrid
+
+        !! **********************************************************************
+        !! 3.4 Calculation of the tendencies of the vertical transport
+        DO ig=1,ngrid
+         IF (storm(ig)) THEN
+           DO l=1,nlayer
+             dqvl_stormdust_mass(ig,l) = (zq_stormdust_mass(ig,l)- &
+                                   zq(ig,l,igcm_stormdust_mass)) /ptimestep
+             dqvl_stormdust_number(ig,l) = (zq_stormdust_number(ig,l)- &
+                                     zq(ig,l,igcm_stormdust_number)) /ptimestep
+             dqvl_dust_mass(ig,l) = (zq_dust_mass(ig,l)-zq(ig,l,igcm_dust_mass)) /ptimestep
+             dqvl_dust_number(ig,l) = (zq_dust_number(ig,l)-zq(ig,l,igcm_dust_number)) /ptimestep
+           ENDDO
+         ENDIF ! IF storm(ig)
+        ENDDO ! DO ig=1,ngrid           
+
+      ! **********************************************************************
+      ! 4. Detrainment: stormdust is converted to background dust
+      ! **********************************************************************
+        !! **********************************************************************
+        !! 4.1 Compute the coefficient of detrainmen
+        DO ig=1,ngrid
+          DO l=1,nlayer
+            IF ((max(abs(wrad(ig,l)),abs(wrad(ig,l+1))) .lt.  & 
+                          wmin) .or. (zq_dust_mass(ig,l) .gt.  &
+                                     10000.*zq_stormdust_mass(ig,l))) THEN
+               coefdetrain=1.
+            ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))   &
+                                                       .le. wmax) THEN
+               coefdetrain=1.*(((1-coefmin)/(wmin-wmax)**2)*     &
+                   (max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))-wmax)**2 &
+                                                               +coefmin)
+            ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1))).gt. wmax ) THEN
+               coefdetrain=coefmin
+            ELSE
+               coefdetrain=coefmin
+            ENDIF
+          ENDDO ! DO l=1,nlayer
+        ENDDO ! DO ig=1,ngrid
         
-        ! 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
+        !! **********************************************************************
+        !! 4.2 Calculation of the tendencies of the detrainment
+        DO ig=1,ngrid
+          DO l=1,nlayer
+            dqdet_stormdust_mass(ig,l)=-coefdetrain*zq_stormdust_mass(ig,l) &
+                                                        /ptimestep
+            dqdet_stormdust_number(ig,l)=-coefdetrain*zq_stormdust_number(ig,l) &
+                                                        /ptimestep
+          ENDDO ! DO l=1,nlayer
+        ENDDO ! DO ig=1,ngrid
+           
+      ! **********************************************************************
+      ! 5. Final tendencies: vertical transport + detrainment
+      ! **********************************************************************
+        DO ig=1,ngrid
+          DO l=1,nlayer      
+          pdqrds(ig,l,igcm_stormdust_mass)=dqdet_stormdust_mass(ig,l) &
+                                                 +dqvl_stormdust_mass(ig,l)
+          pdqrds(ig,l,igcm_stormdust_number)=dqdet_stormdust_number(ig,l) &
+                                                 +dqvl_stormdust_number(ig,l)
+          pdqrds(ig,l,igcm_dust_mass)= -dqdet_stormdust_mass(ig,l) &
+                                       +dqvl_dust_mass(ig,l) 
+          pdqrds(ig,l,igcm_dust_number)= -dqdet_stormdust_number(ig,l) &
+                                       +dqvl_dust_number(ig,l) 
+          ENDDO ! DO l=1,nlayer
+        ENDDO ! DO ig=1,ngrid
+
+      ! **********************************************************************
+      ! 6. To prevent from negative values
+      ! **********************************************************************
+        DO ig=1,ngrid
+          DO l=1,nlayer
+            IF ((pq(ig,l,igcm_stormdust_mass) &
+               +pdqrds(ig,l,igcm_stormdust_mass)*ptimestep .le. 0.) .or. &
+              (pq(ig,l,igcm_stormdust_number) &
+               +pdqrds(ig,l,igcm_stormdust_number)*ptimestep .le. 0.)) THEN
+               pdqrds(ig,l,igcm_stormdust_mass)=-pq(ig,l,igcm_stormdust_mass)/ptimestep
+               pdqrds(ig,l,igcm_stormdust_number)=-pq(ig,l,igcm_stormdust_number)/ptimestep
+            ENDIF
+          ENDDO ! nlayer
+        ENDDO ! DO ig=1,ngrid
+
+        DO ig=1,ngrid
+          DO l=1,nlayer           
+            IF ((pq(ig,l,igcm_dust_mass) &
+               +pdqrds(ig,l,igcm_dust_mass)*ptimestep .le. 0.) .or. &
+              (pq(ig,l,igcm_dust_number) &
+               +pdqrds(ig,l,igcm_dust_number)*ptimestep .le. 0.)) THEN
+               pdqrds(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)/ptimestep
+               pdqrds(ig,l,igcm_dust_number)=-pq(ig,l,igcm_dust_number)/ptimestep
+            ENDIF
+          ENDDO ! nlayer
+        ENDDO ! DO ig=1,ngrid
         
-          !! **********************************************************************
-          !! 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)
+!=======================================================================
+! WRITEDIAGFI
+
         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
 
-!********************************************************************************
-!********************************************************************************
+!=======================================================================
+! VAN LEER
+
       SUBROUTINE vl_storm(nlay,q,pente_max,masse,w,wq)
 !
@@ -801,6 +508,6 @@
 !   Arguments:
 !   ----------
-      integer,intent(in) :: nlay ! number of atmospheric layers
-      real masse(nlay),pente_max
+      INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
+      REAL masse(nlay),pente_max
       REAL q(nlay)
       REAL w(nlay)
@@ -813,11 +520,11 @@
 !
 
-      real dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
-      real newmasse
-      real sigw, Mtot, MQtot
-      integer m
+      REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
+      REAL newmasse
+      REAL sigw, Mtot, MQtot
+      INTEGER m
 
       REAL      SSUM,CVMGP,CVMGT
-      integer ismax,ismin
+      INTEGER ismax,ismin
 
 
@@ -847,7 +554,5 @@
       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  .......
@@ -891,16 +596,9 @@
              if (m.gt.0) then
                 sigw=(w(l+1)+Mtot)/masse(m)
-                wq(l+1)= (MQtot + (-w(l+1)-Mtot)*         &
+                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
@@ -913,8 +611,4 @@
 
        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)     
@@ -967,105 +661,25 @@
 
        enddo
-      end subroutine vl_storm
-!********************************************************************************
-      SUBROUTINE intep_vtemp(nlayer,hm,temp,zlay,envtemp,slpb)
-
-       USE comcstfi_h, only: g,cpp
-
+       
+      END SUBROUTINE vl_storm
+
+!=======================================================================
+! Initialization of the module variables
+
+       subroutine ini_rocketduststorm_mod(ngrid)
+       
        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)
+       
+       integer, intent(in) :: ngrid
+       
+       allocate(dustliftday(ngrid))
+       
+       end subroutine ini_rocketduststorm_mod
+       
+       subroutine end_rocketduststorm_mod
        
        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       
Index: trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 2077)
+++ trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 2079)
@@ -20,5 +20,5 @@
   real,save :: iceradius(2) , dtemisice(2)
   real,save,allocatable :: zmea(:),zstd(:),zsig(:),zgam(:),zthe(:)
-  real,save,allocatable :: hmons(:)
+  real,save,allocatable :: hmons(:),summit(:),base(:)
   real,save,allocatable :: z0(:) ! surface roughness length (m)
   real,save :: z0_default ! default (constant over planet) surface roughness (m)
@@ -57,4 +57,6 @@
     allocate(fluxgrd(ngrid))
     allocate(hmons(ngrid))
+    allocate(summit(ngrid))
+    allocate(base(ngrid))
     
   end subroutine ini_surfdat_h
@@ -82,5 +84,7 @@
     if (allocated(fluxgrd))     deallocate(fluxgrd)
     if (allocated(hmons))       deallocate(hmons)
-
+    if (allocated(summit))      deallocate(summit)
+    if (allocated(base))        deallocate(base)
+    
   end subroutine end_surfdat_h
 
