Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/vdifc.F
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/vdifc.F	(revision 166)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/vdifc.F	(revision 167)
@@ -1,1 +1,1 @@
-link ../../../mars_lmd_new/libf/phymars/vdifc.F
+link vdifc.F.modif
Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/vdifc.F.modif
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/vdifc.F.modif	(revision 167)
+++ trunk/MESOSCALE/LMD_MM_MARS/SRC/WRFV2/mars_lmd_new_storm/libf/phymars/vdifc.F.modif	(revision 167)
@@ -0,0 +1,652 @@
+      SUBROUTINE vdifc(ngrid,nlay,nq,co2ice,ppopsk,
+     $                ptimestep,pcapcal,lecrit,
+     $                pplay,pplev,pzlay,pzlev,pz0,
+     $                pu,pv,ph,pq,ptsrf,pemis,pqsurf,
+     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
+     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
+     $                pdqdif,pdqsdif)
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   subject:
+c   --------
+c   Turbulent diffusion (mixing) for potential T, U, V and tracer
+c
+c   Shema implicite
+c   On commence par rajouter au variables x la tendance physique
+c   et on resoult en fait:
+c      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
+c
+c   author:
+c   ------
+c      Hourdin/Forget/Fournier
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "comcstfi.h"
+#include "callkeys.h"
+#include "surfdat.h"
+#include "comgeomfi.h"
+#include "tracer.h"
+
+#include "watercap.h"
+c
+c   arguments:
+c   ----------
+
+      INTEGER ngrid,nlay
+      REAL ptimestep
+      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
+      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+      REAL ptsrf(ngrid),pemis(ngrid)
+      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
+      REAL pfluxsrf(ngrid)
+      REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
+      REAL pdtsrf(ngrid),pcapcal(ngrid)
+      REAL pq2(ngrid,nlay+1)
+
+c    Argument added for condensation:
+      REAL co2ice (ngrid), ppopsk(ngrid,nlay)
+      logical lecrit
+      REAL pz0
+
+c    Traceurs :
+      integer nq 
+      REAL pqsurf(ngrid,nq)
+      real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) 
+      real pdqdif(ngrid,nlay,nq) 
+      real pdqsdif(ngrid,nq) 
+      
+c   local:
+c   ------
+
+      INTEGER ilev,ig,ilay,nlev
+
+      REAL z4st,zdplanck(ngridmx)
+      REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
+      REAL zcdv(ngridmx),zcdh(ngridmx)
+      REAL zcdv_true(ngridmx),zcdh_true(ngridmx)
+      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
+      REAL zh(ngridmx,nlayermx)
+      REAL ztsrf2(ngridmx)
+      REAL z1(ngridmx),z2(ngridmx)
+      REAL za(ngridmx,nlayermx),zb(ngridmx,nlayermx)
+      REAL zb0(ngridmx,nlayermx)
+      REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
+      REAL zcst1
+      REAL zu2
+
+      EXTERNAL SSUM,SCOPY
+      REAL SSUM
+      LOGICAL firstcall
+      SAVE firstcall
+
+c     variable added for CO2 condensation:
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+      REAL hh , zhcond(ngridmx,nlayermx)
+      REAL latcond,tcond1mb
+      REAL acond,bcond
+      SAVE acond,bcond
+      DATA latcond,tcond1mb/5.9e5,136.27/
+
+c    Tracers :
+c    ~~~~~~~ 
+      INTEGER iq
+      REAL zq(ngridmx,nlayermx,nqmx)
+      REAL zq1temp(ngridmx)
+      REAL rho(ngridmx) ! near surface air density
+      REAL qsat(ngridmx)
+      DATA firstcall/.true./
+
+      REAL kmixmin
+
+c    ** un petit test de coherence
+c       --------------------------
+
+      IF (firstcall) THEN
+         IF(ngrid.NE.ngridmx) THEN
+            PRINT*,'STOP dans vdifc'
+            PRINT*,'probleme de dimensions :'
+            PRINT*,'ngrid  =',ngrid
+            PRINT*,'ngridmx  =',ngridmx
+            STOP
+         ENDIF
+c        To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal)
+         bcond=1./tcond1mb
+         acond=r/latcond
+         PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
+         PRINT*,'          acond,bcond',acond,bcond
+
+        firstcall=.false.
+      ENDIF
+
+
+
+
+
+c-----------------------------------------------------------------------
+c    1. initialisation
+c    -----------------
+
+      nlev=nlay+1
+
+c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
+c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
+c       ----------------------------------------
+
+      DO ilay=1,nlay
+         DO ig=1,ngrid
+            za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g
+         ENDDO
+      ENDDO
+
+      zcst1=4.*g*ptimestep/(r*r)
+      DO ilev=2,nlev-1
+         DO ig=1,ngrid
+            zb0(ig,ilev)=pplev(ig,ilev)*
+     s      (pplev(ig,1)/pplev(ig,ilev))**rcp /
+     s      (ph(ig,ilev-1)+ph(ig,ilev))
+            zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/
+     s      (pplay(ig,ilev-1)-pplay(ig,ilev))
+         ENDDO
+      ENDDO
+      DO ig=1,ngrid
+		 zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
+      ENDDO
+
+c    ** diagnostique pour l'initialisation
+c       ----------------------------------
+
+      IF(lecrit) THEN
+         ig=ngrid/2+1
+         PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
+         DO ilay=1,nlay
+            WRITE(*,'(6f11.5)')
+     s      .01*pplay(ig,ilay),.001*pzlay(ig,ilay),
+     s      pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
+         ENDDO
+         PRINT*,'Pression (mbar) ,altitude (km),zb'
+         DO ilev=1,nlay
+            WRITE(*,'(3f15.7)')
+     s      .01*pplev(ig,ilev),.001*pzlev(ig,ilev),
+     s      zb0(ig,ilev)
+         ENDDO
+      ENDIF
+
+c     Potential Condensation temperature:
+c     -----------------------------------
+
+c     if (callcond) then 
+c       DO ilev=1,nlay
+c         DO ig=1,ngrid
+c           zhcond(ig,ilev) =
+c    &  (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev)
+c         END DO
+c       END DO
+c     else
+        call zerophys(ngrid*nlay,zhcond)
+c     end if
+
+
+c-----------------------------------------------------------------------
+c   2. ajout des tendances physiques
+c      -----------------------------
+
+      DO ilev=1,nlay
+         DO ig=1,ngrid
+            zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
+            zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
+            zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
+            zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev))
+         ENDDO
+      ENDDO
+      if(tracer) then
+        DO iq =1, nq
+         DO ilev=1,nlay
+           DO ig=1,ngrid
+              zq(ig,ilev,iq)=pq(ig,ilev,iq)+pdqfi(ig,ilev,iq)*ptimestep
+           ENDDO
+         ENDDO
+        ENDDO
+      end if
+
+c-----------------------------------------------------------------------
+c   3. schema de turbulence
+c      --------------------
+
+c    ** source d'energie cinetique turbulente a la surface
+c       (condition aux limites du schema de diffusion turbulente
+c       dans la couche limite
+c       ---------------------
+
+      CALL vdif_cd( ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph
+     &             ,zcdv_true,zcdh_true)
+      DO ig=1,ngrid
+        zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
+        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
+        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
+      ENDDO
+
+c    ** schema de diffusion turbulente dans la couche limite
+c       ---------------------------------------------------- 
+
+        CALL vdif_kc(ptimestep,g,pzlev,pzlay
+     &              ,pu,pv,ph,zcdv_true
+     &              ,pq2,zkv,zkh)
+
+      if ((doubleq).and.(ngrid.eq.1)) then
+        kmixmin = 80. !80.! minimum eddy mix coeff in 1D
+        do ilev=1,nlay
+          do ig=1,ngrid
+           zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
+           zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev))
+          end do
+        end do
+      end if
+
+c    ** diagnostique pour le schema de turbulence
+c       -----------------------------------------
+
+      IF(lecrit) THEN
+         PRINT*
+         PRINT*,'Diagnostic for the vertical turbulent mixing'
+         PRINT*,'Cd for momentum and potential temperature'
+
+         PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
+         PRINT*,'Mixing coefficient for momentum and pot.temp.'
+         DO ilev=1,nlay
+            PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
+         ENDDO
+      ENDIF
+
+
+
+
+c-----------------------------------------------------------------------
+c   4. inversion pour l'implicite sur u
+c      --------------------------------
+
+c    ** l'equation est 
+c       u(t+1) =  u(t) + dt * {(du/dt)phys}(t)  +  dt * {(du/dt)difv}(t+1)
+c       avec
+c       /zu/ = u(t) + dt * {(du/dt)phys}(t)   (voir paragraphe 2.)
+c       et
+c       dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1)
+c       donc les entrees sont /zcdv/ pour la condition a la limite sol
+c       et /zkv/ = Ku
+ 
+      CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
+      CALL multipl(ngrid,zcdv,zb0,zb)
+
+      DO ig=1,ngrid
+         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+         zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
+         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+      ENDDO
+
+      DO ilay=nlay-1,1,-1
+         DO ig=1,ngrid
+            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+            zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+
+     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+         ENDDO
+      ENDDO
+
+      DO ig=1,ngrid
+         zu(ig,1)=zc(ig,1)
+      ENDDO
+      DO ilay=2,nlay
+         DO ig=1,ngrid
+            zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
+         ENDDO
+      ENDDO
+
+
+
+
+
+c-----------------------------------------------------------------------
+c   5. inversion pour l'implicite sur v
+c      --------------------------------
+
+c    ** l'equation est 
+c       v(t+1) =  v(t) + dt * {(dv/dt)phys}(t)  +  dt * {(dv/dt)difv}(t+1)
+c       avec
+c       /zv/ = v(t) + dt * {(dv/dt)phys}(t)   (voir paragraphe 2.)
+c       et
+c       dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1)
+c       donc les entrees sont /zcdv/ pour la condition a la limite sol
+c       et /zkv/ = Kv
+
+      DO ig=1,ngrid
+         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+         zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
+         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+      ENDDO
+
+      DO ilay=nlay-1,1,-1
+         DO ig=1,ngrid
+            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+            zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+
+     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+         ENDDO
+      ENDDO
+
+      DO ig=1,ngrid
+         zv(ig,1)=zc(ig,1)
+      ENDDO
+      DO ilay=2,nlay
+         DO ig=1,ngrid
+            zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
+         ENDDO
+      ENDDO
+
+
+
+
+
+c-----------------------------------------------------------------------
+c   6. inversion pour l'implicite sur h sans oublier le couplage
+c      avec le sol (conduction)
+c      ------------------------
+
+c    ** l'equation est 
+c       h(t+1) =  h(t) + dt * {(dh/dt)phys}(t)  +  dt * {(dh/dt)difv}(t+1)
+c       avec
+c       /zh/ = h(t) + dt * {(dh/dt)phys}(t)   (voir paragraphe 2.)
+c       et
+c       dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1)
+c       donc les entrees sont /zcdh/ pour la condition de raccord au sol
+c       et /zkh/ = Kh
+c       -------------
+
+      CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
+      CALL multipl(ngrid,zcdh,zb0,zb)
+
+      DO ig=1,ngrid
+         z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+         zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
+         zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+      ENDDO
+
+      DO ilay=nlay-1,1,-1
+         DO ig=1,ngrid
+            z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $         zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+            zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+
+     $         zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+            zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+         ENDDO
+      ENDDO
+
+c    ** calcul de (d Planck / dT) a la temperature d'interface
+c       ------------------------------------------------------
+
+      z4st=4.*5.67e-8*ptimestep
+      DO ig=1,ngrid
+         zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
+      ENDDO
+
+c    ** calcul de la temperature_d'interface et de sa tendance.
+c       on ecrit que la somme des flux est nulle a l'interface
+c       a t + \delta t,
+c       c'est a dire le flux radiatif a {t + \delta t}
+c       + le flux turbulent a {t + \delta t} 
+c            qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1
+c            (notation K dt = /cpp*b/)        
+c       + le flux dans le sol a t
+c       + l'evolution du flux dans le sol lorsque la temperature d'interface
+c       passe de sa valeur a t a sa valeur a {t + \delta t}.
+c       ----------------------------------------------------
+
+      DO ig=1,ngrid
+         z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1)
+     s     +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
+         z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
+         ztsrf2(ig)=z1(ig)/z2(ig)
+         pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep
+
+c        Modif speciale CO2 condensation:
+c        tconds = 1./(bcond-acond*log(.0095*pplev(ig,1)))
+c        if ((callcond).and.
+c    &      ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then
+c           zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds
+c        else
+            zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig)
+c        end if
+      ENDDO
+
+c    ** et a partir de la temperature au sol on remonte 
+c       -----------------------------------------------
+
+      DO ilay=2,nlay
+         DO ig=1,ngrid
+            hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond
+            zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh
+         ENDDO
+      ENDDO
+
+
+c-----------------------------------------------------------------------
+c   TRACERS
+c   -------
+
+      if(tracer) then
+           PRINT*, 'alphavdifc', alpha_lift(igcm_dust_mass)
+c     Using the wind modified by friction for lifting and  sublimation
+c     ----------------------------------------------------------------
+        DO ig=1,ngrid
+          zu2=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1)
+          zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
+          zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
+        ENDDO
+
+c       Calcul du flux vertical au bas de la premiere couche (dust) :
+c       -----------------------------------------------------------
+        do ig=1,ngridmx  
+          rho(ig) = zb0(ig,1) /ptimestep
+c          zb(ig,1) = 0.
+        end do
+c       Dust lifting:
+        if (lifting) then
+           if (doubleq.AND.submicron) then
+             do ig=1,ngrid
+c              if(co2ice(ig).lt.1) then
+                 pdqsdif(ig,igcm_dust_mass) =
+     &             -alpha_lift(igcm_dust_mass)  
+                 pdqsdif(ig,igcm_dust_number) = 
+     &             -alpha_lift(igcm_dust_number)  
+                 pdqsdif(ig,igcm_dust_submicron) =
+     &             -alpha_lift(igcm_dust_submicron)
+c              end if
+             end do
+           else if (doubleq) then
+            call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
+     &                    pdqsdif)
+!!             do ig=1,ngrid
+           !!! soulevement constant
+!!                 pdqsdif(ig,igcm_dust_mass) =
+!!     &             -alpha_lift(igcm_dust_mass)  
+!!                 pdqsdif(ig,igcm_dust_number) = 
+!!     &             -alpha_lift(igcm_dust_number)  
+!!             end do
+           else if (submicron) then
+             do ig=1,ngrid
+                 pdqsdif(ig,igcm_dust_submicron) =
+     &             -alpha_lift(igcm_dust_submicron)
+             end do
+           else
+            call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,co2ice,
+     &                    pdqsdif)
+           endif !doubleq.AND.submicron
+        else
+           pdqsdif(1:ngrid,1:nq) = 0.
+        end if
+
+c       OU calcul de la valeur de q a la surface (water)  :
+c       ----------------------------------------
+        if (water) then 
+            call watersat(ngridmx,ptsrf,pplev(1,1),qsat)
+        end if
+
+c      Inversion pour l'implicite sur q 
+c       --------------------------------
+        do iq=1,nq
+          CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
+
+          if ((water).and.(iq.eq.igcm_h2o_vap)) then 
+c            This line is required to account for turbulent transport 
+c            from surface (e.g. ice) to mid-layer of atmosphere:
+             CALL multipl(ngrid,zcdv,zb0,zb(1,1))
+             CALL multipl(ngrid,dryness,zb(1,1),zb(1,1)) 
+          else ! (re)-initialize zb(:,1)
+             zb(1:ngrid,1)=0
+          end if
+
+          DO ig=1,ngrid
+               z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+               zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig)
+               zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+          ENDDO
+  
+          DO ilay=nlay-1,2,-1
+               DO ig=1,ngrid
+                z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+
+     $           zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+                zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+
+     $           zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+                zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+               ENDDO
+          ENDDO
+
+          if (water.and.(iq.eq.igcm_h2o_ice)) then
+            ! special case for water ice tracer: do not include
+            ! h2o ice tracer from surface (which is set when handling
+            ! h2o vapour case (see further down).
+            DO ig=1,ngrid
+                z1(ig)=1./(za(ig,1)+zb(ig,1)+
+     $           zb(ig,2)*(1.-zd(ig,2)))
+                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
+     $         zb(ig,2)*zc(ig,2)) *z1(ig)
+            ENDDO
+          else ! general case
+            DO ig=1,ngrid
+                z1(ig)=1./(za(ig,1)+zb(ig,1)+
+     $           zb(ig,2)*(1.-zd(ig,2)))
+                zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+
+     $         zb(ig,2)*zc(ig,2) +
+     $        (-pdqsdif(ig,iq)) *ptimestep) *z1(ig)  !tracer flux from surface
+            ENDDO
+          endif ! of if (water.and.(iq.eq.igcm_h2o_ice))
+  
+          IF ((water).and.(iq.eq.igcm_h2o_vap)) then 
+c           Calculation for turbulent exchange with the surface (for ice)
+            DO ig=1,ngrid
+              zd(ig,1)=zb(ig,1)*z1(ig)
+              zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
+
+              pdqsdif(ig,igcm_h2o_ice)=rho(ig)*dryness(ig)*zcdv(ig)
+     &                       *(zq1temp(ig)-qsat(ig))
+c             write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
+            END DO
+
+            DO ig=1,ngrid
+              if(.not.watercaptag(ig)) then
+                if ((-pdqsdif(ig,igcm_h2o_ice)*ptimestep)
+     &             .gt.pqsurf(ig,igcm_h2o_ice)) then
+c                 write(*,*)'on sublime plus que qsurf!'
+                  pdqsdif(ig,igcm_h2o_ice)=
+     &                         -pqsurf(ig,igcm_h2o_ice)/ptimestep
+c                 write(*,*)'flux vers le sol=',pdqsdif(ig,nq)
+                  z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))
+                  zc(ig,1)=(za(ig,1)*zq(ig,1,igcm_h2o_vap)+
+     $            zb(ig,2)*zc(ig,2) +
+     $            (-pdqsdif(ig,igcm_h2o_ice)) *ptimestep) *z1(ig)
+                  zq1temp(ig)=zc(ig,1)
+                endif   
+              endif ! if (.not.watercaptag(ig))
+c             Starting upward calculations for water :
+               zq(ig,1,igcm_h2o_vap)=zq1temp(ig)
+            ENDDO ! of DO ig=1,ngrid
+          ELSE
+c           Starting upward calculations for simple mixing of tracer (dust)
+            DO ig=1,ngrid
+               zq(ig,1,iq)=zc(ig,1)
+            ENDDO
+          END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap))
+
+          DO ilay=2,nlay
+             DO ig=1,ngrid
+                zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq)
+             ENDDO
+          ENDDO
+        enddo ! of do iq=1,nq
+      end if ! of if(tracer)
+
+c-----------------------------------------------------------------------
+c   8. calcul final des tendances de la diffusion verticale
+c      ----------------------------------------------------
+
+      DO ilev = 1, nlay
+         DO ig=1,ngrid
+            pdudif(ig,ilev)=(    zu(ig,ilev)-
+     $      (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep)    )/ptimestep
+            pdvdif(ig,ilev)=(    zv(ig,ilev)-
+     $      (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep)    )/ptimestep
+            hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ,
+     $           zhcond(ig,ilev))        ! modif co2cond
+            pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep
+         ENDDO
+      ENDDO
+
+    
+      if (tracer) then 
+        DO iq = 1, nq
+          DO ilev = 1, nlay
+            DO ig=1,ngrid
+              pdqdif(ig,ilev,iq)=(zq(ig,ilev,iq)-
+     $      (pq(ig,ilev,iq) + pdqfi(ig,ilev,iq)*ptimestep))/ptimestep
+            ENDDO
+          ENDDO
+        ENDDO
+      end if
+
+c    ** diagnostique final 
+c       ------------------
+
+      IF(lecrit) THEN
+         PRINT*,'In vdif'
+         PRINT*,'Ts (t) and Ts (t+st)'
+         WRITE(*,'(a10,3a15)')
+     s   'theta(t)','theta(t+dt)','u(t)','u(t+dt)'
+         PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
+         DO ilev=1,nlay
+            WRITE(*,'(4f15.7)')
+     s      ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev),
+     s      pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev)
+
+         ENDDO
+      ENDIF
+
+      if (calltherm .and. outptherm) then
+      if (ngrid .eq. 1) then
+        call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',
+     &                       'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)
+        call WRITEDIAGFI(ngrid,'zkh','zkh',
+     &                       'SI',1,zkh)
+      endif
+      endif 
+
+      RETURN
+      END
Index: trunk/MESOSCALE/LMD_MM_MARS/makemeso
===================================================================
--- trunk/MESOSCALE/LMD_MM_MARS/makemeso	(revision 166)
+++ trunk/MESOSCALE/LMD_MM_MARS/makemeso	(revision 167)
@@ -237,6 +237,11 @@
           if [[ "${config}" == "les" ]]
           then 
-            ## voir ci-dessus pour le folder
-            cp ../../SRC/LES$phys/correcfft ./
+            if [[ "${phys}" == "nophys_" ]]
+            then
+              ## voir ci-dessus pour le folder
+              cp ../../SRC/LES$phys/correcfft ./
+            else
+              cp ../../SRC/LES/correcfft ./
+            fi  
             ./correcfft
             \rm correcfft
Index: trunk/MESOSCALE/PLOT/MINIMAL/map_latlon.pro
===================================================================
--- trunk/MESOSCALE/PLOT/MINIMAL/map_latlon.pro	(revision 166)
+++ trunk/MESOSCALE/PLOT/MINIMAL/map_latlon.pro	(revision 167)
@@ -391,5 +391,5 @@
 	max_value=max_contour, $
 	min_value=min_contour, $	
-	c_labels=findgen(n_elements(lev))*1., $
+	c_labels=findgen(n_elements(lev))*0., $
 ;        c_charsize=!P.charsize/2, $
 	color=colcont;75;0;50;75;100tropclair;150;0;, $
