Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F	(revision 608)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F	(revision 609)
@@ -1,5 +1,8 @@
       SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
      &    pq,ccn,tauref,tau,aerosol,reffrad,nueffrad,
-     &    QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
+     &                      QVISsQREF3d,omegaVIS3d,gVIS3d,
+     &                      QIRsQREF3d,omegaIR3d,gIR3d,
+     &                      QREFvis3d,QREFir3d,
+     &                      omegaREFvis3d,omegaREFir3d,
      &    zdqnorm,dsodust,pt)
                                                    
@@ -80,10 +83,25 @@
       REAL reffrad(ngrid,nlayer,naerkind)
       REAL nueffrad(ngrid,nlayer,naerkind)
-      REAL QREFvis3d(ngridmx,nlayermx,naerkind)
-      REAL QREFir3d(ngridmx,nlayermx,naerkind)
-      REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
-      REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
+!      REAL QREFvis3d(ngridmx,nlayermx,naerkind)
+!      REAL QREFir3d(ngridmx,nlayermx,naerkind)
+!      REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
+!      REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
       REAL zdqnorm(ngridmx,nlayermx,2)                       !mass mixing ratio perturbation due to the dust storm. Output for meso_physiq.F
       REAL pt(ngrid,nlayer)                                 !input: temperature. usefull to compute precisely the l_top parameter
+
+c     Aerosol optical properties
+      REAL :: QVISsQREF3d(ngridmx,nlayermx,nsun,naerkind)
+      REAL :: omegaVIS3d(ngridmx,nlayermx,nsun,naerkind)
+      REAL :: gVIS3d(ngridmx,nlayermx,nsun,naerkind)
+
+      REAL :: QIRsQREF3d(ngridmx,nlayermx,nir,naerkind)
+      REAL :: omegaIR3d(ngridmx,nlayermx,nir,naerkind)
+      REAL :: gIR3d(ngridmx,nlayermx,nir,naerkind)
+
+      REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
+      REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
+
+      REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
+      REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind)
 
 
@@ -142,4 +160,5 @@
          logical localstorm        ! =true to create a local dust storm
          real taulocref,ztoploc,radloc,lonloc,latloc  !local dust storm parameters
+         real reffstorm, yeah
          REAL ray(ngridmx)                            !distance from dust storm center
          REAL tauuser(ngridmx)                        !opacity perturbation due to dust storm
@@ -488,5 +507,4 @@
 
       IF (firstcall) THEN
-      WRITE(*,*) " RENORMALISATION !!! "
 c--------------------------------------------------
 c  Parameters of the opacity perturbation
@@ -494,11 +512,59 @@
 
       iaer=1  !!!! PROVISOIRE !!!!
