Index: LMDZ5/trunk/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/conf_phys.F90	(revision 1479)
+++ LMDZ5/trunk/libf/phylmd/conf_phys.F90	(revision 1492)
@@ -13,4 +13,5 @@
   subroutine conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, &
                        ok_LES,&
+                       callstats,&
                        solarlong0,seuil_inversion, &
                        fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
@@ -66,4 +67,5 @@
   logical              :: ok_journe, ok_mensuel, ok_instan, ok_hf
   logical              :: ok_LES
+  LOGICAL              :: callstats
   LOGICAL              :: ok_ade, ok_aie, aerosol_couple
   INTEGER              :: flag_aerosol
@@ -79,4 +81,5 @@
   logical,SAVE        :: ok_journe_omp, ok_mensuel_omp, ok_instan_omp, ok_hf_omp        
   logical,SAVE        :: ok_LES_omp   
+  LOGICAL,SAVE        :: callstats_omp
   LOGICAL,SAVE        :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp
   INTEGER, SAVE       :: flag_aerosol_omp
@@ -1418,4 +1421,13 @@
   ok_LES_omp = .false.                                              
   call getin('OK_LES', ok_LES_omp)                                  
+
+!Config Key  = callstats                                               
+!Config Desc = Pour des sorties callstats                                 
+!Config Def  = .false.                                              
+!Config Help = Pour creer le fichier stats contenant les sorties  
+!              stats                                                  
+!                                                                   
+  callstats_omp = .false.                                              
+  call getin('callstats', callstats_omp)                                  
 !
 !Config Key  = ecrit_LES
@@ -1581,4 +1593,5 @@
     ok_hines = ok_hines_omp
     ok_LES = ok_LES_omp
+    callstats = callstats_omp
     ecrit_LES = ecrit_LES_omp
     carbon_cycle_tr = carbon_cycle_tr_omp
