Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 1912)
+++ trunk/LMDZ.MARS/README	(revision 1913)
@@ -2544,2 +2544,10 @@
 drag_noro.F -> drag_noro_mod.F
 gwstress.F -> gwstress_mod.F
+
+== 04/04/2018 == EM
+- Forgotten in previous commit:  gwprofil.F -> gwprofil_mod.F (here also the
+size of an argument, rho, was incorrect in caller orodrag).
+- Turned newsedim.F into a module newsedim_mod.F
+- Adapted co2cloud.F and improvedCO2clouds.F to not use "newunit" to open file
+(it is perfectly legitimate F2008 Fortran, but older compiler such as gfortran
+on local LMD machines are not there yet).
Index: trunk/LMDZ.MARS/libf/phymars/callsedim.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callsedim.F	(revision 1912)
+++ trunk/LMDZ.MARS/libf/phymars/callsedim.F	(revision 1913)
@@ -13,5 +13,6 @@
      &                      igcm_ccnco2_mass, igcm_ccnco2_number,
      &                      igcm_co2_ice
-      USE comcstfi_h
+      USE newsedim_mod, ONLY: newsedim
+      USE comcstfi_h, ONLY: g
       IMPLICIT NONE
 
@@ -37,5 +38,5 @@
 c   -------------
       
-#include "callkeys.h"
+      include "callkeys.h"
 
 c
@@ -371,5 +372,5 @@
           
                call newsedim(ngrid,nlay,1,1,ptimestep,
-     &         pplev,masse,epaisseur,zt,rd(ir),rho_dust,qr(1,1,ir),
+     &         pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
      &         wq,0.5)
 
Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1912)
+++ trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1913)
@@ -18,5 +18,5 @@
      &     rho_dust, nuiceco2_sed, nuiceco2_ref,
      &     rho_ice_co2,r3n_q,rho_ice,nuice_sed
-     
+      USE newsedim_mod, ONLY: newsedim
       IMPLICIT NONE
 
@@ -156,5 +156,6 @@
       double precision :: Qtemp
       double precision :: ltemp1(10000),ltemp2(10000)
-      integer :: nelem,lebon1,lebon2,uQext
+      integer :: nelem,lebon1,lebon2
+      integer,parameter :: uQext=555
       DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
       DOUBLE PRECISION Qext1bins2(ngrid,nlay)   
@@ -254,5 +255,6 @@
            STOP
         endif
