Index: trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2199)
+++ trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2201)
@@ -11,5 +11,5 @@
 !=======================================================================
 ! calculation of the vertical flux
-! call of vl_storm : Van Leer transport scheme of the dust tracers
+! call of van_leer : Van Leer transport scheme of the dust tracers
 ! detrainement of stormdust in the background dust
 ! -----------------------------------------------------------------------
@@ -97,5 +97,5 @@
 ! Local variables 
 !--------------------------------------------------------
-      INTEGER l,ig,tsub,iq,ll
+      INTEGER l,ig,iq,ll
 !     local variables from callradite.F       
       REAL zdtlw1(ngrid,nlayer)    ! (K/s) storm
@@ -117,7 +117,4 @@
       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)
@@ -140,7 +137,14 @@
             
       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:: wmin =0.25 ! stormdust detrainment if wrad < wmin  
       REAL,PARAMETER:: wmax =10.   ! maximum vertical velocity of the rocket dust storms (m/s)
-      
+
+!     subtimestep
+      INTEGER tsub
+      INTEGER nsubtimestep    !number of subtimestep when calling van_leer 
+      REAL subtimestep        !ptimestep/nsubtimestep
+      REAL dtmax              !considered time needed for dust to cross one layer
+      REAL,PARAMETER:: secu=3.!3.      !coefficient on wspeed to avoid dust crossing many layers during subtimestep
+
 !     diagnostics 
       REAL lapserate(ngrid,nlayer) 
@@ -265,10 +269,4 @@
 
             scheme(ig)=2
-            !! This is the env. lapse rate 
-            zdtvert(ig,1)=0.
-            DO l=1,nlayer-1
-              zdtvert(ig,l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
-              lapserate(ig,l+1)=zdtvert(ig,l+1)
-            ENDDO
       
             !! computing heating rates gradient at boundraies of each layer
@@ -301,5 +299,12 @@
                            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 l=1,nlayer-1
+
+            !! This is the env. lapse rate 
+            zdtvert(ig,1)=0.
+            DO l=1,nlayer-1
+              zdtvert(ig,l+1)=(ztlev(ig,l+1)-ztlev(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
+              lapserate(ig,l+1)=zdtvert(ig,l+1)
+            ENDDO
 
             !! **********************************************************************
@@ -338,39 +343,52 @@
         ENDDO ! DO ig=1,ngrid
 
-        !! **********************************************************************
-        !! 3.2 Vertical transport by a Van Leer scheme
         DO ig=1,ngrid
           IF (storm(ig)) THEN
-           
+        !! **********************************************************************
+        !! 3.2 Compute the subtimestep to conserve the mass in the Van Leer transport
+            dtmax=ptimestep
+            DO l=2,nlayer
+              IF (wrad(ig,l).lt.0.) THEN
+                 dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/  &
+                                   (secu*abs(wrad(ig,l))))
+              ELSE IF (wrad(ig,l).gt.0.) then
+                 dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/  &
+                                   (secu*abs(wrad(ig,l))))
+              ENDIF
+            ENDDO
+            nsubtimestep= int(ptimestep/dtmax) 
+            subtimestep=ptimestep/float(nsubtimestep)
+            !! Mass flux generated by wup in kg/m2
+            DO l=1,nlayer
+               w(ig,l)=wrad(ig,l)*pplev(ig,l)/(r*ztlev(ig,l)) &
+                        *subtimestep
+            ENDDO ! l=1,nlayer
+
+        !! **********************************************************************
+        !! 3.3 Vertical transport by a Van Leer scheme
             !! Mass of atmosphere in the layer
             DO l=1,nlayer
                masse_col(l)=(pplev(ig,l)-pplev(ig,l+1))/g
             ENDDO
-            
-            !! Mass flux in kg/m2
-            DO l=1,nlayer
-               w(ig,l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(ig,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
-            CALL vl_storm(nlayer,zq_vl_col,2.,   &
-                  masse_col,w(ig,:),wqmass(ig,:))
-            CALL vl_storm(nlayer,zn_vl_col,2.,  &
-                  masse_col,w(ig,:),wqnumber(ig,:))
-            !! Mass mixing ratio after transport      
-            mr_stormdust_mass(ig,:) = zq_vl_col(:)
-            mr_stormdust_number(ig,:) = zn_vl_col(:)
-                        
+            !! Mass flux in kg/m2 if you are not using the subtimestep
+            !DO l=1,nlayer
+            !   w(ig,l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(ig,l)))*ptimestep
+            !ENDDO
+            !! Loop over the subtimestep
+            DO tsub=1,nsubtimestep
+              !! Van Leer scheme
+              wqmass(ig,:)=0.
+              wqnumber(ig,:)=0.
+              CALL van_leer(nlayer,mr_stormdust_mass(ig,:),2.,   &
+                    masse_col,w(ig,:),wqmass(ig,:))
+              CALL van_leer(nlayer,mr_stormdust_number(ig,:),2.,  &
+                    masse_col,w(ig,:),wqnumber(ig,:))
+            ENDDO !tsub=...           
+             
           ENDIF ! IF storm(ig)
         ENDDO ! DO ig=1,ngrid  
 
         !! **********************************************************************
-        !! 3.3 Re-calculation of the mixing ratios after vertical transport
+        !! 3.4 Re-calculation of the mixing ratios after vertical transport
         DO ig=1,ngrid
          IF (storm(ig)) THEN
@@ -396,5 +414,5 @@
 
         !! **********************************************************************
-        !! 3.4 Calculation of the tendencies of the vertical transport
+        !! 3.5 Calculation of the tendencies of the vertical transport
         DO ig=1,ngrid
          IF (storm(ig)) THEN
@@ -462,30 +480,30 @@
         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
+!      ! **********************************************************************
+!      ! 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
         
 !=======================================================================
@@ -501,39 +519,28 @@
         END SUBROUTINE rocketduststorm
 
-!=======================================================================
-! VAN LEER
-
-      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
-!
-!
-!   --------------------------------------------------------------------
+!=======================================================================  
+! **********************************************************************
+! Subroutine to determine the vertical transport with 
+! a Van Leer advection scheme (copied from the sedimentation scheme --> see vlz_fi.F)
+!***********************************************************************
+      SUBROUTINE van_leer(nlay,q,pente_max,masse,w,wq)
+
       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 
-!   ---------
-!
+
+!--------------------------------------------------------
+! Input/Output Variables
+!--------------------------------------------------------
+      INTEGER,INTENT(IN) :: nlay       ! number of atmospheric layers
+      REAL,INTENT(IN) ::  masse(nlay)  ! mass of the layer Dp/g
+      REAL,INTENT(IN) :: pente_max     != 2 conseillee
+      REAL,INTENT(INOUT) :: q(nlay)    ! mixing ratio (kg/kg)
+      REAL,INTENT(INOUT) :: w(nlay)    ! atmospheric mass "transferred" at each timestep (kg.m-2)
+      REAL,INTENT(INOUT) :: wq(nlay+1)
+
+!--------------------------------------------------------
+! Local Variables
+!--------------------------------------------------------
+
       INTEGER i,l,j,ii
-!
-
       REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
       REAL newmasse
@@ -541,11 +548,7 @@
       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
-
+! **********************************************************************
+!  Mixing ratio vertical gradient at the levels
+! **********************************************************************
       do l=2,nlay
             dzqw(l)=q(l-1)-q(l)
@@ -554,8 +557,4 @@
 
       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))