-      taulocref = 4.25 !10  ! ref optical depth of the local dust storm
-      ztoploc = 10      ! target pseudo-altitude of local storm (km)
-      radloc = 0.5      ! radius of dust storm (degree)
-      lonloc = 25.      ! center longitude of storm (deg)
-      latloc = -3.      ! center latitude of storm (deg)
+
+        write(*,*) "Add a local storm ?"
+        localstorm=.true. ! default value
+        call getin("localstorm",localstorm)
+        write(*,*) " localstorm = ",localstorm
+
+        IF (localstorm) THEN
+          WRITE(*,*) "********************"
+          WRITE(*,*) "ADDING A LOCAL STORM"
+          WRITE(*,*) "********************"
+
+          !!! DEFAULT CASE
+          !4.25   ! ref optical depth of the local dust storm 
+          !10     ! target pseudo-altitude of local storm (km)
+          !0.5    ! radius of dust storm (degree) (actual radius is twice this)
+          !25.    ! center longitude of storm (deg)
+          !-2.5   ! center latitude of storm (deg)
+          !0.0    ! reff in storm (microns) 0. for background
+
+          write(*,*) "ref opacity of local dust storm"
+              taulocref = 4.25 ! default value
+              call getin("taulocref",taulocref)
+              write(*,*) " taulocref = ",taulocref
+
+          write(*,*) "target altitude of local storm (km)"
+              ztoploc = 10.0 ! default value
+              call getin("ztoploc",ztoploc)
+              write(*,*) " ztoploc = ",ztoploc
+
+          write(*,*) "radius of dust storm (degree)"
+              radloc = 0.5 ! default value
+              call getin("radloc",radloc)
+              write(*,*) " radloc = ",radloc
+
+          write(*,*) "center longitude of storm (deg)"
+              lonloc = 25.0 ! default value
+              call getin("lonloc",lonloc)
+              write(*,*) " lonloc = ",lonloc
+
+          write(*,*) "center latitude of storm (deg)"
+              latloc = -2.5 ! default value
+              call getin("latloc",latloc)
+              write(*,*) " latloc = ",latloc
+
+          write(*,*) "reff storm (mic) 0. for background"
+              reffstorm = 0.0 ! default value
+              call getin("reffstorm",reffstorm)
+              write(*,*) " reffstorm = ",reffstorm
+        ELSE
+          write(*,*) 'I need localstorm = T'
+          STOP
+        ENDIF
 
       DO ig=1,ngrid
+
 
 c---------------------------------------
@@ -509,10 +575,14 @@
      &          (long(ig)*180./pi -lonloc)**2)
 
+      !! transition factor for storm
+      !! increase factor ray diff for steepness
+      yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
+
 c-------------------------------------------------
 c           Tau's new map:
 c------------------------------------------
 
-      tauuser(ig)=max(tauref(ig) * pplev(ig,1) / 700.E0 , taulocref
-     &          * (TANH(2.+(radloc-ray(ig))*2)+1.)/2.)
+      tauuser(ig)=max(tauref(ig) * pplev(ig,1) / 700.E0 , 
+     &          taulocref * yeah)
 
 c---------------------------------------------------------
@@ -529,10 +599,41 @@
 
           ENDDO
-      ENDDO
+
+c---------------------------------------------------------
+c           change reffrad if ever needed
+c----------------------------------------------------------
+
+      IF (reffstorm .gt. 0.) THEN
+          DO l=1,nlayer
+             IF (l .lt. l_top+1) THEN
+
+                WRITE(*,*) "OLD REFFRAD", reffrad(ig,l,iaer)
+
+                reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
+     &          * 1.e-6 * yeah )
+
+                WRITE(*,*) "NEW REFFRAD", reffrad(ig,l,iaer), reffstorm
+     &          * 1.e-6 * yeah
+
+             ENDIF
+          ENDDO
+      ENDIF
+
+c---------------------------------------------------------------------------
+      ENDDO !! boucle sur ngridmx
+c---------------------------------------------------------------------------
+
+      !!! ReComputing 3D scattering parameters: if needed.
+      IF (reffstorm .gt. 0.) THEN
+      CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
+     &                      QVISsQREF3d,omegaVIS3d,gVIS3d,
+     &                      QIRsQREF3d,omegaIR3d,gIR3d,
+     &                      QREFvis3d,QREFir3d,
+     &                      omegaREFvis3d,omegaREFir3d)
+      ENDIF
 
 c----------------------------------------------------------------------------
 c    compute int_factor
 c---------------------------------------------------------------------------
-
 
       DO ig=1,ngrid
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/callradite_tachemobile_z.F
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/callradite_tachemobile_z.F	(revision 608)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/callradite_tachemobile_z.F	(revision 609)
@@ -391,5 +391,8 @@
       CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
      &      pq,ccn,tauref,tau,aerosol,reffrad,nueffrad,
