Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 309)
+++ /trunk/LMDZ.MARS/README	(revision 310)
@@ -1023,2 +1023,9 @@
  M             299   libf/phymars/calltherm_mars.F90
  ^-----------------> Changed r_aspect_thermals from 2. to 1.5 for the GCM version to better match buoyancy profiles.
+
+== 29/09/11 == AS
+--> To easily explore sensitivity to lifting thresholds: in dustlift.F, ustar_seuil=sqrt(stress/rho) 
+    and alpha_lift[dust_mass] can be prescribed through an external stress.def parameter file.
+    --- alpha_lift[dust_number] is computed from alpha_lift[dust_mass] as in initracer.F
+    --- ustar_seuil is more user-friendly than stress because direct comparison with ustar from model 
+--> For the moment this is MESOSCALE only, but potentially useful to everyone
Index: /trunk/LMDZ.MARS/libf/phymars/dustlift.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dustlift.F	(revision 309)
+++ /trunk/LMDZ.MARS/libf/phymars/dustlift.F	(revision 310)
@@ -50,14 +50,22 @@
 !!!! AS: you have to compile with -DMESOSCALE to do so
       REAL alpha
+      REAL r0_lift
       INTEGER ierr
+      REAL ulim
         OPEN(99,file='stress.def',status='old',form='formatted'
      .   ,iostat=ierr)
         !!! no file => default values
         IF(ierr.EQ.0) THEN
-          READ(99,*) stress_seuil
+          READ(99,*) ulim !ulim = sqrt(stress_seuil/rho) avec rho = 0.02. 
+                          !prendre ulim = 1.061 m/s pour avoir stress_seuil = 0.0225 
           READ(99,*) alpha
+          stress_seuil = 0.02 * ulim * ulim
           write(*,*) 'USER-DEFINED threshold: ', stress_seuil, alpha
           CLOSE(99)
-          alpha_lift(1:nq) = alpha
+          alpha_lift(igcm_dust_mass) = alpha
+          r0_lift = radius(igcm_dust_mass) / ref_r0
+          alpha_lift(igcm_dust_number)=r3n_q*
+     &                        alpha_lift(igcm_dust_mass)/r0_lift**3
+          write(*,*) 'set dust number: ', alpha_lift(igcm_dust_number)
         ENDIF
 #endif
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 309)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 310)
@@ -494,4 +494,5 @@
 c       Dust lifting:
         if (lifting) then
+#ifndef MESOSCALE
            if (doubleq.AND.submicron) then
              do ig=1,ngrid
@@ -522,4 +523,8 @@
      &                    pdqsdif)
            endif !doubleq.AND.submicron
+#else
+            call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
+     &                    pdqsdif)
+#endif
         else
            pdqsdif(1:ngrid,1:nq) = 0.
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/2dplot.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/2dplot.py	(revision 309)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/2dplot.py	(revision 310)
@@ -26,7 +26,17 @@
 #title = "Atmospheric temperature"
 
+name = "wrfout_d01_2024-03-22_01:00:00_z"
+itime = 12
+ndiv = 10
+zey = 120
+var = "VMR_ICE"
+title = "Volume mixing ratio of water ice [ppm]"
+vmin = -0.5
+vmax = 400.
+
 nc = Dataset(name)
 
 what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d2=zey )
+
 
 y = nc.variables["vert"][:]
@@ -34,4 +44,10 @@
 horinp = len(what_I_plot[0,:])
 x = np.linspace(0.,horinp*500.,horinp) / 1000.
+xlabel("Horizontal coordinate (km)")
+
+horinp = len(what_I_plot[0,:])
+x = np.linspace(0.,horinp,horinp) 
+xlabel("x grid points")
+
 
 zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
@@ -46,5 +62,4 @@
                        extend='neither',spacing='proportional')
 ptitle(title)
-xlabel("Horizontal coordinate (km)")
 ylabel("Altitude (m)")
 makeplotres(var+str(itime),res=200.,disp=False)
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py	(revision 309)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py	(revision 310)
@@ -23,5 +23,7 @@
     from scipy.ndimage.measurements import minimum_position
     from netCDF4 import Dataset
-    
+    import matplotlib.pyplot as plt
+    import myplot as myp
+   
     ### LOAD NETCDF DATA
     filename = "psfc.nc"
@@ -34,24 +36,51 @@
     allsizesx = []
     allsizesy = []
-    for i in range(0,shape[0]):
-    #for i in range(0,2):
-    
-        print i, ' on ', shape[0]
+    stride = 1 #5
+    for i in range(0,shape[0],stride):
+
         psfc2d = np.array ( psfc [ i, : , : ] )
     
         ############### CRITERION