Index: LMDZ5/trunk/libf/phylmd/fisrtilp.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/fisrtilp.F	(revision 1479)
+++ 	(revision )
@@ -1,613 +1,0 @@
-!
-! $Id$
-!
-c
-      SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
-     s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
-     s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
-     s                   frac_impa, frac_nucl,
-     s                   prfl, psfl, rhcl, zqta, fraca,
-     s                   ztv, zpspsk, ztla, zthl, iflag_cldcon)
-
-c
-      USE dimphy
-      IMPLICIT none
-c======================================================================
-c Auteur(s): Z.X. Li (LMD/CNRS)
-c Date: le 20 mars 1995
-c Objet: condensation et precipitation stratiforme.
-c        schema de nuage
-c======================================================================
-c======================================================================
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-#include "YOMCST.h"
-#include "tracstoke.h"
-#include "fisrtilp.h"
-#include "iniprint.h"
-c
-c Arguments:
-c
-      REAL dtime ! intervalle du temps (s)
-      REAL paprs(klon,klev+1) ! pression a inter-couche
-      REAL pplay(klon,klev) ! pression au milieu de couche
-      REAL t(klon,klev) ! temperature (K)
-      REAL q(klon,klev) ! humidite specifique (kg/kg)
-      REAL d_t(klon,klev) ! incrementation de la temperature (K)
-      REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
-      REAL d_ql(klon,klev) ! incrementation de l'eau liquide
-      REAL rneb(klon,klev) ! fraction nuageuse
-      REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
-      REAL rhcl(klon,klev) ! humidite relative en ciel clair
-      REAL rain(klon) ! pluies (mm/s)
-      REAL snow(klon) ! neige (mm/s)
-      REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
-      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) 
-      REAL ztv(klon,klev)
-      REAL zqta(klon,klev),fraca(klon,klev) 
-      REAL sigma1(klon,klev),sigma2(klon,klev)
-      REAL qltot(klon,klev),ctot(klon,klev)
-      REAL zpspsk(klon,klev),ztla(klon,klev)
-      REAL zthl(klon,klev)
-
-      logical lognormale(klon)
-
-cAA
-c Coeffients de fraction lessivee : pour OFF-LINE
-c
-      REAL pfrac_nucl(klon,klev)
-      REAL pfrac_1nucl(klon,klev)
-      REAL pfrac_impa(klon,klev)
-c
-c Fraction d'aerosols lessivee par impaction et par nucleation
-c POur ON-LINE
-c
-      REAL frac_impa(klon,klev)
-      REAL frac_nucl(klon,klev)
-      real zct      ,zcl
-cAA
-c
-c Options du programme:
-c
-      REAL seuil_neb ! un nuage existe vraiment au-dela
-      PARAMETER (seuil_neb=0.001)
-
-      INTEGER ninter ! sous-intervals pour la precipitation
-      INTEGER ncoreczq  
-      INTEGER iflag_cldcon
-      PARAMETER (ninter=5)
-      LOGICAL evap_prec ! evaporation de la pluie
-      PARAMETER (evap_prec=.TRUE.)
-      REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
-      logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
-
-      real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
-      real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
-      real erf   
-      REAL qcloud(klon)
-c
-      LOGICAL cpartiel ! condensation partielle
-      PARAMETER (cpartiel=.TRUE.)
-      REAL t_coup
-      PARAMETER (t_coup=234.0)
-c
-c Variables locales:
-c
-      INTEGER i, k, n, kk
-      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
-      REAL zrfl(klon), zrfln(klon), zqev, zqevt
-      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
-      REAL ztglace, zt(klon)
-      INTEGER nexpo ! exponentiel pour glace/eau
-      REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
-      REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
-c
-      LOGICAL appel1er
-      SAVE appel1er
-c$OMP THREADPRIVATE(appel1er)
-c
-c---------------------------------------------------------------
-c
-cAA Variables traceurs:
-cAA  Provisoire !!! Parametres alpha du lessivage
-cAA  A priori on a 4 scavenging # possibles
-c
-      REAL a_tr_sca(4)
-      save a_tr_sca
-c$OMP THREADPRIVATE(a_tr_sca)
-c
-c Variables intermediaires
-c
-      REAL zalpha_tr
-      REAL zfrac_lessi
-      REAL zprec_cond(klon)
-cAA
-      REAL zmair, zcpair, zcpeau
-C     Pour la conversion eau-neige
-      REAL zlh_solid(klon), zm_solid
-cIM 
-cym      INTEGER klevm1
-c---------------------------------------------------------------
-c
-c Fonctions en ligne:
-c
-      REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
-      REAL zzz
-#include "YOETHF.h"
-#include "FCTTRE.h"
-      fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
-      fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
-c
-      DATA appel1er /.TRUE./
-cym
-      zdelq=0.0
-      
-      if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
-      IF (appel1er) THEN
-c
-         PRINT*, 'fisrtilp, ninter:', ninter
-         PRINT*, 'fisrtilp, evap_prec:', evap_prec
-         PRINT*, 'fisrtilp, cpartiel:', cpartiel
-         IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
-          PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
-          PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
-c         CALL abort
-         ENDIF
-         appel1er = .FALSE.
-c
-cAA initialiation provisoire
-       a_tr_sca(1) = -0.5
-       a_tr_sca(2) = -0.5
-       a_tr_sca(3) = -0.5
-       a_tr_sca(4) = -0.5
-c
-cAA Initialisation a 1 des coefs des fractions lessivees 
-c
-!cdir collapse
-      DO k = 1, klev
-       DO i = 1, klon
-          pfrac_nucl(i,k)=1.
-          pfrac_1nucl(i,k)=1.
-          pfrac_impa(i,k)=1.
-       ENDDO 
-      ENDDO 
-
-      ENDIF          !  test sur appel1er
-c
-cMAf Initialisation a 0 de zoliq
-c      DO i = 1, klon
-c         zoliq(i)=0.
-c      ENDDO 
-c Determiner les nuages froids par leur temperature
-c  nexpo regle la raideur de la transition eau liquide / eau glace.
-c
-      ztglace = RTT - 15.0
-      nexpo = 6
-ccc      nexpo = 1
-c
-c Initialiser les sorties:
-c
-!cdir collapse
-      DO k = 1, klev+1
-      DO i = 1, klon
-         prfl(i,k) = 0.0
-         psfl(i,k) = 0.0
-      ENDDO
-      ENDDO
-
-!cdir collapse
-      DO k = 1, klev
-      DO i = 1, klon
-         d_t(i,k) = 0.0
-         d_q(i,k) = 0.0
-         d_ql(i,k) = 0.0
-         rneb(i,k) = 0.0
-         radliq(i,k) = 0.0
-         frac_nucl(i,k) = 1. 
-         frac_impa(i,k) = 1. 
-      ENDDO
-      ENDDO
-      DO i = 1, klon
-         rain(i) = 0.0
-         snow(i) = 0.0
-         zoliq(i)=0.
-c     ENDDO
-c
-c Initialiser le flux de precipitation a zero
-c
-c     DO i = 1, klon
-         zrfl(i) = 0.0
-         zneb(i) = seuil_neb
-      ENDDO
-c
-c
-cAA Pour plus de securite 
-
-      zalpha_tr   = 0.
-      zfrac_lessi = 0.
-
-cAA----------------------------------------------------------
-c
-      ncoreczq=0
-c Boucle verticale (du haut vers le bas)
-c
-cIM : klevm1
-cym      klevm1=klev-1
-      DO 9999 k = klev, 1, -1
-c
-cAA----------------------------------------------------------
-c
-      DO i = 1, klon
-         zt(i)=t(i,k)
-         zq(i)=q(i,k)
-      ENDDO
-c
-c Calculer la varition de temp. de l'air du a la chaleur sensible
-C transporter par la pluie.
-C Il resterait a rajouter cet effet de la chaleur sensible sur les
-C flux de surface, du a la diff. de temp. entre le 1er niveau et la
-C surface.
-C
-      IF(k.LE.klevm1) THEN         
-         DO i = 1, klon
-cIM
-            zmair=(paprs(i,k)-paprs(i,k+1))/RG
-            zcpair=RCPD*(1.0+RVTMP2*zq(i))
-            zcpeau=RCPD*RVTMP2
-            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau
-     $           + zmair*zcpair*zt(i) )
-     $           / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
-C     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
-         ENDDO
-      ENDIF
-c
-c
-c Calculer l'evaporation de la precipitation
-c
-
-
-      IF (evap_prec) THEN
-      DO i = 1, klon
-      IF (zrfl(i) .GT.0.) THEN
-         IF (thermcep) THEN
-           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
-           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
-           zqs(i)=MIN(0.5,zqs(i))
-           zcor=1./(1.-RETV*zqs(i))
-           zqs(i)=zqs(i)*zcor
-         ELSE
-           IF (zt(i) .LT. t_coup) THEN
-              zqs(i) = qsats(zt(i)) / pplay(i,k)
-           ELSE
-              zqs(i) = qsatl(zt(i)) / pplay(i,k)
-           ENDIF
-         ENDIF
-         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
-         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
-     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
-     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
-         zqev = MIN (zqev, zqevt)
-         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
-     .                            /RG/dtime
-
-c pour la glace, on rï¿½vapore toute la prï¿½ip dans la couche du dessous
-c la glace venant de la couche du dessus est simplement dans la couche
-c du dessous.
-
-         IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
-
-         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
-     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
-     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-         zrfl(i) = zrfln(i)
-      ENDIF
-      ENDDO
-      ENDIF
-c
-c Calculer Qs et L/Cp*dQs/dT:
-c
-      IF (thermcep) THEN
-         DO i = 1, klon
-           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
-           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
-           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
-           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
-           zqs(i) = MIN(0.5,zqs(i))
-           zcor = 1./(1.-RETV*zqs(i))
-           zqs(i) = zqs(i)*zcor
-           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
-         ENDDO
-      ELSE
-         DO i = 1, klon
-            IF (zt(i).LT.t_coup) THEN
-               zqs(i) = qsats(zt(i))/pplay(i,k)
-               zdqs(i) = dqsats(zt(i),zqs(i))
-            ELSE
-               zqs(i) = qsatl(zt(i))/pplay(i,k)
-               zdqs(i) = dqsatl(zt(i),zqs(i))
-            ENDIF
-         ENDDO
-      ENDIF
-c
-c Determiner la condensation partielle et calculer la quantite
-c de l'eau condensee:
-c
-
-      IF (cpartiel) THEN
-
-c        print*,'Dans partiel k=',k
-c
-c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
-c   nuageuse a partir des PDF de Sandrine Bony.
-c   rneb  : fraction nuageuse
-c   zqn   : eau totale dans le nuage
-c   zcond : eau condensee moyenne dans la maille.
-c           on prend en compte le rï¿½hauffement qui diminue la partie condensee
-c
-c   Version avec les raqts
-
-         if (iflag_pdf.eq.0) then
-
-           do i=1,klon
-            zdelq = min(ratqs(i,k),0.99) * zq(i)
-            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
-            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
-           enddo
-
-         else
-c
-c   Version avec les nouvelles PDFs.
-           do i=1,klon
-              if(zq(i).lt.1.e-15) then
-                ncoreczq=ncoreczq+1
-                zq(i)=1.e-15
-              endif
-           enddo 
-
-              if (iflag_cldcon>=5) then
-
-                 call cloudth(klon,klev,k,ztv,
-     .           zq,zqta,fraca,
-     .           qcloud,ctot,zpspsk,paprs,ztla,zthl,
-     .           ratqs,zqs,t)
-
-                 do i=1,klon
-                 rneb(i,k)=ctot(i,k)
-                 zqn(i)=qcloud(i)
-                 enddo
-
-              endif
-
-! Pour iflag_cldcon<=4, on prend toujours la lognormale
-! Dans le cas iflag_cldcon=5, on prend systÃ©matiquement la bi-gaussienne
-! Dans le cas iflagÃ_cldcon=6, on prend la lognormale en absence des thermiques
-
-            lognormale(:)=
-     .      iflag_cldcon<=4.or.(iflag_cldcon==6.and.fraca(:,k)<1.e-10)
-            do i=1,klon
-            if (lognormale(i)) then
-            zpdf_sig(i)=ratqs(i,k)*zq(i)
-            zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
-            zpdf_delta(i)=log(zq(i)/zqs(i))
-            zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
-            zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
-            zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
-            zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
-            zpdf_e1(i)=1.-erf(zpdf_e1(i))
-            zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
-            zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
-            zpdf_e2(i)=1.-erf(zpdf_e2(i))
-            endif
-            enddo
-
-            do i=1,klon
-            if (lognormale(i)) then
-            if (zpdf_e1(i).lt.1.e-10) then
-               rneb(i,k)=0.
-               zqn(i)=zqs(i)
-            else
-               rneb(i,k)=0.5*zpdf_e1(i)
-               zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
-            endif
-            endif
-            
-           enddo 
-
-
-        endif ! iflag_pdf
-
-        DO i=1,klon
-           IF (rneb(i,k) .LE. 0.0) THEN
-              zqn(i) = 0.0
-              rneb(i,k) = 0.0
-              zcond(i) = 0.0
-              rhcl(i,k)=zq(i)/zqs(i)
-           ELSE IF (rneb(i,k) .GE. 1.0) THEN
-              zqn(i) = zq(i)
-              rneb(i,k) = 1.0                  
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
-              rhcl(i,k)=1.0
-           ELSE
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
-              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
-           ENDIF
-        ENDDO
-!         do i=1,klon
-!            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
-!            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
-!            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
-!c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
-!c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
-!c  la convection.
-!c  ATTENTION !!! Il va falloir verifier tout ca.
-!            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
-!c           print*,'ZDQS ',zdqs(i)
-!c--Olivier
-!            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
-!            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
-!            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
-!c--fin
-!           ENDDO
-      ELSE
-         DO i = 1, klon
-            IF (zq(i).GT.zqs(i)) THEN
-               rneb(i,k) = 1.0
-            ELSE
-               rneb(i,k) = 0.0
-            ENDIF
-            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
-         ENDDO
-      ENDIF
-c
-      DO i = 1, klon
-         zq(i) = zq(i) - zcond(i)
-c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
-         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-      ENDDO
-c
-c Partager l'eau condensee en precipitation et eau liquide nuageuse
-c
-      DO i = 1, klon
-      IF (rneb(i,k).GT.0.0) THEN
-         zoliq(i) = zcond(i)
-         zrho(i) = pplay(i,k) / zt(i) / RD
-         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
-         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
-         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-         zfice(i) = zfice(i)**nexpo
-         zneb(i) = MAX(rneb(i,k), seuil_neb)
-         radliq(i,k) = zoliq(i)/REAL(ninter+1)
-      ENDIF
-      ENDDO
-c
-      DO n = 1, ninter
-      DO i = 1, klon
-      IF (rneb(i,k).GT.0.0) THEN
-         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
-
-         IF (zneb(i).EQ.seuil_neb) THEN
-             ztot = 0.0
-         ELSE
-c  quantite d'eau a eliminer: zchau
-c  meme chose pour la glace: zfroi
-             if (ptconv(i,k)) then
-                zcl   =cld_lc_con
-                zct   =1./cld_tau_con
-                zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
-     .              *fallvc(zrhol(i)) * zfice(i)
-             else
-                zcl   =cld_lc_lsc
-                zct   =1./cld_tau_lsc
-                zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
-     .              *fallvs(zrhol(i)) * zfice(i)
-             endif
-             zchau    = zct   *dtime/REAL(ninter) * zoliq(i)
-     .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
-             ztot    = zchau    + zfroi
-             ztot    = MAX(ztot   ,0.0)
-         ENDIF
-         ztot    = MIN(ztot,zoliq(i))
-         zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
-         radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
-      ENDIF
-      ENDDO
-      ENDDO
-c
-      DO i = 1, klon
-      IF (rneb(i,k).GT.0.0) THEN
-         d_ql(i,k) = zoliq(i)
-         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
-     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
-      ENDIF
-      IF (zt(i).LT.RTT) THEN
-        psfl(i,k)=zrfl(i)
-      ELSE
-        prfl(i,k)=zrfl(i)
-      ENDIF
-      ENDDO
-c
-c Calculer les tendances de q et de t:
-c
-      DO i = 1, klon
-         d_q(i,k) = zq(i) - q(i,k)
-         d_t(i,k) = zt(i) - t(i,k)
-      ENDDO
-c
-cAA--------------- Calcul du lessivage stratiforme  -------------
-
-      DO i = 1,klon
-c
-         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
-     .                * (paprs(i,k)-paprs(i,k+1))/RG
-         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
-cAA lessivage nucleation LMD5 dans la couche elle-meme
-            if (t(i,k) .GE. ztglace) THEN
-               zalpha_tr = a_tr_sca(3)
-            else
-               zalpha_tr = a_tr_sca(4)
-            endif
-            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
-            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
-            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
-c
-c nucleation avec un facteur -1 au lieu de -0.5
-            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
-            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
-         ENDIF
-c
-      ENDDO      ! boucle sur i
-c
-cAA Lessivage par impaction dans les couches en-dessous
-      DO kk = k-1, 1, -1
-        DO i = 1, klon
-          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
-            if (t(i,kk) .GE. ztglace) THEN
-              zalpha_tr = a_tr_sca(1)
-            else
-              zalpha_tr = a_tr_sca(2)
-            endif
-            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
-            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
-            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
-          ENDIF
-        ENDDO
-      ENDDO
-c
-cAA----------------------------------------------------------
-c                     FIN DE BOUCLE SUR K   
- 9999 CONTINUE
-c
-cAA-----------------------------------------------------------
-c
-c Pluie ou neige au sol selon la temperature de la 1ere couche
-c
-      DO i = 1, klon
-      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
-         snow(i) = zrfl(i)
-         zlh_solid(i) = RLSTT-RLVTT
-      ELSE
-         rain(i) = zrfl(i)
-         zlh_solid(i) = 0.
-      ENDIF
-      ENDDO
-C
-C For energy conservation : when snow is present, the solification
-c latent heat is considered.
-      DO k = 1, klev
-        DO i = 1, klon
-          zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
-          zmair=(paprs(i,k)-paprs(i,k+1))/RG
-          zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
-          d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
-        END DO 
-      END DO
-c
-
-      if (ncoreczq>0) then
-         print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
-      endif
-      RETURN
-      END
Index: LMDZ5/trunk/libf/phylmd/fisrtilp.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 1492)
+++ LMDZ5/trunk/libf/phylmd/fisrtilp.F90	(revision 1492)
@@ -0,0 +1,619 @@
+!
+! $Id$
+!
+!
+SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
+     d_t, d_q, d_ql, rneb, radliq, rain, snow, &
+     pfrac_impa, pfrac_nucl, pfrac_1nucl, &
+     frac_impa, frac_nucl, &
+     prfl, psfl, rhcl, zqta, fraca, &
+     ztv, zpspsk, ztla, zthl, iflag_cldcon)
+
+  !
+  USE dimphy
+  IMPLICIT none
+  !======================================================================
+  ! Auteur(s): Z.X. Li (LMD/CNRS)
+  ! Date: le 20 mars 1995
+  ! Objet: condensation et precipitation stratiforme.
+  !        schema de nuage
+  !======================================================================
+  !======================================================================
+  !ym include "dimensions.h"
+  !ym include "dimphy.h"
+  include "YOMCST.h"
+  include "tracstoke.h"
+  include "fisrtilp.h"
+  !
+  ! Arguments:
+  !
+  REAL dtime ! intervalle du temps (s)
+  REAL paprs(klon,klev+1) ! pression a inter-couche
+  REAL pplay(klon,klev) ! pression au milieu de couche
+  REAL t(klon,klev) ! temperature (K)
+  REAL q(klon,klev) ! humidite specifique (kg/kg)
+  REAL d_t(klon,klev) ! incrementation de la temperature (K)
+  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
+  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
+  REAL rneb(klon,klev) ! fraction nuageuse
+  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
+  REAL rhcl(klon,klev) ! humidite relative en ciel clair
+  REAL rain(klon) ! pluies (mm/s)
+  REAL snow(klon) ! neige (mm/s)
+  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
+  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) 
+  REAL ztv(klon,klev)
+  REAL zqta(klon,klev),fraca(klon,klev) 
+  REAL sigma1(klon,klev),sigma2(klon,klev)
+  REAL qltot(klon,klev),ctot(klon,klev)
+  REAL zpspsk(klon,klev),ztla(klon,klev)
+  REAL zthl(klon,klev)
+
+  logical lognormale(klon)
+
+  !AA
+  ! Coeffients de fraction lessivee : pour OFF-LINE
+  !
+  REAL pfrac_nucl(klon,klev)
+  REAL pfrac_1nucl(klon,klev)
+  REAL pfrac_impa(klon,klev)
+  !
+  ! Fraction d'aerosols lessivee par impaction et par nucleation
+  ! POur ON-LINE
+  !
+  REAL frac_impa(klon,klev)
+  REAL frac_nucl(klon,klev)
+  real zct      ,zcl
+  !AA
+  !
+  ! Options du programme:
+  !
+  REAL seuil_neb ! un nuage existe vraiment au-dela
+  PARAMETER (seuil_neb=0.001)
+
+  INTEGER ninter ! sous-intervals pour la precipitation
+  INTEGER ncoreczq  
+  INTEGER iflag_cldcon
+  PARAMETER (ninter=5)
+  LOGICAL evap_prec ! evaporation de la pluie
+  PARAMETER (evap_prec=.TRUE.)
+  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
+  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
+
+  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
+  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
+  real erf   
+  REAL qcloud(klon)
+  !
+  LOGICAL cpartiel ! condensation partielle
+  PARAMETER (cpartiel=.TRUE.)
+  REAL t_coup
+  PARAMETER (t_coup=234.0)
+  !
+  ! Variables locales:
+  !
+  INTEGER i, k, n, kk
+  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
+  REAL zrfl(klon), zrfln(klon), zqev, zqevt
+  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
+  REAL ztglace, zt(klon)
+  INTEGER nexpo ! exponentiel pour glace/eau
+  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
+  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
+  !
+  LOGICAL appel1er
+  SAVE appel1er
+  !$OMP THREADPRIVATE(appel1er)
+  !
+  !---------------------------------------------------------------
+  !
+  !AA Variables traceurs:
+  !AA  Provisoire !!! Parametres alpha du lessivage
+  !AA  A priori on a 4 scavenging # possibles
+  !
+  REAL a_tr_sca(4)
+  save a_tr_sca
+  !$OMP THREADPRIVATE(a_tr_sca)
+  !
+  ! Variables intermediaires
+  !
+  REAL zalpha_tr
+  REAL zfrac_lessi
+  REAL zprec_cond(klon)
+  !AA
+  REAL zmair, zcpair, zcpeau
+  !     Pour la conversion eau-neige
+  REAL zlh_solid(klon), zm_solid
+  !IM 
+  !ym      INTEGER klevm1
+  !---------------------------------------------------------------
+  !
+  ! Fonctions en ligne:
+  !
+  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
+  REAL zzz
+  include "YOETHF.h"
+  include "FCTTRE.h"
+  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
+  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
+  !
+  DATA appel1er /.TRUE./
+  !ym
+  zdelq=0.0
+
+  print*,'NUAGES4 A. JAM'
+  IF (appel1er) THEN
+     !
+     PRINT*, 'fisrtilp, ninter:', ninter
+     PRINT*, 'fisrtilp, evap_prec:', evap_prec
+     PRINT*, 'fisrtilp, cpartiel:', cpartiel
+     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
+        PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
+        PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
+        !         CALL abort
+     ENDIF
+     appel1er = .FALSE.
+     !
+     !AA initialiation provisoire
+     a_tr_sca(1) = -0.5
+     a_tr_sca(2) = -0.5
+     a_tr_sca(3) = -0.5
+     a_tr_sca(4) = -0.5
+     !
+     !AA Initialisation a 1 des coefs des fractions lessivees 
+     !
+     !cdir collapse
+     DO k = 1, klev
+        DO i = 1, klon
+           pfrac_nucl(i,k)=1.
+           pfrac_1nucl(i,k)=1.
+           pfrac_impa(i,k)=1.
+        ENDDO
+     ENDDO
+
+  ENDIF          !  test sur appel1er
+  !
+  !MAf Initialisation a 0 de zoliq
+  !      DO i = 1, klon
+  !         zoliq(i)=0.
+  !      ENDDO 
+  ! Determiner les nuages froids par leur temperature
+  !  nexpo regle la raideur de la transition eau liquide / eau glace.
+  !
+  ztglace = RTT - 15.0
+  nexpo = 6
+  !cc      nexpo = 1
+  !
+  ! Initialiser les sorties:
+  !
+  !cdir collapse
+  DO k = 1, klev+1
+     DO i = 1, klon
+        prfl(i,k) = 0.0
+        psfl(i,k) = 0.0
+     ENDDO
+  ENDDO
+
+  !cdir collapse
+  DO k = 1, klev
+     DO i = 1, klon
+        d_t(i,k) = 0.0
+        d_q(i,k) = 0.0
+        d_ql(i,k) = 0.0
+        rneb(i,k) = 0.0
+        radliq(i,k) = 0.0
+        frac_nucl(i,k) = 1. 
+        frac_impa(i,k) = 1. 
+     ENDDO
+  ENDDO
+  DO i = 1, klon
+     rain(i) = 0.0
+     snow(i) = 0.0
+     zoliq(i)=0.
+     !     ENDDO
+     !
+     ! Initialiser le flux de precipitation a zero
+     !
+     !     DO i = 1, klon
+     zrfl(i) = 0.0
+     zneb(i) = seuil_neb
+  ENDDO
+  !
+  !
+  !AA Pour plus de securite 
+
+  zalpha_tr   = 0.
+  zfrac_lessi = 0.
+
+  !AA----------------------------------------------------------
+  !
+  ncoreczq=0
+  ! Boucle verticale (du haut vers le bas)
+  !
+  !IM : klevm1
+  !ym      klevm1=klev-1
+  DO k = klev, 1, -1
+     !
+     !AA----------------------------------------------------------
+     !
+     DO i = 1, klon
+        zt(i)=t(i,k)
+        zq(i)=q(i,k)
+     ENDDO
+     !
+     ! Calculer la varition de temp. de l'air du a la chaleur sensible
+     ! transporter par la pluie.
+     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
+     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
+     ! surface.
+     !
+     IF(k.LE.klevm1) THEN         
+        DO i = 1, klon
+           !IM
+           zmair=(paprs(i,k)-paprs(i,k+1))/RG
+           zcpair=RCPD*(1.0+RVTMP2*zq(i))
+           zcpeau=RCPD*RVTMP2
+           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
+                + zmair*zcpair*zt(i) ) &
+                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
+           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
+        ENDDO
+     ENDIF
+     !
+     !
+     ! Calculer l'evaporation de la precipitation
+     !
+
+
+     IF (evap_prec) THEN
+        DO i = 1, klon
+           IF (zrfl(i) .GT.0.) THEN
+              IF (thermcep) THEN
+                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
+                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
+                 zqs(i)=MIN(0.5,zqs(i))
+                 zcor=1./(1.-RETV*zqs(i))
+                 zqs(i)=zqs(i)*zcor
+              ELSE
+                 IF (zt(i) .LT. t_coup) THEN
+                    zqs(i) = qsats(zt(i)) / pplay(i,k)
+                 ELSE
+                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
+                 ENDIF
+              ENDIF
+              zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+              zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
+                   * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+              zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
+                   * RG*dtime/(paprs(i,k)-paprs(i,k+1))
+              zqev = MIN (zqev, zqevt)
+              zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
+                   /RG/dtime
+
+              ! pour la glace, on ré-évapore toute la précip dans la
+              ! couche du dessous
+              ! la glace venant de la couche du dessus est simplement
+              ! dans la couche du dessous.
+
+              IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
+
+              zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
+                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
+                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                   * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+              zrfl(i) = zrfln(i)
+           ENDIF
+        ENDDO
+     ENDIF
+     !
+     ! Calculer Qs et L/Cp*dQs/dT:
+     !
+     IF (thermcep) THEN
+        DO i = 1, klon
+           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
+           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
+           zqs(i) = MIN(0.5,zqs(i))
+           zcor = 1./(1.-RETV*zqs(i))
+           zqs(i) = zqs(i)*zcor
+           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
+        ENDDO
+     ELSE
+        DO i = 1, klon
+           IF (zt(i).LT.t_coup) THEN
+              zqs(i) = qsats(zt(i))/pplay(i,k)
+              zdqs(i) = dqsats(zt(i),zqs(i))
+           ELSE
+              zqs(i) = qsatl(zt(i))/pplay(i,k)
+              zdqs(i) = dqsatl(zt(i),zqs(i))
+           ENDIF
+        ENDDO
+     ENDIF
+     !
+     ! Determiner la condensation partielle et calculer la quantite
+     ! de l'eau condensee:
+     !
+
+     IF (cpartiel) THEN
+
+        !        print*,'Dans partiel k=',k
+        !
+        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
+        !   nuageuse a partir des PDF de Sandrine Bony.
+        !   rneb  : fraction nuageuse
+        !   zqn   : eau totale dans le nuage
+        !   zcond : eau condensee moyenne dans la maille.
+        !  on prend en compte le réchauffement qui diminue la partie
+        ! condensee
+        !
+        !   Version avec les raqts
+
+        if (iflag_pdf.eq.0) then
+
+           do i=1,klon
+              zdelq = min(ratqs(i,k),0.99) * zq(i)
+              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
+              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
+           enddo
+
+        else
+           !
+           !   Version avec les nouvelles PDFs.
+           do i=1,klon
+              if(zq(i).lt.1.e-15) then
+                 ncoreczq=ncoreczq+1
+                 zq(i)=1.e-15
+              endif
+           enddo
+
+           if (iflag_cldcon>=5) then
+
+              call cloudth(klon,klev,k,ztv, &
+                   zq,zqta,fraca, &
+                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+                   ratqs,zqs,t)
+
+              do i=1,klon
+                 rneb(i,k)=ctot(i,k)
+                 zqn(i)=qcloud(i)
+              enddo
+
+           endif
+
+           if (iflag_cldcon <= 4) then
+              lognormale = .true.
+           elseif (iflag_cldcon == 6) then
+              ! lognormale en l'absence des thermiques
+              lognormale = fraca(:,k) < 1e-10
+           else
+              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
+              ! bi-gaussienne
+              lognormale = .false.
+           end if
+
+           do i=1,klon
+              if (lognormale(i)) then
+                 zpdf_sig(i)=ratqs(i,k)*zq(i)
+                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
+                 zpdf_delta(i)=log(zq(i)/zqs(i))
+                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
+                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
+                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
+                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
+                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
+                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
+                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
+                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
+              endif
+           enddo
+
+           do i=1,klon
+              if (lognormale(i)) then
+                 if (zpdf_e1(i).lt.1.e-10) then
+                    rneb(i,k)=0.
+                    zqn(i)=zqs(i)
+                 else
+                    rneb(i,k)=0.5*zpdf_e1(i)
+                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+                 endif
+              endif
+
+           enddo
+
+
+        endif ! iflag_pdf
+
+        DO i=1,klon
+           IF (rneb(i,k) .LE. 0.0) THEN
+              zqn(i) = 0.0
+              rneb(i,k) = 0.0
+              zcond(i) = 0.0
+              rhcl(i,k)=zq(i)/zqs(i)
+           ELSE IF (rneb(i,k) .GE. 1.0) THEN
+              zqn(i) = zq(i)
+              rneb(i,k) = 1.0                  
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+              rhcl(i,k)=1.0
+           ELSE
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+           ENDIF
+        ENDDO
+        !         do i=1,klon
+        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
+        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
+        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
+        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
+        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
+        !c  la convection.
+        !c  ATTENTION !!! Il va falloir verifier tout ca.
+        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+        !c           print*,'ZDQS ',zdqs(i)
+        !c--Olivier
+        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
+        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
+        !c--fin
+        !           ENDDO
+     ELSE
+        DO i = 1, klon
+           IF (zq(i).GT.zqs(i)) THEN
+              rneb(i,k) = 1.0
+           ELSE
+              rneb(i,k) = 0.0
+           ENDIF
+           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
+        ENDDO
+     ENDIF
+     !
+     DO i = 1, klon
+        zq(i) = zq(i) - zcond(i)
+        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
+        zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+     ENDDO
+     !
+     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
+     !
+     DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
+           zoliq(i) = zcond(i)
+           zrho(i) = pplay(i,k) / zt(i) / RD
+           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
+           zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
+           zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+           zfice(i) = zfice(i)**nexpo
+           zneb(i) = MAX(rneb(i,k), seuil_neb)
+           radliq(i,k) = zoliq(i)/REAL(ninter+1)
+        ENDIF
+     ENDDO
+     !
+     DO n = 1, ninter
+        DO i = 1, klon
+           IF (rneb(i,k).GT.0.0) THEN
+              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+
+              IF (zneb(i).EQ.seuil_neb) THEN
+                 ztot = 0.0
+              ELSE
+                 !  quantite d'eau a eliminer: zchau
+                 !  meme chose pour la glace: zfroi
+                 if (ptconv(i,k)) then
+                    zcl   =cld_lc_con
+                    zct   =1./cld_tau_con
+                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
+                         *fallvc(zrhol(i)) * zfice(i)
+                 else
+                    zcl   =cld_lc_lsc
+                    zct   =1./cld_tau_lsc
+                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
+                         *fallvs(zrhol(i)) * zfice(i)
+                 endif
+                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
+                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
+                 ztot    = zchau    + zfroi
+                 ztot    = MAX(ztot   ,0.0)
+              ENDIF
+              ztot    = MIN(ztot,zoliq(i))
+              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
+              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
+           ENDIF
+        ENDDO
+     ENDDO
+     !
+     DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
+           d_ql(i,k) = zoliq(i)
+           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
+                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+        ENDIF
+        IF (zt(i).LT.RTT) THEN
+           psfl(i,k)=zrfl(i)
+        ELSE
+           prfl(i,k)=zrfl(i)
+        ENDIF
+     ENDDO
+     !
+     ! Calculer les tendances de q et de t:
+     !
+     DO i = 1, klon
+        d_q(i,k) = zq(i) - q(i,k)
+        d_t(i,k) = zt(i) - t(i,k)
+     ENDDO
+     !
+     !AA--------------- Calcul du lessivage stratiforme  -------------
+
+     DO i = 1,klon
+        !
+        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
+             * (paprs(i,k)-paprs(i,k+1))/RG
+        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+           !AA lessivage nucleation LMD5 dans la couche elle-meme
+           if (t(i,k) .GE. ztglace) THEN
+              zalpha_tr = a_tr_sca(3)
+           else
+              zalpha_tr = a_tr_sca(4)
+           endif
+           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
+           !
+           ! nucleation avec un facteur -1 au lieu de -0.5
+           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
+           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+        ENDIF
+        !
+     ENDDO      ! boucle sur i
+     !
+     !AA Lessivage par impaction dans les couches en-dessous
+     DO kk = k-1, 1, -1
+        DO i = 1, klon
+           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+              if (t(i,kk) .GE. ztglace) THEN
+                 zalpha_tr = a_tr_sca(1)
+              else
+                 zalpha_tr = a_tr_sca(2)
+              endif
+              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
+              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
+           ENDIF
+        ENDDO
+     ENDDO
+     !
+     !AA----------------------------------------------------------
+     !                     FIN DE BOUCLE SUR K   
+  end DO
+  !
+  !AA-----------------------------------------------------------
+  !
+  ! Pluie ou neige au sol selon la temperature de la 1ere couche
+  !
+  DO i = 1, klon
+     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
+        snow(i) = zrfl(i)
+        zlh_solid(i) = RLSTT-RLVTT
+     ELSE
+        rain(i) = zrfl(i)
+        zlh_solid(i) = 0.
+     ENDIF
+  ENDDO
+  !
+  ! For energy conservation : when snow is present, the solification
+  ! latent heat is considered.
+  DO k = 1, klev
+     DO i = 1, klon
+        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
+        zmair=(paprs(i,k)-paprs(i,k+1))/RG
+        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
+        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
+     END DO
+  END DO
+  !
+
+  if (ncoreczq>0) then
+     print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
+  endif
+
+END SUBROUTINE fisrtilp
Index: LMDZ5/trunk/libf/phylmd/inistats.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/inistats.F	(revision 1492)
+++ LMDZ5/trunk/libf/phylmd/inistats.F	(revision 1492)
@@ -0,0 +1,133 @@
+      subroutine inistats(ierr)
+
+      implicit none
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+#include "comvert.h"
+#include "comconst.h"
+#include "statto.h"
+#include "netcdf.inc"
+
+      integer,intent(out) :: ierr
+      integer :: nid
+      integer :: l,nsteppd
+      real, dimension(llm) ::  sig_s
+      integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
+      real, dimension(istime) :: lt
+      integer :: nvarid
+      real, dimension(llm) :: pseudoalt
+
+      write (*,*) 
+      write (*,*) '                        || STATS ||'
+      write (*,*) 
+      write (*,*) 'daysec',daysec
+      write (*,*) 'dtphys',dtphys
+      nsteppd=nint(daysec/dtphys)
+      write (*,*) 'nsteppd=',nsteppd
+      if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec)
+     &   stop'Dans Instat:  1jour .ne. n pas physiques'
+
+      if(mod(nsteppd,istime).ne.0)
+     &   stop'Dans Instat:  1jour .ne. n*istime pas physiques'
+
+      istats=nsteppd/istime
+      write (*,*) 'istats=',istats
+      write (*,*) 'Storing ',istime,'times per day'
+      write (*,*) 'thus every ',istats,'physical timestep '
+      write (*,*) 
+
+      do l= 1, llm
+         sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
+         pseudoalt(l)=log(preff/presnivs(l))*8.  
+      enddo
+
+      ierr = NF_CREATE("stats.nc",NF_CLOBBER,nid)
+      if (ierr.ne.NF_NOERR) then
+         write (*,*) NF_STRERROR(ierr)
+         stop ""
+      endif
+
+      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_lat)
+      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_lon)
+      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
+      ierr = NF_DEF_DIM (nid, "llmp1", llm+1, idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
+
+      ierr = NF_ENDDEF(nid)
+      call def_var_stats(nid,"Time","Time",
+     &            "hours since 0000-00-0 00:00:00",1,
+     &            idim_time,nvarid,ierr)
+! Time is initialised later by mkstats subroutine
+
+      call def_var_stats(nid,"latitude","latitude",
+     &            "degrees_north",1,idim_lat,nvarid,ierr)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
+#endif
+      call def_var_stats(nid,"longitude","East longitude",
+     &            "degrees_east",1,idim_lon,nvarid,ierr)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+#endif
+
+! Niveaux verticaux, aps et bps
+      ierr = NF_REDEF (nid)
+! presnivs
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"presnivs", NF_DOUBLE, 1,idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"presnivs", NF_FLOAT, 1,idim_llm,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",15,
+     &                        "Vertical levels")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"Pa")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',4,"down")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid, 
+     &                          presnivs(1:llm))
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid, 
+     &                        presnivs(1:llm))
+#endif 
+! Pseudo alts
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,idim_llm,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,idim_llm,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
+#endif 
+!      call def_var_stats(nid,"aps","hybrid pressure at midlayers"," ", 
+!     &            1,idim_llm,nvarid,ierr)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
+!#endif
+
+!      call def_var_stats(nid,"bps","hybrid sigma at midlayers"," ", 
+!     &            1,idim_llm,nvarid,ierr)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
+!#else 
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
+!#endif
+
+      ierr=NF_CLOSE(nid)
+
+      end subroutine inistats
+
Index: LMDZ5/trunk/libf/phylmd/mkstat.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/mkstat.F90	(revision 1492)
+++ LMDZ5/trunk/libf/phylmd/mkstat.F90	(revision 1492)
@@ -0,0 +1,164 @@
+subroutine mkstats(ierr)
+
+ 
+!
+!  This program writes a stats.nc file from sums and sums of squares
+!  to means and standard deviations and also writes netcdf style
+!  file so that the data can be viewed easily.  The data file is
+!  overwritten in place.  
+!  SRL  21 May 1996
+!  Yann W. july 2003
+
+
+implicit none
+
+#include "dimensions.h"
+#include "statto.h"
+#include "netcdf.inc"
+
+integer,parameter :: iip1=iim+1
+integer,parameter :: jjp1=jjm+1
+integer :: ierr,nid,nbvar,i,ndims,lt,nvarid
+integer, dimension(4) :: id,varid,start,size
+integer, dimension(5) :: dimids
+character (len=50) :: name,nameout,units,title
+real, dimension(iip1,jjp1,llm) :: sum3d,square3d,mean3d,sd3d
+real, dimension(iip1,jjp1) :: sum2d,square2d,mean2d,sd2d
+real, dimension(istime) :: time
+real, dimension(jjp1) :: lat
+real, dimension(iip1) :: lon
+real, dimension(llm) :: alt
+logical :: lcopy=.true.
+!integer :: latid,lonid,altid,timeid
+integer :: meanid,sdid
+!integer, dimension(4) :: dimout
+
+! Incrementation of count for the last step, which is not done in wstats
+count(istime)=count(istime)+1
+
+ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
+
+! We catch the id of dimensions of the stats file
+
+ierr= NF_INQ_DIMID(nid,"latitude",id(1))
+ierr= NF_INQ_DIMID(nid,"longitude",id(2))
+ierr= NF_INQ_DIMID(nid,"altitude",id(3))
+ierr= NF_INQ_DIMID(nid,"Time",id(4))
+
+ierr= NF_INQ_VARID(nid,"latitude",varid(1))
+ierr= NF_INQ_VARID(nid,"longitude",varid(2))
+ierr= NF_INQ_VARID(nid,"altitude",varid(3))
+ierr= NF_INQ_VARID(nid,"Time",varid(4))
+
+! Time initialisation
+
+do i=1,istime
+   time(i)=i*24./istime
+#ifdef NC_DOUBLE
+   ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),i,1,time(i))
+#else
+   ierr= NF_PUT_VARA_REAL(nid,varid(4),i,1,time(i))
+#endif
+enddo
+
+! We catche the values of the variables
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(1),lat)
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(2),lon)
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(3),alt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,varid(1),lat)
+         ierr = NF_GET_VAR_REAL(nid,varid(2),lon)
+         ierr = NF_GET_VAR_REAL(nid,varid(3),alt)
+#endif
+
+! We catch the number of variables in the stats file
+ierr = NF_INQ_NVARS(nid,nbvar)
+
+! to catche the "real" number of variables (without the "additionnal variables")
+nbvar=(nbvar-4)/2 
+
+do i=1,nbvar
+   varid=(i-1)*2+5
+
+   ! What's the variable's name?
+   ierr=NF_INQ_VARNAME(nid,varid,name)
+   write(*,*) "OK variable ",name
+   ! Its units?
+   units=" "
+   ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
+   ! Its title?
+   title=" "
+   ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
+   ! Its number of dimensions?   
+   ierr=NF_INQ_VARNDIMS(nid,varid,ndims)
+   ! Its values?
+
+   if(ndims==4) then ! lat, lon, alt & time
+
+!      dimout(1)=lonid
+!      dimout(2)=latid
+!      dimout(3)=altid
+!      dimout(4)=timeid
+
+      size=(/iip1,jjp1,llm,1/)
+      do lt=1,istime
+         start=(/1,1,1,lt/)
+         ! Extraction of the "source" variables
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum3d)
+         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square3d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum3d)
+         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square3d)
+#endif
+         ! Calculation of these variables
+         mean3d=sum3d/count(lt)
+         sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2))
+         ! Writing of the variables
+#ifdef NC_DOUBLE
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean3d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd3d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean3d)
+         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd3d)
+#endif
+      enddo
+
+    else if (ndims.eq.3) then
+
+!      dimout(1)=lonid
+!      dimout(2)=latid
+!      dimout(3)=timeid
+
+      size=(/iip1,jjp1,1,0/)
+      do lt=1,istime
+         start=(/1,1,lt,0/)
+         ! Extraction of the "source" variables
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum2d)
+         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square2d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum2d)
+         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square2d)
+#endif
+         ! Calculation of these variables
+         mean2d=sum2d/count(lt)
+         sd2d=sqrt(max(0.,square2d/count(lt)-mean2d**2))
+         ! Writing of the variables
+#ifdef NC_DOUBLE
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean2d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd2d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean2d)
+         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd2d)
+#endif
+      enddo
+
+    endif 
+enddo
+
+ierr= NF_CLOSE(nid)
+
+end
Index: LMDZ5/trunk/libf/phylmd/orografi_strato.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/orografi_strato.F	(revision 1479)
+++ LMDZ5/trunk/libf/phylmd/orografi_strato.F	(revision 1492)
@@ -2004,9 +2004,9 @@
 
       DO 110 JK=1,NLEV