-     &      QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
+     &                      QVISsQREF3d,omegaVIS3d,gVIS3d,
+     &                      QIRsQREF3d,omegaIR3d,gIR3d,
+     &                      QREFvis3d,QREFir3d,
+     &                      omegaREFvis3d,omegaREFir3d,
      &      zdqnorm,dsodust,pt)
 
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F	(revision 608)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F	(revision 609)
@@ -1,1 +1,1 @@
-link ../../../mars_lmd_new/libf/phymars/callsedim.F
+link callsedim.F.modified
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F.modified
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F.modified	(revision 609)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/callsedim.F.modified	(revision 609)
@@ -0,0 +1,370 @@
+      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
+     $                pplev,zlev, pt, rdust, rice,
+     &                pq, pdqfi, pdqsed,pdqs_sed,nq)
+
+! to use  'getin'
+      USE ioipsl_getincom
+
+      IMPLICIT NONE
+
+c=======================================================================
+c      Sedimentation of the  Martian aerosols
+c      depending on their density and radius
+c
+c      F.Forget 1999
+c
+c      Modified by J.-B. Madeleine 2010: Now includes the doubleq
+c        technique in order to have only one call to callsedim in
+c        physiq.F.
+c
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "tracer.h"
+#include "callkeys.h"
+
+c
+c   arguments:
+c   ----------
+
+      INTEGER ngrid              ! number of horizontal grid points
+      INTEGER nlay               ! number of atmospheric layers
+      REAL ptimestep             ! physics time step (s)
+      REAL pplev(ngrid,nlay+1)   ! pressure at inter-layers (Pa)
+      REAL pt(ngrid,nlay)        ! temperature at mid-layer (K)
+      REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
+c    Aerosol radius provided by the water ice microphysical scheme:
+      REAL rdust(ngrid,nlay)     ! Dust geometric mean radius (m)
+      REAL rice(ngrid,nlay)      ! Ice geometric mean radius (m)
+
+c    Traceurs :
+      integer nq             ! number of tracers
+      real pq(ngrid,nlay,nq)  ! tracers (kg/kg)
+      real pdqfi(ngrid,nlay,nq)  ! tendency before sedimentation (kg/kg.s-1)
+      real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
+      real pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
+      
+c   local:
+c   ------
+
+      REAL CBRT
+      EXTERNAL CBRT
+
+      INTEGER l,ig, iq
+
+      real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
+      real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
+      real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
+      real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
+      real r0(ngridmx,nlayermx) ! geometric mean doubleq radius (m)
+c     Sedimentation radius of water ice
+      real rfall(ngridmx,nlayermx)
+c     Sedimentation effective variance of water ice
+      REAL, PARAMETER :: nuice_sed = 0.45  !! TESTS_JB  !! 0.1 avant
+
+c     Discrete size distributions (doubleq)
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c       1) Parameters used to represent the changes in fall
+c          velocity as a function of particle size;
+      integer nr,ir
+      parameter (nr=12) !(nr=7) ! number of bins
+      real rd(nr),qr(ngridmx,nlayermx,nr)
+      real rdi(nr+1)    ! extreme and intermediate radii
+      real Sq(ngridmx,nlayermx)
+      real rdmin,rdmax,rdimin,rdimax
+      data rdmin/1.e-8/ !/1.e-7/
+      data rdmax/30.e-6/
+      data rdimin/1.e-10/
+      data rdimax/1e-4/
+      save rd, rdi
+
+c       2) Second size distribution for the log-normal integration
+c          (the mass mixing ratio is computed for each radius)
+
+      integer ninter, iint
+      parameter (ninter=4)  ! nombre de point entre chaque rayon rdi
+      real rr(ninter,nr)
+      save rr
+      integer radpower
+
+c       3) Other local variables used in doubleq
+
+      real reff(ngridmx,nlayermx,2)  ! for diagnostic only
+      INTEGER idust_mass   ! index of tracer containing dust mass
+                           !   mix. ratio
+      INTEGER idust_number ! index of tracer containing dust number
+                           !   mix. ratio
+      SAVE idust_mass,idust_number
+
+c     Firstcall:
+
+      LOGICAL firstcall
+      SAVE firstcall
+      DATA firstcall/.true./
+
+      REAL betaset
+      SAVE betaset
+
+
+c    ** un petit test de coherence
+c       --------------------------
+
+      IF (firstcall) THEN
+         IF(ngrid.NE.ngridmx) THEN
+            PRINT*,'STOP dans callsedim'
+            PRINT*,'probleme de dimensions :'
+            PRINT*,'ngrid  =',ngrid
+            PRINT*,'ngridmx  =',ngridmx
+            STOP
+         ENDIF
+
+         write(*,*) "beta coefficient for sedimentation"
+         betaset = 0.5 ! default value
+         call getin("betaset",betaset)
+         write(*,*) " betaset = ",betaset
+
+c       Doubleq: initialization
+        IF (doubleq) THEN
+         do ir=1,nr
+             rd(ir)= rdmin*(rdmax/rdmin)**(float(ir-1)/float(nr-1))
+         end do
+         rdi(1)=rdimin
+         do ir=2,nr
+           rdi(ir)= sqrt(rd(ir-1)*rd(ir))
+         end do
+         rdi(nr+1)=rdimax
+
+         do ir=1,nr
+           do iint=1,ninter
+             rr(iint,ir)=
+     &        rdi(ir)*
+     &        (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
+c             write(*,*) rr(iint,ir)
+           end do
+         end do
+
+      ! identify tracers corresponding to mass mixing ratio and
+      ! number mixing ratio
+        idust_mass=0      ! dummy initialization
+        idust_number=0    ! dummy initialization
+
+        do iq=1,nq
+          if (noms(iq).eq."dust_mass") then
+            idust_mass=iq
+          endif
+          if (noms(iq).eq."dust_number") then
+            idust_number=iq
+          endif
+        enddo
+
+        ! check that we did find the tracers
+        if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
+          write(*,*) 'callsedim: error! could not identify'
+          write(*,*) ' tracers for dust mass and number mixing'
+          write(*,*) ' ratio and doubleq is activated!'
+          stop
+        endif
+        ENDIF !of if (doubleq)
+
+        IF (water) THEN
+          write(*,*) "water_param nueff Sedimentation:", nuice_sed
+          IF (activice) THEN
+            write(*,*) "water_param nueff Radiative:", nuice_ref
+          ENDIF
+        ENDIF
+      
+        firstcall=.false.
+      ENDIF ! of IF (firstcall)
+
+c-----------------------------------------------------------------------
+c    1. Initialization
+c    -----------------
+
+      zqi(1:ngrid,1:nlay,1:nqmx) = 0.
+c     Updating the mass mixing ratio with the tendencies coming
+c       from other parameterizations:
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      do iq=1,nq
+        do l=1,nlay
+          do ig=1,ngrid
+            zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
+          enddo
+        enddo 
+      enddo
+
+c    Computing the different layer properties
+c    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c    Mass (kg.m-2), thickness(m), crossing time (s)  etc.
+
+      do  l=1,nlay
+        do ig=1, ngrid
+          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
+          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
+        end do
+      end do
+
+c =================================================================
+      do iq=1,nq
+        if(radius(iq).gt.1.e-9) then   ! no sedim for gaz
+
+c -----------------------------------------------------------------
+c         DOUBLEQ CASE
+c -----------------------------------------------------------------
+          if (doubleq.and.
+     &        ((iq.eq.idust_mass).or.
+     &         (iq.eq.idust_number))) then
+
+c           Computing size distribution:
+c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            do  l=1,nlay
+              do ig=1, ngrid
+                r0(ig,l)=
+     &            CBRT(r3n_q*zqi(ig,l,idust_mass)/
+     &            max(zqi(ig,l,idust_number),0.01))
+                r0(ig,l)=min(max(r0(ig,l),1.e-10),500.e-6)
+              end do
+            end do
+
+c        Computing mass mixing ratio for each particle size
+c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+          IF (iq.EQ.idust_mass) then
+            radpower = 2
+          ELSE
+            radpower = -1
+          ENDIF
+          Sq(1:ngrid,1:nlay) = 0.
+          do ir=1,nr
+            do l=1,nlay
+              do ig=1,ngrid
+c                ****************
+c                Size distribution integration
+c                (Trapezoid Integration Method)
+                 qr(ig,l,ir)=0.5*(rr(2,ir)-rr(1,ir))*
+     &             (rr(1,ir)**radpower)*
+     &             exp(-(log(rr(1,ir)/r0(ig,l)))**2/(2*varian**2))
+                 do iint=2,ninter-1
+                   qr(ig,l,ir)=qr(ig,l,ir) +
+     &             0.5*(rr(iint+1,ir)-rr(iint-1,ir))*
+     &             (rr(iint,ir)**radpower)*
+     &             exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
+     &             (2*varian**2))
+                 end do
+                 qr(ig,l,ir)=qr(ig,l,ir) +
+     &             0.5*(rr(ninter,ir)-rr(ninter-1,ir))*
+     &             (rr(ninter,ir)**radpower)*
+     &             exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
+     &             (2*varian**2))
+
+c                **************** old method (not recommended!)
+c                qr(ig,l,ir)=(rd(ir)**(5-3*iq))*
+c    &           exp( -(log(rd(ir)/r0(ig,l)))**2 / (2*varian**2) )
+c                ******************************
+
+                 Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
+              enddo
+            enddo
+          enddo
+
+          do ir=1,nr
+            do l=1,nlay
+              do ig=1,ngrid
+                 qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
+              enddo
+            enddo
+          enddo
+
+c         Computing sedimentation for each tracer
+c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+c         call zerophys(ngridmx*nlayermx,zqi(1,1,iq))
+          zqi(1:ngrid,1:nlay,iq) = 0.
+c         call zerophys(ngridmx,pdqs_sed(1,iq))
+          pdqs_sed(1:ngrid,iq) = 0.
+
+          do ir=1,nr
+             call newsedim(ngrid,nlay,1,ptimestep,
+     &       pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
+     &       wq,betaset)
+
+c            Tendencies
+c            ~~~~~~~~~~
+             do ig=1,ngrid
+               pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
+     &                                + wq(ig,1)/ptimestep
+             end do
+             DO l = 1, nlay
+               DO ig=1,ngrid
+                 zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
+               ENDDO
+             ENDDO
+          enddo ! of do ir=1,nr
+c -----------------------------------------------------------------
+c         WATER CYCLE CASE
+c -----------------------------------------------------------------
+          else if (water.and.(iq.eq.igcm_h2o_ice)) then
+c           if (doubleq) then
+c             do l=1,nlay
+c               do ig=1,ngrid
+c                 rfall(ig,l)=max( rice(ig,l),rdust(ig,l) )
+c                 rfall(ig,l)=min(rfall(ig,l),1.e-4)
+c               enddo
+c             enddo ! of do l=1,nlay
+c           else
+              do l=1,nlay
+                do ig=1,ngrid
+c               For water cycle, a typical dust size is assumed:
+c               r(z)=r0*exp(-z/H) with r0=0.8 micron and H=18 km.
+c               rfall(ig,l)=max( rice(ig,l)*1.5,rdust(ig,l) )
+                rfall(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3.,
+     &                           rdust(ig,l) )
+c modif FranckMM pour ameliorer cycle H2O: rfall= 20 microns
+c mars commente pour l'instant               rfall(ig,l)=20.e-6
+                rfall(ig,l)=min(rfall(ig,l),1.e-4)
+                enddo
+              enddo ! of do l=1,nlay
+c           endif
+            call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
+     &      pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi(1,1,iq),
+     &      wq,1.0)
+c           Tendencies
+c           ~~~~~~~~~~
+            do ig=1,ngrid 
+              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
+            end do
+c -----------------------------------------------------------------
+c         GENERAL CASE
+c -----------------------------------------------------------------
+          else
+            call newsedim(ngrid,nlay,1,ptimestep,
+     &      pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
+     &      zqi(1,1,iq),wq,1.0)
+c           Tendencies
+c           ~~~~~~~~~~
+            do ig=1,ngrid 
+              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
+            end do
+          endif ! of if doubleq and if water
+c -----------------------------------------------------------------
+
+          DO l = 1, nlay
+            DO ig=1,ngrid
+              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
+     $        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
+            ENDDO
+          ENDDO
+
+        endif ! of if(radius(iq).gt.1.e-9)
+c =================================================================
+      enddo ! of do iq=1,nq
+ 
+      RETURN
+      END
+
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/localstorm.def
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/localstorm.def	(revision 609)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/localstorm.def	(revision 609)
@@ -0,0 +1,39 @@
+### "Add a local storm ?" 
+###      default is T
+#########################
+localstorm        = .true. 
+
+### "ref opacity of local dust storm" 
+###      default is 4.25
+#####################################
+taulocref         = 4.25 
+
+### "target altitude of local storm (km)" 
+###      default is 10.0
+#########################################
+ztoploc           = 10.0 
+
+### "radius of dust storm (degree)" 
+###      default is 0.5
+###################################
+radloc            = 0.5 
+
+### "center longitude of storm (deg)" 
+###      default is 25.0
+#####################################
+lonloc            = 25.0 
+
+### "center latitude of storm (deg)" 
+###      default is -2.5
+####################################
+latloc            = -2.5 
+
+### "reff storm (mic) 0. for background" 
+###      default is 0.0 (use background reff)
+#############################################
+reffstorm         = 0.0 
+
+### "beta coefficient for sedimentation"
+###      default is 0.5
+########################################
+betaset           = 0.5 
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/physiq.F
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/physiq.F	(revision 608)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/physiq.F	(revision 609)
@@ -296,4 +296,5 @@
       SAVE ccn  !! in case iradia != 1
       real rdust(ngridmx,nlayermx) ! dust geometric mean radius (m)