@@ -563,5 +562,4 @@
                 dzq(l)=0.
             endif
-#endif
             dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
             dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
@@ -571,23 +569,20 @@
       dzq(nlay)=0.
 
-! ---------------------------------------------------------------
-!   .... 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:
+! **********************************************************************
+!  Vertical advection
+! ********************************************************************** 
+
+       !! No flux at the model top:
        wq(nlay+1)=0.
 
-!      Surface flux up:
+       !! 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
-
-!      1) Compute wq where w < 0 (up)
+       do l = 1,nlay-1
+
+!      1) Compute wq where w < 0 (up) (UPWARD TRANSPORT)
 !      ===============================
 
          if(w(l+1).le.0)then
-
 !         Regular scheme (transfered mass < 1 layer)
 !         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -595,4 +590,8 @@
             sigw=w(l+1)/masse(l)
             wq(l+1)=w(l+1)*(q(l)-0.5*(1.+sigw)*dzq(l))
+!!-------------------------------------------------------
+!          The following part should not be needed in the
+!          case of an integration over subtimesteps
+!!-------------------------------------------------------
 !         Extended scheme (transfered mass > 1 layer)
 !         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -621,5 +620,5 @@
           endif ! (-w(l+1).le.masse(l))
      
-!      2) Compute wq where w > 0 (down)    
+!      2) Compute wq where w > 0 (down) (DOWNWARD TRANSPORT)     
 !      ===============================
 
@@ -630,7 +629,9 @@
           if(w(l).le.masse(l))then
             sigw=w(l)/masse(l)
-            wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l))
-            
-
+            wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l))            
+!!-------------------------------------------------------
+!          The following part should not be needed in the
+!          case of an integration over subtimesteps
+!!-------------------------------------------------------
 !         Extended scheme (transfered mass > 1 layer)
 !         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -661,5 +662,5 @@
        enddo ! l = 1,nlay-1
        
-       do l = 1,nlay          ! loop different than when w<0
+       do l = 1,nlay
 
 !         it cannot entrain more than available mass !
@@ -672,5 +673,5 @@
        enddo
        
-      END SUBROUTINE vl_storm
+      END SUBROUTINE van_leer
 
 !=======================================================================