-      ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2,1) 
+      ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2+1,1) 
       IF(ZPM1R.GE.ZSIGT)THEN
          nktopg=JK
       ENDIF
-      ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2,1) 
+      ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2+1,1) 
       IF(ZPM1R.GE.ZTOP)THEN
          nstra=JK
Index: LMDZ5/trunk/libf/phylmd/phys_output_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 1479)
+++ LMDZ5/trunk/libf/phylmd/phys_output_mod.F90	(revision 1492)
@@ -427,4 +427,6 @@
   type(ctrl_out),save :: o_pres         = ctrl_out((/ 2, 3, 10, 10, 1 /),'pres')
   type(ctrl_out),save :: o_paprs        = ctrl_out((/ 2, 3, 10, 10, 1 /),'paprs')
+  type(ctrl_out),save :: o_mass        = ctrl_out((/ 2, 3, 10, 10, 1 /),'mass')
+
   type(ctrl_out),save :: o_rneb         = ctrl_out((/ 2, 5, 10, 10, 1 /),'rneb')
   type(ctrl_out),save :: o_rnebcon      = ctrl_out((/ 2, 5, 10, 10, 1 /),'rnebcon')
@@ -1057,5 +1059,5 @@
 ! Couplage conv-CL
  IF (iflag_con.GE.3) THEN