+      real reffdust(ngridmx,nlayermx) ! effective radius (m)
       real qtot1,qtot2 ! total aerosol mass
       integer igmin, lmin
@@ -1046,4 +1047,9 @@
              ENDDO
            ENDDO
+
+             DO ig=1,ngrid
+              zdqssed_diag(ig)= zdqssed(ig,igcm_dust_mass)
+             ENDDO
+
          END IF   ! of IF (sedimentation)
 
@@ -1456,6 +1462,8 @@
      .           * mugaz / mmol(igcm_h2o_ice)
       ENDIF
+
       !! Dust quantity integration along the vertical axe
       dustot(:)=0
+      reffdust(:,:)=0
       do ig=1,ngrid
        do l=1,nlayermx
@@ -1463,4 +1471,5 @@
      &               zq(ig,l,igcm_dust_mass)
      &               *  (pplev(ig,l) - pplev(ig,l+1)) / g
+        reffdust(ig,l) = rdust(ig,l)*ref_r0
        enddo
       enddo
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/run.def
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/run.def	(revision 609)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/run.def	(revision 609)
@@ -0,0 +1,4 @@
+# some definitions for the physics, in file 'callphys.def'
+INCLUDEDEF=callphys.def
+
+INCLUDEDEF=localstorm.def
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/stress.def
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/stress.def	(revision 609)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/stress.def	(revision 609)
@@ -0,0 +1,29 @@
+10.
+0.000001
+
+0.5     !! ulim
+0.002   !! alpha
+
+10.
+0.000001
+
+0.5
+0.002
+
+0.5
+0.005
+
+0.5
+0.001
+
+0.5
+0.05
+
+0.6
+0.1
+
+0.6
+1.
+
+1.061
+1.