-        open(newunit=uQext,file=trim(datafile)//
+!        open(newunit=uQext,file=trim(datafile)//
+        open(unit=uQext,file=trim(datafile)//
      &       '/optprop_co2ice_1mic.dat'
      &       ,FORM='formatted')
Index: trunk/LMDZ.MARS/libf/phymars/gwprofil.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/gwprofil.F	(revision 1912)
+++ 	(revision )
@@ -1,214 +1,0 @@
-      SUBROUTINE GWPROFIL
-     *         ( klon, klev
-     *         , kgwd ,kdx  , ktest
-     *         , KKCRIT, KKCRITH, KCRIT ,  kkenvh, kknu,kknu2
-     *         , PAPHM1, PRHO   , PSTAB , PTFR , PVPH , PRI , PTAU
-     *         , ptauf ,pdmod   , pnu   , psig ,pgamma, pvar      )
-
-C**** *GWPROFIL*
-C
-C     PURPOSE.
-C     --------
-C
-C**   INTERFACE.
-C     ----------
-C          FROM *GWDRAG*
-C
-C        EXPLICIT ARGUMENTS :
-C        --------------------
-C     ==== INPUTS ===
-C     ==== OUTPUTS ===
-C
-C        IMPLICIT ARGUMENTS :   NONE
-C        --------------------
-C
-C     METHOD:
-C     -------
-C     THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
-C     IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
-C     AND THE TOP OF THE BLOCKED LAYER (KKENVH).
-C     IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE 
-C     BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR 
-C     NONLINEAR GRAVITY WAVE BREAKING.
-C     ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
-C     LEVEL (KCRIT) OR WHEN IT BREAKS.
-C     
-C
-C
-C     EXTERNALS.
-C     ----------
-C
-C
-C     REFERENCE.
-C     ----------
-C
-C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
-C
-C     AUTHOR.
-C     -------
-C
-C     MODIFICATIONS.
-C     --------------
-C     PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
-C-----------------------------------------------------------------------
-      use dimradmars_mod, only: ndlo2
-      implicit none
-C
-
-C
-
-      integer klon,klev,kidia,kfdia
-#include "yoegwd.h"
-
-C-----------------------------------------------------------------------
-C
-C*       0.1   ARGUMENTS
-C              ---------
-C
-      integer kgwd
-      INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2)
-     *       ,kdx(NDLO2),ktest(NDLO2)
-     *       ,kkenvh(NDLO2),kknu(NDLO2),kknu2(NDLO2)
-C
-      REAL PAPHM1(NDLO2,klev+1), PSTAB(NDLO2,klev+1),
-     *     PRHO  (NDLO2,klev+1), PVPH (NDLO2,klev+1),
-     *     PRI   (NDLO2,klev+1), PTFR (NDLO2), PTAU(NDLO2,klev+1),
-     *     ptauf (NDLO2,klev+1)
-     
-      REAL pdmod (NDLO2) , pnu (NDLO2) , psig(NDLO2),
-     *     pgamma(NDLO2) , pvar(NDLO2)
-     
-C-----------------------------------------------------------------------
-C
-C*       0.2   LOCAL ARRAYS
-C              ------------
-C
-c   declarations pour 'implicit none"
-      real zsqr,zalfa,zriw,zalpha,zb,zdel,zdz2n,zdelp,zdelpt
-
-      integer ji,jk,jl,ilevh
-      REAL ZDZ2 (NDLO2,klev) , ZNORM(NDLO2) , zoro(NDLO2)
-      REAL ZTAU (NDLO2,klev+1)
-C
-C-----------------------------------------------------------------------
-C
-C*         1.    INITIALIZATION
-C                --------------
-
-
-      kidia=1
-      kfdia=klon
-
- 100  CONTINUE
-C
-C
-C*    COMPUTATIONAL CONSTANTS.
-C     ------------- ----------
-C
-      ilevh=KLEV/3
-C
-      DO 400 ji=1,kgwd
-      jl=kdx(ji)
-      Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0)
-      ZTAU(JL,KKNU(JL)+1)=PTAU(JL,KKNU(JL)+1)
-      ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1)
-  400 CONTINUE
-C
-      DO 430 JK=KLEV,2,-1
-C
-C
-C*         4.1    CONSTANT WAVE STRESS UNTIL TOP OF THE
-C                 BLOCKING LAYER.
-  410 CONTINUE
-C
-      DO 411 ji=1,kgwd
-      jl=kdx(ji)
-           IF(JK.GE.KKNU2(JL)) THEN
-           PTAU(JL,JK)=ZTAU(JL,KLEV+1)
-           ENDIF
- 411  CONTINUE             
-C
-C*         4.15   CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
-C                 LOW LEVEL FLOW LAYER.
- 415  CONTINUE
-C        
-C
-C*         4.2    WAVE DISPLACEMENT AT NEXT LEVEL.
-C
-  420 CONTINUE
-C
-      DO 421 ji=1,kgwd
-      jl=kdx(ji)
-      IF(JK.LT.KKNU2(JL)) THEN
-      ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK)
-     *                                                    *zoro(jl)
-      ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec)
-      ENDIF
-  421 CONTINUE
-C
-C*         4.3    WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
-C*                AND STRESS:  BREAKING EVALUATION AND CRITICAL 
-C                 LEVEL
-C
-      DO 431 ji=1,kgwd
-      jl=kdx(ji)
-          IF(JK.LT.KKNU2(JL)) THEN
-          IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN
-            PTAU(JL,JK)=0.0
-          ELSE
-               ZSQR=SQRT(PRI(JL,JK))
-               ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK)
-               ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2
-               IF(ZRIW.LT.GRCRIT) THEN
-                 ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
-                 ZB=1./GRCRIT+2./ZSQR
-                 ZALPHA=0.5*(-ZB+SQRT(ZDEL))
-                 ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK)
-                 PTAU(JL,JK)=ZNORM(JL)*ZDZ2N
-               ELSE
-                 PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
-               ENDIF
-            PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1))
-          ENDIF
-          ENDIF
-  431 CONTINUE
-  
-  430 CONTINUE
-  440 CONTINUE
-  
-c     write(*,*) 'ptau'
-c     write(*,99) ((ji,ilevh,ptau(ji,ilevh),ji=1,NDLO2),
-c    .                  ilevh=1,klev+1)
- 99   FORMAT(i3,i3,f15.5)
-
-
-C  REORGANISATION OF THE STRESS PROFILE
-C  IF BREAKING OCCURS AT LOW LEVEL:
-
-      DO 530 ji=1,kgwd
-      jl=kdx(ji)
-      ZTAU(JL,KKENVH(JL))=PTAU(JL,KKENVH(JL))
-      ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL))
- 530  CONTINUE      
-
-      DO 531 JK=1,KLEV
-      
-      DO 532 ji=1,kgwd
-      jl=kdx(ji)
-                
-         IF(JK.GT.KKCRITH(JL).AND.JK.LT.KKENVH(JL))THEN
-
-          ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKENVH(JL))
-          ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KKENVH(JL))
-          PTAU(JL,JK)=ZTAU(JL,KKENVH(JL)) +
-     .                (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KKENVH(JL)) )*
-     .                ZDELP/ZDELPT
-     
-        ENDIF
-            
- 532  CONTINUE    
- 
- 531  CONTINUE        
-
-      RETURN
-      END
Index: trunk/LMDZ.MARS/libf/phymars/gwprofil_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/gwprofil_mod.F	(revision 1913)
+++ trunk/LMDZ.MARS/libf/phymars/gwprofil_mod.F	(revision 1913)
@@ -0,0 +1,221 @@
+      MODULE gwprofil_mod
+      
+      IMPLICIT NONE
+      
+      CONTAINS
+      
+      SUBROUTINE GWPROFIL
+     *         ( klon, klev
+     *         , kgwd ,kdx  , ktest
+     *         , KKCRIT, KKCRITH, KCRIT ,  kkenvh, kknu,kknu2
+     *         , PAPHM1, PRHO   , PSTAB , PTFR , PVPH , PRI , PTAU
+     *         , ptauf ,pdmod   , pnu   , psig ,pgamma, pvar      )
+
+C**** *GWPROFIL*
+C
+C     PURPOSE.
+C     --------
+C
+C**   INTERFACE.
+C     ----------
+C          FROM *GWDRAG*
+C
+C        EXPLICIT ARGUMENTS :
+C        --------------------
+C     ==== INPUTS ===
+C     ==== OUTPUTS ===
+C
+C        IMPLICIT ARGUMENTS :   NONE
+C        --------------------
+C
+C     METHOD:
+C     -------
+C     THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
+C     IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
+C     AND THE TOP OF THE BLOCKED LAYER (KKENVH).
+C     IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE 
+C     BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR 
+C     NONLINEAR GRAVITY WAVE BREAKING.
+C     ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
+C     LEVEL (KCRIT) OR WHEN IT BREAKS.
+C     
+C
+C
+C     EXTERNALS.
+C     ----------
+C
+C
+C     REFERENCE.
+C     ----------
+C
+C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
+C
+C     AUTHOR.
+C     -------
+C
+C     MODIFICATIONS.
+C     --------------
+C     PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
+C-----------------------------------------------------------------------
+      use dimradmars_mod, only: ndlo2
+      implicit none
+C
+
+C
+
+      integer klon,klev,kidia,kfdia
+#include "yoegwd.h"
+
+C-----------------------------------------------------------------------
+C
+C*       0.1   ARGUMENTS
+C              ---------
+C
+      integer kgwd
+      INTEGER KKCRIT(NDLO2),KKCRITH(NDLO2),KCRIT(NDLO2)
+     *       ,kdx(NDLO2),ktest(NDLO2)
+     *       ,kkenvh(NDLO2),kknu(NDLO2),kknu2(NDLO2)
+C
+      REAL PAPHM1(NDLO2,klev+1), PSTAB(NDLO2,klev+1),
+     *     PRHO  (NDLO2,klev+1), PVPH (NDLO2,klev+1),
+     *     PRI   (NDLO2,klev+1), PTFR (NDLO2), PTAU(NDLO2,klev+1),
+     *     ptauf (NDLO2,klev+1)
+     
+      REAL pdmod (NDLO2) , pnu (NDLO2) , psig(NDLO2),
+     *     pgamma(NDLO2) , pvar(NDLO2)
+     
+C-----------------------------------------------------------------------
+C
+C*       0.2   LOCAL ARRAYS
+C              ------------
+C
+c   declarations pour 'implicit none"
+      real zsqr,zalfa,zriw,zalpha,zb,zdel,zdz2n,zdelp,zdelpt
+
+      integer ji,jk,jl,ilevh
+      REAL ZDZ2 (NDLO2,klev) , ZNORM(NDLO2) , zoro(NDLO2)
+      REAL ZTAU (NDLO2,klev+1)
+C
+C-----------------------------------------------------------------------
+C
+C*         1.    INITIALIZATION
+C                --------------
+
+
+      kidia=1
+      kfdia=klon
+
+ 100  CONTINUE
+C
+C
+C*    COMPUTATIONAL CONSTANTS.
+C     ------------- ----------
+C
+      ilevh=KLEV/3
+C
+      DO 400 ji=1,kgwd
+      jl=kdx(ji)
+      Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0)
+      ZTAU(JL,KKNU(JL)+1)=PTAU(JL,KKNU(JL)+1)
+      ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1)
+  400 CONTINUE
+C
+      DO 430 JK=KLEV,2,-1
+C
+C
+C*         4.1    CONSTANT WAVE STRESS UNTIL TOP OF THE
+C                 BLOCKING LAYER.
+  410 CONTINUE
+C
+      DO 411 ji=1,kgwd
+      jl=kdx(ji)
+           IF(JK.GE.KKNU2(JL)) THEN
+           PTAU(JL,JK)=ZTAU(JL,KLEV+1)
+           ENDIF
+ 411  CONTINUE             
+C
+C*         4.15   CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
+C                 LOW LEVEL FLOW LAYER.
+ 415  CONTINUE
+C        
+C
+C*         4.2    WAVE DISPLACEMENT AT NEXT LEVEL.
+C
+  420 CONTINUE
+C
+      DO 421 ji=1,kgwd
+      jl=kdx(ji)
+      IF(JK.LT.KKNU2(JL)) THEN
+      ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK)
+     *                                                    *zoro(jl)
+      ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec)
+      ENDIF
+  421 CONTINUE
+C
+C*         4.3    WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
+C*                AND STRESS:  BREAKING EVALUATION AND CRITICAL 
+C                 LEVEL
+C
+      DO 431 ji=1,kgwd
+      jl=kdx(ji)
+          IF(JK.LT.KKNU2(JL)) THEN
+          IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN
+            PTAU(JL,JK)=0.0
+          ELSE
+               ZSQR=SQRT(PRI(JL,JK))
+               ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK)
+               ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2
+               IF(ZRIW.LT.GRCRIT) THEN
+                 ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
+                 ZB=1./GRCRIT+2./ZSQR
+                 ZALPHA=0.5*(-ZB+SQRT(ZDEL))
+                 ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK)
+                 PTAU(JL,JK)=ZNORM(JL)*ZDZ2N
+               ELSE
+                 PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
+               ENDIF
+            PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1))
+          ENDIF
+          ENDIF
+  431 CONTINUE
+  
+  430 CONTINUE
+  440 CONTINUE
+  
+c     write(*,*) 'ptau'
+c     write(*,99) ((ji,ilevh,ptau(ji,ilevh),ji=1,NDLO2),
+c    .                  ilevh=1,klev+1)
+ 99   FORMAT(i3,i3,f15.5)
+
+
+C  REORGANISATION OF THE STRESS PROFILE
+C  IF BREAKING OCCURS AT LOW LEVEL:
+
+      DO 530 ji=1,kgwd
+      jl=kdx(ji)
+      ZTAU(JL,KKENVH(JL))=PTAU(JL,KKENVH(JL))
+      ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL))
+ 530  CONTINUE      
+
+      DO 531 JK=1,KLEV
+      
+      DO 532 ji=1,kgwd
+      jl=kdx(ji)
+                
+         IF(JK.GT.KKCRITH(JL).AND.JK.LT.KKENVH(JL))THEN
+
+          ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKENVH(JL))
+          ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KKENVH(JL))
+          PTAU(JL,JK)=ZTAU(JL,KKENVH(JL)) +
+     .                (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KKENVH(JL)) )*
+     .                ZDELP/ZDELPT
+     
+        ENDIF
+            
+ 532  CONTINUE    
+ 
+ 531  CONTINUE        
+
+      END SUBROUTINE GWPROFIL
+
+      END MODULE gwprofil_mod
Index: trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F	(revision 1912)
+++ trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F	(revision 1913)
@@ -167,5 +167,6 @@
       integer nelem,lebon1,lebon2
       double precision :: ltemp1(130),ltemp2(130)