-    IF (iflag_coupl.EQ.1) THEN
+    IF (iflag_coupl>=1) THEN
  CALL histdef2d(iff,o_ale_bl%flag,o_ale_bl%name, "ALE BL", "m2/s2") 
  CALL histdef2d(iff,o_alp_bl%flag,o_alp_bl%name, "ALP BL", "m2/s2") 
@@ -1108,4 +1110,5 @@
  CALL histdef3d(iff,o_pres%flag,o_pres%name, "Air pressure", "Pa" )
  CALL histdef3d(iff,o_paprs%flag,o_paprs%name, "Air pressure Inter-Couches", "Pa" )
+ CALL histdef3d(iff,o_mass%flag,o_mass%name, "Masse Couches", "kg/m2" )
  CALL histdef3d(iff,o_rneb%flag,o_rneb%name, "Cloud fraction", "-")
  CALL histdef3d(iff,o_rnebcon%flag,o_rnebcon%name, "Convective Cloud Fraction", "-")
Index: LMDZ5/trunk/libf/phylmd/phys_output_write.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_write.h	(revision 1479)
+++ LMDZ5/trunk/libf/phylmd/phys_output_write.h	(revision 1492)
@@ -104,4 +104,9 @@
      s                   o_psol%name,itau_w,zx_tmp_fi2d)
        ENDIF