-        ave = np.mean(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy
-        where = np.where(psfc2d - ave < -0.3)  ## comme le papier Phoenix
+        ave = np.mean(psfc2d,dtype=np.float64)  ## dtype otherwise inaccuracy
+
+        #limdp = -0.2 ## on en loupe pas mal
+        #where = np.where(psfc2d - ave < limdp)  ## comme le papier Phoenix
     
-        std = np.std(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy
-        fac = 2.5  ## not too low, otherwise vortices are not caught. 4 good choice.
+        std = np.std(psfc2d,dtype=np.float64)   ## dtype otherwise inaccuracy
+        fac = 4.   ## how many sigmas. not too low, otherwise vortices are not caught. 4 good choice.
+                   ## 2.5 clearly too low, 3.5 not too bad, 4 probably good
+        fac = 3.5
+        #fac = 2.5
         lim = ave - fac*std
-        print lim, ave, std
         where = np.where(psfc2d < lim)
         ############### END CRITERION
-    
+   
         lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
         lab[where] = 1.  ## do not treat points close to 'mean' (background) pressure
-    
+
+        draw = False
+        if draw:
+        ##################################################################################
+            vmin = -0.3
+            vmax =  0.0
+            ndiv = 3
+            palette = plt.get_cmap(name="YlGnBu") 
+            what_I_plot = psfc2d-ave
+            zevmin, zevmax = myp.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
+            what_I_plot = myp.bounds(what_I_plot,zevmin,zevmax)
+            zelevels = np.linspace(zevmin,zevmax)
+            fig = plt.figure() 
+            sub = myp.definesubplot(2,fig)
+            plt.subplot(sub) 
+            plt.contourf(what_I_plot,zelevels,cmap=palette) 
+            plt.colorbar(fraction=0.05,pad=0.03,format="%.1f",\
+                           ticks=np.linspace(zevmin,zevmax,ndiv+1),\
+                           extend='both',spacing='proportional')
+            plt.subplot(sub+1)
+            palette = plt.get_cmap(name="binary")
+            plt.contourf(lab,2,cmap=palette)
+            plt.show() 
+        ##################################################################################
+  
         xx = []
         yy = []
@@ -69,9 +98,18 @@
         sizex = detsize( xx, res = 10, loga=False, thres=2 )
         sizey = detsize( yy, res = 10, loga=False, thres=2 )
+        #sizex = detsize( xx, res = 10, loga=False, thres=3 )
+        #sizey = detsize( yy, res = 10, loga=False, thres=3 )
         ###
+        #if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex  ### un peu limite dans certains cas
+        if (len(sizex) > len(sizey))     : sizey = sizex        ### plus fidele mais petit souci lorsque PBC
+        elif (len(sizex) == len(sizey))  : 
+          if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex
+          else                                  : sizex = sizey
+        else                            : sizex = sizey
         allsizesx = np.append(allsizesx,sizex)
         allsizesy = np.append(allsizesy,sizey)
+        print i, ' on ', shape[0], ' caught ', len(sizex), ' vortices ', sizex
         ########
-    
+  
     allsizesx.sort()
     allsizesy.sort()
@@ -79,11 +117,6 @@
     return allsizesx, allsizesy
 
-def writeascii ( tab, filename ):
-    mydata = tab
-    myfile = open(filename, 'w')
-    for line in mydata:
-        myfile.write(str(line) + '\n')
-    myfile.close()
-    return
+#########################################################################
+#########################################################################
 
 import matplotlib.pyplot as plt
@@ -91,21 +124,23 @@
 import numpy as np
 import matplotlib.mlab as mlab
+import mymath as mym
 
 save = True
-#save = False
+save = False
+
 if save:
     allsizesx, allsizesy = getsize("psfc.nc")
-    writeascii(allsizesx,'allsizex.txt')
-    writeascii(allsizesy,'allsizey.txt')
-    
+    ### sauvegarde texte pour inspection
+    mym.writeascii(allsizesx,'allsizex.txt')
+    mym.writeascii(allsizesy,'allsizey.txt')
+    ### sauvegarde binaire pour utilisation python
     myfile = open('allsizex.bin', 'wb')
     pickle.dump(allsizesx, myfile)
     myfile.close()
-    
     myfile = open('allsizey.bin', 'wb')
     pickle.dump(allsizesy, myfile)
     myfile.close()
 
-
+### load files
 myfile = open('allsizex.bin', 'r')
 allsizesx = pickle.load(myfile)
@@ -113,16 +148,33 @@
 allsizesy = pickle.load(myfile)
 
+### append sizes
 plothist = np.append(allsizesx,allsizesy)
+plothist = allsizesx
+#plothist = allsizesy
 plothist.sort()
 
-nbins=100
+### make bins
+nbins = 100
 zebins = [2.0]
+#nbins = 8
+#zebins = [19.0]
+#nbins = 15
+#zebins = [11.] ##12 non mais donne un peu la meme chose
+#zebins = [20.] ##20 non car trop pres du premier
+nbins = 100
+zebins = [2./np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
+nbins = 200
+zebins = [0.2/np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
+
 for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
 zebins = np.array(zebins)
-print zebins
-
+### select reasonable bins for DD
 zebins = zebins [ zebins > 10. ] 
 zebins = zebins [ zebins < 1000. ]
-print zebins
+#print 'bins ',zebins
+print 'corrected bins ',zebins
+print 'max value ',mym.max(plothist),mym.max(allsizesx),mym.max(allsizesy) 
+
+#plothist = [20.,30.,40.,50.,70.,100.,130.,190.,260.,370.,520.]  ## pour calibrer
 
 #### HISTOGRAM
@@ -130,8 +182,7 @@
              log=True,\
              bins=zebins,\
-             #cumulative=-1,\
+#             cumulative=-1,\
         )
 plt.xscale('log')
-
 plt.show()
 
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py	(revision 309)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py	(revision 310)
@@ -14,2 +14,10 @@
         return u'\u00b0'
 
+def writeascii ( tab, filename ):
+    mydata = tab
+    myfile = open(filename, 'w')
+    for line in mydata:
+        myfile.write(str(line) + '\n')
+    myfile.close()
+    return
+
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py	(revision 309)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py	(revision 310)
@@ -450,9 +450,12 @@
     from mymath import max,min,mean
     ### might be convenient to add the missing value in arguments
-    what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
+    #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
+    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
+    else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
     print "new min ", min(what_I_plot)
     what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
-    what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
+    what_I_plot[ what_I_plot > zevmax ] = zevmax
     print "new max ", max(what_I_plot)
+    
     return what_I_plot
 
@@ -508,5 +511,5 @@
              "USTM":         "YlOrRd",\
              "HFX":          "RdYlBu",\
-             "ICETOT":       "YlGnBu",\
+             "ICETOT":       "YlGnBu_r",\
              "MTOT":         "PuBu",\
              "TAU_ICE":      "Blues",\
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 309)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/aeropacity_tachemobile_z.F	(revision 310)
@@ -497,6 +497,6 @@
       ztoploc = 10      ! target pseudo-altitude of local storm (km)
       radloc = 0.5      ! radius of dust storm (degree)
-      lonloc = 80       ! center longitude of storm (deg)
-      latloc = -25      ! center latitude of storm (deg)
+      lonloc = 25.      ! center longitude of storm (deg)
+      latloc = -3.      ! center latitude of storm (deg)
 
       DO ig=1,ngrid
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/initracer.F
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/initracer.F	(revision 310)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/STORM_JULIEN_LAST/initracer.F	(revision 310)
@@ -0,0 +1,749 @@
+      SUBROUTINE initracer(qsurf,co2ice)
+
+       IMPLICIT NONE
+c=======================================================================
+c   subject:
+c   --------
+c   Initialization related to tracer 
+c   (transported dust, water, chemical species, ice...)
+c
+c   Name of the tracer
+c
+c   Test of dimension :
+c   Initialize COMMON tracer in tracer.h, using tracer names provided
+c   by the dynamics in "advtrac.h"
+c
+c   Old conventions: (not used any more)
+c
+c   If water=T : q(iq=nqmx) is the water mass mixing ratio
+c     and q(iq=nqmx-1) is the ice mass mixing ratio
+
+c   If there is transported dust, it uses iq=1 to iq=dustbin
+c   If there is no transported dust : dustbin=0
+c   If doubleq=T : q(iq=1) is the dust mass mixing ratio
+c                  q(iq=2) is the dust number mixing ratio 
+
+c   If (photochem.or.thermochem) there is "ncomp" chemical species (ncomp
+c   is set in aeronomars/chimiedata.h) using the ncomp iq values starting at
+c      iq=nqchem_min = dustbin+1   (nqchem_min is defined in inifis.F)
+c   
+c
+c   author: F.Forget
+c   ------
+c    Modifs: Franck Montmessin, Sebastien Lebonnois (june 2003)
+c            Ehouarn Millour (oct. 2008) identify tracers by their names
+c=======================================================================
+
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "callkeys.h"
+#include "tracer.h"
+#include "advtrac.h"
+#include "comgeomfi.h"
+#include "chimiedata.h"
+
+#include "surfdat.h"
+
+      real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g.  kg.m-2)
+      real co2ice(ngridmx)           ! co2 ice mass on surface (e.g.  kg.m-2)
+      integer iq,ig,count
+      real r0_lift , reff_lift, nueff_lift
+c     Ratio of small over large dust particles (used when both 
+c       doubleq and the submicron mode are active); In Montmessin
+c       et al. (2002), a value of 25 has been deduced;
+      real, parameter :: popratio = 25.
+      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
+      character(len=20) :: txt ! to store some text
+
+c-----------------------------------------------------------------------
+c  radius(nqmx)      ! aerosol particle radius (m)
+c  rho_q(nqmx)       ! tracer densities (kg.m-3)
+c  alpha_lift(nqmx)  ! saltation vertical flux/horiz flux ratio (m-1)
+c  alpha_devil(nqmx) ! lifting coeeficient by dust devil
+c  rho_dust          ! Mars dust density
+c  rho_ice           ! Water ice density
+c  nuice_ref         ! Effective variance nueff of the
+c                    !   water-ice size distributions
+c  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
+c  varian            ! Characteristic variance of log-normal distribution
+c-----------------------------------------------------------------------
+
+      integer :: nqchem_start
+
+! Initialization: get tracer names from the dynamics and check if we are
+!                 using 'old' tracer convention ('q01',q02',...)
+!                 or new convention (full tracer names)
+      ! check if tracers have 'old' names
+      oldnames=.false.
+      count=0
+      do iq=1,nqmx
+        txt=" "
+        write(txt,'(a1,i2.2)') 'q',iq
+        if (txt.eq.tnom(iq)) then
+          count=count+1
+        endif
+      enddo ! of do iq=1,nqmx
+      
+      if (count.eq.nqmx) then
+        write(*,*) "initracer: tracers seem to follow old naming ",
+     &             "convention (q01,q02,...)"
+        write(*,*) "   => will work for now ... "
+        write(*,*) "      but you should run newstart to rename them"
+        oldnames=.true.
+      endif
+
+      ! copy/set tracer names
+      if (oldnames) then ! old names (derived from iq & option values)
+        ! 1. dust:
+        if (dustbin.ne.0) then ! transported dust
+          do iq=1,dustbin
+            txt=" "
+            write(txt,'(a4,i2.2)') 'dust',iq
+            noms(iq)=txt
+            mmol(iq)=100.
+          enddo
+          ! special case if doubleq
+          if (doubleq) then
+            if (dustbin.ne.2) then
+              write(*,*) 'initracer: error expected dustbin=2'
+            else
+!              noms(1)='dust_mass'   ! dust mass mixing ratio
+!              noms(2)='dust_number' ! dust number mixing ratio
+! same as above, but avoid explict possible out-of-bounds on array
+               noms(1)='dust_mass'   ! dust mass mixing ratio
+               do iq=2,2
+               noms(iq)='dust_number' ! dust number mixing ratio
+               enddo
+            endif
+          endif
+        endif
+        ! 2. water & ice
+        if (water) then
+          noms(nqmx)='h2o_vap'
+          mmol(nqmx)=18.
+!            noms(nqmx-1)='h2o_ice'
+!            mmol(nqmx-1)=18.
+          !"loop" to avoid potential out-of-bounds in array
+          do iq=nqmx-1,nqmx-1
+            noms(iq)='h2o_ice'
+            mmol(iq)=18.
+          enddo
+        endif
+        ! 3. Chemistry
+        if (photochem .or. callthermos) then
+          if (doubleq) then
+            nqchem_start=3
+          else
+            nqchem_start=dustbin+1
+          end if
+          do iq=nqchem_start,nqchem_start+ncomp-1
+            noms(iq)=nomchem(iq-nqchem_start+1)
+            mmol(iq)=mmolchem(iq-nqchem_start+1)
+          enddo
+        endif ! of if (photochem .or. callthermos)
+        ! 4. Other tracers ????
+        if ((dustbin.eq.0).and.(.not.water)) then
+          noms(1)='co2'
+          mmol(1)=44
+          if (nqmx.eq.2) then
+            noms(nqmx)='Ar_N2'
+            mmol(nqmx)=30
+          endif
+        endif
+      else
+        ! copy tracer names from dynamics
+        do iq=1,nqmx
+          noms(iq)=tnom(iq)
+        enddo
+        ! mmol(:) array is initialized later (see below)
+      endif ! of if (oldnames)
+
+      ! special modification when using 'old' tracers:
+      ! qsurf(nqmx) was h2o ice whereas q(nqmx) was water vapour
+      ! and (if iceparty) q(nqmx-1) was null whereas q(nqmx-1) was water ice
+      if (oldnames.and.water) then
+        write(*,*)'initracer: moving surface water ice to index ',nqmx-1
+        ! "loop" to avoid potential out-of-bounds on arrays
+        do iq=nqmx-1,nqmx-1
+          qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1)
+        enddo
+        qsurf(1:ngridmx,nqmx)=0
+      endif 
+      
+c------------------------------------------------------------
+c     Test Dimensions tracers
+c------------------------------------------------------------
+
+! Ehouarn: testing number of tracers is obsolete...
+!      if(photochem.or.thermochem) then
+!          if (water) then
+!              if ((dustbin+ncomp+2).ne.nqmx) then
+!                 print*,'initracer: tracer dimension problem:'
+!                 print*,"(dustbin+ncomp+2).ne.nqmx"
+!                 print*,"ncomp: ",ncomp
+!                 print*,"dustbin: ",dustbin
+!                 print*,"nqmx: ",nqmx
+!                 print*,'Change ncomp in chimiedata.h'
+!               endif
+!          else
+!              if ((dustbin+ncomp+1).ne.nqmx) then
+!                 print*,'initracer: tracer dimension problem:'
+!                 print*,"(dustbin+ncomp+1).ne.nqmx"
+!                 print*,"ncomp: ",ncomp
+!                 print*,"dustbin: ",dustbin
+!                 print*,"nqmx: ",nqmx
+!                 print*,'Change ncomp in chimiedata.h'
+!                 STOP
+!               endif
+!            endif
+!      endif
+
+c------------------------------------------------------------
+c         NAME and molar mass of the tracer 
+c------------------------------------------------------------
+
+   
+! Identify tracers by their names: (and set corresponding values of mmol)
+      ! 0. initialize tracer indexes to zero:
+      do iq=1,nqmx
+        igcm_dustbin(iq)=0
+      enddo
+      igcm_dust_mass=0
+      igcm_dust_number=0
+      igcm_dust_submicron=0
+      igcm_h2o_vap=0
+      igcm_h2o_ice=0
+      igcm_co2=0
+      igcm_co=0
+      igcm_o=0
+      igcm_o1d=0
+      igcm_o2=0
+      igcm_o3=0
+      igcm_h=0
+      igcm_h2=0
+      igcm_oh=0
+      igcm_ho2=0
+      igcm_h2o2=0
+      igcm_n2=0
+      igcm_ar=0
+      igcm_ar_n2=0
+      igcm_n=0
+      igcm_no=0
+      igcm_no2=0
+      igcm_n2d=0
+      igcm_co2plus=0
+      igcm_oplus=0
+      igcm_o2plus=0
+      igcm_coplus=0
+      igcm_cplus=0
+      igcm_nplus=0
+      igcm_noplus=0
+      igcm_n2plus=0
+      igcm_hplus=0
+      igcm_elec=0
+
+
+      ! 1. find dust tracers
+      count=0
+      if (dustbin.gt.0) then
+        do iq=1,nqmx
+          txt=" "
+          write(txt,'(a4,i2.2)')'dust',count+1
+          if (noms(iq).eq.txt) then
+            count=count+1
+            igcm_dustbin(count)=iq
+            mmol(iq)=100.
+          endif
+        enddo !do iq=1,nqmx
+      endif ! of if (dustbin.gt.0)
+      if (doubleq) then
+        do iq=1,nqmx
+          if (noms(iq).eq."dust_mass") then
+            igcm_dust_mass=iq
+            count=count+1
+          endif
+          if (noms(iq).eq."dust_number") then
+            igcm_dust_number=iq
+            count=count+1
+          endif
+        enddo
+      endif ! of if (doubleq)
+      if (submicron) then
+        do iq=1,nqmx
+          if (noms(iq).eq."dust_submicron") then
+            igcm_dust_submicron=iq
+            mmol(iq)=100.
+            count=count+1
+          endif
+        enddo
+      endif ! of if (submicron)
+      ! 2. find chemistry and water tracers
+      do iq=1,nqmx
+        if (noms(iq).eq."co2") then
+          igcm_co2=iq
+          mmol(igcm_co2)=44.
+          count=count+1
+        endif
+        if (noms(iq).eq."co") then
+          igcm_co=iq
+          mmol(igcm_co)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."o") then
+          igcm_o=iq
+          mmol(igcm_o)=16.
+          count=count+1
+        endif
+        if (noms(iq).eq."o1d") then
+          igcm_o1d=iq
+          mmol(igcm_o1d)=16.
+          count=count+1
+        endif
+        if (noms(iq).eq."o2") then
+          igcm_o2=iq
+          mmol(igcm_o2)=32.
+          count=count+1
+        endif
+        if (noms(iq).eq."o3") then
+          igcm_o3=iq
+          mmol(igcm_o3)=48.
+          count=count+1
+        endif
+        if (noms(iq).eq."h") then
+          igcm_h=iq
+          mmol(igcm_h)=1.
+          count=count+1
+        endif
+        if (noms(iq).eq."h2") then
+          igcm_h2=iq
+          mmol(igcm_h2)=2.
+          count=count+1
+        endif
+        if (noms(iq).eq."oh") then
+          igcm_oh=iq
+          mmol(igcm_oh)=17.
+          count=count+1
+        endif
+        if (noms(iq).eq."ho2") then
+          igcm_ho2=iq
+          mmol(igcm_ho2)=33.
+          count=count+1
+        endif
+        if (noms(iq).eq."h2o2") then
+          igcm_h2o2=iq
+          mmol(igcm_h2o2)=34.
+          count=count+1
+        endif
+        if (noms(iq).eq."n2") then
+          igcm_n2=iq
+          mmol(igcm_n2)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."ar") then
+          igcm_ar=iq
+          mmol(igcm_ar)=40.
+          count=count+1
+        endif
+	if (noms(iq).eq."n") then
+          igcm_n=iq
+          mmol(igcm_n)=14.
+          count=count+1
+        endif
+        if (noms(iq).eq."no") then
+          igcm_no=iq
+          mmol(igcm_no)=30.
+          count=count+1
+        endif
+        if (noms(iq).eq."no2") then
+          igcm_no2=iq
+          mmol(igcm_no2)=46.
+          count=count+1
+        endif
+        if (noms(iq).eq."n2d") then
+          igcm_n2d=iq
+          mmol(igcm_n2d)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."co2plus") then
+          igcm_co2plus=iq
+          mmol(igcm_co2plus)=44.
+          count=count+1
+        endif
+        if (noms(iq).eq."oplus") then
+          igcm_oplus=iq
+          mmol(igcm_oplus)=16.
+          count=count+1
+        endif
+        if (noms(iq).eq."o2plus") then
+          igcm_o2plus=iq
+          mmol(igcm_o2plus)=32.
+          count=count+1
+        endif
+        if (noms(iq).eq."coplus") then
+          igcm_coplus=iq
+          mmol(igcm_coplus)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."cplus") then
+          igcm_cplus=iq
+          mmol(igcm_cplus)=12.
+          count=count+1
+        endif
+        if (noms(iq).eq."nplus") then
+          igcm_nplus=iq
+          mmol(igcm_nplus)=14.
+          count=count+1
+        endif
+        if (noms(iq).eq."noplus") then
+          igcm_noplus=iq
+          mmol(igcm_noplus)=30.
+          count=count+1
+        endif
+	if (noms(iq).eq."n2plus") then
+          igcm_n2plus=iq
+          mmol(igcm_n2plus)=28.
+          count=count+1
+        endif
+        if (noms(iq).eq."hplus") then
+          igcm_hplus=iq
+          mmol(igcm_hplus)=1.
+          count=count+1
+        endif
+        if (noms(iq).eq."elec") then
+          igcm_elec=iq
+          mmol(igcm_elec)=1./1822.89
+          count=count+1
+        endif
+        if (noms(iq).eq."h2o_vap") then
+          igcm_h2o_vap=iq
+          mmol(igcm_h2o_vap)=18.
+          count=count+1
+        endif
+        if (noms(iq).eq."h2o_ice") then
+          igcm_h2o_ice=iq
+          mmol(igcm_h2o_ice)=18.
+          count=count+1
+        endif
+        ! Other stuff: e.g. for simulations using co2 + neutral gaz
+        if (noms(iq).eq."Ar_N2") then
+          igcm_ar_n2=iq
+          mmol(igcm_ar_n2)=30.
+          count=count+1
+        endif
+	
+
+      enddo ! of do iq=1,nqmx
+!      count=count+nbqchem
+      
+      ! check that we identified all tracers:
+      if (count.ne.nqmx) then
+        write(*,*) "initracer: found only ",count," tracers"
+        write(*,*) "               expected ",nqmx
+        do iq=1,count
+          write(*,*)'      ',iq,' ',trim(noms(iq))
+        enddo
+        stop
+      else
+        write(*,*) "initracer: found all expected tracers, namely:"
+        do iq=1,nqmx
+          write(*,*)'      ',iq,' ',trim(noms(iq))
+        enddo
+      endif
+
+      ! if water cycle but iceparty=.false., there will nevertheless be
+      ! water ice at the surface (iceparty is not used anymore, but this
+      ! part is still relevant, as we want to stay compatible with the
+      ! older versions).
+      if (water.and.(igcm_h2o_ice.eq.0)) then
+        igcm_h2o_ice=igcm_h2o_vap ! so that qsurf(i_h2o_ice) is identified
+                                  ! even though there is no q(i_h2o_ice)
+      else
+       ! surface ice qsurf(i_h2o_ice) was loaded twice by phyetat0,
+       ! as qsurf(i_h2o_vap) & as qsurf(i_h2o_ice), so to be clean:
+       if (igcm_h2o_vap.ne.0) then
+         qsurf(1:ngridmx,igcm_h2o_vap)=0
+       endif
+      endif
+
+c------------------------------------------------------------
+c     Initialisation tracers ....
+c------------------------------------------------------------
+      call zerophys(nqmx,rho_q)
+
+      rho_dust=2500.  ! Mars dust density (kg.m-3)
+      rho_ice=920.    ! Water ice density (kg.m-3)
+      nuice_ref=0.1   ! Effective variance nueff of the
+                      ! water-ice size distributions
+
+      if (doubleq) then
+c       "doubleq" technique 
+c       -------------------
+c      (transport of mass and number mixing ratio)
+c       iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
+
+        if( (nqmx.lt.2).or.(water.and.(nqmx.lt.4)) ) then
+          write(*,*)'initracer: nqmx is too low : nqmx=', nqmx
+          write(*,*)'water= ',water,' doubleq= ',doubleq   
+        end if
+
+        nueff_lift = 0.5
+        varian=sqrt(log(1.+nueff_lift))
+
+        rho_q(igcm_dust_mass)=rho_dust
+        rho_q(igcm_dust_number)=rho_dust
+
+c       Intermediate calcul for computing geometric mean radius r0
+c       as a function of mass and number mixing ratio Q and N
+c       (r0 = (r3n_q * Q/ N)^(1/3))
+        r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust)
+
+c       Intermediate calcul for computing effective radius reff
+c       from geometric mean radius r0
+c       (reff = ref_r0 * r0)
+        ref_r0 = exp(2.5*varian**2)
+        
+c       lifted dust :
+c       '''''''''''
+        reff_lift = 3.0e-6 !3.e-6 !Effective radius of lifted dust (m)
+        alpha_devil(igcm_dust_mass)=9.e-9   !  dust devil lift mass coeff
+c       alpha_lift(igcm_dust_mass)=3.0e-15  !  Lifted mass coeff
+        alpha_lift(igcm_dust_mass)= 3.3e-10!1.e-6 !Lifted mass coeff
+
+        r0_lift = reff_lift/ref_r0
+        alpha_devil(igcm_dust_number)=r3n_q*
+     &                        alpha_devil(igcm_dust_mass)/r0_lift**3
+        alpha_lift(igcm_dust_number)=r3n_q*
+     &                        alpha_lift(igcm_dust_mass)/r0_lift**3
+
+        radius(igcm_dust_mass) = reff_lift
+        radius(igcm_dust_number) = reff_lift
+
+        write(*,*) "initracer: doubleq_param reff_lift:", reff_lift
+        write(*,*) "initracer: doubleq_param nueff_lift:", nueff_lift
+        write(*,*) "initracer: doubleq_param alpha_lift:",
+     &    alpha_lift(igcm_dust_mass)
+
+      else
+
+       if (dustbin.gt.1) then
+        print*,'initracer: STOP!',
+     $   ' properties of dust need to be set in initracer !!!'
+        stop
+
+       else if (dustbin.eq.1) then
+
+c       This will be used for 1 dust particle size:
+c       ------------------------------------------
+        radius(igcm_dustbin(1))=3.e-6
+        alpha_lift(igcm_dustbin(1))=0.0e-6
+        alpha_devil(igcm_dustbin(1))=7.65e-9
+        rho_q(igcm_dustbin(1))=rho_dust
+
+       endif
+      end if    ! (doubleq)
+
+c     Submicron dust mode:
+c     --------------------
+
+      if (submicron) then
+        radius(igcm_dust_submicron)=0.1e-6
+        rho_q(igcm_dust_submicron)=rho_dust
+        if (doubleq) then
+c         If doubleq is also active, we use the population ratio:
+          alpha_lift(igcm_dust_submicron) = 
+     &      alpha_lift(igcm_dust_number)*popratio*
+     &      rho_q(igcm_dust_submicron)*4./3.*pi*
+     &      radius(igcm_dust_submicron)**3.
+          alpha_devil(igcm_dust_submicron)=1.e-30
+        else
+          alpha_lift(igcm_dust_submicron)=1e-6
+          alpha_devil(igcm_dust_submicron)=1.e-30
+        endif ! (doubleq)
+      end if  ! (submicron)
+
+c     Initialization for photochemistry:
+c     ---------------------------------
+      if (photochem) then
+      ! initialize chemistry+water (water will be correctly initialized below)
+      ! by initializing everything which is not dust ...
+        do iq=1,nqmx
+          txt=noms(iq)
+          if (txt(1:4).ne."dust") then
+            radius(iq)=0.
+            alpha_lift(iq) =0.
+            alpha_devil(iq)=0.
+          endif
+        enddo ! do iq=1,nqmx
+      endif
+
+c     Initialization for water vapor
+c     ------------------------------
+      if(water) then
+         radius(igcm_h2o_vap)=0.
+         alpha_lift(igcm_h2o_vap) =0.
+         alpha_devil(igcm_h2o_vap)=0.
+         if(water.and.(nqmx.ge.2)) then
+           radius(igcm_h2o_ice)=3.e-6
+           rho_q(igcm_h2o_ice)=rho_ice
+           alpha_lift(igcm_h2o_ice) =0.
+           alpha_devil(igcm_h2o_ice)=0.
+         elseif(water.and.(nqmx.lt.2)) then
+            write(*,*) 'nqmx is too low : nqmx=', nqmx
+            write(*,*) 'water= ',water
+         endif
+
+      end if  ! (water)
+
+c     Output for records:
+c     ~~~~~~~~~~~~~~~~~~
+      write(*,*)
+      Write(*,*) '******** initracer : dust transport parameters :'
+      write(*,*) 'alpha_lift = ', alpha_lift
+      write(*,*) 'alpha_devil = ', alpha_devil
+      write(*,*) 'radius  = ', radius
+      if(doubleq) then
+        write(*,*) 'reff_lift (um) =  ', reff_lift
+        write(*,*) 'size distribution variance  = ', varian
+        write(*,*) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
+      end if
+
+!
+!     some extra (possibly redundant) sanity checks for tracers: 
+!     ---------------------------------------------------------
+
+       if (doubleq) then 
+       ! verify that we indeed have dust_mass and dust_number tracers 
+         if (igcm_dust_mass.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use doubleq option without ",
+     &                "a dust_mass tracer !"
+           stop
+         endif
+         if (igcm_dust_number.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use doubleq option without ",
+     &                "a dust_number tracer !"
+           stop
+         endif
+       endif
+
+       if ((.not.doubleq).and.(dustbin.gt.0)) then
+       ! verify that we indeed have 'dustbin' dust tracers
+         count=0
+         do iq=1,dustbin
+           if (igcm_dustbin(iq).ne.0) then
+             count=count+1
+           endif
+         enddo
+         if (count.ne.dustbin) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  dusbin is set to ",dustbin,
+     &                " but we only have the following dust tracers:"
+           do iq=1,count
+             write(*,*)"   ",trim(noms(igcm_dustbin(iq)))
+           enddo
+           stop
+         endif
+       endif
+
+       if (water) then
+       ! verify that we indeed have h2o_vap and h2o_ice tracers
+         if (igcm_h2o_vap.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use water option without ",
+     &                "an h2o_vap tracer !"
+           stop
+         endif
+         if (igcm_h2o_ice.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use water option without ",
+     &                "an h2o_ice tracer !"
+           stop
+         endif
+       endif
+
+       if (photochem .or. callthermos) then
+       ! verify that we indeed have the chemistry tracers
+         if (igcm_co2.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a co2 tracer !"
+         stop
+         endif
+         if (igcm_co.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a co tracer !"
+         stop
+         endif
+         if (igcm_o.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a o tracer !"
+         stop
+         endif
+         if (igcm_o1d.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a o1d tracer !"
+         stop
+         endif
+         if (igcm_o2.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "an o2 tracer !"
+         stop
+         endif
+         if (igcm_o3.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "an o3 tracer !"
+         stop
+         endif
+         if (igcm_h.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a h tracer !"
+         stop
+         endif
+         if (igcm_h2.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a h2 tracer !"
+         stop
+         endif
+         if (igcm_oh.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "an oh tracer !"
+         stop
+         endif
+         if (igcm_ho2.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a ho2 tracer !"
+         stop
+         endif
+         if (igcm_h2o2.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a h2o2 tracer !"
+         stop
+         endif
+         if (igcm_n2.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "a n2 tracer !"
+         stop
+         endif
+         if (igcm_ar.eq.0) then
+           write(*,*) "initracer: error !!"
+           write(*,*) "  cannot use chemistry option without ",
+     &                "an ar tracer !"
+         stop
+         endif
+       endif ! of if (photochem .or. callthermos)
+
+      end
Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/initracer.F
===================================================================
--- /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/initracer.F	(revision 309)
+++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/initracer.F	(revision 310)
@@ -1,1 +1,1 @@
-link ../../../mars_lmd_new/libf/phymars/initracer.F
+link STORM_JULIEN_LAST/initracer.F