-      integer ibin,uMeteor,j
+      integer ibin,j
+      integer,parameter :: uMeteor=666
 
       IF (firstcall) THEN
@@ -254,5 +255,6 @@
            endif
 !used Variables
-           open(newunit=uMeteor,file=trim(datafile)//
+!           open(newunit=uMeteor,file=trim(datafile)//
+           open(unit=uMeteor,file=trim(datafile)//
      &          '/Meteo_flux_Plane.dat'
      &          ,FORM='formatted')
Index: trunk/LMDZ.MARS/libf/phymars/newsedim.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/newsedim.F	(revision 1912)
+++ 	(revision )
@@ -1,213 +1,0 @@
-      SUBROUTINE newsedim(ngrid,nlay,naersize,nrhosize,ptimestep,
-     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq,beta)
-      USE comcstfi_h
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c      Compute sedimentation of 1 tracer 
-c      of radius rd (m) and density rho (kg.m-3) 
-c
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-c
-c   arguments:
-c   ----------
-
-      INTEGER,INTENT(IN) :: ngrid,nlay,naersize,nrhosize
-      REAL,INTENT(IN) :: ptimestep            ! pas de temps physique (s)
-      REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
-      REAL,INTENT(IN) :: pt(ngrid,nlay) ! temperature au centre des couches (K)
-      real,intent(in) :: masse (ngrid,nlay) ! masse d'une couche (kg)
-      real,intent(in) :: epaisseur (ngrid,nlay)  ! epaisseur d'une couche (m)
-      real,intent(in) :: rd(naersize)             ! particle radius (m)
-      real,intent(in) :: rho(nrhosize)             ! particle density (kg.m-3)
-
-
-c    Traceurs :
-      real,intent(inout) :: pqi(ngrid,nlay)  ! traceur   (e.g. ?/kg)
-      real,intent(out) :: wq(ngrid,nlay+1)  ! flux de traceur durant timestep (?/m-2)
-      real,intent(in) :: beta ! correction for the shape of the particles
-                !   (see Murphy et al. JGR 1990 vol.95)
-                !   beta=1 for spheres
-                !   beta=0.85 for irregular particles
-                !   beta=0.5 for disk shaped particles
-      
-c   local:
-c   ------
-
-      INTEGER l,ig, k, i
-      REAL rfall,rhofall
-
-      LOGICAL,SAVE :: firstcall=.true.
-
-c    Traceurs :
-c    ~~~~~~~~ 
-      real traversee (ngrid,nlay)
-      real vstokes(ngrid,nlay)
-      real w(ngrid,nlay)
-      real ptop, dztop, Ep, Stra
-
-
-c    Physical constant
-c    ~~~~~~~~~~~~~~~~~
-c     Gas molecular viscosity (N.s.m-2)
-      real,parameter :: visc=1.e-5       ! CO2
-c     Effective gas molecular radius (m)
-      real,parameter :: molrad=2.2e-10   ! CO2
-
-c     local and saved variable
-      real,save :: a,b
-
-
-c    ** un petit test de coherence
-c       --------------------------
-
-      ! AS: OK firstcall absolute
-      IF (firstcall) THEN
-
-         firstcall=.false.
-
-
-c       Preliminary calculations for sedimenation velocity :
-c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-c       - Constant to compute stokes speed simple formulae
-c        (Vstokes =  b * rho* r**2   avec   b= (2/9) * rho * g / visc
-         b = 2./9. * g / visc
-      ENDIF ! of IF(firstcall)
-      
-c       - Constant  to compute gas mean free path
-c        l= (T/P)*a, avec a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
-         a = 0.707*8.31/(4*3.1416* molrad*molrad * 6.023e23)
-
-c       - Correction to account for non-spherical shape (Murphy et al.  1990)
-         a = a * beta
-
-
-
-c-----------------------------------------------------------------------
-c    1. initialisation
-c    -----------------
-
-c     Sedimentation velocity (m/s)
-c     ~~~~~~~~~~~~~~~~~~~~~~
-c     (stokes law corrected for low pressure by the Cunningham
-c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
-
-        do  l=1,nlay
-          do ig=1, ngrid
-c           radius
-            if (naersize.eq.1) then 
-              rfall=rd(1)
-            else
-              i=ngrid*(l-1)+ig
-              rfall=rd(i)
-            endif  
-c           density
-            if (nrhosize.eq.1) then 
-              rhofall=rho(1)
-            else
-              i=ngrid*(l-1)+ig
-              rhofall=rho(i)
-            endif  
-c           vstokes
-            vstokes(ig,l) = b * rhofall * rfall*rfall *
-     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
-
-c           Layer crossing time (s) :
-            traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l)
-          end do
-        end do
-
-
-c     Calcul de la masse d'atmosphere correspondant a q transferee
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c     (e.g. on recherche le niveau  en dessous de laquelle le traceur
-c      va traverser le niveau intercouche l : "dztop" est sa hauteur
-c      au dessus de l (m), "ptop" est sa pression (Pa))
-
-      do  l=1,nlay
-        do ig=1, ngrid
-             
-             dztop = vstokes(ig,l)*  ptimestep 
-             Ep=0
-             k=0
-
-            w(ig,l) = 0. !! JF+AS ajout initialisation
-c **************************************************************
-c            Simple Method
-
-cc             w(ig,l) =
-cc     &       (1.- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
-cccc           write(*,*) 'OK simple method l,w =', l, w(ig,l)
-cccc           write(*,*) 'OK simple method dztop =', dztop
-
-             w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
-             !!! Diagnostic: JF. Fix: AS. Date: 05/11 
-             !!! Probleme arrondi avec la quantite ci-dessus
-             !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit
-             !!! ---> dans ce cas on utilise le developpement limite !
-             !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2            
-             IF ( w(ig,l) .eq. 0. ) THEN
-                w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
-             ELSE
-                w(ig,l) = w(ig,l) * pplev(ig,l) / g  
-             ENDIF
-
-
-c **************************************************************
-cccc         Complex method :
-            if (dztop.gt.epaisseur(ig,l)) then                !!!if on traverse plus d'une couche
-cccc            Cas ou on "epuise" la couche l : On calcule le flux
-cccc            Venant de dessus en tenant compte de la variation de Vstokes
-c **************************************************************
-               Ep= epaisseur(ig,l)
-               Stra= traversee(ig,l)
-               do while(dztop.gt.Ep.and.l+k+1.le.nlay)
-                 k=k+1
-                 dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra)
-                 Ep = Ep + epaisseur(ig,l+k)
-                 Stra = Stra + traversee(ig,l+k)
-               enddo 
-               Ep = Ep - epaisseur(ig,l+k)
-             !ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
-
-             !!! JF+AS 05/11 Probleme arrondi potentiel, meme solution que ci-dessus
-             ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
-             IF ( ptop .eq. 1. ) THEN
-                !PRINT*, 'newsedim: exposant trop petit ', ig, l
-                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
-             ELSE
-                ptop=pplev(ig,l+k) * ptop 
-             ENDIF
-             w(ig,l) = (pplev(ig,l) - Ptop)/g
-
-            endif                !!!!!if complex method
-
-
-cc           write(*,*) 'OK new    method l,w =', l, w(ig,l)
-cc           write(*,*) 'OK new    method dztop =', dztop
-cc       if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop
-cc       if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep
-cc            if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l)
-cc            if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l)
-cc            if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l)
-c **************************************************************
-
-
-        end do
-      end do
-
-      call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)
-c         write(*,*) ' newsed: wq(6), wq(7), q(6)',
-c    &                wq(1,6),wq(1,7),pqi(1,6)
-
-
-      RETURN
-      END
-
Index: trunk/LMDZ.MARS/libf/phymars/newsedim_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/newsedim_mod.F	(revision 1913)
+++ trunk/LMDZ.MARS/libf/phymars/newsedim_mod.F	(revision 1913)
@@ -0,0 +1,220 @@
+      MODULE newsedim_mod
+      
+      IMPLICIT NONE
+      
+      CONTAINS
+
+      SUBROUTINE newsedim(ngrid,nlay,naersize,nrhosize,ptimestep,
+     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq,beta)
+      USE comcstfi_h, ONLY: r,g
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c      Compute sedimentation of 1 tracer 
+c      of radius rd (m) and density rho (kg.m-3) 
+c
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+c
+c   arguments:
+c   ----------
+
+      INTEGER,INTENT(IN) :: ngrid,nlay,naersize,nrhosize
+      REAL,INTENT(IN) :: ptimestep            ! pas de temps physique (s)
+      REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
+      REAL,INTENT(IN) :: pt(ngrid,nlay) ! temperature au centre des couches (K)
+      real,intent(in) :: masse (ngrid,nlay) ! masse d'une couche (kg)
+      real,intent(in) :: epaisseur (ngrid,nlay)  ! epaisseur d'une couche (m)
+      real,intent(in) :: rd(naersize)             ! particle radius (m)
+      real,intent(in) :: rho(nrhosize)             ! particle density (kg.m-3)
+
+
+c    Traceurs :
+      real,intent(inout) :: pqi(ngrid,nlay)  ! traceur   (e.g. ?/kg)
+      real,intent(out) :: wq(ngrid,nlay+1)  ! flux de traceur durant timestep (?/m-2)
+      real,intent(in) :: beta ! correction for the shape of the particles
+                !   (see Murphy et al. JGR 1990 vol.95)
+                !   beta=1 for spheres
+                !   beta=0.85 for irregular particles
+                !   beta=0.5 for disk shaped particles
+      
+c   local:
+c   ------
+
+      INTEGER l,ig, k, i
+      REAL rfall,rhofall
+
+      LOGICAL,SAVE :: firstcall=.true.
+
+c    Traceurs :
+c    ~~~~~~~~ 
+      real traversee (ngrid,nlay)
+      real vstokes(ngrid,nlay)
+      real w(ngrid,nlay)
+      real ptop, dztop, Ep, Stra
+
+
+c    Physical constant
+c    ~~~~~~~~~~~~~~~~~
+c     Gas molecular viscosity (N.s.m-2)
+      real,parameter :: visc=1.e-5       ! CO2
+c     Effective gas molecular radius (m)
+      real,parameter :: molrad=2.2e-10   ! CO2
+
+c     local and saved variable
+      real,save :: a,b
+
+
+c    ** un petit test de coherence
+c       --------------------------
+
+      ! AS: OK firstcall absolute
+      IF (firstcall) THEN
+
+         firstcall=.false.
+
+
+c       Preliminary calculations for sedimenation velocity :
+c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+c       - Constant to compute stokes speed simple formulae
+c        (Vstokes =  b * rho* r**2   avec   b= (2/9) * rho * g / visc
+         b = 2./9. * g / visc
+      ENDIF ! of IF(firstcall)
+      
+c       - Constant  to compute gas mean free path
+c        l= (T/P)*a, avec a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
+         a = 0.707*8.31/(4*3.1416* molrad*molrad * 6.023e23)
+
+c       - Correction to account for non-spherical shape (Murphy et al.  1990)
+         a = a * beta
+
+
+
+c-----------------------------------------------------------------------
+c    1. initialisation
+c    -----------------
+
+c     Sedimentation velocity (m/s)
+c     ~~~~~~~~~~~~~~~~~~~~~~
+c     (stokes law corrected for low pressure by the Cunningham
+c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
+
+        do  l=1,nlay
+          do ig=1, ngrid
+c           radius
+            if (naersize.eq.1) then 
+              rfall=rd(1)
+            else
+              i=ngrid*(l-1)+ig
+              rfall=rd(i)
+            endif  
+c           density
+            if (nrhosize.eq.1) then 
+              rhofall=rho(1)
+            else
+              i=ngrid*(l-1)+ig
+              rhofall=rho(i)
+            endif  
+c           vstokes
+            vstokes(ig,l) = b * rhofall * rfall*rfall *
+     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
+
+c           Layer crossing time (s) :
+            traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l)
+          end do
+        end do
+
+
+c     Calcul de la masse d'atmosphere correspondant a q transferee
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c     (e.g. on recherche le niveau  en dessous de laquelle le traceur
+c      va traverser le niveau intercouche l : "dztop" est sa hauteur
+c      au dessus de l (m), "ptop" est sa pression (Pa))
+
+      do  l=1,nlay
+        do ig=1, ngrid
+             
+             dztop = vstokes(ig,l)*  ptimestep 
+             Ep=0
+             k=0
+
+            w(ig,l) = 0. !! JF+AS ajout initialisation
+c **************************************************************
+c            Simple Method
+
+cc             w(ig,l) =
+cc     &       (1.- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
+cccc           write(*,*) 'OK simple method l,w =', l, w(ig,l)
+cccc           write(*,*) 'OK simple method dztop =', dztop
+
+             w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
+             !!! Diagnostic: JF. Fix: AS. Date: 05/11 
+             !!! Probleme arrondi avec la quantite ci-dessus
+             !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit
+             !!! ---> dans ce cas on utilise le developpement limite !
+             !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2            
+             IF ( w(ig,l) .eq. 0. ) THEN
+                w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
+             ELSE
+                w(ig,l) = w(ig,l) * pplev(ig,l) / g  
+             ENDIF
+
+
+c **************************************************************
+cccc         Complex method :
+            if (dztop.gt.epaisseur(ig,l)) then                !!!if on traverse plus d'une couche
+cccc            Cas ou on "epuise" la couche l : On calcule le flux
+cccc            Venant de dessus en tenant compte de la variation de Vstokes
+c **************************************************************
+               Ep= epaisseur(ig,l)
+               Stra= traversee(ig,l)
+               do while(dztop.gt.Ep.and.l+k+1.le.nlay)
+                 k=k+1
+                 dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra)
+                 Ep = Ep + epaisseur(ig,l+k)
+                 Stra = Stra + traversee(ig,l+k)
+               enddo 
+               Ep = Ep - epaisseur(ig,l+k)
+             !ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
+
+             !!! JF+AS 05/11 Probleme arrondi potentiel, meme solution que ci-dessus
+             ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
+             IF ( ptop .eq. 1. ) THEN
+                !PRINT*, 'newsedim: exposant trop petit ', ig, l
+                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
+             ELSE
+                ptop=pplev(ig,l+k) * ptop 
+             ENDIF
+             w(ig,l) = (pplev(ig,l) - Ptop)/g
+
+            endif                !!!!!if complex method
+
+
+cc           write(*,*) 'OK new    method l,w =', l, w(ig,l)
+cc           write(*,*) 'OK new    method dztop =', dztop
+cc       if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop
+cc       if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep
+cc            if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l)
+cc            if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l)
+cc            if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l)
+c **************************************************************
+
+
+        end do
+      end do
+
+      call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)
+c         write(*,*) ' newsed: wq(6), wq(7), q(6)',
+c    &                wq(1,6),wq(1,7),pqi(1,6)
+
+
+      END SUBROUTINE newsedim
+      
+      END MODULE newsedim_mod
+
Index: trunk/LMDZ.MARS/libf/phymars/orodrag_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/orodrag_mod.F	(revision 1912)
+++ trunk/LMDZ.MARS/libf/phymars/orodrag_mod.F	(revision 1913)
@@ -78,4 +78,5 @@
       use dimradmars_mod, only: ndlo2
       USE gwstress_mod, ONLY: gwstress
+      USE gwprofil_mod, ONLY: gwprofil
       USE comcstfi_h, ONLY: g, cpp
       implicit none
@@ -123,5 +124,6 @@
       integer ji,jk,jl,klevm1,ilevp1
 C      real gkwake
-      real ztmst,pvar(NDLO2,4),ztauf,zrtmst,zdelp,zb,zc,zbet
+      real ztmst,pvar(NDLO2,4),ztauf(NDLO2,klev+1)
+      real zrtmst,zdelp,zb,zc,zbet
       real zconb,zabsv,zzd1,ratio,zust,zvst,zdis,ztemp
 C