+
+       IF (o_mass%flag(iff)<=lev_files(iff)) THEN
+      CALL histwrite_phy(nid_files(iff),o_mass%name,itau_w,zmasse)
+        ENDIF
+
 
        IF (o_qsurf%flag(iff)<=lev_files(iff)) THEN
@@ -691,5 +696,5 @@
 ! Couplage convection-couche limite
       IF (iflag_con.GE.3) THEN
-      IF (iflag_coupl.EQ.1) THEN
+      IF (iflag_coupl>=1) THEN
        IF (o_ale_bl%flag(iff)<=lev_files(iff)) THEN
        CALL histwrite_phy(nid_files(iff),o_ale_bl%name,itau_w,ale_bl)
@@ -698,5 +703,5 @@
        CALL histwrite_phy(nid_files(iff),o_alp_bl%name,itau_w,alp_bl)
        ENDIF
-      ENDIF !iflag_coupl.EQ.1
+      ENDIF !iflag_coupl>=1
       ENDIF !(iflag_con.GE.3)
 
Index: LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1479)
+++ LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1492)
@@ -158,4 +158,8 @@
       save ok_LES                            
 c$OMP THREADPRIVATE(ok_LES)                  
+c
+      LOGICAL callstats ! sortir le fichier stats 
+      save callstats                            
+c$OMP THREADPRIVATE(callstats)                  
 c
       LOGICAL ok_region ! sortir le fichier regional
