Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/atm2geo.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/atm2geo.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/atm2geo.F	(revision 222)
@@ -0,0 +1,44 @@
+C
+      SUBROUTINE atm2geo ( im, jm, pte, ptn, plon, plat, pxx, pyy, pzz )
+cc
+cc Change wind local atmospheric coordinates to
+cc geocentric
+cc
+c$$$      INCLUDE 'param.h'
+c
+      INTEGER, INTENT (in)              :: im, jm
+      REAL, DIMENSION (im,jm), INTENT (in) :: pte, ptn
+      REAL, DIMENSION (im,jm), INTENT (in) :: plon, plat
+      REAL, DIMENSION (im,jm), INTENT(out) :: pxx, pyy, pzz
+c
+      REAL, PARAMETER :: rpi = 3.141592653E0
+      REAL, PARAMETER :: rad = rpi / 180.0E0
+c
+      REAL, DIMENSION (im,jm) :: zsinlon, zcoslon
+      REAL, DIMENSION (im,jm) :: zsinlat, zcoslat
+c
+      LOGICAL, SAVE :: linit = .FALSE.
+c
+c$$$      IF ( .NOT. linit ) THEN 
+          zsinlon = SIN (rad * plon)
+          zcoslon = COS (rad * plon)
+          zsinlat = SIN (rad * plat)
+          zcoslat = COS (rad * plat)
+          linit = .TRUE.
+c$$$      ENDIF 
+c
+      pxx = - zsinlon * pte - zsinlat * zcoslon * ptn
+      pyy =   zcoslon * pte - zsinlat * zsinlon * ptn
+      pzz =   zcoslat * ptn
+c
+c Value at North Pole
+      pxx ( :,  1) = - ptn ( 1, 1)
+      pyy ( :,  1) = - pte ( 1, 1)
+      pzz ( :,  1) = 0.0
+c Value at South Pole
+      pxx ( :, jm) = + ptn ( 1, jm)
+      pyy ( :, jm) = + pte ( 1, jm)
+      pzz ( :, jm) = 0.0
+c
+      RETURN 
+      END SUBROUTINE atm2geo
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/calcratqs.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/calcratqs.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/calcratqs.F	(revision 222)
@@ -0,0 +1,161 @@
+c $Header$
+      SUBROUTINE calcratqs (
+     I            paprs,pplay,q_seri,d_t_con,d_t_ajs
+     O           ,ratqs,zpt_conv)
+      IMPLICIT none
+c======================================================================
+c
+c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
+c
+c Objet: Moniteur general de la physique du modele
+cAA      Modifications quant aux traceurs :
+cAA                  -  uniformisation des parametrisations ds phytrac
+cAA                  -  stockage des moyennes des champs necessaires
+cAA                     en mode traceur off-line 
+c======================================================================
+c    modif   ( P. Le Van ,  12/10/98 )
+c
+c  Arguments:
+c
+c paprs---input-R-pression pour chaque inter-couche (en Pa)
+c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
+c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
+#include "dimensions.h"
+#include "dimphy.h"
+      REAL paprs(klon,klev+1)
+      REAL pplay(klon,klev)
+      REAL d_t_con(klon,klev)
+      REAL d_t_ajs(klon,klev)
+      REAL ratqs(klon,klev)
+      LOGICAL pt_conv(klon,klev)
+      REAL q_seri(klon,klev)
+
+      logical firstcall
+      save firstcall
+      data firstcall/.true./
+
+
+      REAL ratqsmin,ratqsmax,zx,epmax
+      REAL ratqs1,ratqs2,ratqs3,ratqs4
+      REAL ratqsc1,ratqsc2,ratqsc3,ratqsc4
+      INTEGER i,k
+      INTEGER flag_ratqs
+      save ratqsmin,ratqsmax,flag_ratqs,epmax
+      save ratqs1,ratqs2,ratqs3,ratqs4
+      save ratqsc1,ratqsc2,ratqsc3,ratqsc4
+      real zpt_conv(klon,klev)
+
+      REAL zx_min
+      PARAMETER (zx_min=1.0)
+      REAL zx_max
+      PARAMETER (zx_max=0.1)
+
+      flag_ratqs=1
+      flag_ratqs=0
+c
+c Appeler le processus de condensation a grande echelle
+c et le processus de precipitation
+c
+      if (flag_ratqs.eq.0) then
+
+         ratqsmax=0.01
+         ratqsmin=0.3
+
+         if (firstcall) print*,'RATQS ANCIEN '
+         do k=1,klev
+         do i=1,klon
+            zx = pplay(i,k)/paprs(i,1)
+            zx = (zx_max-zx)/(zx_max-zx_min)
+            zx = MIN(MAX(zx,0.0),1.0)
+            zx = zx * zx * zx
+            ratqs(i,k)= zx * (ratqsmax-ratqsmin) + ratqsmin
+         enddo
+         enddo
+
+      else
+
+c  On aplique un ratqs "interactif" a toutes les mailles affectees
+c  par la convection ou se trouvant "sous" une maille affectee.
+         do i=1,klon
+            pt_conv(i,klev)=.false.
+         enddo
+         do k=klev-1,1,-1
+            do i=1,klon
+               pt_conv(i,k)=pt_conv(i,k+1).or.
+     s               (abs(d_t_con(i,k))+abs(d_t_ajs(i,k))).gt.1.e-8
+               if(pt_conv(i,k)) then
+                  zpt_conv(i,k)=1.
+               else
+                  zpt_conv(i,k)=0.
+               endif
+            enddo
+         enddo
+
+         if (flag_ratqs.eq.1) then
+
+            ratqsmin=0.4
+            ratqsmax=0.99
+            if (firstcall) print*,'RATQS INTERACTIF '
+            do k=1,klev
+                do i=1,klon
+                   if (pt_conv(i,k)) then
+                      ratqs(i,k)=0.01
+     s                +0.25*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
+                      ratqs(i,k)=min(ratqs(i,k),ratqsmax)
+                      ratqs(i,k)=max(ratqs(i,k),0.1)
+                   else
+                      ratqs(i,k)=0.01+(ratqsmin-0.01)*
+     s             min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
+                   endif
+                enddo
+            enddo
+         else if (flag_ratqs.eq.2) then
+            do k=1,klev
+                do i=1,klon
+                   ratqs(i,k)=0.001+
+     s             (q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
+                   if (pt_conv(i,k)) then
+                      ratqs(i,k)=min(ratqs(i,k),ratqsmax)
+                   else
+                      ratqs(i,k)=min(ratqs(i,k),ratqsmin)
+                   endif
+                enddo
+            enddo
+         else
+            do k=1,klev
+               do i=1,klon
+                  if (pplay(i,k).ge.95000.) then
+                     if (pt_conv(i,k)) then
+                        ratqs(i,k)=ratqsc1
+                     else
+                        ratqs(i,k)=ratqs1
+                     endif
+                  else if (pplay(i,k).ge.75000.) then
+                     if (pt_conv(i,k)) then
+                        ratqs(i,k)=ratqsc2
+                     else
+                        ratqs(i,k)=ratqs2
+                     endif
+                  else if (pplay(i,k).ge.50000.) then
+                     if (pt_conv(i,k)) then
+                        ratqs(i,k)=ratqsc3
+                     else
+                        ratqs(i,k)=ratqs3
+                     endif
+                  else
+                     if (pt_conv(i,k)) then
+                        ratqs(i,k)=ratqsc4
+                     else
+                        ratqs(i,k)=ratqs4
+                     endif
+                  endif
+               enddo
+            enddo
+         endif
+
+      endif
+
+      firstcall=.false.
+
+      return
+      end
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clift.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clift.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clift.F	(revision 222)
@@ -0,0 +1,70 @@
+c $Header$
+        SUBROUTINE CLIFT (P,T,RR,RS,PLCL,DPLCLDT,DPLCLDQ)
+C***************************************************************
+C*                                                             *
+C* CLIFT : COMPUTE LIFTING CONDENSATION LEVEL AND ITS          *
+C*         DERIVATIVES RELATIVE TO T AND R                     *
+C*   (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980)   *
+C*                                                             *
+C* written by   : GRANDPEIX Jean-Yves, 17/11/98, 12.39.01      *
+C* modified by :                                               *
+C***************************************************************
+C*
+C*Arguments :
+C*
+C* Input :  P = pressure of level from wich lifting is performed
+C*          T = temperature of level P
+C*          RR = vapour mixing ratio at level P
+C*          RS = vapour saturation mixing ratio at level P
+C*
+C* Output : PLCL = lifting condensation level
+C*          DPLCLDT = derivative of PLCL relative to T
+C*          DPLCLDQ = derivative of PLCL relative to R
+C*
+ccccccccccccccccccccccc
+c constantes coherentes avec le modele du Centre Europeen
+c      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
+c      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
+c      CPD = 3.5 * RD
+c      CPV = 4.0 * RV
+c      CL = 4218.0
+c      CI=2090.0
+c      CPVMCL=CL-CPV
+c      CLMCI=CL-CI
+c      EPS=RD/RV
+c      ALV0=2.5008E+06
+c      ALF0=3.34E+05
+c
+c on utilise les constantes thermo du Centre Europeen: (sb)
+c
+#include "YOMCST.h"
+c
+       CPD = RCPD
+       CPV = RCPV
+       CL = RCW
+       CPVMCL = CL-CPV
+       EPS = RD/RV
+       ALV0 = RLVTT
+c
+c
+c      Bolton formula coefficients :
+      A = 1669.0
+      B = 122.0
+c
+      RH=RR/RS
+      CHI=T/(A-B*RH-T)
+      PLCL=P*(RH**CHI)
+c
+      ALV = ALV0 - CPVMCL*(T-273.15)
+c
+c -- sb: correction:
+c       DPLCLDQ = PLCL*CHI*( 1./RR - B*CHI/T/RS*ALOG(RH) )
+      DPLCLDQ = PLCL*CHI*( 1./RR + B*CHI/T/RS*ALOG(RH) )
+c sb --
+c
+      DPLCLDT = PLCL*CHI*((A-B*RH*(1.+ALV/RV/T))/T**2*CHI*ALOG(RH)
+     $                    - ALV/RV/T**2 )
+c
+c
+      RETURN
+      END
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/conema3.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/conema3.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/conema3.F	(revision 222)
@@ -0,0 +1,281 @@
+c $Header$
+      SUBROUTINE conema (dtime,paprs,pplay,t,q,u,v,tra,ntra,
+     .             work1,work2,d_t,d_q,d_u,d_v,d_tra,
+     .             rain, snow, kbas, ktop,
+     .             upwd,dnwd,dnwdbis,bas,top,Ma,cape,tvp,rflag,
+     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)
+
+      IMPLICIT none
+c======================================================================
+c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
+c Objet: schema de convection de Emanuel (1991) interface
+c Mai 1998: Interface modifiee pour implementation dans LMDZ
+c======================================================================
+c Arguments:
+c dtime---input-R-pas d'integration (s)
+c paprs---input-R-pression inter-couches (Pa)
+c pplay---input-R-pression au milieu des couches (Pa)
+c t-------input-R-temperature (K)
+c q-------input-R-humidite specifique (kg/kg)
+c u-------input-R-vitesse du vent zonal (m/s)
+c v-------input-R-vitesse duvent meridien (m/s)
+c tra-----input-R-tableau de rapport de melange des traceurs
+c work*: input et output: deux variables de travail,
+c                            on peut les mettre a 0 au debut
+c
+C d_t-----output-R-increment de la temperature
+c d_q-----output-R-increment de la vapeur d'eau
+c d_u-----output-R-increment de la vitesse zonale
+c d_v-----output-R-increment de la vitesse meridienne
+c d_tra---output-R-increment du contenu en traceurs
+c rain----output-R-la pluie (mm/s)
+c snow----output-R-la neige (mm/s)
+c kbas----output-R-bas du nuage (integer)
+c ktop----output-R-haut du nuage (integer)
+c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
+c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
+c dnwdbis-output-R-unsaturated downdraft mass flux (kg/m**2/s)
+c bas-----output-R-bas du nuage (real)
+c top-----output-R-haut du nuage (real)
+c Ma------output-R-flux ascendant non dilue (kg/m**2/s)
+c cape----output-R-CAPE
+c tvp-----output-R-virtual temperature of the lifted parcel
+c rflag---output-R-flag sur le fonctionnement de convect
+c pbase---output-R-pression a la base du nuage (Pa)
+c bbase---output-R-buoyancy a la base du nuage (K)
+c dtvpdt1-output-R-derivative of parcel virtual temp wrt T1 
+c dtvpdq1-output-R-derivative of parcel virtual temp wrt Q1 
+c dplcldt-output-R-derivative of the PCP pressure wrt T1
+c dplcldr-output-R-derivative of the PCP pressure wrt Q1
+c======================================================================
+c
+#include "dimensions.h"
+#include "dimphy.h"
+      INTEGER i, l,m,itra
+      INTEGER ntra,ntrac !number of tracers; if no tracer transport
+                         ! is needed, set ntra = 1 (or 0)
+      PARAMETER (ntrac=nqmx-2)
+      REAL dtime
+c
+
+c 
+      REAL paprs(klon,klev+1), pplay(klon,klev)
+      REAL t(klon,klev), q(klon,klev), d_t(klon,klev), d_q(klon,klev)
+      REAL u(klon,klev), v(klon,klev), tra(klon,klev,ntra)
+      REAL d_u(klon,klev), d_v(klon,klev), d_tra(klon,klev,ntra)
+      REAL work1(klon,klev), work2(klon,klev)
+      REAL upwd(klon,klev), dnwd(klon,klev), dnwdbis(klon,klev)
+      REAL rain(klon)
+      REAL snow(klon)
+      REAL cape(klon), tvp(klon,klev), rflag(klon)
+      REAL pbase(klon), bbase(klon)
+      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
+      REAL dplcldt(klon), dplcldr(klon)
+      INTEGER kbas(klon), ktop(klon)
+c
+      REAL em_t(klev)
+      REAL em_q(klev)
+      REAL em_qs(klev)
+      REAL em_u(klev), em_v(klev), em_tra(klev,ntrac)
+      REAL em_ph(klev+1), em_p(klev)
+      REAL em_work1(klev), em_work2(klev)
+      REAL em_precip, em_d_t(klev), em_d_q(klev)
+      REAL em_d_u(klev), em_d_v(klev), em_d_tra(klev,ntrac)
+      REAL em_upwd(klev), em_dnwd(klev), em_dnwdbis(klev)
+      REAL em_dtvpdt1(klev), em_dtvpdq1(klev)
+      REAL em_dplcldt, em_dplcldr
+      SAVE em_t,em_q, em_qs, em_ph, em_p, em_work1, em_work2
+      SAVE em_u,em_v, em_tra
+      SAVE em_d_u,em_d_v, em_d_tra
+      SAVE em_precip, em_d_t, em_d_q, em_upwd, em_dnwd, em_dnwdbis
+      INTEGER em_bas, em_top
+      SAVE em_bas, em_top
+c
+      REAL zx_t, zx_qs, zdelta, zcor
+      INTEGER iflag
+      REAL sigsum
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     VARIABLES A SORTIR
+cccccccccccccccccccccccccccccccccccccccccccccccccc
+ 
+      REAL emmip(klev) !variation de flux ascnon dilue i et i+1
+      SAVE emmip
+      real emMke(klev)
+      save emMke
+      real top
+      real bas
+      real emMa(klev)
+      save emMa
+      real Ma(klon,klev)
+      real Ment(klev,klev)
+      real Qent(klev,klev)
+      real TPS(klev),TLS(klev)
+      real SIJ(klev,klev)
+      real em_CAPE, em_TVP(klev)
+      real em_pbase, em_bbase
+      integer iw,j,k,ix,iy
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+ 
+c
+c$$$      print*,'debut conema'
+
+      DO 999 i = 1, klon
+      DO l = 1, klev+1
+         em_ph(l) = paprs(i,l) / 100.0
+      ENDDO
+c
+      DO l = 1, klev
+         em_p(l) = pplay(i,l) / 100.0
+         em_t(l) = t(i,l)
+         em_q(l) = q(i,l)
+         em_u(l) = u(i,l)
+         em_v(l) = v(i,l)
+         do itra = 1, ntra
+          em_tra(l,itra) = tra(i,l,itra)
+         enddo
+c$$$      print*,'em_t',em_t
+c$$$      print*,'em_q',em_q
+c$$$      print*,'em_qs',em_qs
+c$$$      print*,'em_u',em_u
+c$$$      print*,'em_v',em_v
+c$$$      print*,'em_tra',em_tra
+c$$$      print*,'em_p',em_p
+
+ 
+c
+         zx_t = em_t(l)
+         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
+         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(l)/100.0
+         zx_qs=MIN(0.5,zx_qs)
+c$$$       print*,'zx_qs',zx_qs
+         zcor=1./(1.-retv*zx_qs) 
+         zx_qs=zx_qs*zcor
+         em_qs(l) = zx_qs
+c$$$      print*,'em_qs',em_qs
+c
+         em_work1(l) = work1(i,l)
+         em_work2(l) = work2(i,l)
+         emMke(l)=0
+c        emMa(l)=0
+c        Ma(i,l)=0
+     
+         em_dtvpdt1(l) = 0.
+         em_dtvpdq1(l) = 0.
+         dtvpdt1(i,l) = 0.
+         dtvpdq1(i,l) = 0.
+      ENDDO
+c
+      em_dplcldt = 0.
+      em_dplcldr = 0.
+      rain(i) = 0.0
+      snow(i) = 0.0
+      kbas(i) = 1
+      ktop(i) = 1
+c ajout SB:
+      bas = 1
+      top = 1
+ 
+ 
+c sb3d      write(*,1792) (em_work1(m),m=1,klev)
+1792  format('sig avant convect ',/,10(1X,E13.5))
+c
+c sb d      write(*,1793) (em_work2(m),m=1,klev)
+1793  format('w avant convect ',/,10(1X,E13.5))
+ 
+c$$$      print*,'avant convect' 
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c 
+      CALL convect(dtime, em_t, em_q, em_qs,em_u ,em_v ,
+     .              em_tra, em_p, em_ph,
+     .              klev, klev+1, klev-1,ntra, dtime, iflag,
+     .              em_d_t, em_d_q,em_d_u,em_d_v,
+     .              em_d_tra, em_precip,
+     .              em_bas, em_top,em_upwd, em_dnwd, em_dnwdbis,
+     .              em_work1, em_work2,emmip,emMke,emMa,Ment,
+c SB 11sept98     .              Qent,TPS,TLS,SIJ)
+c 19oct98     .              Qent,TPS,TLS,SIJ,em_CAPE,em_TVP)
+     .  Qent,TPS,TLS,SIJ,em_CAPE,em_TVP,em_pbase,em_bbase,
+     .  em_dtvpdt1,em_dtvpdq1,em_dplcldt,em_dplcldr)
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c SB:
+      if (iflag.ne.1 .and. iflag.ne.4) then
+         em_CAPE = 0.
+      do l = 1, klev
+         em_upwd(l) = 0.
+         em_dnwd(l) = 0.
+         em_dnwdbis(l) = 0.
+         emMa(l) = 0.
+         em_TVP(l) = 0.
+      enddo
+      endif
+c fin SB
+c
+c  If sig has been set to zero, then set Ma to zero
+c
+      sigsum = 0.
+      do k = 1,klev
+        sigsum = sigsum + em_work1(k)
+      enddo
+      if (sigsum .eq. 0.0) then
+        do k = 1,klev
+          emMa(k) = 0.
+        enddo
+      endif
+c
+c sb3d       print*,'i, iflag=',i,iflag
+c 
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c       SORTIE DES ICB ET INB
+c       en fait inb et icb correspondent au niveau ou se trouve
+c       le nuage,le numero d'interface
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ 
+c modif SB:
+      if (iflag.EQ.1 .or. iflag.EQ.4) then
+       top=em_top
+       bas=em_bas
+       kbas(i) = em_bas
+       ktop(i) = em_top
+      endif
+ 
+      pbase(i) = em_pbase
+      bbase(i) = em_bbase
+      rain(i) = em_precip/ 86400.0
+      snow(i) = 0.0
+      cape(i) = em_CAPE
+      rflag(i) = float(iflag)
+c SB      kbas(i) = em_bas
+c SB      ktop(i) = em_top
+      dplcldt(i) = em_dplcldt
+      dplcldr(i) = em_dplcldr
+      DO l = 1, klev
+         d_t(i,l) = dtime * em_d_t(l) 
+         d_q(i,l) = dtime * em_d_q(l)
+         d_u(i,l) = dtime * em_d_u(l)
+         d_v(i,l) = dtime * em_d_v(l)
+         do itra = 1, ntra
+         d_tra(i,l,itra) = dtime * em_d_tra(l,itra)
+         enddo
+         upwd(i,l) = em_upwd(l)
+         dnwd(i,l) = em_dnwd(l)
+         dnwdbis(i,l) = em_dnwdbis(l)
+         work1(i,l) = em_work1(l)
+         work2(i,l) = em_work2(l)
+         Ma(i,l)=emMa(l)
+         tvp(i,l)=em_TVP(l)
+         dtvpdt1(i,l) = em_dtvpdt1(l)
+         dtvpdq1(i,l) = em_dtvpdq1(l)
+      ENDDO
+  999 CONTINUE
+ 
+
+      RETURN
+      END
+
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/conf_phys.F90	(revision 222)
@@ -0,0 +1,102 @@
+!
+! $Header$
+!
+
+!module conf_phys
+!
+!  use IOIPSL
+!  implicit none
+!
+!  public conf_phys
+!
+!contains
+
+  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan)
+
+!
+! Configuration du module "physique" de LMDZ a l'aide de la fonction
+! GETIN de IOIPSL
+!
+! LF 05/2001
+!
+
+  use IOIPSL     
+
+!
+! ocean:      type d'ocean (force, slab, couple)
+! ok_veget:   type de modele de vegetation
+! ok_journe:  sorties journalieres
+! ok_mensuel: sorties mensuelles
+! ok_instan:  sorties instantanees
+
+
+
+! Sortie:
+  character (len = 6)  :: ocean
+  logical              :: ok_veget
+  logical              :: ok_journe, ok_mensuel, ok_instan
+
+!
+!
+!
+
+
+!Config Key  = OCEAN 
+!Config Desc = Type d'ocean
+!Config Def  = force
+!Config Help = Type d'ocean utilise: force, slab,couple
+!
+  ocean = 'force '
+  call getin('OCEAN', ocean)
+!
+!Config Key  = VEGET 
+!Config Desc = Type de modele de vegetation
+!Config Def  = .false.
+!Config Help = Type de modele de vegetation utilise
+!
+  ok_veget = .false.
+  call getin('VEGET', ok_veget)
+!
+!Config Key  = OK_journe
+!Config Desc = Pour des sorties journalieres 
+!Config Def  = .false.
+!Config Help = Pour creer le fichier histday contenant les sorties
+!              journalieres 
+!
+  ok_journe = .false.
+  call getin('OK_journe', ok_journe)
+!
+!Config Key  = OK_mensuel
+!Config Desc = Pour des sorties mensuelles 
+!Config Def  = .true.
+!Config Help = Pour creer le fichier histmth contenant les sorties
+!              mensuelles 
+!
+  ok_mensuel = .true.
+  call getin('OK_mensuel', ok_mensuel)
+!
+!Config Key  = OK_instan
+!Config Desc = Pour des sorties instantanees 
+!Config Def  = .false.
+!Config Help = Pour creer le fichier histins contenant les sorties
+!              instantanees 
+!
+  ok_instan = .false.
+  call getin('OK_instan', ok_instan)
+
+!
+!
+!
+
+    write(numout,*)' Config ocean = ', ocean
+    write(numout,*)' Config veget = ', ok_veget
+    write(numout,*)' Sortie journaliere = ', ok_journe
+    write(numout,*)' Sortie mensuelle = ', ok_mensuel
+    write(numout,*)' Sortie instantanee = ', ok_instan
+
+
+  return
+
+  end subroutine conf_phys
+
+!end module conf_phys
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/convect3.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/convect3.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/convect3.F	(revision 222)
@@ -0,0 +1,1252 @@
+c $Header$
+      SUBROUTINE CONVECT
+     *    (DTIME,  T1,   R1,   RS,    U,  V,  TRA,   P,     PH,
+     *     ND,       NDP1,     NL, NTRA,  DELT,  IFLAG,
+     *     FT, FR, FU,  FV,  FTRA,  PRECIP,
+     *     icb,inb,   upwd,dnwd,dnwd0,SIG, W0,Mike,Mke,
+     *     Ma,MENTS,QENTS,TPS,TLS,SIGIJ,CAPE,TVP,PBASE,BUOYBASE,
+     *     DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR)
+C
+C    ***  THE PARAMETER NA SHOULD IN GENERAL EQUAL ND   ***
+C
+c#################################################################
+cFleur       Introduction des traceurs dans convect3 le 6 juin 200
+c#################################################################
+#include "dimensions.h"
+#include "dimphy.h"
+      PARAMETER (NA=60)
+
+      integer NTRAC
+      PARAMETER (NTRAC=nqmx-2)
+
+      INTEGER NENT(NA)
+      REAL T1(ND),R1(ND),RS(ND),U(ND),V(ND),TRA(ND,NTRA)
+      REAL P(ND),PH(NDP1)
+      REAL FT(ND),FR(ND),FU(ND),FV(ND),FTRA(ND,NTRA)
+      REAL SIG(ND),W0(ND)
+      REAL UENT(NA,NA),VENT(NA,NA),TRAENT(NA,NA,NTRAC),TRATM(NA)
+      REAL UP(NA),VP(NA),TRAP(NA,NTRAC)
+      REAL M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA)
+      REAL SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA)
+      REAL RP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA)
+      REAL SIGP(NA),B(NA),TP(NA),CPN(NA)
+      REAL LV(NA),LVCP(NA),H(NA),HP(NA),GZ(NA)
+      REAL T(NA),RR(NA)
+C
+      REAL BUOY(NA)     !  Lifted parcel buoyancy
+      REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
+C                                      temperature wrt T1 and Q1
+      REAL CLW_NEW(NA),QI(NA)
+C
+      LOGICAL ICE_CONV
+      PARAMETER (ICE_CONV=.TRUE.)
+ 
+cccccccccccccccccccccccccccccccccccccccccccccc
+c     declaration des variables a sortir
+ccccccccccccccccccccccccccccccccccccccccccccc
+      real Mke(nd)
+      real Mike(nd)
+      real Ma(nd)
+      real TPS(ND) !temperature dans les ascendances non diluees
+      real TLS(ND) !temperature potentielle
+      real MENTS(nd,nd)
+      real QENTS(nd,nd)
+      REAL SIGIJ(KLEV,KLEV)
+      REAL PBASE ! pressure at the cloud base level
+      REAL BUOYBASE ! buoyancy at the cloud base level
+cccccccccccccccccccccccccccccccccccccccccccccc
+ 
+ 
+ 
+c
+      real dnwd0(nd)  !  precipitation driven unsaturated downdraft flux
+      real dnwd(nd), dn1  ! in-cloud saturated downdraft mass flux
+      real upwd(nd), up1  ! in-cloud saturated updraft mass flux
+C
+C   ***         ASSIGN VALUES OF THERMODYNAMIC CONSTANTS        ***
+C   ***             THESE SHOULD BE CONSISTENT WITH             ***
+C   ***              THOSE USED IN CALLING PROGRAM              ***
+C   ***     NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT  ***
+C
+c sb      CPD=1005.7
+c sb      CPV=1870.0
+c sb      CL=4190.0
+c sb      CPVMCL=CL-CPV
+c sb      RV=461.5
+c sb      RD=287.04
+c sb      EPS=RD/RV
+c sb      ALV0=2.501E6
+ccccccccccccccccccccccc
+c constantes coherentes avec le modele du Centre Europeen
+c sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
+c sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
+c sb      CPD = 3.5 * RD
+c sb      CPV = 4.0 * RV
+c sb      CL = 4218.0
+c sb      CPVMCL=CL-CPV
+c sb      EPS=RD/RV
+c sb      ALV0=2.5008E+06
+cccccccccccccccccccccc
+c on utilise les constantes thermo du Centre Europeen: (SB)
+c
+#include "YOMCST.h"
+c
+       CPD = RCPD
+       CPV = RCPV
+       CL = RCW
+       CPVMCL = CL-CPV
+       EPS = RD/RV
+       ALV0 = RLVTT
+c
+       NK = 1 ! origin level of the lifted parcel
+c
+cccccccccccccccccccccc
+C
+C           ***  INITIALIZE OUTPUT ARRAYS AND PARAMETERS  ***
+C
+      DO 5 I=1,ND
+         FT(I)=0.0
+         FR(I)=0.0
+         FU(I)=0.0
+         FV(I)=0.0
+         DO 4 J=1,NTRA
+          FTRA(I,J)=0.0
+    4    CONTINUE
+         T(I)=T1(I)
+         RR(I)=R1(I)
+    5 CONTINUE
+      DO 7 I=1,NL
+         RDCP=(RD*(1.-RR(I))+RR(I)*RV)/ (CPD*(1.-RR(I))+RR(I)*CPV)
+         TH(I)=T(I)*(1000.0/P(I))**RDCP
+    7 CONTINUE
+C
+*************************************************************
+**    CALCUL DES TEMPERATURES POTENTIELLES A SORTIR
+*************************************************************
+      do i=1,ND
+      RDCP=(RD*(1.-RR(I))+RR(I)*RV)/ (CPD*(1.-RR(I))+RR(I)*CPV)
+ 
+      TLS(i)=T(I)*(1000.0/P(I))**RDCP
+      enddo
+ 
+ 
+ 
+ 
+************************************************************
+ 
+ 
+      PRECIP=0.0
+      IFLAG=1
+C
+C   ***                    SPECIFY PARAMETERS                        ***
+C   ***  PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE   ***
+C   ***       PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO         ***
+C   ***  PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP.      ***
+C   ***            EFFICIENCY IS ASSUMED TO BE UNITY                 ***
+C   ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
+C   ***  SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE      ***
+C   ***                        OF CLOUD                              ***
+C   ***    ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF    ***
+C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
+C   ***    (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY)    ***
+C   ***           (BETA MUST BE LESS THAN OR EQUAL TO 1)             ***
+C   ***    DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE    ***
+C   ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
+C   ***                     IT MUST BE LESS THAN 0                   ***
+C
+      PBCRIT=150.0
+      PTCRIT=500.0
+      SIGD=0.01
+      SPFAC=0.15
+c sb:
+      EPMAX=0.993 ! precip efficiency less than unity
+c      EPMAX=1. ! precip efficiency less than unity
+C
+Cjyg
+CCC      BETA=0.96
+C  Beta is now expressed as a function of the characteristic time
+C  of the convective process.
+CCC        Old value : TAU = 15000.   !(for dtime = 600.s)
+CCC        Other value (inducing little change) :TAU = 8000.
+      TAU = 8000.
+      BETA = 1.-DTIME/TAU
+Cjyg
+CCC      ALPHA=1.0
+      ALPHA=1.5E-3*DTIME/TAU
+C        Increase alpha in order to compensate W decrease
+      ALPHA = ALPHA*1.5
+C
+Cjyg (voir CONVECT 3)
+CCC      DTCRIT=-0.2
+      DTCRIT=-2.
+Cgf&jyg
+CCC     DT pour l'overshoot.
+      DTOVSH = -0.2
+ 
+C
+C           ***        INCREMENT THE COUNTER       ***
+C
+      SIG(ND)=SIG(ND)+1.0
+      SIG(ND)=AMIN1(SIG(ND),12.1)
+C
+C           ***    IF NOPT IS AN INTEGER OTHER THAN 0, CONVECT     ***
+C           ***     RETURNS ARRAYS T AND R THAT MAY HAVE BEEN      ***
+C           ***  ALTERED BY DRY ADIABATIC ADJUSTMENT; OTHERWISE    ***
+C           ***        THE RETURNED ARRAYS ARE UNALTERED.          ***
+C
+      NOPT=0
+C
+C     ***            PERFORM DRY ADIABATIC ADJUSTMENT            ***
+C
+C     ***  DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A  ***
+C     ***                BOUNDARY LAYER SCHEME !!!               ***
+C
+c test sb      DO 30 I=NL-1,1,-1
+c test sb         JN=0
+c test sb         DO 10 J=I+1,NL
+c test sb   10    IF(TH(J).LT.TH(I))JN=J
+c test sb         IF(JN.EQ.0)GOTO 30
+c test sb         AHM=0.0
+c test sb         RM=0.0
+c test sb         UM=0.0
+c test sb         VM=0.0
+c test sb         DO K=1,NTRA
+c test sb          TRATM(K)=0.0
+c test sb         END DO
+c test sb         DO 15 J=I,JN
+c test sb          AHM=AHM+(CPD*(1.-RR(J))+RR(J)*CPV)*T(J)*(PH(J)-PH(J+1))
+c test sb          RM=RM+RR(J)*(PH(J)-PH(J+1))
+c test sb          UM=UM+U(J)*(PH(J)-PH(J+1))
+c test sb          VM=VM+V(J)*(PH(J)-PH(J+1))
+c test sb          DO K=1,NTRA
+c test sb           TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1))
+c test sb          END DO
+c test sb   15    CONTINUE
+c test sb         DPHINV=1./(PH(I)-PH(JN+1))
+c test sb         RM=RM*DPHINV
+c test sb         UM=UM*DPHINV
+c test sb         VM=VM*DPHINV
+c test sb         DO K=1,NTRA
+c test sb          TRATM(K)=TRATM(K)*DPHINV
+c test sb         END DO
+c test sb         A2=0.0
+c test sb         DO 20 J=I,JN
+c test sb            RR(J)=RM
+c test sb          U(J)=UM
+c test sb          V(J)=VM
+c test sb          DO K=1,NTRA
+c test sb           TRA(J,K)=TRATM(K)
+c test sb          END DO
+c test sb            RDCP=(RD*(1.-RR(J))+RR(J)*RV)/ (CPD*(1.-RR(J))+RR(J)*CPV)
+c test sb            X=(0.001*P(J))**RDCP
+c test sb            T(J)=X
+c test sb            A2=A2+(CPD*(1.-RR(J))+RR(J)*CPV)*X*(PH(J)-PH(J+1))
+c test sb   20    CONTINUE
+c test sb         DO 25 J=I,JN
+c test sb            TH(J)=AHM/A2
+c test sb            T(J)=T(J)*TH(J)
+c test sb   25    CONTINUE
+c test sb   30 CONTINUE
+C
+C   ***   RESET INPUT ARRAYS IF NOPT .NE. 0   ***
+C
+      IF (NOPT.NE.0)THEN
+         DO 35 I=1,ND
+            T1(I)=T(I)
+            R1(I)=RR(I)
+   35    CONTINUE
+      END IF
+C
+C  *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY
+C
+      GZ(1)=0.0
+      CPN(1)=CPD*(1.-RR(1))+RR(1)*CPV
+      H(1)=T(1)*CPN(1)
+      DO 40 I=2,NL
+        TVX=T(I)*(1.+RR(I)/EPS-RR(I))
+        TVY=T(I-1)*(1.+RR(I-1)/EPS-RR(I-1))
+        GZ(I)=GZ(I-1)+0.5*RD*(TVX+TVY)*(P(I-1)-P(I))/PH(I)
+        CPN(I)=CPD*(1.-RR(I))+CPV*RR(I)
+        H(I)=T(I)*CPN(I)+GZ(I)
+   40 CONTINUE
+C
+C   ***  CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT LOWEST MODEL LEVEL ***
+C   ***       (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980)     ***
+C
+      IF (T(1).LT.250.0.OR.RR(1).LE.0.0)THEN
+         IFLAG=0
+c sb3d         print*,'je suis passe par 366'
+         RETURN
+      END IF
+
+cjyg1 Utilisation de la subroutine CLIFT
+CC      RH=RR(1)/RS(1)
+CC      CHI=T(1)/(1669.0-122.0*RH-T(1))
+CC      PLCL=P(1)*(RH**CHI)
+      CALL CLIFT(P(1),T(1),RR(1),RS(1),PLCL,DPLCLDT,DPLCLDR)
+cjyg2
+c sb3d      PRINT *,' em_plcl,p1,t1,r1,rs1,rh '
+c sb3d     $        ,PLCL,P(1),T(1),RR(1),RS(1),RH
+c
+      IF (PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN
+         IFLAG=2
+         RETURN
+      END IF
+Cjyg1
+C     Essais de modification de ICB
+C
+C   ***  CALCULATE FIRST LEVEL ABOVE LCL (=ICB)  ***
+C
+CC      ICB=NL-1
+CC      DO 50 I=2,NL-1
+CC         IF(P(I).LT.PLCL)THEN
+CC            ICB=MIN(ICB,I)   ! ICB sup ou egal a 2
+CC         END IF
+CC   50 CONTINUE
+CC      IF(ICB.EQ.(NL-1))THEN
+CC         IFLAG=3
+CC         RETURN
+CC      END IF
+C
+C   *** CALCULATE LAYER CONTAINING LCL (=ICB)   ***
+C
+      ICB=NL-1
+c sb      DO 50 I=2,NL-1
+      DO 50 I=3,NL-1 ! modif sb pour que ICB soit sup/egal a 2
+C   la modification consiste a comparer PLCL a PH et non a P:
+C   ICB est defini par :  PH(ICB)<PLCL<PH(ICB-!)
+         IF(PH(I).LT.PLCL)THEN
+            ICB=MIN(ICB,I)
+         END IF
+   50 CONTINUE
+      IF(ICB.EQ.(NL-1))THEN
+         IFLAG=3
+         RETURN
+      END IF
+      ICB = ICB-1 ! ICB sup ou egal a 2 
+Cjyg2
+C
+C
+ 
+C   *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL      ***
+C   ***  TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC             ***
+C   ***                   LIQUID WATER CONTENT                             ***
+C
+ 
+cjyg1
+c make sure that "Cloud base" seen by TLIFT is actually the 
+c fisrt level where adiabatic ascent is saturated 
+       IF (PLCL .GT. P(ICB)) THEN
+c sb        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,TVP,TP,CLW,ND,NL)
+        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK,TVP,TP,CLW,ND,NL
+     :            ,DTVPDT1,DTVPDQ1)
+       ELSE
+c sb        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,TVP,TP,CLW,ND,NL)
+        CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,NK,TVP,TP,CLW,ND,NL
+     :            ,DTVPDT1,DTVPDQ1)
+       ENDIF
+cjyg2
+ 
+******************************************************************************
+****     SORTIE DE LA TEMPERATURE DE L ASCENDANCE NON DILUE
+******************************************************************************
+        do i=1,ND
+        TPS(i)=TP(i)
+        enddo
+ 
+ 
+******************************************************************************
+ 
+C
+C   ***  SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF   ***
+C   ***          PRECIPITATION FALLING OUTSIDE OF CLOUD           ***
+C   ***      THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)     ***
+C
+      DO 55 I=1,NL
+         PDEN=PTCRIT-PBCRIT
+c
+cjyg
+ccc         EP(I)=(P(ICB)-P(I)-PBCRIT)/PDEN
+c sb         EP(I)=(PLCL-P(I)-PBCRIT)/PDEN
+         EP(I)=(PLCL-P(I)-PBCRIT)/PDEN * EPMAX ! sb
+c
+         EP(I)=AMAX1(EP(I),0.0)
+c sb         EP(I)=AMIN1(EP(I),1.0)
+         EP(I)=AMIN1(EP(I),EPMAX) ! sb
+         SIGP(I)=SPFAC
+C
+C   ***       CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL     ***
+C   ***                    VIRTUAL TEMPERATURE                    ***
+C
+         TV(I)=T(I)*(1.+RR(I)/EPS-RR(I))
+Ccd1
+C    . Keep all liquid water in lifted parcel (-> adiabatic CAPE)
+C
+ccc    TVP(I)=TVP(I)-TP(I)*(RR(1)-EP(I)*CLW(I))
+c!!! sb         TVP(I)=TVP(I)-TP(I)*RR(1) ! calcule dans tlift
+Ccd2
+C
+C   ***       Calculate first estimate of buoyancy
+C
+         BUOY(I) = TVP(I) - TV(I)
+   55 CONTINUE
+C
+C   ***   Set Cloud Base Buoyancy at (Plcl+DPbase) level buoyancy
+C
+      DPBASE = -40.   !That is 400m above LCL
+      PBASE = PLCL + DPBASE
+      TVPBASE = TVP(ICB  )*(PBASE -P(ICB+1))/(P(ICB)-P(ICB+1))
+     $         +TVP(ICB+1)*(P(ICB)-PBASE   )/(P(ICB)-P(ICB+1))
+      TVBASE = TV(ICB  )*(PBASE -P(ICB+1))/(P(ICB)-P(ICB+1))
+     $        +TV(ICB+1)*(P(ICB)-PBASE   )/(P(ICB)-P(ICB+1))
+C
+      BUOYBASE = TVPBASE - TVBASE
+C
+CC       Set buoyancy = BUOYBASE for all levels below BASE.
+CC       For safety, set : BUOY(ICB) = BUOYBASE
+      DO I = ICB,NL
+        IF (P(I) .GE. PBASE) THEN
+          BUOY(I) = BUOYBASE
+        ENDIF
+      ENDDO
+      BUOY(ICB) = BUOYBASE
+C
+c sb3d      print *,'buoybase,tvp_tv,tvpbase,tvbase,pbase,plcl'
+c sb3d     $,        buoybase,tvp(icb)-tv(icb),tvpbase,tvbase,pbase,plcl
+c sb3d      print *,'TVP ',(tvp(i),i=1,nl)
+c sb3d      print *,'TV ',(tv(i),i=1,nl)
+c sb3d      print *, 'P ',(p(i),i=1,nl)
+c sb3d      print *,'ICB ',icb
+C
+C   ***   MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE  ***
+C   ***    AND CLOUD BASE, AND THAT LIFTED AIR IS POSITIVELY BUOYANT  ***
+C   ***                         AT CLOUD BASE                         ***
+C   ***       IF NOT, RETURN TO CALLING PROGRAM AFTER RESETTING       ***
+C   ***                        SIG(I) AND W0(I)                       ***
+C
+Cjyg
+CCC      TDIF=TVP(ICB)-TV(ICB)
+      TDIF = BUOY(ICB)
+      ATH1=TH(1)
+Cjyg
+CCC      ATH=TH(ICB-1)-1.0
+      ATH=TH(ICB-1)-5.0
+c      ATH=0.                          ! ajout sb
+c      IF (ICB.GT.1) ATH=TH(ICB-1)-5.0 ! modif sb
+      IF(TDIF.LT.DTCRIT.OR.ATH.GT.ATH1)THEN
+         DO 60 I=1,NL
+            SIG(I)=BETA*SIG(I)-2.*ALPHA*TDIF*TDIF
+            SIG(I)=AMAX1(SIG(I),0.0)
+            W0(I)=BETA*W0(I)
+   60    CONTINUE
+         IFLAG=0
+         RETURN
+      END IF
+C
+ 
+ 
+C   ***  IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY ***
+C   ***        NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS       ***
+C
+      DO 70 I=1,NL
+         HP(I)=H(I)
+         WT(I)=0.001
+         RP(I)=RR(I)
+         UP(I)=U(I)
+         VP(I)=V(I)
+         DO 71 J=1,NTRA
+          TRAP(I,J)=TRA(I,J)
+   71    CONTINUE
+         NENT(I)=0
+         WATER(I)=0.0
+         EVAP(I)=0.0
+         B(I)=0.0
+         MP(I)=0.0
+         M(I)=0.0
+         LV(I)=ALV0-CPVMCL*(T(I)-273.15)
+         LVCP(I)=LV(I)/CPN(I)
+         DO 70 J=1,NL
+            QENT(I,J)=RR(J)
+            ELIJ(I,J)=0.0
+            MENT(I,J)=0.0
+            SIJ(I,J)=0.0
+          UENT(I,J)=U(J)
+          VENT(I,J)=V(J)
+          DO 70 K=1,NTRA
+           TRAENT(I,J,K)=TRA(J,K)
+   70 CONTINUE
+ 
+      DELTI=1.0/DELT
+C
+C  ***  FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S       ***
+C  ***                LEVEL OF NEUTRAL BUOYANCY                   ***
+C
+      INB=NL-1
+      DO 80 I=ICB,NL-1
+Cjyg
+CCC         IF((TVP(I)-TV(I)).LT.DTCRIT)THEN
+         IF(BUOY(I).LT.DTOVSH)THEN
+            INB=MIN(INB,I)
+         END IF
+   80 CONTINUE
+ 
+ 
+ 
+ 
+C
+C   ***          RESET SIG(I) AND W0(I) FOR I>INB AND I<ICB       ***
+C
+      IF(INB.LT.(NL-1))THEN
+         DO 85 I=INB+1,NL-1
+Cjyg
+CCC            SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(INB)-TVP(INB))*
+CCC     1              ABS(TV(INB)-TVP(INB))
+            SIG(I)=BETA*SIG(I)+2.*ALPHA*BUOY(INB)*
+     1              ABS(BUOY(INB))
+            SIG(I)=AMAX1(SIG(I),0.0)
+            W0(I)=BETA*W0(I)
+   85    CONTINUE
+      END IF
+      DO 87 I=1,ICB
+Cjyg
+CCC         SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(ICB)-TVP(ICB))*
+CCC     1           (TV(ICB)-TVP(ICB))
+         SIG(I)=BETA*SIG(I)-2.*ALPHA*BUOY(ICB)*BUOY(ICB)
+         SIG(I)=AMAX1(SIG(I),0.0)
+         W0(I)=BETA*W0(I)
+   87 CONTINUE
+C
+C   ***    RESET FRACTIONAL AREAS OF UPDRAFTS AND W0 AT INITIAL TIME    ***
+C   ***           AND AFTER 10 TIME STEPS OF NO CONVECTION              ***
+C
+ 
+      IF(SIG(ND).LT.1.5.OR.SIG(ND).GT.12.0)THEN
+         DO 90 I=1,NL-1
+            SIG(I)=0.0
+            W0(I)=0.0
+   90    CONTINUE
+      END IF
+C
+C   ***   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL   ***
+C
+      DO 95 I=ICB,INB
+         HP(I)=H(1)+(LV(I)+(CPD-CPV)*T(I))*EP(I)*CLW(I)
+   95 CONTINUE
+C
+C   ***  CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY (CAPE),  ***
+C   ***     VERTICAL VELOCITY (W), FRACTIONAL AREA COVERED BY     ***
+C   ***     UNDILUTE UPDRAFT (SIG),  AND UPDRAFT MASS FLUX (M)    ***
+C
+      CAPE=0.0
+C
+      DO 98 I=ICB+1,INB
+Cjyg1
+CCC         CAPE=CAPE+RD*(TVP(I-1)-TV(I-1))*(PH(I-1)-PH(I))/P(I-1)
+CCC         DCAPE=RD*BUOY(I-1)*(PH(I-1)-PH(I))/P(I-1)
+CCC         DLNP=(PH(I-1)-PH(I))/P(I-1)
+C          The interval on which CAPE is computed starts at PBASE :
+         DELTAP = MIN(PBASE,PH(I-1))-MIN(PBASE,PH(I))
+         CAPE=CAPE+RD*BUOY(I-1)*DELTAP/P(I-1)
+         DCAPE=RD*BUOY(I-1)*DELTAP/P(I-1)
+         DLNP=DELTAP/P(I-1)
+Cjyg2
+c sb3d         print *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape
+         CAPE=AMAX1(0.0,CAPE)
+C
+         SIGOLD=SIG(I)
+         DTMIN=100.0
+         DO 97 J=ICB,I-1
+Cjyg
+CCC            DTMIN=AMIN1(DTMIN,(TVP(J)-TV(J)))
+            DTMIN=AMIN1(DTMIN,BUOY(J))
+   97    CONTINUE
+c sb3d     print *, 'DTMIN, BETA, ALPHA, SIG = ',DTMIN,BETA,ALPHA,SIG(I)
+         SIG(I)=BETA*SIG(I)+ALPHA*DTMIN*ABS(DTMIN)
+         SIG(I)=AMAX1(SIG(I),0.0)
+         SIG(I)=AMIN1(SIG(I),0.01)
+         FAC=AMIN1(((DTCRIT-DTMIN)/DTCRIT),1.0)
+Cjyg
+CC    Essais de reduction de la vitesse
+CC         FAC = FAC*.5
+C
+         W=(1.-BETA)*FAC*SQRT(CAPE)+BETA*W0(I)
+         AMU=0.5*(SIG(I)+SIGOLD)*W
+         M(I)=AMU*0.007*P(I)*(PH(I)-PH(I+1))/TV(I)
+         W0(I)=W
+   98 CONTINUE
+      W0(ICB)=0.5*W0(ICB+1)
+      M(ICB)=0.5*M(ICB+1)*(PH(ICB)-PH(ICB+1))/(PH(ICB+1)-PH(ICB+2))
+      SIG(ICB)=SIG(ICB+1)
+      SIG(ICB-1)=SIG(ICB)
+cjyg1
+c sb3d      print *, 'Cloud base, c. top, CAPE',ICB,INB,cape
+c sb3d      print *, 'SIG ',(sig(i),i=1,inb)
+c sb3d      print *, 'W ',(w0(i),i=1,inb)
+c sb3d      print *, 'M ',(m(i), i=1,inb)
+c sb3d      print *, 'Dt1 ',(tvp(i)-tv(i),i=1,inb)
+c sb3d      print *, 'Dt_vrai ',(buoy(i),i=1,inb)
+Cjyg2
+C
+C   ***  CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING  ***
+C   ***     RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING     ***
+C   ***                        FRACTION (SIJ)                          ***
+C
+ 
+ 
+      DO 170 I=ICB,INB
+         RTI=RR(1)-EP(I)*CLW(I)
+         DO 160 J=ICB-1,INB
+            BF2=1.+LV(J)*LV(J)*RS(J)/(RV*T(J)*T(J)*CPD)
+            ANUM=H(J)-HP(I)+(CPV-CPD)*T(J)*(RTI-RR(J))
+            DENOM=H(I)-HP(I)+(CPD-CPV)*(RR(I)-RTI)*T(J)
+            DEI=DENOM
+            IF(ABS(DEI).LT.0.01)DEI=0.01
+            SIJ(I,J)=ANUM/DEI
+            SIJ(I,I)=1.0
+            ALTEM=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI-RS(J)
+            ALTEM=ALTEM/BF2
+            CWAT=CLW(J)*(1.-EP(J))
+            STEMP=SIJ(I,J)
+            IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR.
+     1      ALTEM.GT.CWAT).AND.J.GT.I)THEN
+            ANUM=ANUM-LV(J)*(RTI-RS(J)-CWAT*BF2)
+            DENOM=DENOM+LV(J)*(RR(I)-RTI)
+            IF(ABS(DENOM).LT.0.01)DENOM=0.01
+            SIJ(I,J)=ANUM/DENOM
+            ALTEM=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI-RS(J)
+            ALTEM=ALTEM-(BF2-1.)*CWAT
+            END IF
+ 
+ 
+            IF(SIJ(I,J).GT.0.0.AND.SIJ(I,J).LT.0.95)THEN
+               QENT(I,J)=SIJ(I,J)*RR(I)+(1.-SIJ(I,J))*RTI
+               UENT(I,J)=SIJ(I,J)*U(I)+(1.-SIJ(I,J))*U(NK)
+               VENT(I,J)=SIJ(I,J)*V(I)+(1.-SIJ(I,J))*V(NK)
+               DO K=1,NTRA
+               TRAENT(I,J,K)=SIJ(I,J)*TRA(I,K)+(1.-SIJ(I,J))*
+     1         TRA(NK,K)
+               END DO
+               ELIJ(I,J)=ALTEM
+               ELIJ(I,J)=AMAX1(0.0,ELIJ(I,J))
+               MENT(I,J)=M(I)/(1.-SIJ(I,J))
+               NENT(I)=NENT(I)+1
+            END IF
+            SIJ(I,J)=AMAX1(0.0,SIJ(I,J))
+            SIJ(I,J)=AMIN1(1.0,SIJ(I,J))
+  160    CONTINUE
+C
+C   ***   IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS  ***
+C   ***   AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES  ***
+C
+         IF(NENT(I).EQ.0)THEN
+            MENT(I,I)=M(I)
+            QENT(I,I)=RR(NK)-EP(I)*CLW(I)
+           UENT(I,I)=U(NK)
+           VENT(I,I)=V(NK)
+           DO J=1,NTRA
+            TRAENT(I,I,J)=TRA(NK,J)
+           END DO
+            ELIJ(I,I)=CLW(I)
+            SIJ(I,I)=1.0
+         END IF
+C
+         DO J = ICB-1,INB
+           SIGIJ(I,J) = SIJ(I,J)
+         ENDDO
+C	
+  170 CONTINUE
+C
+C   ***  NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL  ***
+C   ***              PROBABILITIES OF MIXING                     ***
+C
+ 
+      DO 200 I=ICB,INB
+      IF(NENT(I).NE.0)THEN
+       QP=RR(1)-EP(I)*CLW(I)
+       ANUM=H(I)-HP(I)-LV(I)*(QP-RS(I))+(CPV-CPD)*T(I)*
+     1    (QP-RR(I))
+       DENOM=H(I)-HP(I)+LV(I)*(RR(I)-QP)+
+     1    (CPD-CPV)*T(I)*(RR(I)-QP)
+       IF(ABS(DENOM).LT.0.01)DENOM=0.01
+       SCRIT=ANUM/DENOM
+       ALT=QP-RS(I)+SCRIT*(RR(I)-QP)
+       IF(SCRIT.LE.0.0.OR.ALT.LE.0.0)SCRIT=1.0
+       SMAX=0.0
+       ASIJ=0.0
+        DO 175 J=INB,ICB-1,-1
+        IF(SIJ(I,J).GT.1.0E-16.AND.SIJ(I,J).LT.0.95)THEN
+         WGH=1.0
+         IF(J.GT.I)THEN
+          SJMAX=AMAX1(SIJ(I,J+1),SMAX)
+          SJMAX=AMIN1(SJMAX,SCRIT)
+          SMAX=AMAX1(SIJ(I,J),SMAX)
+          SJMIN=AMAX1(SIJ(I,J-1),SMAX)
+          SJMIN=AMIN1(SJMIN,SCRIT)
+          IF(SIJ(I,J).LT.(SMAX-1.0E-16))WGH=0.0
+          SMID=AMIN1(SIJ(I,J),SCRIT)
+         ELSE
+          SJMAX=AMAX1(SIJ(I,J+1),SCRIT)
+          SMID=AMAX1(SIJ(I,J),SCRIT)
+          SJMIN=0.0
+          IF(J.GT.1)SJMIN=SIJ(I,J-1)
+          SJMIN=AMAX1(SJMIN,SCRIT)
+         END IF
+         DELP=ABS(SJMAX-SMID)
+         DELM=ABS(SJMIN-SMID)
+         ASIJ=ASIJ+WGH*(DELP+DELM)
+         MENT(I,J)=MENT(I,J)*(DELP+DELM)*WGH
+        END IF
+  175       CONTINUE
+       ASIJ=AMAX1(1.0E-16,ASIJ)
+       ASIJ=1.0/ASIJ
+       DO 180 J=ICB-1,INB
+        MENT(I,J)=MENT(I,J)*ASIJ
+  180    CONTINUE
+       ASUM=0.0
+       BSUM=0.0
+       DO 190 J=ICB-1,INB
+        ASUM=ASUM+MENT(I,J)
+        MENT(I,J)=MENT(I,J)*SIG(J)
+        BSUM=BSUM+MENT(I,J)
+  190       CONTINUE
+       BSUM=AMAX1(BSUM,1.0E-16)
+       BSUM=1.0/BSUM
+       DO 195 J=ICB-1,INB
+        MENT(I,J)=MENT(I,J)*ASUM*BSUM	
+  195       CONTINUE
+       CSUM=0.0
+       DO 197 J=ICB-1,INB
+        CSUM=CSUM+MENT(I,J)
+  197       CONTINUE
+ 
+       IF(CSUM.LT.M(I))THEN
+        NENT(I)=0
+        MENT(I,I)=M(I)
+        QENT(I,I)=RR(1)-EP(I)*CLW(I)
+          UENT(I,I)=U(NK)
+          VENT(I,I)=V(NK)
+          DO J=1,NTRA
+           TRAENT(I,I,J)=TRA(NK,J)
+          END DO
+        ELIJ(I,I)=CLW(I)
+        SIJ(I,I)=1.0
+       END IF
+      END IF
+  200      CONTINUE
+ 
+ 
+ 
+***************************************************************
+**       CALCUL DES MENTS(I,J) ET DES QENTS(I,J)
+**************************************************************
+ 
+         DO im=1,nd
+         do jm=1,nd
+ 
+         QENTS(im,jm)=QENT(im,jm)
+         MENTS(im,jm)=MENT(im,jm)
+         enddo
+         enddo
+ 
+***********************************************************
+ 
+ 
+ 
+C
+C   ***  CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING    ***
+C   ***             DOWNDRAFT CALCULATION                      ***
+C
+        IF(EP(INB).LT.0.0001)GOTO 405
+C
+C   ***  INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER   ***
+C   ***                AND CONDENSED WATER FLUX                    ***
+C
+        WFLUX=0.0
+        TINV=1./3.
+C
+C    ***                    BEGIN DOWNDRAFT LOOP                    ***
+C
+        DO 400 I=INB,1,-1
+C
+C    ***              CALCULATE DETRAINED PRECIPITATION             ***
+C
+ 
+ 
+        WDTRAIN=10.0*EP(I)*M(I)*CLW(I)
+        IF(I.GT.1)THEN
+         DO 320 J=1,I-1
+       AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I)
+       AWAT=AMAX1(AWAT,0.0)
+  320    WDTRAIN=WDTRAIN+10.0*AWAT*MENT(J,I)
+        END IF
+C
+C    ***    FIND RAIN WATER AND EVAPORATION USING PROVISIONAL   ***
+C    ***              ESTIMATES OF RP(I)AND RP(I-1)             ***
+C
+ 
+ 
+        WT(I)=45.0
+      IF(I.LT.INB)THEN
+       RP(I)=RP(I+1)+(CPD*(T(I+1)-T(I))+GZ(I+1)-GZ(I))/LV(I)
+       RP(I)=0.5*(RP(I)+RR(I))
+      END IF
+      RP(I)=AMAX1(RP(I),0.0)
+      RP(I)=AMIN1(RP(I),RS(I))
+      RP(INB)=RR(INB)
+      IF(I.EQ.1)THEN
+       AFAC=P(1)*(RS(1)-RP(1))/(1.0E4+2000.0*P(1)*RS(1))
+      ELSE
+       RP(I-1)=RP(I)+(CPD*(T(I)-T(I-1))+GZ(I)-GZ(I-1))/LV(I)
+       RP(I-1)=0.5*(RP(I-1)+RR(I-1))
+       RP(I-1)=AMIN1(RP(I-1),RS(I-1))
+       RP(I-1)=AMAX1(RP(I-1),0.0)
+       AFAC1=P(I)*(RS(I)-RP(I))/(1.0E4+2000.0*P(I)*RS(I))
+       AFAC2=P(I-1)*(RS(I-1)-RP(I-1))/(1.0E4+
+     1    2000.0*P(I-1)*RS(I-1))
+       AFAC=0.5*(AFAC1+AFAC2)
+      END IF
+      IF(I.EQ.INB)AFAC=0.0
+        AFAC=AMAX1(AFAC,0.0)
+        BFAC=1./(SIGD*WT(I))
+C
+Cjyg1
+CCC        SIGT=1.0
+CCC        IF(I.GE.ICB)SIGT=SIGP(I)
+C Prise en compte de la variation progressive de SIGT dans
+C les couches ICB et ICB-1:
+C 	Pour PLCL<PH(I+1), PR1=0 & PR2=1
+C 	Pour PLCL>PH(I),   PR1=1 & PR2=0
+C 	Pour PH(I+1)<PLCL<PH(I), PR1 est la proportion a cheval
+C    sur le nuage, et PR2 est la proportion sous la base du
+C    nuage.
+         PR1 =(PLCL-PH(I+1))/(PH(I)-PH(I+1))
+         PR1 = MAX(0.,MIN(1.,PR1))
+         PR2 = (PH(I)-PLCL)/(PH(I)-PH(I+1))
+         PR2 = MAX(0.,MIN(1.,PR2))
+         SIGT = SIGP(I)*PR1 + PR2
+c sb3d         print *,'i,sigt,pr1,pr2', i,sigt,pr1,pr2
+Cjyg2
+C
+        B6=BFAC*50.*SIGD*(PH(I)-PH(I+1))*SIGT*AFAC
+        C6=WATER(I+1)+BFAC*WDTRAIN-50.*SIGD*BFAC*
+     1   (PH(I)-PH(I+1))*EVAP(I+1)
+      IF(C6.GT.0.0)THEN
+         REVAP=0.5*(-B6+SQRT(B6*B6+4.*C6))
+         EVAP(I)=SIGT*AFAC*REVAP
+         WATER(I)=REVAP*REVAP
+      ELSE
+       EVAP(I)=-EVAP(I+1)+0.02*(WDTRAIN+SIGD*WT(I)*
+     1    WATER(I+1))/(SIGD*(PH(I)-PH(I+1)))
+      END IF
+ 
+ 
+C
+C    ***  CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER     ***
+C    ***              HYDROSTATIC APPROXIMATION                 ***
+C
+        IF(I.EQ.1)GOTO 360
+      TEVAP=AMAX1(0.0,EVAP(I))
+      DELTH=AMAX1(0.001,(TH(I)-TH(I-1)))
+      MP(I)=10.*LVCP(I)*SIGD*TEVAP*(P(I-1)-P(I))/DELTH
+C
+C    ***           IF HYDROSTATIC ASSUMPTION FAILS,             ***
+C    ***   SOLVE CUBIC DIFFERENCE EQUATION FOR DOWNDRAFT THETA  ***
+C    ***  AND MASS FLUX FROM TWO SIMULTANEOUS DIFFERENTIAL EQNS ***
+C
+      AMFAC=SIGD*SIGD*70.0*PH(I)*(P(I-1)-P(I))*
+     1   (TH(I)-TH(I-1))/(TV(I)*TH(I))
+      AMP2=ABS(MP(I+1)*MP(I+1)-MP(I)*MP(I))
+      IF(AMP2.GT.(0.1*AMFAC))THEN
+         XF=100.0*SIGD*SIGD*SIGD*(PH(I)-PH(I+1))
+         TF=B(I)-5.0*(TH(I)-TH(I-1))*T(I)/(LVCP(I)*SIGD*TH(I))
+         AF=XF*TF+MP(I+1)*MP(I+1)*TINV
+         BF=2.*(TINV*MP(I+1))**3+TINV*MP(I+1)*XF*TF+50.*
+     1    (P(I-1)-P(I))*XF*TEVAP
+         FAC2=1.0
+         IF(BF.LT.0.0)FAC2=-1.0
+         BF=ABS(BF)
+         UR=0.25*BF*BF-AF*AF*AF*TINV*TINV*TINV
+         IF(UR.GE.0.0)THEN
+          SRU=SQRT(UR)
+          FAC=1.0
+          IF((0.5*BF-SRU).LT.0.0)FAC=-1.0
+          MP(I)=MP(I+1)*TINV+(0.5*BF+SRU)**TINV+
+     1     FAC*(ABS(0.5*BF-SRU))**TINV
+         ELSE
+          D=ATAN(2.*SQRT(-UR)/(BF+1.0E-28))
+          IF(FAC2.LT.0.0)D=3.14159-D
+          MP(I)=MP(I+1)*TINV+2.*SQRT(AF*TINV)*COS(D*TINV)
+         END IF
+         MP(I)=AMAX1(0.0,MP(I))
+         B(I-1)=B(I)+100.0*(P(I-1)-P(I))*TEVAP/(MP(I)+SIGD*0.1)-
+     1    10.0*(TH(I)-TH(I-1))*T(I)/(LVCP(I)*SIGD*TH(I))
+         B(I-1)=AMAX1(B(I-1),0.0)
+      END IF
+ 
+ 
+C
+C   ***         LIMIT MAGNITUDE OF MP(I) TO MEET CFL CONDITION      ***
+C
+      AMPMAX=2.0*(PH(I)-PH(I+1))*DELTI
+      AMP2=2.0*(PH(I-1)-PH(I))*DELTI
+      AMPMAX=AMIN1(AMPMAX,AMP2)
+      MP(I)=AMIN1(MP(I),AMPMAX)
+C
+C    ***      FORCE MP TO DECREASE LINEARLY TO ZERO                 ***
+C    ***       BETWEEN CLOUD BASE AND THE SURFACE                   ***
+C
+          IF(P(I).GT.P(ICB))THEN
+           MP(I)=MP(ICB)*(P(1)-P(I))/(P(1)-P(ICB))
+          END IF
+  360   CONTINUE
+C
+C    ***       FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT     ***
+C
+        IF(I.EQ.INB)GOTO 400
+      RP(I)=RR(I)
+        IF(MP(I).GT.MP(I+1))THEN
+        RP(I)=RP(I+1)*MP(I+1)+RR(I)*(MP(I)-MP(I+1))+
+     1       5.*SIGD*(PH(I)-PH(I+1))*(EVAP(I+1)+EVAP(I))
+        RP(I)=RP(I)/MP(I)
+          UP(I)=UP(I+1)*MP(I+1)+U(I)*(MP(I)-MP(I+1))
+         UP(I)=UP(I)/MP(I)
+          VP(I)=VP(I+1)*MP(I+1)+V(I)*(MP(I)-MP(I+1))
+         VP(I)=VP(I)/MP(I)
+          DO J=1,NTRA
+           TRAP(I,J)=TRAP(I+1,J)*MP(I+1)+
+     s     TRAP(I,J)*(MP(I)-MP(I+1))
+           TRAP(I,J)=TRAP(I,J)/MP(I)
+          END DO
+        ELSE
+        IF(MP(I+1).GT.1.0E-16)THEN
+           RP(I)=RP(I+1)+5.0*SIGD*(PH(I)-PH(I+1))*(EVAP(I+1)+
+     1      EVAP(I))/MP(I+1)
+            UP(I)=UP(I+1)
+            VP(I)=VP(I+1)
+            DO J=1,NTRA
+             TRAP(I,J)=TRAP(I+1,J)
+            END DO
+        END IF
+        END IF
+      RP(I)=AMIN1(RP(I),RS(I))
+      RP(I)=AMAX1(RP(I),0.0)
+  400   CONTINUE
+C
+C   ***  CALCULATE SURFACE PRECIPITATION IN MM/DAY     ***
+C
+        PRECIP=WT(1)*SIGD*WATER(1)*8640.0
+  405   CONTINUE
+C
+C   ***  CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE  ***
+C   ***                      AND MIXING RATIO                        ***
+C
+      DPINV=1.0/(PH(1)-PH(2))
+        AM=0.0
+        DO 410 K=2,INB
+  410   AM=AM+M(K)
+      IF((0.1*DPINV*AM).GE.DELTI)IFLAG=4
+      FT(1)=0.1*DPINV*AM*(T(2)-T(1)+(GZ(2)-GZ(1))/CPN(1))
+        FT(1)=FT(1)-0.5*LVCP(1)*SIGD*(EVAP(1)+EVAP(2))
+        FT(1)=FT(1)-0.09*SIGD*MP(2)*T(1)*B(1)*DPINV
+      FT(1)=FT(1)+0.01*SIGD*WT(1)*(CL-CPD)*WATER(2)*(T(2)-
+     1   T(1))*DPINV/CPN(1)
+        FR(1)=0.1*MP(2)*(RP(2)-RR(1))*
+     1    DPINV+SIGD*0.5*(EVAP(1)+EVAP(2))
+        FR(1)=FR(1)+0.1*AM*(RR(2)-RR(1))*DPINV
+        FU(1)=FU(1)+0.1*DPINV*(MP(2)*(UP(2)-U(1))+AM*(U(2)-U(1)))
+        FV(1)=FV(1)+0.1*DPINV*(MP(2)*(VP(2)-V(1))+AM*(V(2)-V(1)))
+        DO J=1,NTRA
+         FTRA(1,J)=FTRA(1,J)+0.1*DPINV*(MP(2)*(TRAP(2,J)-TRA(1,J))+
+     1    AM*(TRA(2,J)-TRA(1,J)))
+        END DO
+        AMDE=0.0
+        DO 415 J=2,INB
+         FR(1)=FR(1)+0.1*DPINV*MENT(J,1)*(QENT(J,1)-RR(1))
+         FU(1)=FU(1)+0.1*DPINV*MENT(J,1)*(UENT(J,1)-U(1))
+         FV(1)=FV(1)+0.1*DPINV*MENT(J,1)*(VENT(J,1)-V(1))
+         DO K=1,NTRA
+          FTRA(1,K)=FTRA(1,K)+0.1*DPINV*MENT(J,1)*(TRAENT(J,1,K)-
+     1     TRA(1,K))
+         END DO
+  415      CONTINUE
+C
+C   ***  CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO  ***
+C   ***               AT LEVELS ABOVE THE LOWEST LEVEL                   ***
+C
+C   ***  FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES  ***
+C   ***                      THROUGH EACH LEVEL                          ***
+C
+ 
+ 
+        DO 500 I=2,INB
+        DPINV=1.0/(PH(I)-PH(I+1))
+      CPINV=1.0/CPN(I)
+        AMP1=0.0
+        DO 440 K=I+1,INB+1
+  440   AMP1=AMP1+M(K)
+        DO 450 K=1,I
+        DO 450 J=I+1,INB+1
+         AMP1=AMP1+MENT(K,J)
+  450   CONTINUE
+      IF((0.1*DPINV*AMP1).GE.DELTI)IFLAG=4
+        AD=0.0
+        DO 470 K=1,I-1
+        DO 470 J=I,INB
+  470   AD=AD+MENT(J,K)
+      FT(I)=0.1*DPINV*(AMP1*(T(I+1)-T(I)+(GZ(I+1)-GZ(I))*
+     1   CPINV)-AD*(T(I)-T(I-1)+(GZ(I)-GZ(I-1))*CPINV))
+     2   -0.5*SIGD*LVCP(I)*(EVAP(I)+EVAP(I+1))
+      RAT=CPN(I-1)*CPINV
+        FT(I)=FT(I)-0.09*SIGD*(MP(I+1)*T(I)*
+     1    B(I)-MP(I)*T(I-1)*RAT*B(I-1))*DPINV
+      FT(I)=FT(I)+0.1*DPINV*MENT(I,I)*(HP(I)-H(I)+
+     1    T(I)*(CPV-CPD)*(RR(I)-QENT(I,I)))*CPINV
+      FT(I)=FT(I)+0.01*SIGD*WT(I)*(CL-CPD)*WATER(I+1)*
+     1    (T(I+1)-T(I))*DPINV*CPINV
+        FR(I)=0.1*DPINV*(AMP1*(RR(I+1)-RR(I))-
+     1    AD*(RR(I)-RR(I-1)))
+        FU(I)=FU(I)+0.1*DPINV*(AMP1*(U(I+1)-U(I))-
+     1    AD*(U(I)-U(I-1)))
+        FV(I)=FV(I)+0.1*DPINV*(AMP1*(V(I+1)-V(I))-
+     1    AD*(V(I)-V(I-1)))
+        DO K=1,NTRA
+         FTRA(I,K)=FTRA(I,K)+0.1*DPINV*(AMP1*(TRA(I+1,K)-
+     1    TRA(I,K))-AD*(TRA(I,K)-TRA(I-1,K)))
+        END DO
+        DO 480 K=1,I-1
+       AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I)
+       AWAT=AMAX1(AWAT,0.0)
+         FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-AWAT
+     1    -RR(I))
+         FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
+         FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
+         DO J=1,NTRA
+          FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
+     1     TRA(I,J))
+         END DO
+  480   CONTINUE
+      DO 490 K=I,INB
+       FR(I)=FR(I)+0.1*DPINV*MENT(K,I)*(QENT(K,I)-RR(I))
+         FU(I)=FU(I)+0.1*DPINV*MENT(K,I)*(UENT(K,I)-U(I))
+         FV(I)=FV(I)+0.1*DPINV*MENT(K,I)*(VENT(K,I)-V(I))
+         DO J=1,NTRA
+          FTRA(I,J)=FTRA(I,J)+0.1*DPINV*MENT(K,I)*(TRAENT(K,I,J)-
+     1     TRA(I,J))
+         END DO
+  490      CONTINUE
+        FR(I)=FR(I)+0.5*SIGD*(EVAP(I)+EVAP(I+1))+0.1*(MP(I+1)*
+     1    (RP(I+1)-RR(I))-MP(I)*(RP(I)-RR(I-1)))*DPINV
+        FU(I)=FU(I)+0.1*(MP(I+1)*(UP(I+1)-U(I))-MP(I)*
+     1    (UP(I)-U(I-1)))*DPINV
+        FV(I)=FV(I)+0.1*(MP(I+1)*(VP(I+1)-V(I))-MP(I)*
+     1    (VP(I)-V(I-1)))*DPINV
+        DO J=1,NTRA
+         FTRA(I,J)=FTRA(I,J)+0.1*DPINV*(MP(I+1)*(TRAP(I+1,J)-TRA(I,J))-
+     1    MP(I)*(TRAP(I,J)-TRAP(I-1,J)))
+        END DO
+  500   CONTINUE
+ 
+ 
+ 
+C
+C   ***   MOVE THE DETRAINMENT AT LEVEL INB DOWN TO LEVEL INB-1   ***
+C   ***        IN SUCH A WAY AS TO PRESERVE THE VERTICALLY        ***
+C   ***          INTEGRATED ENTHALPY AND WATER TENDENCIES         ***
+C
+      AX=0.1*MENT(INB,INB)*(HP(INB)-H(INB)+T(INB)*
+     1    (CPV-CPD)*(RR(INB)-QENT(INB,INB)))/(CPN(INB)*
+     2    (PH(INB)-PH(INB+1)))
+      FT(INB)=FT(INB)-AX
+      FT(INB-1)=FT(INB-1)+AX*CPN(INB)*(PH(INB)-PH(INB+1))/
+     1    (CPN(INB-1)*(PH(INB-1)-PH(INB)))
+      BX=0.1*MENT(INB,INB)*(QENT(INB,INB)-RR(INB))/
+     1    (PH(INB)-PH(INB+1))
+      FR(INB)=FR(INB)-BX
+      FR(INB-1)=FR(INB-1)+BX*(PH(INB)-PH(INB+1))/
+     1    (PH(INB-1)-PH(INB))
+      CX=0.1*MENT(INB,INB)*(UENT(INB,INB)-U(INB))/
+     1    (PH(INB)-PH(INB+1))
+      FU(INB)=FU(INB)-CX
+      FU(INB-1)=FU(INB-1)+CX*(PH(INB)-PH(INB+1))/
+     1    (PH(INB-1)-PH(INB))
+      DX=0.1*MENT(INB,INB)*(VENT(INB,INB)-V(INB))/
+     1    (PH(INB)-PH(INB+1))
+      FV(INB)=FV(INB)-DX
+      FV(INB-1)=FV(INB-1)+DX*(PH(INB)-PH(INB+1))/
+     1    (PH(INB-1)-PH(INB))
+      DO J=1,NTRA
+      EX=0.1*MENT(INB,INB)*(TRAENT(INB,INB,J)
+     1    -TRA(INB,J))/(PH(INB)-PH(INB+1))
+      FTRA(INB,J)=FTRA(INB,J)-EX
+      FTRA(INB-1,J)=FTRA(INB-1,J)+EX*
+     1     (PH(INB)-PH(INB+1))/(PH(INB-1)-PH(INB))
+      ENDDO   
+C
+C   ***    HOMOGINIZE TENDENCIES BELOW CLOUD BASE    ***
+C
+      ASUM=0.0
+      BSUM=0.0
+      CSUM=0.0
+        DSUM=0.0
+      DO 650 I=1,ICB-1
+       ASUM=ASUM+FT(I)*(PH(I)-PH(I+1))
+         BSUM=BSUM+FR(I)*(LV(I)+(CL-CPD)*(T(I)-T(1)))*
+     1    (PH(I)-PH(I+1))
+       CSUM=CSUM+(LV(I)+(CL-CPD)*(T(I)-T(1)))*(PH(I)-PH(I+1))
+       DSUM=DSUM+T(I)*(PH(I)-PH(I+1))/TH(I)
+  650      CONTINUE
+      DO 700 I=1,ICB-1
+       FT(I)=ASUM*T(I)/(TH(I)*DSUM)
+       FR(I)=BSUM/CSUM
+  700      CONTINUE
+C
+C   ***           RESET COUNTER AND RETURN           ***
+C
+      SIG(ND)=2.0
+c
+c
+      do i = 1, nd
+         upwd(i) = 0.0
+         dnwd(i) = 0.0
+c sb       dnwd0(i) = - mp(i)
+      enddo
+c
+      do i = 1, nl
+       dnwd0(i) = - mp(i)
+      enddo
+      do i = nl+1, nd
+       dnwd0(i) = 0.
+      enddo
+c
+      do i = icb, inb
+         upwd(i) = 0.0
+         dnwd(i) = 0.0
+
+         do k =i, inb
+            up1=0.0
+            dn1=0.0
+            do n = 1, i-1
+               up1 = up1 + ment(n,k)
+               dn1 = dn1 - ment(k,n)
+            enddo
+            upwd(i) = upwd(i)+ m(k) + up1
+            dnwd(i) = dnwd(i) + dn1
+         enddo
+        enddo
+ 
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c        DETERMINATION DE LA VARIATION DE FLUX ASCENDANT ENTRE
+C        DEUX NIVEAU NON DILUE Mike
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ 
+ 
+c sb      do i=1,ND
+c sb      Mike(i)=M(i)
+c sb      enddo
+ 
+      do i = 1, NL
+       Mike(i) = M(i)
+      enddo
+      do i = NL+1, ND
+       Mike(i) = 0.
+      enddo
+ 
+      do i=1,nd
+      Ma(i)=0
+      enddo
+ 
+c sb      do i=1,nd
+c sb      do j=i,nd
+c sb      Ma(i)=Ma(i)+M(j)
+c sb      enddo
+c sb      enddo
+
+      do i = 1, NL
+      do j = i, NL
+       Ma(i) = Ma(i) + M(j)
+      enddo
+      enddo
+c
+      do i = NL+1, ND
+       Ma(i) = 0.
+      enddo
+c 
+      do i=1,ICB-1
+      Ma(i)=0
+      enddo
+ 
+ 
+ 
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+C        ICB REPRESENTE DE NIVEAU OU SE TROUVE LA
+c        BASE DU NUAGE , ET INB LE TOP DU NUAGE
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ 
+ 
+       do i=1,ND
+       Mke(i)=upwd(i)+dnwd(i)
+       enddo
+c$$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c$$$         call writeg1d(1,klev,ma,'ma  ','ma  ')
+c$$$          call writeg1d(1,klev,upwd,'upwd  ','upwd  ')
+c$$$          call writeg1d(1,klev,dnwd,'dnwd  ','dnwd  ')
+c$$$          call writeg1d(1,klev,dnwd0,'dnwd0  ','dnwd0  ')
+c$$$          call writeg1d(1,klev,tvp,'tvp  ','tvp  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,3),'tra3  ','tra3  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,4),'tra4  ','tra4  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,5),'tra5  ','tra5  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,6),'tra6  ','tra6  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,7),'tra7  ','tra7  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,8),'tra8  ','tra8  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,9),'tra9  ','tra9  ')
+c$$$          call writeg1d(1,klev,tra(1:klev,10),'tra10','tra10')
+c$$$          call writeg1d(1,klev,tra(1:klev,11),'tra11','tra11')
+c$$$          call writeg1d(1,klev,tra(1:klev,12),'tra12','tra12')
+c$$$          call writeg1d(1,klev,tra(1:klev,13),'tra13','tra13')
+c$$$          call writeg1d(1,klev,tra(1:klev,14),'tra14','tra14')
+c$$$          call writeg1d(1,klev,tra(1:klev,15),'tra15','tra15')
+c$$$          call writeg1d(1,klev,tra(1:klev,16),'tra16','tra16')
+c$$$          call writeg1d(1,klev,tra(1:klev,17),'tra17','tra17')
+c$$$          call writeg1d(1,klev,tra(1:klev,18),'tra18','tra18')
+c$$$          call writeg1d(1,klev,tra(1:klev,19),'tra19','tra19')
+c$$$          call writeg1d(1,klev,tra(1:klev,20),'tra20','tra20 ')
+c$$$          call writeg1d(1,klev,trap(1:klev,1),'trp1','trp1')
+c$$$          call writeg1d(1,klev,trap(1:klev,2),'trp2','trp2')
+c$$$          call writeg1d(1,klev,trap(1:klev,3),'trp3','trp3')
+c$$$          call writeg1d(1,klev,trap(1:klev,4),'trp4','trp4')
+c$$$          call writeg1d(1,klev,trap(1:klev,5),'trp5','trp5')
+c$$$          call writeg1d(1,klev,trap(1:klev,10),'trp10','trp10')
+c$$$          call writeg1d(1,klev,trap(1:klev,12),'trp12','trp12')
+c$$$          call writeg1d(1,klev,trap(1:klev,15),'trp15','trp15')
+c$$$          call writeg1d(1,klev,trap(1:klev,20),'trp20','trp20')
+c$$$          call writeg1d(1,klev,ftra(1:klev,1),'ftr1  ','ftr1  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,2),'ftr2  ','ftr2  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,3),'ftr3  ','ftr3  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,4),'ftr4  ','ftr4  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,5),'ftr5  ','ftr5  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,6),'ftr6  ','ftr6  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,7),'ftr7  ','ftr7  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,8),'ftr8  ','ftr8  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,9),'ftr9  ','ftr9  ')
+c$$$          call writeg1d(1,klev,ftra(1:klev,10),'ftr10','ftr10')
+c$$$          call writeg1d(1,klev,ftra(1:klev,11),'ftr11','ftr11')
+c$$$          call writeg1d(1,klev,ftra(1:klev,12),'ftr12','ftr12')
+c$$$          call writeg1d(1,klev,ftra(1:klev,13),'ftr13','ftr13')
+c$$$          call writeg1d(1,klev,ftra(1:klev,14),'ftr14','ftr14')
+c$$$          call writeg1d(1,klev,ftra(1:klev,15),'ftr15','ftr15')
+c$$$          call writeg1d(1,klev,ftra(1:klev,16),'ftr16','ftr16')
+c$$$          call writeg1d(1,klev,ftra(1:klev,17),'ftr17','ftr17')
+c$$$          call writeg1d(1,klev,ftra(1:klev,18),'ftr18','ftr18')
+c$$$          call writeg1d(1,klev,ftra(1:klev,19),'ftr19','ftr19')
+c$$$          call writeg1d(1,klev,ftra(1:klev,20),'ftr20','ftr20 ')
+c$$$          call writeg1d(1,klev,mp,'mp  ','mp ')
+c$$$          call writeg1d(1,klev,Mke,'Mke  ','Mke ')
+
+ 
+ 
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c
+        RETURN
+        END
+C ---------------------------------------------------------------------------
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/read_pstoke.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/read_pstoke.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/read_pstoke.F	(revision 222)
@@ -0,0 +1,431 @@
+c
+c $Header$
+c
+	subroutine read_pstoke(irec,
+     .   zrec,zklono,zklevo,airefi,phisfi,
+     .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,
+     .   frac_impa,frac_nucl,pyu1,pyv1,ftsol,psrf)
+
+
+       IMPLICIT NONE
+
+#include "netcdf.inc"
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comgeom.h"
+#include "temps.h"
+#include "ener.h"
+#include "logic.h"
+#include "description.h"
+#include "serre.h"
+#include "indicesol.h"
+#include "control.h"
+#include "dimphy.h"
+	
+	  integer*4 klono,klevo,imo,jmo
+	  parameter (imo=iim/2,jmo=(jjm+1)/2)
+	  parameter(klono=(jmo-1)*imo+2,klevo=llm)
+	  REAL*4 phisfi(klono)
+          REAL*4 phisfi2(imo,jmo+1),airefi2(imo,jmo+1)
+
+          REAL*4 mfu(klono,klevo), mfd(klono,klevo)
+          REAL*4 en_u(klono,klevo), de_u(klono,klevo)
+          REAL*4 en_d(klono,klevo), de_d(klono,klevo)
+          REAL*4 coefh(klono,klevo)
+
+          REAL*4 mfu2(imo,jmo+1,klevo), mfd2(imo,jmo+1,klevo)
+          REAL*4 en_u2(imo,jmo+1,klevo), de_u2(imo,jmo+1,klevo)
+          REAL*4 en_d2(imo,jmo+1,klevo), de_d2(imo,jmo+1,klevo)
+          REAL*4 coefh2(imo,jmo+1,klevo)
+
+          REAL*4 pl(klevo)
+          integer irec
+          integer*4 xid,yid,zid,tid
+          real zrec,zklono,zklevo,zim,zjm
+          integer*4 ncrec,ncklono,ncklevo,ncim,ncjm
+
+          real*4 airefi(klono)
+          character namedim
+
+c  !! attention !!
+c attention il y a aussi le pb de def klono
+c dim de phis??
+	  
+	 
+          REAL*4 frac_impa(klono,klevo), frac_nucl(klono,klevo)
+          REAL*4 frac_impa2(imo,jmo+1,klevo), 
+     .     frac_nucl2(imo,jmo+1,klevo)
+          REAL*4 pyu1(klono), pyv1(klono)
+          REAL*4 pyu12(imo,jmo+1), pyv12(imo,jmo+1)
+          REAL*4 ftsol(klono,nbsrf)
+          REAL*4 psrf(klono,nbsrf) 
+	  REAL*4 ftsol1(klono),ftsol2(klono),ftsol3(klono),ftsol4(klono)
+          REAL*4 psrf1(klono),psrf2(klono),psrf3(klono),psrf4(klono)
+          REAL*4 ftsol12(imo,jmo+1),ftsol22(imo,jmo+1),
+     .     ftsol32(imo,jmo+1),
+     .     ftsol42(imo,jmo+1)
+          REAL*4 psrf12(imo,jmo+1),psrf22(imo,jmo+1),psrf32(imo,jmo+1),
+     .     psrf42(imo,jmo+1)
+		REAL*4 t(klono,klevo)
+		REAL*4 t2(imo,jmo+1)	
+	  integer ncidp
+          save ncidp
+		integer varidt
+          integer varidmfu, varidmfd, varidps, varidenu, variddeu	
+          integer varidend,varidded,varidch,varidfi,varidfn
+          integer varidyu1,varidyv1,varidpl,varidai,varididvt
+          integer varidfts1,varidfts2,varidfts3,varidfts4
+          integer varidpsr1,varidpsr2,varidpsr3,varidpsr4
+          save varidmfu, varidmfd, varidps, varidenu, variddeu
+          save varidend,varidded,varidch,varidfi,varidfn
+          save varidyu1,varidyv1,varidpl,varidai,varididvt
+          save varidfts1,varidfts2,varidfts3,varidfts4
+          save varidpsr1,varidpsr2,varidpsr3,varidpsr4
+		save varidt
+
+          integer l, i
+          integer start(4),count(4),status
+          real rcode
+          logical first
+          save first
+          data first/.true./
+
+
+
+c ---------------------------------------------
+c   Initialisation de la lecture des fichiers
+c ---------------------------------------------
+
+      if (irec .eq. 0) then
+
+            ncidp=NCOPN('phystoke.nc',NCNOWRIT,rcode)
+
+            varidps=NCVID(ncidp,'phis',rcode)
+            print*,'ncidp,varidps',ncidp,varidps
+
+            varidpl=NCVID(ncidp,'sig_s',rcode)
+            print*,'ncidp,varidpl',ncidp,varidpl
+
+            varidai=NCVID(ncidp,'aire',rcode)
+            print*,'ncidp,varidai',ncidp,varidai
+
+	        varidt=NCVID(ncidp,'t',rcode)
+                print*,'ncidp,varidt',ncidp,varidt
+            varidmfu=NCVID(ncidp,'mfu',rcode)
+            print*,'ncidp,varidmfu',ncidp,varidmfu
+
+            varidmfd=NCVID(ncidp,'mfd',rcode)
+            print*,'ncidp,varidmfd',ncidp,varidmfd
+
+            varidenu=NCVID(ncidp,'en_u',rcode)
+            print*,'ncidp,varidenu',ncidp,varidenu
+
+            variddeu=NCVID(ncidp,'de_u',rcode)
+            print*,'ncidp,variddeu',ncidp,variddeu
+
+            varidend=NCVID(ncidp,'en_d',rcode)
+            print*,'ncidp,varidend',ncidp,varidend
+	
+            varidded=NCVID(ncidp,'de_d',rcode)
+            print*,'ncidp,varidded',ncidp,varidded
+	
+            varidch=NCVID(ncidp,'coefh',rcode)
+            print*,'ncidp,varidch',ncidp,varidch
+	
+	    varidfi=NCVID(ncidp,'frac_impa',rcode)
+            print*,'ncidp,varidfi',ncidp,varidfi
+	
+	    varidfn=NCVID(ncidp,'frac_nucl',rcode)
+            print*,'ncidp,varidfn',ncidp,varidfn
+	
+            varidyu1=NCVID(ncidp,'pyu1',rcode)
+            print*,'ncidp,varidyu1',ncidp,varidyu1
+	
+            varidyv1=NCVID(ncidp,'pyv1',rcode)
+            print*,'ncidp,varidyv1',ncidp,varidyv1
+	
+            varidfts1=NCVID(ncidp,'ftsol1',rcode)
+            print*,'ncidp,varidfts1',ncidp,varidfts1
+	
+            varidfts2=NCVID(ncidp,'ftsol2',rcode)
+            print*,'ncidp,varidfts2',ncidp,varidfts2
+         
+            varidfts3=NCVID(ncidp,'ftsol3',rcode)
+            print*,'ncidp,varidfts3',ncidp,varidfts3
+  
+            varidfts4=NCVID(ncidp,'ftsol4',rcode)
+            print*,'ncidp,varidfts4',ncidp,varidfts4
+	
+            varidpsr1=NCVID(ncidp,'psrf1',rcode)
+            print*,'ncidp,varidpsr1',ncidp,varidpsr1
+	
+            varidpsr2=NCVID(ncidp,'psrf2',rcode)
+            print*,'ncidp,varidpsr2',ncidp,varidpsr2
+	
+	    varidpsr3=NCVID(ncidp,'psrf3',rcode)
+            print*,'ncidp,varidpsr3',ncidp,varidpsr3
+
+            varidpsr4=NCVID(ncidp,'psrf4',rcode)
+            print*,'ncidp,varidpsr4',ncidp,varidpsr4
+	
+c ID pour les dimensions
+
+            status = nf_inq_dimid(ncidp,'y',yid)
+            status = nf_inq_dimid(ncidp,'x',xid)
+            status = nf_inq_dimid(ncidp,'sig_s',zid)
+            status = nf_inq_dimid(ncidp,'time_counter',tid)
+
+c lecture des dimensions
+
+            status = nf_inq_dim(ncidp,yid,namedim,ncjm)
+            status = nf_inq_dim(ncidp,xid,namedim,ncim)
+            status = nf_inq_dim(ncidp,zid,namedim,ncklevo)
+            status = nf_inq_dim(ncidp,tid,namedim,ncrec)
+	
+            zrec=ncrec
+            zklevo=ncklevo
+            zim=ncim
+            zjm=ncjm
+	
+	    zklono=zim*(zjm-2)+2
+	
+	    write(*,*) 'read_pstoke : zrec = ', zrec
+            write(*,*) 'read_pstoke : zklevo = ', zklevo
+            write(*,*) 'read_pstoke : zim = ', zim 
+            write(*,*) 'read_pstoke : zjm = ', zjm
+            write(*,*) 'read_pstoke : zklono = ', zklono
+
+c niveaux de pression
+
+            status=NF_GET_VARA_REAL(ncidp,varidpl,1,zklevo,pl)
+
+c lecture de aire et phis
+	
+      start(1)=1
+      start(2)=1
+      start(3)=1
+      start(4)=0
+
+      count(1)=zim
+      count(2)=zjm
+      count(3)=1
+      count(4)=0
+
+c phis
+      status=NF_GET_VARA_REAL(ncidp,varidps,start,count,phisfi2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,phisfi2)
+c      call dump2d(iip1-1,jjp1,phisfi2,'PHISNC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,phisfi2,phisfi)
+
+c aire
+      status=NF_GET_VARA_REAL(ncidp,varidai,start,count,airefi2)
+c      call correctbid(iim,jjp1*1,airefi2)
+c       call dump2d(iip1-1,jjp1,airefi2,'AIRENC')
+       call gr_ecrit_fi(1,klono,imo,jmo+1,airefi2,airefi)
+      else
+
+      print*,'ok1'
+
+c ---------------------
+c   lecture des champs
+c ---------------------
+	
+	print*,'WARNING!!! Il n y a pas de test de coherence'
+        print*,'sur le nombre de niveaux verticaux dans le fichier nc'
+
+      start(1)=1
+      start(2)=1
+      start(3)=1
+      start(4)=irec
+
+      count(1)=zim
+      count(2)=zjm
+      count(3)=zklevo
+      count(4)=1
+
+c frac_impa 
+
+      status=NF_GET_VARA_REAL(ncidp,varidfi,start,count,frac_impa2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,frac_impa2)
+c      call dump2d(iip1-1,jjp1,frac_impa2,'FINC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_impa2,frac_impa)
+
+c frac_nucl 
+
+      status=NF_GET_VARA_REAL(ncidp,varidfn,start,count,frac_nucl2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,frac_nucl2)
+c      call dump2d(iip1-1,jjp1,frac_nucl2,'FINC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_nucl2,frac_nucl)
+
+c abder t
+      status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,t2,t)
+
+c mfu
+      status=NF_GET_VARA_REAL(ncidp,varidmfu,start,count,mfu2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,mfu2)
+c      call dump2d(iip1-1,jjp1,mfu2,'MFUNC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfu2,mfu)
+
+c mfd
+      status=NF_GET_VARA_REAL(ncidp,varidmfd,start,count,mfd2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,mfd2)
+c      call dump2d(iip1-1,jjp1,mfd2,'MFDNC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfd2,mfd)
+
+c en_u 
+      status=NF_GET_VARA_REAL(ncidp,varidenu,start,count,en_u2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,en_u2)
+c      call dump2d(iip1-1,jjp1,en_u2,'ENUNC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_u2,en_u)
+
+c de_u 
+      status=NF_GET_VARA_REAL(ncidp,variddeu,start,count,de_u2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,de_u2)
+c      call dump2d(iip1-1,jjp1,de_u2,'DEUNC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_u2,de_u)
+
+c en_d 
+      status=NF_GET_VARA_REAL(ncidp,varidend,start,count,en_d2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,en_d2)
+c      call dump2d(iip1-1,jjp1,en_d2,'ENDNC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_d2,en_d)
+
+c de_d 
+      status=NF_GET_VARA_REAL(ncidp,varidded,start,count,de_d2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,de_d2)
+c      call dump2d(iip1-1,jjp1,de_d2,'DEDNC COUCHE 1')
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_d2,de_d)
+
+c coefh 
+      status=NF_GET_VARA_REAL(ncidp,varidch,start,count,coefh2)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*klevo,coefh2)
+c      call dump2d(iip1-1,jjp1,coefh2,'CHNC COUCHE 1')
+       call gr_ecrit_fi(klevo,klono,imo,jmo+1,coefh2,coefh)
+
+      start(3)=irec
+      start(4)=0
+      count(3)=1
+      count(4)=0
+
+c pyu1
+      status=NF_GET_VARA_REAL(ncidp,varidyu1,start,count,pyu12)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,pyu12)
+c      call dump2d(iip1-1,jjp1,pyu12,'PYU1NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,pyu12,pyu1)
+
+c pyv1
+      status=NF_GET_VARA_REAL(ncidp,varidyv1,start,count,pyv12)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,pyv12)
+c      call dump2d(iip1-1,jjp1,pyv12,'PYV1NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,pyv12,pyv1)
+
+c ftsol1
+      status=NF_GET_VARA_REAL(ncidp,varidfts1,start,count,ftsol12)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,ftsol12)
+c      call dump2d(iip1-1,jjp1,ftsol12,'FTS1NC')
+       call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol12,ftsol1)
+
+c ftsol2
+      status=NF_GET_VARA_REAL(ncidp,varidfts2,start,count,ftsol22)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,ftsol22)
+c      call dump2d(iip1-1,jjp1,ftsol22,'FTS2NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol22,ftsol2)
+
+c ftsol3
+      status=NF_GET_VARA_REAL(ncidp,varidfts3,start,count,ftsol32)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,ftsol32)
+c      call dump2d(iip1-1,jjp1,ftsol32,'FTS3NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol32,ftsol3)
+
+c ftsol4
+      status=NF_GET_VARA_REAL(ncidp,varidfts4,start,count,ftsol42)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,ftsol42)
+c      call dump2d(iip1-1,jjp1,ftsol42,'FTS4NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol42,ftsol4)
+
+c psrf1 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr1,start,count,psrf12)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,psrf12)
+c      call dump2d(iip1-1,jjp1,psrf12,'PSRF1NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf12,psrf1)
+
+c psrf2 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr2,start,count,psrf22)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,psrf22)
+c      call dump2d(iip1-1,jjp1,psrf22,'PSRF2NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf22,psrf2)
+
+c psrf3 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr3,start,count,psrf32)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,psrf32)
+c      call dump2d(iip1-1,jjp1,psrf32,'PSRF3NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf32,psrf3)
+
+c psrf4 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr4,start,count,psrf42)
+c      print*,'WARNING!!! Correction bidon pour palier a un '
+c      print*,'probleme dans la creation des fichiers nc'
+c      call correctbid(iim,jjp1*1,psrf42)
+c      call dump2d(iip1-1,jjp1,psrf42,'PSRF4NC')
+      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf42,psrf4)
+	
+          do i = 1,klono
+	
+	psrf(i,1) = psrf1(i)
+        psrf(i,2) = psrf2(i)
+        psrf(i,3) = psrf3(i)
+        psrf(i,4) = psrf4(i)
+  
+        ftsol(i,1) = ftsol1(i)
+        ftsol(i,2) = ftsol2(i)
+        ftsol(i,3) = ftsol3(i)
+        ftsol(i,4) = ftsol4(i)
+	
+          enddo
+	
+	endif
+	
+	return
+	
+	end
+
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/read_pstoke0.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/read_pstoke0.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/read_pstoke0.F	(revision 222)
@@ -0,0 +1,358 @@
+c
+c $Header$
+c
+	subroutine read_pstoke0(irec,
+     .   zrec,zkon,zkev,airefi,phisfi,
+     .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,
+     .   frac_impa,frac_nucl,pyu1,pyv1,ftsol,psrf)
+
+
+       IMPLICIT NONE
+
+#include "netcdf.inc"
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comgeom.h"
+#include "temps.h"
+#include "ener.h"
+#include "logic.h"
+#include "description.h"
+#include "serre.h"
+#include "indicesol.h"
+#include "control.h"
+#include "dimphy.h"
+	  
+	  integer*4 kon,kev,zkon,zkev
+	  parameter(kon=iim*(jjm-1)+2,kev=llm)
+	  REAL*4 phisfi(kon)
+          REAL*4 phisfi2(iim,jjm+1),airefi2(iim,jjm+1)
+
+          REAL*4 mfu(kon,kev), mfd(kon,kev)
+          REAL*4 en_u(kon,kev), de_u(kon,kev)
+          REAL*4 en_d(kon,kev), de_d(kon,kev)
+          REAL*4 coefh(kon,kev)
+		REAL*4 t(kon,kev)
+
+          REAL*4 mfu2(iim,jjm+1,kev), mfd2(iim,jjm+1,kev)
+          REAL*4 en_u2(iim,jjm+1,kev), de_u2(iim,jjm+1,kev)
+          REAL*4 en_d2(iim,jjm+1,kev), de_d2(iim,jjm+1,kev)
+          REAL*4 coefh2(iim,jjm+1,kev)
+		REAL*4 t2(iim,jjm+1,kev)
+
+          REAL*4 pl(kev)
+          integer irec
+          integer*4 xid,yid,zid,tid
+          integer*4 zrec,zim,zjm
+          integer*4 ncrec,nckon,nckev,ncim,ncjm
+
+          real*4 airefi(kon)
+          character namedim
+
+c  !! attention !!
+c attention il y a aussi le pb de def kon
+c dim de phis??
+
+          REAL*4 frac_impa(kon,kev), frac_nucl(kon,kev)
+          REAL*4 frac_impa2(iim,jjm+1,kev), 
+     .     frac_nucl2(iim,jjm+1,kev)
+          REAL*4 pyu1(kon), pyv1(kon)
+          REAL*4 pyu12(iim,jjm+1), pyv12(iim,jjm+1)
+          REAL*4 ftsol(kon,nbsrf)
+          REAL*4 psrf(kon,nbsrf) 
+	  REAL*4 ftsol1(kon),ftsol2(kon),ftsol3(kon),ftsol4(kon)
+          REAL*4 psrf1(kon),psrf2(kon),psrf3(kon),psrf4(kon)
+          REAL*4 ftsol12(iim,jjm+1),ftsol22(iim,jjm+1),
+     .     ftsol32(iim,jjm+1),
+     .     ftsol42(iim,jjm+1)
+          REAL*4 psrf12(iim,jjm+1),psrf22(iim,jjm+1),psrf32(iim,jjm+1),
+     .     psrf42(iim,jjm+1)
+	
+	  integer ncidp
+          save ncidp
+          integer varidmfu, varidmfd, varidps, varidenu, variddeu	
+		integer varidt
+          integer varidend,varidded,varidch,varidfi,varidfn
+          integer varidyu1,varidyv1,varidpl,varidai,varididvt
+          integer varidfts1,varidfts2,varidfts3,varidfts4
+          integer varidpsr1,varidpsr2,varidpsr3,varidpsr4
+          save varidmfu, varidmfd, varidps, varidenu, variddeu
+		save varidt
+          save varidend,varidded,varidch,varidfi,varidfn
+          save varidyu1,varidyv1,varidpl,varidai,varididvt
+          save varidfts1,varidfts2,varidfts3,varidfts4
+          save varidpsr1,varidpsr2,varidpsr3,varidpsr4
+
+          integer l, i
+          integer start(4),count(4),status
+          real rcode
+          logical first
+          save first
+          data first/.true./
+
+
+
+c ---------------------------------------------
+c   Initialisation de la lecture des fichiers
+c ---------------------------------------------
+
+      if (irec .eq. 0) then
+
+            ncidp=NCOPN('phystoke.nc',NCNOWRIT,rcode)
+
+            varidps=NCVID(ncidp,'phis',rcode)
+            print*,'ncidp,varidps',ncidp,varidps
+
+            varidpl=NCVID(ncidp,'sig_s',rcode)
+            print*,'ncidp,varidpl',ncidp,varidpl
+
+            varidai=NCVID(ncidp,'aire',rcode)
+            print*,'ncidp,varidai',ncidp,varidai
+
+            varidmfu=NCVID(ncidp,'mfu',rcode)
+            print*,'ncidp,varidmfu',ncidp,varidmfu
+
+		varidt=NCVID(ncidp,'t',rcode)
+                print*,'ncidp,varidt',ncidp,varidt
+
+            varidmfd=NCVID(ncidp,'mfd',rcode)
+            print*,'ncidp,varidmfd',ncidp,varidmfd
+
+            varidenu=NCVID(ncidp,'en_u',rcode)
+            print*,'ncidp,varidenu',ncidp,varidenu
+
+            variddeu=NCVID(ncidp,'de_u',rcode)
+            print*,'ncidp,variddeu',ncidp,variddeu
+
+            varidend=NCVID(ncidp,'en_d',rcode)
+            print*,'ncidp,varidend',ncidp,varidend
+	
+            varidded=NCVID(ncidp,'de_d',rcode)
+            print*,'ncidp,varidded',ncidp,varidded
+	
+            varidch=NCVID(ncidp,'coefh',rcode)
+            print*,'ncidp,varidch',ncidp,varidch
+	
+	    varidfi=NCVID(ncidp,'frac_impa',rcode)
+            print*,'ncidp,varidfi',ncidp,varidfi
+	
+	    varidfn=NCVID(ncidp,'frac_nucl',rcode)
+            print*,'ncidp,varidfn',ncidp,varidfn
+	
+            varidyu1=NCVID(ncidp,'pyu1',rcode)
+            print*,'ncidp,varidyu1',ncidp,varidyu1
+	
+            varidyv1=NCVID(ncidp,'pyv1',rcode)
+            print*,'ncidp,varidyv1',ncidp,varidyv1
+	
+            varidfts1=NCVID(ncidp,'ftsol1',rcode)
+            print*,'ncidp,varidfts1',ncidp,varidfts1
+	
+            varidfts2=NCVID(ncidp,'ftsol2',rcode)
+            print*,'ncidp,varidfts2',ncidp,varidfts2
+         
+            varidfts3=NCVID(ncidp,'ftsol3',rcode)
+            print*,'ncidp,varidfts3',ncidp,varidfts3
+  
+            varidfts4=NCVID(ncidp,'ftsol4',rcode)
+            print*,'ncidp,varidfts4',ncidp,varidfts4
+	
+            varidpsr1=NCVID(ncidp,'psrf1',rcode)
+            print*,'ncidp,varidpsr1',ncidp,varidpsr1
+	
+            varidpsr2=NCVID(ncidp,'psrf2',rcode)
+            print*,'ncidp,varidpsr2',ncidp,varidpsr2
+	
+	    varidpsr3=NCVID(ncidp,'psrf3',rcode)
+            print*,'ncidp,varidpsr3',ncidp,varidpsr3
+
+            varidpsr4=NCVID(ncidp,'psrf4',rcode)
+            print*,'ncidp,varidpsr4',ncidp,varidpsr4
+	
+c ID pour les dimensions
+
+            status = nf_inq_dimid(ncidp,'y',yid)
+            status = nf_inq_dimid(ncidp,'x',xid)
+            status = nf_inq_dimid(ncidp,'sig_s',zid)
+            status = nf_inq_dimid(ncidp,'time_counter',tid)
+
+c lecture des dimensions
+
+            status = nf_inq_dim(ncidp,yid,namedim,ncjm)
+            status = nf_inq_dim(ncidp,xid,namedim,ncim)
+            status = nf_inq_dim(ncidp,zid,namedim,nckev)
+            status = nf_inq_dim(ncidp,tid,namedim,ncrec)
+	
+            zrec=ncrec
+            zkev=nckev
+            zim=ncim
+            zjm=ncjm
+	
+	    zkon=zim*(zjm-2)+2
+	
+	    write(*,*) 'read_pstoke : zrec = ', zrec
+            write(*,*) 'read_pstoke : kev = ', zkev
+            write(*,*) 'read_pstoke : zim = ', zim 
+            write(*,*) 'read_pstoke : zjm = ', zjm
+            write(*,*) 'read_pstoke : kon = ', zkon
+
+c niveaux de pression
+
+            status=NF_GET_VARA_REAL(ncidp,varidpl,1,kev,pl)
+
+c lecture de aire et phis
+	
+      start(1)=1
+      start(2)=1
+      start(3)=1
+      start(4)=0
+
+      count(1)=zim
+      count(2)=zjm
+      count(3)=1
+      count(4)=0
+
+c 
+c phis
+      status=NF_GET_VARA_REAL(ncidp,varidps,start,count,phisfi2)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,phisfi2,phisfi)
+
+c aire
+      status=NF_GET_VARA_REAL(ncidp,varidai,start,count,airefi2)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,airefi2,airefi)
+      else
+
+      print*,'ok1'
+
+c ---------------------
+c   lecture des champs
+c ---------------------
+	
+	print*,'WARNING!!! Il n y a pas de test de coherence'
+        print*,'sur le nombre de niveaux verticaux dans le fichier nc'
+
+      start(1)=1
+      start(2)=1
+      start(3)=1
+      start(4)=irec
+
+      count(1)=zim
+      count(2)=zjm
+      count(3)=kev
+      count(4)=1
+
+c frac_impa 
+
+      status=NF_GET_VARA_REAL(ncidp,varidfi,start,count,frac_impa2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,frac_impa2,frac_impa)
+
+c frac_nucl 
+
+      status=NF_GET_VARA_REAL(ncidp,varidfn,start,count,frac_nucl2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,frac_nucl2,frac_nucl)
+
+c abder t
+      status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,t2,t)
+
+c mfu
+      status=NF_GET_VARA_REAL(ncidp,varidmfu,start,count,mfu2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,mfu2,mfu)
+
+c mfd
+      status=NF_GET_VARA_REAL(ncidp,varidmfd,start,count,mfd2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,mfd2,mfd)
+
+c en_u 
+      status=NF_GET_VARA_REAL(ncidp,varidenu,start,count,en_u2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,en_u2,en_u)
+
+c de_u 
+      status=NF_GET_VARA_REAL(ncidp,variddeu,start,count,de_u2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,de_u2,de_u)
+
+c en_d 
+      status=NF_GET_VARA_REAL(ncidp,varidend,start,count,en_d2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,en_d2,en_d)
+
+c de_d 
+      status=NF_GET_VARA_REAL(ncidp,varidded,start,count,de_d2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,de_d2,de_d)
+
+c coefh 
+	print*,'LECTURE de coefh a irec =',irec
+       status=NF_GET_VARA_REAL(ncidp,varidch,start,count,coefh2)
+       call gr_ecrit_fi(kev,kon,iim,jjm+1,coefh2,coefh)
+
+      start(3)=irec
+      start(4)=0
+      count(3)=1
+      count(4)=0
+
+c pyu1
+	print*,'LECTURE de yu1 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidyu1,start,count,pyu12)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,pyu12,pyu1)
+
+c pyv1
+        print*,'LECTURE de yv1 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidyv1,start,count,pyv12)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,pyv12,pyv1)
+
+c ftsol1
+        print*,'LECTURE de ftsol1 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidfts1,start,count,ftsol12)
+       call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol12,ftsol1)
+
+c ftsol2
+        print*,'LECTURE de ftsol2 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidfts2,start,count,ftsol22)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol22,ftsol2)
+
+c ftsol3
+	 print*,'LECTURE de ftsol3 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidfts3,start,count,ftsol32)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol32,ftsol3)
+
+c ftsol4
+      status=NF_GET_VARA_REAL(ncidp,varidfts4,start,count,ftsol42)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol42,ftsol4)
+
+c psrf1 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr1,start,count,psrf12)
+c      call dump2d(iip1-1,jjm+1,psrf12,'PSRF1NC')
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf12,psrf1)
+
+c psrf2 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr2,start,count,psrf22)
+c      call dump2d(iip1-1,jjm+1,psrf22,'PSRF2NC')
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf22,psrf2)
+
+c psrf3 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr3,start,count,psrf32)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf32,psrf3)
+
+c psrf4 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr4,start,count,psrf42)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf42,psrf4)
+	
+          do i = 1,kon
+	
+	psrf(i,1) = psrf1(i)
+        psrf(i,2) = psrf2(i)
+        psrf(i,3) = psrf3(i)
+        psrf(i,4) = psrf4(i)
+  
+        ftsol(i,1) = ftsol1(i)
+        ftsol(i,2) = ftsol2(i)
+        ftsol(i,3) = ftsol3(i)
+        ftsol(i,4) = ftsol4(i)
+	
+          enddo
+	
+	endif
+	
+	return
+	
+	end
+
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/tilft43.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/tilft43.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/tilft43.F	(revision 222)
@@ -0,0 +1,93 @@
+c $Header$
+        SUBROUTINE TLIFT43(P,T,Q,QS,GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK)
+        REAL GZ(ND),TPK(ND),CLW(ND),P(ND)
+        REAL T(ND),Q(ND),QS(ND),TVP(ND),LV0
+C
+C   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
+C
+c -- sb:
+c!      CPD=1005.7
+c!      CPV=1870.0
+c!      CL=4190.0
+c!      RV=461.5
+c!      RD=287.04
+c!      LV0=2.501E6
+c!      G=9.8
+c!      ROWL=1000.0
+c ajouts:
+#include "YOMCST.h"
+        CPD = RCPD
+        CPV = RCPV
+        CL = RCW
+        LV0 = RLVTT
+        G = RG
+        ROWL= RATM/100.
+        GRAVITY = RG !sb: Pr que gravite ne devienne pas humidite!
+C sb --
+C
+        CPVMCL=CL-CPV
+        EPS=RD/RV
+        EPSI=1./EPS
+C
+C   ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
+C
+        AH0=(CPD*(1.-Q(NK))+CL*Q(NK))*T(NK)+Q(NK)*(LV0-CPVMCL*(
+     1   T(NK)-273.15))+GZ(NK)
+        CPP=CPD*(1.-Q(NK))+Q(NK)*CPV
+        CPINV=1./CPP
+C
+        IF(KK.EQ.1)THEN
+C
+C   ***   CALCULATE LIFTED PARCEL QUANTITIES BELOW CLOUD BASE   ***
+C
+        DO 50 I=1,ICB-1
+         CLW(I)=0.0
+   50   CONTINUE
+        DO 100 I=NK,ICB-1
+         TPK(I)=T(NK)-(GZ(I)-GZ(NK))*CPINV
+         TVP(I)=TPK(I)*(1.+Q(NK)*EPSI)
+  100   CONTINUE
+        END IF
+C
+C    ***  FIND LIFTED PARCEL QUANTITIES ABOVE CLOUD BASE    ***
+C
+        NST=ICB
+        NSB=ICB
+        IF(KK.EQ.2)THEN  
+         NST=NL
+         NSB=ICB+1
+        END IF
+        DO 300 I=NSB,NST
+         TG=T(I)
+         QG=QS(I)
+         ALV=LV0-CPVMCL*(T(I)-273.15)
+         DO 200 J=1,2
+          S=CPD+ALV*ALV*QG/(RV*T(I)*T(I))
+          S=1./S
+          AHG=CPD*TG+(CL-CPD)*Q(NK)*T(I)+ALV*QG+GZ(I)
+          TG=TG+S*(AH0-AHG)
+          TG=MAX(TG,35.0)
+          TC=TG-273.15
+          DENOM=243.5+TC
+          IF(TC.GE.0.0)THEN  
+           ES=6.112*EXP(17.67*TC/DENOM)
+          ELSE  
+           ES=EXP(23.33086-6111.72784/TG+0.15215*LOG(TG))
+          END IF  
+          QG=EPS*ES/(P(I)-ES*(1.-EPS))
+  200    CONTINUE
+         ALV=LV0-CPVMCL*(T(I)-273.15)
+         TPK(I)=(AH0-(CL-CPD)*Q(NK)*T(I)-GZ(I)-ALV*QG)/CPD
+         CLW(I)=Q(NK)-QG
+         CLW(I)=MAX(0.0,CLW(I))
+         RG=QG/(1.-Q(NK))
+         TVP(I)=TPK(I)*(1.+RG*EPSI)
+  300   CONTINUE
+
+c -- sb:
+        RG = GRAVITY  ! RG redevient la gravite de YOMCST (sb)
+c sb --
+
+        RETURN
+        END
+
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/tlift.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/tlift.F	(revision 222)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/tlift.F	(revision 222)
@@ -0,0 +1,243 @@
+c $Header$
+        SUBROUTINE TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK,
+     .                  TVP,TPK,CLW,ND,NL,
+     .                  DTVPDT1,DTVPDQ1)
+C
+C     Argument NK ajoute (jyg) = Niveau de depart de la
+C     convection
+C
+        PARAMETER (NA=60)
+        REAL GZ(ND),TPK(ND),CLW(ND)
+        REAL T(ND),RR(ND),RS(ND),TVP(ND),P(ND)
+        REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
+C                                   temperature wrt T1 and Q1
+C
+        REAL CLW_NEW(NA),QI(NA)
+        REAL DTPDT1(NA),DTPDQ1(NA)      ! Derivatives of parcel temperature
+C                                   wrt T1 and Q1
+ 
+C
+        LOGICAL ICE_CONV
+C
+C   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
+C
+c sb        CPD=1005.7
+c sb      CPV=1870.0
+c sb        CL=4190.0
+c sb        CPVMCL=2320.0
+c sb        RV=461.5
+c sb        RD=287.04
+c sb        EPS=RD/RV
+c sb        ALV0=2.501E6
+ccccccccccccccccccccccc
+c constantes coherentes avec le modele du Centre Europeen
+c sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
+c sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
+c sb      CPD = 3.5 * RD
+c sb      CPV = 4.0 * RV
+c sb      CL = 4218.0
+c sb      CI=2090.0
+c sb      CPVMCL=CL-CPV
+c sb      CLMCI=CL-CI
+c sb      EPS=RD/RV
+c sb      ALV0=2.5008E+06
+c sb      ALF0=3.34E+05
+ 
+cccccccccccc
+c on utilise les constantes thermo du Centre Europeen: (SB)
+c
+#include "YOMCST.h"
+       GRAVITY = RG !sb: Pr que gravite ne devienne pas humidite!
+c
+       CPD = RCPD
+       CPV = RCPV
+       CL = RCW
+       CI = RCS
+       CPVMCL = CL-CPV
+       CLMCI = CL-CI
+       EPS = RD/RV
+       ALV0 = RLVTT
+       ALF0 = RLMLT ! (ALF0 = RLSTT-RLVTT)
+c 
+cccccccccccccccccccccc
+C
+C   ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
+C
+        ICB1=MAX(ICB,2)
+        ICB1=MIN(ICB,NL)
+C
+Cjyg1
+CC      CPP=CPD*(1.-RR(1))+RR(1)*CPV
+      CPP=CPD*(1.-RR(NK))+RR(NK)*CPV
+Cjyg2
+      CPINV=1./CPP
+Cjyg1
+C         ICB may be below condensation level
+CCC        DO 100 I=1,ICB1-1
+CCC         TPK(I)=T(1)-GZ(I)*CPINV
+CCC         TVP(I)=TPK(I)*(1.+RR(1)/EPS)
+        DO 50 I=1,ICB1
+         CLW(I)=0.0
+50      CONTINUE
+C
+        DO 100 I=NK,ICB1
+         TPK(I)=T(NK)-(GZ(I) - GZ(NK))*CPINV
+Cjyg1
+CCC         TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
+         TVP(I)=TPK(I)*(1.+RR(NK)/EPS-RR(NK))
+Cjyg2
+         DTVPDT1(I) = 1.+RR(NK)/EPS-RR(NK)
+         DTVPDQ1(I) = TPK(I)*(1./EPS-1.)
+C
+Cjyg2
+ 
+  100   CONTINUE
+ 
+C
+C    ***  FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO    ***
+C
+Cjyg1
+CC        AH0=(CPD*(1.-RR(1))+CL*RR(1))*T(1)
+CC     $     +RR(1)*(ALV0-CPVMCL*(T(1)-273.15))
+        AH0=(CPD*(1.-RR(NK))+CL*RR(NK))*T(NK)
+     $     +RR(NK)*(ALV0-CPVMCL*(T(NK)-273.15)) + GZ(NK)
+Cjyg2
+C
+Cjyg1
+        IMIN = ICB1
+C         If ICB is below LCL, start loop at ICB+1
+        IF (PLCL .LT. P(ICB1)) IMIN = MIN(IMIN+1,NL)
+C
+CCC        DO 300 I=ICB1,NL
+        DO 300 I=IMIN,NL
+Cjyg2
+         ALV=ALV0-CPVMCL*(T(I)-273.15)
+         ALF=ALF0+CLMCI*(T(I)-273.15)
+ 
+        RG=RS(I)
+        TG=T(I)
+C       S=CPD+ALV*ALV*RG/(RV*T(I)*T(I))
+Cjyg1
+CC        S=CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*RG/(RV*T(I)*T(I))
+        S=CPD*(1.-RR(NK))+CL*RR(NK)+ALV*ALV*RG/(RV*T(I)*T(I))
+Cjyg2
+        S=1./S
+ 
+        DO 200 J=1,2
+Cjyg1
+CC         AHG=CPD*TG+(CL-CPD)*RR(1)*TG+ALV*RG+GZ(I)
+         AHG=CPD*TG+(CL-CPD)*RR(NK)*TG+ALV*RG+GZ(I)
+Cjyg2
+        TG=TG+S*(AH0-AHG)
+        TC=TG-273.15
+        DENOM=243.5+TC
+        DENOM=MAX(DENOM,1.0)
+C
+C       FORMULE DE BOLTON POUR PSAT
+C
+        ES=6.112*EXP(17.67*TC/DENOM)
+        RG=EPS*ES/(P(I)-ES*(1.-EPS))
+ 
+ 
+  200   CONTINUE
+ 
+Cjyg1
+CC        TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(1))
+        TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(NK))
+Cjyg2
+C       TPK(I)=(AH0-GZ(I)-ALV*RG-(CL-CPD)*T(I)*RR(1))/CPD
+ 
+Cjyg1
+CC        CLW(I)=RR(1)-RG
+        CLW(I)=RR(NK)-RG
+Cjyg2
+        CLW(I)=MAX(0.0,CLW(I))
+Cjyg1
+CCC        TVP(I)=TPK(I)*(1.+RG/EPS)
+        TVP(I)=TPK(I)*(1.+RG/EPS-RR(NK))
+Cjyg2
+C
+Cjyg1       Derivatives
+C
+        DTPDT1(I) = CPD*S
+        DTPDQ1(I) = ALV*S
+C
+         DTVPDT1(I) = DTPDT1(I)*(1. + RG/EPS -
+     .           RR(NK) + ALV*RG/(RD*TPK(I)) )
+        DTVPDQ1(I) = DTPDQ1(I)*(1. + RG/EPS -
+     .           RR(NK) + ALV*RG/(RD*TPK(I)) ) - TPK(I)
+C
+Cjyg2
+ 
+  300   CONTINUE
+C
+      ICE_CONV = .FALSE.
+
+      IF (ICE_CONV) THEN
+C
+CJAM
+C       RAJOUT DE LA PROCEDURE ICEFRAC
+C
+c sb        CALL ICEFRAC(T,CLW,CLW_NEW,QI,ND,NL)
+ 
+        DO 400 I=ICB1,NL
+        IF (T(I).LT.263.15) THEN
+        TG=TPK(I)
+        TC=TPK(I)-273.15
+        DENOM=243.5+TC
+        ES=6.112*EXP(17.67*TC/DENOM)
+        ALV=ALV0-CPVMCL*(T(I)-273.15)
+        ALF=ALF0+CLMCI*(T(I)-273.15)
+ 
+        DO J=1,4
+        ESI=EXP(23.33086-(6111.72784/TPK(I))+0.15215*LOG(TPK(I)))
+        QSAT_NEW=EPS*ESI/(P(I)-ESI*(1.-EPS))
+CCC        SNEW= CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
+        SNEW= CPD*(1.-RR(NK))+CL*RR(NK)
+     .        +ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
+C
+        SNEW=1./SNEW
+        TPK(I)=TG+(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
+c$$$        PRINT*,'################################'
+c$$$        PRINT*,TPK(I)
+c$$$        PRINT*,(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
+        ENDDO
+CCC        CLW(I)=RR(1)-QSAT_NEW
+        CLW(I)=RR(NK)-QSAT_NEW
+        CLW(I)=MAX(0.0,CLW(I))
+Cjyg1
+CCC        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
+        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS-RR(NK))
+Cjyg2
+        ELSE
+        CONTINUE
+        ENDIF
+ 
+  400   CONTINUE
+C
+      ENDIF
+C
+ 
+******************************************************
+** BK :  RAJOUT DE LA TEMPERATURE DES ASCENDANCES
+**   NON DILUES AU  NIVEAU KLEV = ND
+**   POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
+********************************************************
+ 
+      TPK(NL+1)=TPK(NL)
+ 
+*******************************************************
+
+      RG = GRAVITY  ! RG redevient la gravite de YOMCST (sb)
+ 
+ 
+        RETURN
+        END
+ 
+ 
+ 
+ 
+ 
+ 
+ 
+ 