@@ -1150,4 +1154,5 @@
 !     and 360
 
+      INTEGER ierr
 #include "YOMCST.h"
 #include "YOETHF.h"
@@ -1222,4 +1227,5 @@
      .     ok_instan, ok_hf,
      .     ok_LES,
+     .     callstats,
      .     solarlong0,seuil_inversion,
      .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
@@ -2459,4 +2465,14 @@
       endif
 ! ----------------------------------------------------------------------
+!IM/FH: 2011/02/23 
+! Couplage Thermiques/Emanuel seulement si T<0
+      if (iflag_coupl==2) then
+        print*,'Couplage Thermiques/Emanuel seulement si T<0'
+        do i=1,klon
+           if (t_seri(i,lmax_th(i))>273.) then
+              Ale_bl(i)=0.
+           endif
+        enddo
+      endif
 
          endif
@@ -2834,5 +2850,8 @@
 ! de la convection profonde.
 
+!IM/FH: 2011/02/23 
+! definition des points sur lesquels ls thermiques sont actifs
          if (prt_level>9)write(*,*)'TEST SCHEMA DE NUAGES '
+         ptconvth(:,:)=fm_therm(:,:)>0.
          do k=1,klev
             do i=1,klon
@@ -3695,5 +3714,58 @@
 c====================================================================
 c
-      
+
+c        -----------------------------------------------------------------
+c        WSTATS: Saving statistics
+c        -----------------------------------------------------------------
+c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+c        which can later be used to make the statistic files of the run:
+c        "stats")          only possible in 3D runs !
+
+         
+         IF (callstats) THEN
+
+           call wstats(klon,o_psol%name,"Surface pressure","Pa"
+     &                 ,2,paprs(:,1))
+           call wstats(klon,o_tsol%name,"Surface temperature","K",
+     &                 2,zxtsol)
+           zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
+           call wstats(klon,o_precip%name,"Precip Totale liq+sol",
+     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
+           call wstats(klon,o_plul%name,"Large-scale Precip",
+     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
+           call wstats(klon,o_pluc%name,"Convective Precip",
+     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           call wstats(klon,o_sols%name,"Solar rad. at surf.",
+     &                 "W/m2",2,solsw)
+           call wstats(klon,o_soll%name,"IR rad. at surf.",
+     &                 "W/m2",2,sollw)
+          zx_tmp_fi2d(:) = topsw(:)-toplw(:)
+          call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA",
+     &                 "W/m2",2,zx_tmp_fi2d)
+
+
+
+           call wstats(klon,o_temp%name,"Air temperature","K",
+     &                 3,t_seri)
+           call wstats(klon,o_vitu%name,"Zonal wind","m.s-1",
+     &                 3,u_seri)
+           call wstats(klon,o_vitv%name,"Meridional wind",
+     &                "m.s-1",3,v_seri)
+           call wstats(klon,o_vitw%name,"Vertical wind",
+     &                "m.s-1",3,omega)
+           call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg",
+     &                 3,q_seri)
+ 
+
+
+           IF(lafin) THEN
+             write (*,*) "Writing stats..."
+             call mkstats(ierr)
+           ENDIF
+
+         ENDIF !if callstats
+     
 
       IF (lafin) THEN
Index: LMDZ5/trunk/libf/phylmd/readaerosol.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/readaerosol.F90	(revision 1479)
+++ LMDZ5/trunk/libf/phylmd/readaerosol.F90	(revision 1492)
@@ -7,5 +7,5 @@
 CONTAINS
 
-SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
+SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
 
 !****************************************************************************************
@@ -27,5 +27,6 @@
   ! Input arguments
   CHARACTER(len=7), INTENT(IN) :: name_aero
-  CHARACTER(len=*), INTENT(IN) :: type  ! correspond to aer_type in clesphys.h
+  CHARACTER(len=*), INTENT(IN) :: type  ! actuel, annuel, scenario or preind
+  CHARACTER(len=8), INTENT(IN) :: filename
   INTEGER, INTENT(IN)          :: iyr_in
 
@@ -58,5 +59,5 @@
      ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
      ! pt_out has dimensions (klon, klev_src, 12)
-     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
      
 
@@ -67,5 +68,5 @@
      ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
      ! pt_out has dimensions (klon, klev_src, 12)
-     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
      
   ELSE IF (type == 'annuel') THEN
@@ -76,5 +77,5 @@
      ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 
      ! pt_out has dimensions (klon, klev_src, 12)
-     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
      
   ELSE IF (type == 'scenario') THEN
@@ -86,5 +87,5 @@
         ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
         ! pt_out has dimensions (klon, klev_src, 12)
-        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
         
      ELSE IF (iyr_in .GE. 2100) THEN
@@ -93,5 +94,5 @@
         ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
         ! pt_out has dimensions (klon, klev_src, 12)
-        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
         
      ELSE
@@ -113,5 +114,5 @@
         ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
         ! pt_out has dimensions (klon, klev_src, 12)
-        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
         
         ! If to read two decades:
@@ -125,5 +126,5 @@
            ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month 
            ! pt_2 has dimensions (klon, klev_src, 12)
-           CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2)
+           CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2)
            ! Test for same number of vertical levels
            IF (klev_src /= klev_src2) THEN
@@ -160,5 +161,5 @@
 
   ELSE
-     WRITE(lunout,*)'This option is not implemented : aer_type = ', type
+     WRITE(lunout,*)'This option is not implemented : aer_type = ', type,' name_aero=',name_aero
      CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
   END IF ! type
@@ -168,5 +169,5 @@
 
 
-  SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
+  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
 !****************************************************************************************
 ! Read 12 month aerosol from file and distribute to local process on physical grid. 
@@ -200,4 +201,5 @@
     CHARACTER(len=7), INTENT(IN)          :: varname
     CHARACTER(len=4), INTENT(IN)          :: cyr
+    CHARACTER(len=8), INTENT(IN)          :: filename
 
 ! Output arguments
@@ -213,5 +215,4 @@
 ! Local variables
     CHARACTER(len=30)     :: fname
-    CHARACTER(len=8)      :: filename='aerosols'
     CHARACTER(len=30)     :: cvar
     INTEGER               :: ncid, dimid, varid
@@ -242,7 +243,8 @@
 ! 1) Open file 
 !****************************************************************************************
-       fname = filename//cyr//'.nc'
+! Add suffix to filename
+       fname = trim(filename)//cyr//'.nc'
   
-       WRITE(lunout,*) 'reading ', TRIM(fname)
+       WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
        CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) )
 
@@ -283,14 +285,4 @@
           CALL abort_gcm('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
        END IF
-
-! 1.5) Check number of month in file opened
-!
-!**************************************************************************************************
-       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
-       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
-!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
-       IF (nbr_tsteps /= 12 ) THEN
-         CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)',1)
-       ENDIF
 
 
@@ -335,4 +327,12 @@
 
        IF (new_file) THEN
+! ++) Check number of month in file opened
+!**************************************************************************************************
+       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
+       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
+!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
+       IF (nbr_tsteps /= 12 ) THEN
+         CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)',1)
+       ENDIF
 
 ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
Index: LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90	(revision 1479)
+++ LMDZ5/trunk/libf/phylmd/readaerosol_interp.F90	(revision 1492)
@@ -92,4 +92,6 @@
   LOGICAL,SAVE       :: debug=.FALSE.! Debugging in this subroutine
 !$OMP THREADPRIVATE(vert_interp, debug)
+  CHARACTER(len=8)      :: type
+  CHARACTER(len=8)      :: filename
 
 
@@ -173,5 +175,32 @@
      ! Reading values corresponding to the closest year taking into count the choice of aer_type. 
      ! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
-     CALL readaerosol(name_aero(id_aero), aer_type, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
+     ! If aer_type=mix1 or mix2, the run type and file name depends on the aerosol.
+     IF (aer_type=='preind' .OR. aer_type=='actuel' .OR. aer_type=='annuel' .OR. aer_type=='scenario') THEN
+        ! Standard case
+        filename='aerosols'
+        type=aer_type
+     ELSE IF (aer_type == 'mix1') THEN
+        ! Special case using a mix of decenal sulfate file and annual aerosols(all aerosols except sulfate)
+        IF (name_aero(id_aero) == 'SO4') THEN
+           filename='so4.run '
+           type='scenario'
+        ELSE
+           filename='aerosols'
+           type='annuel'
+        END IF
+     ELSE  IF (aer_type == 'mix2') THEN
+        ! Special case using a mix of decenal sulfate file and natrual aerosols
+        IF (name_aero(id_aero) == 'SO4') THEN
+           filename='so4.run '
+           type='scenario'
+        ELSE
+           filename='aerosols'
+           type='preind'
+        END IF
+     ELSE
+        CALL abort_gcm('readaerosol_interp', 'this aer_type not supported',1)
+     END IF
+
+     CALL readaerosol(name_aero(id_aero), type, filename, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
           psurf_year(:,:,id_aero), load_year(:,:,id_aero))
      IF (.NOT. ALLOCATED(var_year)) THEN
@@ -182,5 +211,6 @@
 
      ! Reading values corresponding to the preindustrial concentrations.
-     CALL readaerosol(name_aero(id_aero), 'preind', iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
+     type='preind'
+     CALL readaerosol(name_aero(id_aero), type, filename, iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
           pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))
 
Index: LMDZ5/trunk/libf/phylmd/statto.h
===================================================================
--- LMDZ5/trunk/libf/phylmd/statto.h	(revision 1492)
+++ LMDZ5/trunk/libf/phylmd/statto.h	(revision 1492)
@@ -0,0 +1,34 @@
+
+!  statto:
+!     This include file controls the production of statistics.
+!     Some variables could be set in a namelist, but it is easier to
+!     do it here since arrays can then be dimensioned using parameters
+!     and values shouldn't have to change too often.   SRL
+
+!     Calculate stats every istats physics timesteps, starting at first
+!     call.  If istats=0 then don't do statistics at all.  Check value
+!     if number of physics timesteps changes.
+	integer istats
+
+!     Calculate itime independent sums and sums of squares,
+!     example, istat=1,istime=1 gives a single time mean
+	integer, parameter :: istime=12
+
+!     Number of 2D and 3D variables on which to do statistics.
+	integer n2dvar, n3dvar
+	parameter (n2dvar = 8, n3dvar = 5)
+
+!     Units for writing stats header and data
+	integer usdata
+
+!     count tab to know the variable record
+        integer count(istime)
+
+!     Record of the number of stores made for each time.
+	integer nstore(istime)
+
+! Size of the "controle" array
+        integer, parameter :: cntrlsize=15
+
+!       common /sttcom/ dummy,nstore,istats,usdata
+        common /sttcom/ nstore,istats,usdata,count
Index: LMDZ5/trunk/libf/phylmd/wstats.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/wstats.F90	(revision 1492)
+++ LMDZ5/trunk/libf/phylmd/wstats.F90	(revision 1492)
@@ -0,0 +1,351 @@
+subroutine wstats(ngrid,nom,titre,unite,dim,px)
+
+implicit none
+
+#include "dimensions.h"
+#include "statto.h"
+#include "netcdf.inc"
+
+integer,intent(in) :: ngrid
+character (len=*),intent(in) :: nom,titre,unite
+integer,intent(in) :: dim
+real, dimension(ngrid,llm),intent(in) :: px
+integer,parameter :: iip1=iim+1
+integer,parameter :: jjp1=jjm+1
+real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
+real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
+character (len=50) :: namebis
+character (len=50), save :: firstvar
+integer :: ierr,varid,nbdim,nid
+integer :: meanid,sdid
+integer, dimension(4)  :: id,start,size
+logical, save :: firstcall=.TRUE.
+integer :: l,i,j,ig0
+integer,save :: index
+
+integer, save :: step=0
+
+
+if (firstcall) then
+   firstcall=.false.
+   firstvar=trim((nom))
+   call inistats(ierr)
+endif
+
+if (firstvar==nom) then ! If we're back to the first variable
+      step = step + 1
+endif
+
+if (mod(step,istats).ne.0) then
+   RETURN
+endif
+
+ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
+
+namebis=trim(nom)
+ierr= NF_INQ_VARID(nid,namebis,meanid)
+
+if (ierr.ne.NF_NOERR) then
+
+   if (firstvar==nom) then 
+      index=1
+      count=0
+   endif
+
+
+!declaration de la variable
+
+! choix du nom des coordonnees
+   ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+   ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+   if (dim.eq.3) then
+      ierr= NF_INQ_DIMID(nid,"altitude",id(3))
+      ierr= NF_INQ_DIMID(nid,"Time",id(4))
+      nbdim=4
+   else if (dim==2) then
+      ierr= NF_INQ_DIMID(nid,"Time",id(3))
+      nbdim=3
+   endif
+
+   write (*,*) "====================="
+   write (*,*) "STATS: creation de ",nom
+   namebis=trim(nom)
+   call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
+   call inivar(nid,meanid,ngrid,dim,index,px,ierr)
+   namebis=trim(nom)//"_sd"
+   call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
+   call inivar(nid,sdid,ngrid,dim,index,px,ierr)
+
+   ierr= NF_CLOSE(nid)
+   return
+
+else
+   namebis=trim(nom)//"_sd"
+   ierr= NF_INQ_VARID(nid,namebis,sdid)
+
+endif
+
+if (firstvar==nom) then 
+   count(index)=count(int(index))+1
+   index=index+1
+   if (index>istime) then
+      index=1
+   endif
+endif
+
+if (count(index)==0) then
+   if (dim.eq.3) then
+      start=(/1,1,1,index/)
+      size=(/iip1,jjp1,llm,1/)
+      mean3d=0
+      sd3d=0
+   else if (dim.eq.2) then
+      start=(/1,1,index,0/)
+      size=(/iip1,jjp1,1,0/)
+      mean2d=0
+      sd2d=0
+   endif
+else
+   if (dim.eq.3) then
+      start=(/1,1,1,index/)
+      size=(/iip1,jjp1,llm,1/)
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean3d)
+      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd3d)
+#else
+      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean3d)
+      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd3d)
+#endif
+      if (ierr.ne.NF_NOERR) then
+         write (*,*) NF_STRERROR(ierr)
+         stop ""
+      endif
+
+   else if (dim.eq.2) then
+      start=(/1,1,index,0/)
+      size=(/iip1,jjp1,1,0/)
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean2d)
+      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd2d)
+#else
+      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean2d)
+      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd2d)
+#endif
+      if (ierr.ne.NF_NOERR) then
+         write (*,*) NF_STRERROR(ierr)
+         stop ""
+      endif
+   endif
+endif
+
+if (dim.eq.3) then
+
+!  Passage variable physique -->  variable dynamique
+
+   DO l=1,llm
+      DO i=1,iip1
+         dx3(i,1,l)=px(1,l)
+         dx3(i,jjp1,l)=px(ngrid,l)
+      ENDDO
+      DO j=2,jjm
+         ig0= 1+(j-2)*iim
+         DO i=1,iim
+            dx3(i,j,l)=px(ig0+i,l)
+         ENDDO
+         dx3(iip1,j,l)=dx3(1,j,l)
+      ENDDO
+   ENDDO
+
+   mean3d(:,:,:)= mean3d(:,:,:)+dx3(:,:,:)
+   sd3d(:,:,:)= sd3d(:,:,:)+dx3(:,:,:)**2
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean3d)
+   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd3d)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean3d)
+   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd3d)
+#endif
+   if (ierr.ne.NF_NOERR) then
+     write (*,*) NF_STRERROR(ierr)
+     stop ""
+   endif
+
+else if (dim.eq.2) then
+
+!    Passage variable physique -->  physique dynamique
+
+  DO i=1,iip1
+     dx2(i,1)=px(1,1)
+     dx2(i,jjp1)=px(ngrid,1)
+  ENDDO
+  DO j=2,jjm
+     ig0= 1+(j-2)*iim
+     DO i=1,iim
+        dx2(i,j)=px(ig0+i,1)
+     ENDDO
+     dx2(iip1,j)=dx2(1,j)
+  ENDDO
+
+   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
+   sd2d(:,:)= sd2d(:,:)+dx2(:,:)**2
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean2d)
+   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd2d)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean2d)
+   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd2d)
+#endif
+   if (ierr.ne.NF_NOERR) then
+     write (*,*) NF_STRERROR(ierr)
+     stop ""
+   endif
+
+endif
+
+ierr= NF_CLOSE(nid)
+
+end subroutine wstats
+
+!======================================================
+subroutine inivar(nid,varid,ngrid,dim,index,px,ierr)
+
+implicit none
+
+include "dimensions.h"
+!!!!!!include "dimphys.h"
+include "netcdf.inc"
+
+integer, intent(in) :: nid,varid,dim,index,ngrid
+real, dimension(ngrid,llm), intent(in) :: px
+integer, intent(out) :: ierr
+
+integer,parameter :: iip1=iim+1
+integer,parameter :: jjp1=jjm+1
+
+integer :: l,i,j,ig0
+integer, dimension(4) :: start,size
+real, dimension(iip1,jjp1,llm) :: dx3
+real, dimension(iip1,jjp1) :: dx2
+
+if (dim.eq.3) then
+
+   start=(/1,1,1,index/)
+   size=(/iip1,jjp1,llm,1/)
+
+!  Passage variable physique -->  variable dynamique
+
+   DO l=1,llm
+      DO i=1,iip1
+         dx3(i,1,l)=px(1,l)
+         dx3(i,jjp1,l)=px(ngrid,l)
+      ENDDO
+      DO j=2,jjm
+         ig0= 1+(j-2)*iim
+         DO i=1,iim
+            dx3(i,j,l)=px(ig0+i,l)
+         ENDDO
+         dx3(iip1,j,l)=dx3(1,j,l)
+      ENDDO
+   ENDDO
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx3)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx3)
+#endif
+
+else if (dim.eq.2) then
+
+      start=(/1,1,index,0/)
+      size=(/iip1,jjp1,1,0/)
+
+!    Passage variable physique -->  physique dynamique
+
+  DO i=1,iip1
+     dx2(i,1)=px(1,1)
+     dx2(i,jjp1)=px(ngrid,1)
+  ENDDO
+  DO j=2,jjm
+     ig0= 1+(j-2)*iim
+     DO i=1,iim
+        dx2(i,j)=px(ig0+i,1)
+     ENDDO
+     dx2(iip1,j)=dx2(1,j)
+  ENDDO
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx2)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx2)
+#endif
+
+endif
+
+end subroutine inivar
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
+
+! This subroutine defines variable 'name' in a (pre-existing and opened)
+! NetCDF file (known from its NetCDF ID 'nid').
+! The number of dimensions 'nbdim' of the variable, as well as the IDs of
+! corresponding dimensions must be set (in array 'dimids').
+! Upon successfull definition of the variable, 'nvarid' contains the
+! NetCDF ID of the variable.
+! The variables' attributes 'title' (Note that 'long_name' would be more
+! appropriate) and 'units' are also set. 
+
+implicit none
+
+#include "netcdf.inc"
+
+integer,intent(in) :: nid ! NetCDF file ID
+character(len=*),intent(in) :: name ! the variable's name
+character(len=*),intent(in) :: title ! 'title' attribute of variable
+character(len=*),intent(in) :: units ! 'units' attribute of variable
+integer,intent(in) :: nbdim ! number of dimensions of the variable
+integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
+                                              ! the variable is defined along
+integer,intent(out) :: nvarid ! NetCDF ID of the variable
+integer,intent(out) :: ierr ! returned NetCDF staus code
+
+! 1. Switch to NetCDF define mode 
+ierr=NF_REDEF(nid)
+
+! 2. Define the variable
+#ifdef NC_DOUBLE
+ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
+#else
+ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
+#endif
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+! 3. Write attributes
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
+                     len_trim(adjustl(title)),adjustl(title))
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
+                     len_trim(adjustl(units)),adjustl(units))
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+! 4. Switch out of NetCDF define mode
+ierr = NF_ENDDEF(nid)
+
+end subroutine def_var_stats
+
