Index: /LMDZ4/trunk/libf/dyn3d/bilan_dyn.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/bilan_dyn.F	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/bilan_dyn.F	(revision 541)
@@ -2,7 +2,4 @@
 ! $Header$
 !
-c
-c $Header$
-c
       SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
      s  ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
Index: /LMDZ4/trunk/libf/dyn3d/conf_gcm.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/conf_gcm.F	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/conf_gcm.F	(revision 541)
@@ -574,4 +574,12 @@
       ENDIF
 c
+!Config  Key  = offline
+!Config  Desc = Nouvelle eau liquide
+!Config  Def  = n
+!Config  Help = Permet de mettre en route la
+!Config         nouvelle parametrisation de l'eau liquide !
+       offline = .FALSE.
+       CALL getin('offline',offline)
+
 
       write(lunout,*)' #########################################'
@@ -615,4 +623,5 @@
       write(lunout,*)' tauxx = ', tauxx
       write(lunout,*)' tauyy = ', tauyy
+      write(lunout,*)' offline = ', offline
 
       RETURN
@@ -720,4 +729,13 @@
        CALL getin('ysinus',ysinus)
 c
+!Config  Key  = offline
+!Config  Desc = Nouvelle eau liquide
+!Config  Def  = n
+!Config  Help = Permet de mettre en route la
+!Config         nouvelle parametrisation de l'eau liquide !
+       offline = .FALSE.
+       CALL getin('offline',offline)
+      write(lunout,*)' offline = ', offline
+
 
       write(lunout,*)' #########################################'
@@ -761,4 +779,5 @@
       write(lunout,*)' taux = ', taux
       write(lunout,*)' tauy = ', tauy
+      write(lunout,*)' offline = ', offline
 c
       RETURN
Index: /LMDZ4/trunk/libf/dyn3d/control.h
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/control.h	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/control.h	(revision 541)
@@ -6,11 +6,12 @@
 
       COMMON/control/nday,day_step,
-     $              iperiod,iapp_tracvl,iconser,iecri,idissip,iphysiq ,
-     $              periodav,ecritphy,iecrimoy,dayref,anneeref,
-     $              raz_date
+     .              iperiod,iapp_tracvl,iconser,iecri,idissip,iphysiq ,
+     .              periodav,ecritphy,iecrimoy,dayref,anneeref,
+     .              raz_date,offline
 
       INTEGER   nday,day_step,iperiod,iapp_tracvl,iconser,iecri,
-     $          idissip,iphysiq,iecrimoy,dayref,anneeref, raz_date
+     .          idissip,iphysiq,iecrimoy,dayref,anneeref, raz_date
       REAL periodav, ecritphy
+      logical offline
 
 c-----------------------------------------------------------------------
Index: /LMDZ4/trunk/libf/dyn3d/dynetat0.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/dynetat0.F	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/dynetat0.F	(revision 541)
@@ -106,8 +106,8 @@
       IF ( tab_cntrl(24).EQ.1. )  THEN
         fxyhypb  = . TRUE .
-        dzoomx   = tab_cntrl(25)
-        dzoomy   = tab_cntrl(26)
-        taux     = tab_cntrl(28)
-        tauy     = tab_cntrl(29)
+c        dzoomx   = tab_cntrl(25)
+c        dzoomy   = tab_cntrl(26)
+c        taux     = tab_cntrl(28)
+c        tauy     = tab_cntrl(29)
       ELSE
         fxyhypb = . FALSE .
Index: /LMDZ4/trunk/libf/dyn3d/dynredem.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/dynredem.F	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/dynredem.F	(revision 541)
@@ -394,4 +394,5 @@
 #include "comgeom.h"
 #include "advtrac.h"
+#include "temps.h"
 
       INTEGER nq, l
@@ -406,4 +407,7 @@
       INTEGER ierr
       INTEGER iq
+      INTEGER length
+      PARAMETER (length = 100)
+      REAL tab_cntrl(length) ! tableau des parametres du run
       character*20 modname
       character*80 abort_message
@@ -435,4 +439,25 @@
 #endif
       PRINT*, "Enregistrement pour ", nb, time
+
+c
+c  Re-ecriture du tableau de controle, itaufin n'est plus defini quand
+c  on passe dans dynredem0
+      ierr = NF_INQ_VARID (nid, "controle", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         abort_message="dynredem1: Le champ <controle> est absent"
+         ierr = 1
+         CALL abort_gcm(modname,abort_message,ierr)
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
+#endif
+       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
 
 c  Ecriture des champs
Index: /LMDZ4/trunk/libf/dyn3d/fluxstokenc.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/fluxstokenc.F	(revision 541)
+++ /LMDZ4/trunk/libf/dyn3d/fluxstokenc.F	(revision 541)
@@ -0,0 +1,155 @@
+      SUBROUTINE fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
+     . time_step,itau )
+
+       USE IOIPSL
+c
+c     Auteur :  F. Hourdin
+c
+c
+ccc   ..   Modif. P. Le Van  ( 20/12/97 )  ...
+c
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "tracstoke.h"
+#include "temps.h"
+
+      REAL time_step,t_wrt, t_ops
+      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+      REAL masse(ip1jmp1,llm),teta(ip1jmp1,llm),phi(ip1jmp1,llm)
+      REAL phis(ip1jmp1)
+
+      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
+      REAL massem(ip1jmp1,llm),tetac(ip1jmp1,llm),phic(ip1jmp1,llm)
+
+      REAL pbarug(ip1jmp1,llm),pbarvg(iip1,jjm,llm),wg(ip1jmp1,llm)
+
+      REAL pbarvst(iip1,jjp1,llm),zistdyn
+	real dtcum
+
+      INTEGER iadvtr,ndex(1) 
+      integer nscal
+      real tst(1),ist(1),istp(1)
+      INTEGER ij,l,irec,i,j,itau
+      INTEGER fluxid, fluxvid,fluxdid
+ 
+      SAVE iadvtr, massem,pbaruc,pbarvc,irec
+      SAVE phic,tetac
+      logical first
+      save first
+      data first/.true./
+      DATA iadvtr/0/
+
+      if(first) then
+
+	CALL initfluxsto( 'fluxstoke',
+     .  time_step,istdyn* time_step,istdyn* time_step,
+     . nqmx, fluxid,fluxvid,fluxdid) 
+	
+	ndex(1) = 0
+        call histwrite(fluxid, 'phis', 1, phis, iip1*jjp1, ndex)
+        call histwrite(fluxid, 'aire', 1, aire, iip1*jjp1, ndex)
+	
+	ndex(1) = 0
+        nscal = 1
+        tst(1) = time_step
+        call histwrite(fluxdid, 'dtvr', 1, tst, nscal, ndex)
+        ist(1)=istdyn
+        call histwrite(fluxdid, 'istdyn', 1, ist, nscal, ndex)
+        istp(1)= istphy
+        call histwrite(fluxdid, 'istphy', 1, istp, nscal, ndex)
+	
+	first = .false.
+
+      endif
+
+
+      IF(iadvtr.EQ.0) THEN
+         CALL initial0(ijp1llm,phic)
+         CALL initial0(ijp1llm,tetac)
+         CALL initial0(ijp1llm,pbaruc)
+         CALL initial0(ijmllm,pbarvc)
+      ENDIF
+
+c   accumulation des flux de masse horizontaux
+      DO l=1,llm
+         DO ij = 1,ip1jmp1
+            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
+            tetac(ij,l) = tetac(ij,l) + teta(ij,l)
+            phic(ij,l) = phic(ij,l) + phi(ij,l)
+         ENDDO
+         DO ij = 1,ip1jm
+            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
+         ENDDO
+      ENDDO
+
+c   selection de la masse instantannee des mailles avant le transport.
+      IF(iadvtr.EQ.0) THEN
+         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
+      ENDIF
+
+      iadvtr   = iadvtr+1
+
+
+c   Test pour savoir si on advecte a ce pas de temps
+      IF ( iadvtr.EQ.istdyn ) THEN
+c    normalisation
+      DO l=1,llm
+         DO ij = 1,ip1jmp1
+            pbaruc(ij,l) = pbaruc(ij,l)/float(istdyn)
+            tetac(ij,l) = tetac(ij,l)/float(istdyn)
+            phic(ij,l) = phic(ij,l)/float(istdyn)
+         ENDDO
+         DO ij = 1,ip1jm
+            pbarvc(ij,l) = pbarvc(ij,l)/float(istdyn)
+         ENDDO
+      ENDDO
+
+c   traitement des flux de masse avant advection.
+c     1. calcul de w
+c     2. groupement des mailles pres du pole.
+
+        CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
+
+        do l=1,llm
+           do j=1,jjm
+              do i=1,iip1
+                 pbarvst(i,j,l)=pbarvg(i,j,l)
+              enddo
+           enddo
+           do i=1,iip1
+              pbarvst(i,jjp1,l)=0.
+           enddo
+        enddo
+
+         iadvtr=0
+	Print*,'ITAU auqel on stoke les fluxmasses',itau
+	
+	call histwrite(fluxid, 'masse', itau, massem,
+     .               iip1*jjp1*llm, ndex)
+	
+	call histwrite(fluxid, 'pbaru', itau, pbarug,
+     .               iip1*jjp1*llm, ndex)
+	
+	call histwrite(fluxvid, 'pbarv', itau, pbarvg,
+     .               iip1*jjm*llm, ndex)
+	
+        call histwrite(fluxid, 'w' ,itau, wg, 
+     .             iip1*jjp1*llm, ndex) 
+	
+	call histwrite(fluxid, 'teta' ,itau, tetac, 
+     .             iip1*jjp1*llm, ndex) 
+	
+	call histwrite(fluxid, 'phi' ,itau, phic, 
+     .             iip1*jjp1*llm, ndex) 
+	
+C
+
+      ENDIF ! if iadvtr.EQ.istdyn
+
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3d/gcm.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/gcm.F	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/gcm.F	(revision 541)
@@ -56,6 +56,5 @@
 #include "com_io_dyn.h"
 #include "iniprint.h"
-
-c#include "tracstoke.h"
+#include "tracstoke.h"
 
 
@@ -118,7 +117,4 @@
 c-jld 
 
-
-      LOGICAL offline  ! Controle du stockage ds "fluxmass"
-      PARAMETER (offline=.false.)
 
       character*80 dynhist_file, dynhistave_file
@@ -368,4 +364,11 @@
 #endif
 
+c  Choix des frequences de stokage pour le offline
+c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
+c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
+      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
+      istphy=istdyn/iphysiq     
+
+
 c
 c-----------------------------------------------------------------------
Index: /LMDZ4/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/leapfrog.F	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/leapfrog.F	(revision 541)
@@ -149,7 +149,4 @@
 c-jld 
 
-      LOGICAL offline  ! Controle du stockage ds "fluxmass"
-      PARAMETER (offline=.false.)
-
       character*80 dynhist_file, dynhistave_file
       character*20 modname
@@ -199,6 +196,8 @@
 
 #ifdef CPP_IOIPSL
-      if (ok_guide) then
+      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
         call guide(itau,ucov,vcov,teta,q,masse,ps)
+      else
+        print*,'attention on ne guide pas les 6 dernieres heures'
       endif
 #endif
@@ -289,6 +288,6 @@
 
 #ifdef CPP_IOIPSL
-c           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
-c    .   time_step, itau)
+           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
+     .   dtvr, itau)
 #endif
 
@@ -486,11 +485,11 @@
 
             IF( itau. EQ. itaufinp1 ) then  
-c      write(79,*) 'ucov',ucov
-c      write(80,*) 'vcov',vcov
-c      write(81,*) 'teta',teta
-c      write(82,*) 'ps',ps
-c      write(83,*) 'q',q
-c      WRITE(85,*) 'q1 = ',q(:,:,1)
-c      WRITE(86,*) 'q3 = ',q(:,:,3)
+       write(79,*) 'ucov',ucov
+       write(80,*) 'vcov',vcov
+       write(81,*) 'teta',teta
+       write(82,*) 'ps',ps
+       write(83,*) 'q',q
+       WRITE(85,*) 'q1 = ',q(:,:,1)
+       WRITE(86,*) 'q3 = ',q(:,:,3)
 
               abort_message = 'Simulation finished'
Index: /LMDZ4/trunk/libf/dyn3d/wrgrads.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3d/wrgrads.F	(revision 540)
+++ /LMDZ4/trunk/libf/dyn3d/wrgrads.F	(revision 541)
@@ -29,5 +29,5 @@
       writectl=.false.
 
-c     print*,if,iid(if),jid(if),ifd(if),jfd(if)
+      print*,if,iid(if),jid(if),ifd(if),jfd(if)
       iii=iid(if)
       iji=jid(if)
@@ -38,6 +38,6 @@
       lm=lmd(if)
 
-c      print*,'im,jm,lm,name,firsttime(if)'
-c      print*,im,jm,lm,name,firsttime(if)
+      print*,'im,jm,lm,name,firsttime(if)'
+      print*,im,jm,lm,name,firsttime(if)
 
       if(firsttime(if)) then
@@ -81,6 +81,6 @@
       endif
 
-c     print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
-c     print*,ivar(if),nvar(if),var(ivar(if),if),writectl
+      print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
+      print*,ivar(if),nvar(if),var(ivar(if),if),writectl
       do l=1,nl
          irec(if)=irec(if)+1
Index: /LMDZ4/trunk/libf/phylmd/FCTTRE.inc
===================================================================
--- /LMDZ4/trunk/libf/phylmd/FCTTRE.inc	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/FCTTRE.inc	(revision 541)
@@ -16,5 +16,6 @@
 !
       FOEEW ( PTARG,PDELARG ) = EXP ( &
-     &          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
+!    &          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-TEST_RTT) &
+     &          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-273.16) &
      & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
 !
@@ -31,8 +32,10 @@
      &           + 25.21935 * EXP( - 2999.924 / ptarg))
 !
-      dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg &
+!      dqsats(ptarg,pqsarg) = TEST_RLVTT/TEST_RCPD*pqsarg * (3.56654/ptarg &
+      dqsats(ptarg,pqsarg) = 0.2500800E+07/0.1004709E+04*pqsarg * (3.56654/ptarg &
      &                     +2484.896*LOG(10.)/ptarg**2 &
      &                     -0.00320991*LOG(10.))
-      dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* &
+!     dqsatl(ptarg,pqsarg) = TEST_RLVTT/TEST_RCPD*pqsarg*LOG(10.)* &
+      dqsatl(ptarg,pqsarg) = 0.2500800E+07/0.1004709E+04*pqsarg*LOG(10.)* &
      &                (2948.964/ptarg**2-5.028/LOG(10.)/ptarg &
      &                +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) &
Index: /LMDZ4/trunk/libf/phylmd/calltherm.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/calltherm.F	(revision 541)
+++ /LMDZ4/trunk/libf/phylmd/calltherm.F	(revision 541)
@@ -0,0 +1,116 @@
+!
+! $Header$
+!
+      subroutine calltherm(dtime
+     s      ,pplay,paprs,pphi
+     s      ,u_seri,v_seri,t_seri,q_seri
+     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
+     s      ,fm_therm,entr_therm)
+
+      implicit none
+#include "dimensions.h"
+#include "dimphy.h"
+#include "thermcell.h"
+
+c  A inclure eventuellement dans les fichiers de configuration
+      data r_aspect_thermals,l_mix_thermals,tho_thermals/4.,10.,0./
+      data w2di_thermals/0/
+
+      REAL dtime
+
+      REAL u_seri(klon,klev),v_seri(klon,klev)
+      REAL t_seri(klon,klev),q_seri(klon,klev)
+      REAL paprs(klon,klev+1)
+      REAL pplay(klon,klev)
+      REAL pphi(klon,klev)
+
+CFH Update Thermiques
+      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
+      REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
+      real fm_therm(klon,klev+1),entr_therm(klon,klev)
+
+
+c variables locales
+      REAL d_t_the(klon,klev), d_q_the(klon,klev)
+      REAL d_u_the(klon,klev),d_v_the(klon,klev)
+c
+      real zfm_therm(klon,klev+1),zentr_therm(klon,klev),zdt
+      save zentr_therm,zfm_therm
+
+      integer i,k
+
+*********************************************************
+
+c  Modele du thermique
+c  ===================
+c         print*,'thermiques: WARNING on passe t au lieu de t_seri'
+       print*,'avant isplit ',nsplit_thermals
+
+
+         fm_therm(:,:)=0.
+         entr_therm(:,:)=0.
+
+c   tests sur les valeurs negatives de l'eau
+         do k=1,klev
+            do i=1,klon
+               if (.not.q_seri(i,k).ge.0.) then
+                   print*,'WARN eau<0 avant therm i=',i,'  k=',k
+     s         ,' dq,q',d_q_the(i,k),q_seri(i,k)
+                  q_seri(i,k)=1.e-15
+               endif
+            enddo
+         enddo
+
+
+         zdt=dtime/float(nsplit_thermals)
+         do isplit=1,nsplit_thermals
+
+            CALL thermcell(klon,klev,zdt
+     s      ,pplay,paprs,pphi
+     s      ,u_seri,v_seri,t_seri,q_seri
+     s      ,d_u_the,d_v_the,d_t_the,d_q_the
+     s      ,zfm_therm,zentr_therm
+     s      ,r_aspect_thermals,l_mix_thermals,w2di_thermals
+     s      ,tho_thermals,3)
+
+c  transformation de la derivee en tendance
+            d_t_the(:,:)=d_t_the(:,:)*dtime/float(nsplit_thermals)
+            d_u_the(:,:)=d_u_the(:,:)*dtime/float(nsplit_thermals)
+            d_v_the(:,:)=d_v_the(:,:)*dtime/float(nsplit_thermals)
+            d_q_the(:,:)=d_q_the(:,:)*dtime/float(nsplit_thermals)
+            fm_therm(:,:)=fm_therm(:,:)
+     s      +zfm_therm(:,:)/float(nsplit_thermals)
+            entr_therm(:,:)=entr_therm(:,:)
+     s       +zentr_therm(:,:)/float(nsplit_thermals)
+            fm_therm(:,klev+1)=0.
+
+
+
+c  accumulation de la tendance
+            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
+            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
+            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
+            d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
+
+c  incrementation des variables meteo
+            t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
+            u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
+            v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
+            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
+
+c   tests sur les valeurs negatives de l'eau
+            DO k = 1, klev
+            DO i = 1, klon
+               if (.not.q_seri(i,k).ge.0.) then
+                   print*,'WARN eau<0 apres therm i=',i,'  k=',k
+     s         ,' dq,q',d_q_the(i,k),q_seri(i,k)
+                  q_seri(i,k)=1.e-15
+               endif
+            ENDDO
+            ENDDO
+
+         enddo ! isplit
+
+      return
+
+      end
Index: /LMDZ4/trunk/libf/phylmd/clmain.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/clmain.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/clmain.F	(revision 541)
@@ -17,6 +17,6 @@
      .                  d_t,d_q,d_u,d_v,d_ts,
      .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,
+     .                  q2,
      .                  dflux_t,dflux_q,
-cIM cf JLD    .                  zcoefh,zu1,zv1, t2m, q2m, u10m, v10m)
      .                  zcoefh,zu1,zv1, t2m, q2m, u10m, v10m,
      .                  fqcalving,ffonte, run_off_lic_0)
@@ -89,4 +89,5 @@
 #include "dimsoil.h"
 #include "iniprint.h"
+#include "compbl.h"
 c
       REAL dtime
@@ -183,10 +184,20 @@
 cAA      REAL yflxsrf(klon,nbtr)
 c
-      LOGICAL contreg
-      PARAMETER (contreg=.TRUE.)
 c
       LOGICAL ok_nonloc
       PARAMETER (ok_nonloc=.FALSE.)
       REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
+
+      real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)
+      real ykmm(klon,klev+1),ykmn(klon,klev+1)
+      real ykmq(klon,klev+1)
+      real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)
+      real q2diag(klon,klev+1)
+      real yustar(klon),y_cd_m(klon),y_cd_h(klon)
+
+
+
+
+
 c
 #include "YOMCST.h"
@@ -220,5 +231,5 @@
       LOGICAL first_appel
       SAVE first_appel
-      DATA first_appel/.false./
+      DATA first_appel/.true./
       LOGICAL debugindex
       SAVE debugindex
@@ -254,5 +265,5 @@
       endif
       IF (first_appel) THEN
-          first_appel=.false.
+!          first_appel=.false.
 !
 ! initialisation sorties netcdf
@@ -479,5 +490,12 @@
 c
 c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c calculer Cdrag et les coefficients d'echange
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c  Calcul anciens du LMD. Effectues de toutes facons.
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
       CALL coefkz(nsrf, knon, ypaprs, ypplay,
 cIM 261103
@@ -487,6 +505,25 @@
      .            yqsurf, 
      .            ycoefm, ycoefh)
+
       CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
      .                  ycoefm0, ycoefh0)
+      if (first_appel) then
+        if (prt_level > 9) THEN
+          WRITE(lunout,*)'Apres coefkz2 '
+          WRITE(lunout,*)'nsrf,knon,yts,yrugos,yqsurf',
+     .    nsrf,knon,yts,yrugos,yqsurf
+          WRITE(lunout,*)'ypaprs(1,k),ypplay(1,k),yu,yv,yt'
+          do k=1,klev
+            WRITE(lunout,*)ypaprs(1,k),ypplay(1,k),
+     .      yu(1,k),yv(1,k),yt(1,k)
+          enddo
+          do k=1,klev
+            WRITE(lunout,*)ycoefm(1,k),ycoefh(1,k),
+     .      ycoefm0(1,k),ycoefh0(1,k)
+          enddo
+        ENDIF
+        first_appel=.false.
+      endif
+
       DO k = 1, klev
       DO i = 1, knon
@@ -495,4 +532,5 @@
       ENDDO
       ENDDO
+
 c
 cIM cf JLD : on seuille ycoefm et ycoefh
@@ -528,6 +566,78 @@
 cIM: 261103
 
-c
+
+      IF (iflag_pbl.ge.3) then
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+         yzlay(1:knon,1)=
+     .   RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))
+     .   *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
+         do k=2,klev
+            yzlay(1:knon,k)=
+     .      yzlay(1:knon,k-1)+RD*0.5*(yt(1:knon,k-1)+yt(1:knon,k))
+     .      /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG
+         enddo
+         do k=1,klev
+            yteta(1:knon,k)=
+     .      yt(1:knon,k)*(ypaprs(1:knon,1)/ypplay(1:knon,k))**rkappa
+     .      *(1.+0.61*yq(1:knon,k))
+         enddo
+         yzlev(1:knon,1)=0.
+         yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
+         do k=2,klev
+            yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))
+         enddo
+         DO k = 1, klev+1
+            DO j = 1, knon
+               i = ni(j)
+               yq2(j,k)=q2(i,k,nsrf)
+            enddo
+         enddo
+
+
+c   Bug introduit volontairement pour converger avec les resultats
+c  du papier sur les thermiques.
+         if (1.eq.1) then
+         y_cd_m(1:knon) = ycoefm(1:knon,1)
+         y_cd_h(1:knon) = ycoefh(1:knon,1)
+         else
+         y_cd_h(1:knon) = ycoefm(1:knon,1)
+         y_cd_m(1:knon) = ycoefh(1:knon,1)
+         endif
+         call ustarhb(knon,yu,yv,y_cd_m, yustar)
+
+        if (prt_level > 9) THEN
+          WRITE(lunout,*)'USTAR = ',yustar
+        ENDIF
+
+c   iflag_pbl peut etre utilise comme longuer de melange
+
+         if (iflag_pbl.ge.11) then
+            call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt
+     s      ,yzlev,yzlay,yu,yv,yteta
+     s      ,y_cd_m,yq2,q2diag,ykmm,ykmn,yustar,
+     s      iflag_pbl)
+         else
+            call yamada4(knon,dtime,rg,rd,ypaprs,yt
+     s      ,yzlev,yzlay,yu,yv,yteta
+     s      ,y_cd_m,yq2,ykmm,ykmn,ykmq,yustar,
+     s      iflag_pbl)
+         endif
+
+         ycoefm(1:knon,1)=y_cd_m(1:knon)
+         ycoefh(1:knon,1)=y_cd_h(1:knon)
+         ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
+         ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
+
+
+      ENDIF
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c calculer la diffusion des vitesses "u" et "v"
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
      s            y_d_u,y_flux_u)
@@ -545,5 +655,7 @@
 c$$$      enddo
 
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c calculer la diffusion de "q" et de "h"
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       CALL clqh(dtime, itap, date0,jour, debut,lafin,
      e          rlon, rlat, cufi, cvfi,
@@ -560,5 +672,4 @@
      s          pctsrf_new, yagesno,
      s          y_d_t, y_d_q, y_d_ts, yz0_new,
-cIM cf JLD    s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q)
      s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,
      s          y_fqcalving,y_ffonte,y_run_off_lic_0)
@@ -699,23 +810,4 @@
 c
         qairsol(j) = yqsurf(j)
-c$$$        IF (nsrf.EQ.1) THEN
-c$$$          qairsol(j) = yqsurf(j)
-c$$$        ELSE IF(nsrf.GT.1) THEN
-c$$$         zt = ts(i,nsrf)
-c$$$         IF (thermcep) THEN
-c$$$           zdelta = MAX(0.,SIGN(1.,RTT-zt))
-c$$$           zqs = R2ES * FOEEW(zt,zdelta) / ypplay(j,1)
-c$$$           zqs = MIN(0.5,zqs)
-c$$$           zcor = 1./(1.-RETV*zqs)
-c$$$           zqs = zqs*zcor
-c$$$         ELSE
-c$$$           IF (zt .LT. t_coup) THEN
-c$$$             zqs = qsats(zt) / ypplay(j,1)
-c$$$           ELSE
-c$$$             zqs = qsatl(zt) / ypplay(j,1)
-c$$$           ENDIF
-c$$$         ENDIF   
-c$$$         qairsol(j) = zqs
-c$$$        ENDIF   
       ENDDO
 c
@@ -778,4 +870,13 @@
        ENDDO 
 #endif
+
+      do j=1,knon
+         do k=1,klev+1
+         i=ni(j)
+         q2(i,k,nsrf)=yq2(j,k)
+         enddo
+      enddo
+c
+
 99999 CONTINUE
 c
@@ -802,5 +903,4 @@
      s                pctsrf_new, agesno,
      s                d_t, d_q, d_ts, z0_new, 
-cIM cf JLD    s                flux_t, flux_q,dflux_s,dflux_l)
      s                flux_t, flux_q,dflux_s,dflux_l,
      s                fqcalving,ffonte,run_off_lic_0)
@@ -820,4 +920,5 @@
 #include "indicesol.h"
 #include "dimsoil.h"
+#include "iniprint.h"
 
 c Arguments:
@@ -901,6 +1002,5 @@
       REAL zdelz
 c======================================================================
-      logical contreg
-      parameter (contreg=.true.)
+#include "compbl.h"
 c======================================================================
 c Rajout pour l'interface
@@ -944,12 +1044,6 @@
       ENDIF
 C
-      if (.not. contreg) then
-        do k = 2, klev
-          do i = 1, knon
-            gamq(i,k) = 0.0
-            gamt(i,k) = 0.0
-          enddo
-        enddo
-      else
+C
+      if (iflag_pbl.eq.1) then
         do k = 3, klev
           do i = 1, knon
@@ -961,4 +1055,11 @@
           gamq(i,2) = 0.0
           gamt(i,2) = -2.5e-03
+        enddo
+      else
+        do k = 2, klev
+          do i = 1, knon
+            gamq(i,k) = 0.0
+            gamt(i,k) = 0.0
+          enddo
         enddo
       endif
@@ -1089,5 +1190,4 @@
      s evap, fluxsens, fluxlat, dflux_l, dflux_s,              
      s tsol_rad, tsurf_new, alb_new, alblw, emis_new, z0_new, 
-cIM cf JLD    s pctsrf_new, agesno)
      s pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0)
 
@@ -1162,4 +1262,5 @@
 #include "dimensions.h"
 #include "dimphy.h"
+#include "iniprint.h"
       INTEGER knon
       REAL dtime
@@ -1288,4 +1389,5 @@
 #include "YOMCST.h"
 #include "indicesol.h"
+#include "iniprint.h"
 c
 c Arguments:
@@ -1335,6 +1437,6 @@
       LOGICAL opt_ec ! formule du Centre Europeen dans l'atmosphere
       PARAMETER (opt_ec=.FALSE.)
-      LOGICAL contreg ! utiliser le contre-gradient dans Ri
-      PARAMETER (contreg=.TRUE.)
+
+#include "compbl.h"
 c
 c Variables locales:
@@ -1379,10 +1481,12 @@
 c
       IF (appel1er) THEN
-         PRINT*, 'coefkz, opt_ec:', opt_ec
-         PRINT*, 'coefkz, richum:', richum
-         IF (richum) PRINT*, 'coefkz, ratqs:', ratqs
-         PRINT*, 'coefkz, isommet:', isommet
-         PRINT*, 'coefkz, tvirtu:', tvirtu
-         appel1er = .FALSE.
+        if (prt_level > 9) THEN
+          WRITE(lunout,*)'coefkz, opt_ec:', opt_ec
+          WRITE(lunout,*)'coefkz, richum:', richum
+          IF (richum) WRITE(lunout,*)'coefkz, ratqs:', ratqs
+          WRITE(lunout,*)'coefkz, isommet:', isommet
+          WRITE(lunout,*)'coefkz, tvirtu:', tvirtu
+          appel1er = .FALSE.
+        endif
       ENDIF
 c
@@ -1399,22 +1503,16 @@
       ENDDO
 
-c$$$      IF(nsrf.NE.1) THEN
-c$$$      do i = 1, knon
-c$$$        qsurf(i) = qsatl(ts(i))/paprs(i,1)
-c$$$      enddo
-c$$$      ENDIF
-
 c
 c Prescrire la valeur de contre-gradient
 c
-      IF (.NOT.contreg) THEN
-         DO k = 2, klev
-            gamt(k) = 0.0
-         ENDDO
-      ELSE
+      if (iflag_pbl.eq.1) then
          DO k = 3, klev
             gamt(k) = -1.0E-03
          ENDDO
          gamt(2) = -2.5E-03
+      else
+         DO k = 2, klev
+            gamt(k) = 0.0
+         ENDDO
       ENDIF
 cIM cf JLD/ GKtest
@@ -1469,7 +1567,4 @@
       ENDDO
 
-      IF (check) THEN
-      PRINT*,' isommet=',isommet,' knon=',knon
-      ENDIF
 
       DO k = 2, isommet
@@ -1608,4 +1703,5 @@
 #include "YOMCST.h"
 #include "indicesol.h"
+#include "iniprint.h"
 c
 c Arguments:
@@ -1702,4 +1798,5 @@
 #include "YOMCST.h"
 #include "indicesol.h"
+#include "iniprint.h"
       REAL tau_gl ! temps de relaxation pour la glace de mer
 ccc      PARAMETER (tau_gl=86400.0*30.0)
@@ -1783,4 +1880,5 @@
 #include "dimphy.h"
 #include "YOMCST.h"
+#include "iniprint.h"
 c
 c Arguments:
Index: /LMDZ4/trunk/libf/phylmd/conf_phys.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/conf_phys.F90	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/conf_phys.F90	(revision 541)
@@ -9,5 +9,6 @@
  &                     ratqsbas,ratqshaut,if_ebil, &
  &		       ok_ade, ok_aie, &
- &                     bl95_b0, bl95_b1)
+ &                     bl95_b0, bl95_b1,&
+ &                     iflag_thermals,nsplit_thermals)
 
    use IOIPSL
@@ -19,5 +20,6 @@
 #include "YOMCST.inc"
 !IM : on inclut/initialise les taux de CH4, N2O, CFC11 et CFC12
-#include "clesphys.inc"
+include "clesphys.inc"
+include "compbl.h"
 !
 ! Configuration de la "physique" de LMDZ a l'aide de la fonction
@@ -51,4 +53,5 @@
   real                 :: zzz
 
+  integer :: iflag_thermals,nsplit_thermals
 !
 !
@@ -514,4 +517,35 @@
   ok_kzmin = .true.
   call getin('ok_kzmin',ok_kzmin)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! PARAMETER FOR THE PLANETARY BOUNDARY LAYER
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!Config Key  = iflag_pbl
+!Config Desc =
+!Config Def  = 1
+!Config Help =
+!
+  iflag_pbl = 1
+  call getin('iflag_pbl',iflag_pbl)
+!
+!Config Key  = iflag_thermals
+!Config Desc =
+!Config Def  = 0
+!Config Help =
+!
+  iflag_thermals = 0
+  call getin('iflag_thermals',iflag_thermals)
+!
+!
+!Config Key  = nsplit_thermals
+!Config Desc =
+!Config Def  = 1
+!Config Help =
+!
+  nsplit_thermals = 1
+  call getin('nsplit_thermals',nsplit_thermals)
+
+
 
 !
@@ -606,4 +640,6 @@
   write(numout,*)' lev_histday = ',lev_histday 
   write(numout,*)' lev_histmth = ',lev_histmth 
+  write(numout,*)' iflag_pbl = ', iflag_pbl
+  write(numout,*)' iflag_thermals = ', iflag_thermals
 
   return
Index: /LMDZ4/trunk/libf/phylmd/initphysto.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/initphysto.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/initphysto.F	(revision 541)
@@ -188,7 +188,17 @@
 c coefh frac_impa,frac_nucl
 	
-	call histdef(fileid, 'coefh', ' ', ' ',
-     .             iim, jjm+1, nhoriid, llm, 1, llm, zvertiid,
-     .             32, 'inst(X)', t_ops, t_wrt)
+	call histdef(fileid, "coefh", " ", " ",
+     .             iim, jjm+1, nhoriid, llm, 1, llm, zvertiid,
+     .             32, "inst(X)", t_ops, t_wrt)
+
+c abderrahmane le 16 09 02
+        call histdef(fileid, "fm_th", " ", " ",
+     .             iim, jjm+1, nhoriid, llm, 1, llm, zvertiid,
+     .             32, "inst(X)", t_ops, t_wrt)
+
+        call histdef(fileid, "en_th", " ", " ",
+     .             iim, jjm+1, nhoriid, llm, 1, llm, zvertiid,
+     .             32, "inst(X)", t_ops, t_wrt)
+c fin aj
 	
 	write(*,*) 'apres coefh ds initphysto'	
@@ -283,5 +293,6 @@
 
       call histend(fileid)
-      if (ok_sync) call histsync(fileid)
+c     if (ok_sync) call histsync(fileid)
+      if (ok_sync) call histsync
 
 	
Index: /LMDZ4/trunk/libf/phylmd/initrrnpb.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/initrrnpb.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/initrrnpb.F	(revision 541)
@@ -74,6 +74,4 @@
       DO i = 1,klon
         masktr(i,it) = 0.
-c       IF ( NINT(pctsrf(i,3)) .EQ. 1 ) masktr(i,it) = 1.
-c       fshtr(i,it) = s * masktr(i,it) * pctsrf(i,3)
         IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i,it) = 1.
         fshtr(i,it) = s * masktr(i,it)
Index: /LMDZ4/trunk/libf/phylmd/oasis.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/oasis.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/oasis.F	(revision 541)
@@ -2,5 +2,4 @@
 ! $Header$
 !
-C $Id$
 C****
 C
@@ -43,5 +42,6 @@
 #include "mpiclim.h"
 c
-#include "oasis.h"      ! contains the name of communication technique. Here
+#include "oasis.h"      
+                        ! contains the name of communication technique. Here
                         ! cchan=CLIM only is possible.
 c			! ctype=MPI2
@@ -557,4 +557,19 @@
       END
 
+      SUBROUTINE halte
+      print *, 'Attention dans oasis.F, halte est non defini'
+      RETURN
+      END
+
+      SUBROUTINE locread
+      print *, 'Attention dans oasis.F, locread est non defini'
+      RETURN
+      END
+
+      SUBROUTINE locwrite
+      print *, 'Attention dans oasis.F, locwrite est non defini'
+      RETURN
+      END
+
       SUBROUTINE pipe_model_define
       print*,'Attention dans oasis.F, pipe_model_define est non defini'
@@ -577,4 +592,54 @@
       END
 
+      SUBROUTINE clim_stepi
+      print *, 'Attention dans oasis.F, clim_stepi est non defini'
+      RETURN
+      END
+
+      SUBROUTINE clim_start
+      print *, 'Attention dans oasis.F, clim_start est non defini'
+      RETURN
+      END
+
+      SUBROUTINE clim_import
+      print *, 'Attention dans oasis.F, clim_import est non defini'
+      RETURN
+      END
+
+      SUBROUTINE clim_export
+      print *, 'Attention dans oasis.F, clim_export est non defini'
+      RETURN
+      END
+
+      SUBROUTINE clim_init
+      print *, 'Attention dans oasis.F, clim_init est non defini'
+      RETURN
+      END
+
+      SUBROUTINE clim_define
+      print *, 'Attention dans oasis.F, clim_define est non defini'
+      RETURN
+      END
+
+      SUBROUTINE clim_quit
+      print *, 'Attention dans oasis.F, clim_quit est non defini'
+      RETURN
+      END
+
+      SUBROUTINE svipc_write
+      print *, 'Attention dans oasis.F, svipc_write est non defini'
+      RETURN
+      END
+
+      SUBROUTINE svipc_close
+      print *, 'Attention dans oasis.F, svipc_close est non defini'
+      RETURN
+      END
+
+      SUBROUTINE svipc_read
+      print *, 'Attention dans oasis.F, svipc_read est non defini'
+      RETURN
+      END
+
       SUBROUTINE quitcpl
       print *, 'Attention dans oasis.F, quitcpl est non defini'
Index: /LMDZ4/trunk/libf/phylmd/phyetat0.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/phyetat0.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/phyetat0.F	(revision 541)
@@ -665,31 +665,4 @@
            ENDDO
          ENDDO
-c          IF (nsrf.GT.99) THEN
-c            PRINT*, "Trop de sous-mailles"
-c            CALL abort
-c          ENDIF
-c          WRITE(str2,'(i2.2)') nsrf
-c          ierr = NF_INQ_VARID (nid, "ALBLW"//str2, nvarid)
-c           IF (ierr.NE.NF_NOERR) THEN
-c             PRINT*, "phyetat0: Le champ <ALBLW"//str2//"> est absent"
-c             CALL abort
-c          ENDIF
-c#ifdef NC_DOUBLE
-c           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alblw(1,nsrf))
-c#else
-c           ierr = NF_GET_VAR_REAL(nid, nvarid, alblw(1,nsrf))
-c#endif
-c          IF (ierr.NE.NF_NOERR) THEN
-c            PRINT*, "phyetat0: Lecture echouee pour <ALBLW"//str2//">"
-c            CALL abort
-c          ENDIF
-c          xmin = 1.0E+20
-c          xmax = -1.0E+20
-c          DO i = 1, klon
-c             xmin = MIN(alblw(i,nsrf),xmin)
-c             xmax = MAX(alblw(i,nsrf),xmax)
-c          ENDDO
-c          PRINT*,'Albedo du sol ALBLW**:', nsrf, xmin, xmax
-c        ENDDO
       ELSE
          PRINT*, 'phyetat0: Le champ <ALBLW> est present'
Index: /LMDZ4/trunk/libf/phylmd/phyredem.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/phyredem.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/phyredem.F	(revision 541)
@@ -3,5 +3,4 @@
 !
 c
-cIM   SUBROUTINE phyredem (fichnom,dtime,radpas,co2_ppm,solaire,
       SUBROUTINE phyredem (fichnom,dtime,radpas,
      .           rlat,rlon, pctsrf,tsol,tsoil,deltat,qsurf,qsol,snow,
@@ -30,6 +29,4 @@
       INTEGER radpas
       REAL rlat(klon), rlon(klon)
-cIM   REAL co2_ppm
-cIM   REAL solaire
       REAL tsol(klon,nbsrf)
       REAL tsoil(klon,nsoilmx,nbsrf)
Index: /LMDZ4/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/physiq.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/physiq.F	(revision 541)
@@ -2,6 +2,4 @@
 ! $Header$
 !
-C
-c $Header$
 c
       SUBROUTINE physiq (nlon,nlev,nqmax,
@@ -88,4 +86,5 @@
 #include "advtrac.h"
 #include "iniprint.h"
+#include "thermcell.h"
 c======================================================================
       LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
@@ -157,4 +156,10 @@
       PARAMETER (ok_region=.FALSE.)
 c======================================================================
+c     pour phsystoke avec thermiques
+      REAL fm_therm(klon,klev+1)
+      REAL entr_therm(klon,klev)
+      real q2(klon,klev+1,nbsrf)
+      save q2
+c======================================================================
 c
       INTEGER ivap          ! indice de traceurs pour vapeur d'eau
@@ -594,7 +599,7 @@
       REAL zxffonte(klon), zxfqcalving(klon)
 
-      LOGICAL offline           ! Controle du stockage ds "physique"
-      PARAMETER (offline=.false.)
-      INTEGER physid
+c$$$      LOGICAL offline           ! Controle du stockage ds "physique"
+c$$$      PARAMETER (offline=.false.)
+c$$$      INTEGER physid
       REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
       save pfrac_impa
@@ -807,4 +812,5 @@
       REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
       REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
+      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
       REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
       REAL rneb(klon,klev)
@@ -1034,6 +1040,6 @@
      .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
      .                  ok_ade, ok_aie, 
-     .                  bl95_b0, bl95_b1)
-cIM  .                  , RI0)
+     .                  bl95_b0, bl95_b1,
+     .                  iflag_thermals,nsplit_thermals)
 
 c
@@ -1053,4 +1059,6 @@
      .       run_off_lic_0)
 
+c   ATTENTION : il faudra a terme relire q2 dans l'etat initial
+         q2(:,:,:)=1.e-8
 c
          radpas = NINT( 86400./dtime/nbapp_rad)
@@ -1508,4 +1516,5 @@
      s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
      s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
+     s            q2,
      s            dsens, devap,
      s            ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
@@ -1850,11 +1859,42 @@
 c Appeler l'ajustement sec
 c
-      CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
-      DO k = 1, klev
-      DO i = 1, klon
-         t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
-         q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
-      ENDDO
-      ENDDO
+c===================================================================
+c Convection seche (thermiques ou ajustement)
+c===================================================================
+c
+      d_t_ajs(:,:)=0.
+      d_u_ajs(:,:)=0.
+      d_v_ajs(:,:)=0.
+      d_q_ajs(:,:)=0.
+      fm_therm(:,:)=0.
+      entr_therm(:,:)=0.
+c
+      print*,'AVANT LA CONVECTION SECHE , iflag_thermals='
+     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+      if(iflag_thermals.lt.0) then
+c  Rien
+c  ====
+         print*,'pas de convection'
+      else if(iflag_thermals.eq.0) then
+
+c  Ajustement sec
+c  ==============
+         print*,'ajsec'
+         CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
+         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
+         q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
+      else
+c  Thermiques
+c  ==========
+         print*,'JUSTE AVANT , iflag_thermals='
+     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+         call calltherm(pdtphys
+     s      ,pplay,paprs,pphi
+     s      ,u_seri,v_seri,t_seri,q_seri
+     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
+     s      ,fm_therm,entr_therm)
+      endif
+c
+c===================================================================
 c
       IF (if_ebil.ge.2) THEN 
@@ -2573,6 +2613,6 @@
      I                   u,v,t,paprs,pplay,
      I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
-     I                   ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
-     I                   frac_impa, frac_nucl,
+     I                   ycoefh,fm_therm,entr_therm,yu1,yv1,ftsol,
+     I                   pctsrf,rlat,frac_impa, frac_nucl,
      I                   rlon,presnivs,pphis,pphi,
      I                   albsol,
@@ -2589,7 +2629,9 @@
       IF (offline) THEN
 
+         print*,'Attention on met a 0 les thermiques pour phystoke'
 	 call phystokenc (
      I                   nlon,nlev,pdtphys,rlon,rlat,
      I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
+     I                   fm_therm,entr_therm,
      I                   ycoefh,yu1,yv1,ftsol,pctsrf,
      I                   frac_impa, frac_nucl,
Index: /LMDZ4/trunk/libf/phylmd/phystokenc.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/phystokenc.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/phystokenc.F	(revision 541)
@@ -7,6 +7,7 @@
      I                   nlon,nlev,pdtphys,rlon,rlat,
      I                   pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
+     I                   pfm_therm,pentr_therm,
      I                   pcoefh,yu1,yv1,ftsol,pctsrf,
-     I                   pfrac_impa,pfrac_nucl,
+     I                   frac_impa,frac_nucl,
      I                   pphis,paire,dtime,itap)
       USE ioipsl
@@ -40,5 +41,7 @@
       real pdtphys ! pas d'integration pour la physique (seconde)
 c
-      integer physid, itap,ndex(1)
+      integer physid, itap
+      save physid
+      integer ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
 
 c   convection:
@@ -51,5 +54,5 @@
       REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
       REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
-	REAL pt(klon,klev)
+	real pt(klon,klev),t(klon,klev)
 c
       REAL rlon(klon), rlat(klon), dtime
@@ -62,10 +65,18 @@
       REAL yv1(klon)
       REAL yu1(klon),pphis(klon),paire(klon)
+
+c   Les Thermiques : (Abderr 25 11 02)
+c   ---------------
+      REAL pfm_therm(klon,klev+1)
+	real fm_therm1(klon,klev)
+      REAL pentr_therm(klon,klev)
+      REAL entr_therm(klon,klev)
+      REAL fm_therm(klon,klev)
 c
 c   Lessivage:
 c   ----------
 c
-      REAL pfrac_impa(klon,klev)
-      REAL pfrac_nucl(klon,klev)
+      REAL frac_impa(klon,klev)
+      REAL frac_nucl(klon,klev)
 c
 c Arguments necessaires pour les sources et puits de traceur
@@ -84,8 +95,4 @@
       REAL de_d(klon,klev) ! flux detraine dans le panache descendant
       REAL coefh(klon,klev) ! flux detraine dans le panache descendant
-	REAL t(klon,klev)
-      REAL frac_impa(klon,klev)
-      REAL frac_nucl(klon,klev)
-      REAL rain(klon)
 
       REAL pyu1(klon),pyv1(klon)
@@ -98,8 +105,9 @@
       integer iadvtr,irec
       real zmin,zmax
-
+      logical ok_sync
+ 
       save t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,dtcum
+	save fm_therm,entr_therm
       save iadvtr,irec
-      save frac_impa,frac_nucl,rain
       save pyu1,pyv1,pftsol,ppsrf
 
@@ -109,4 +117,6 @@
 c======================================================================
 
+      ok_sync = .true.
+	print*,'Dans phystokenc.F'
       print*,'iadvtr= ',iadvtr
       print*,'istphy= ',istphy
@@ -120,19 +130,19 @@
 	write(*,*) 'apres initphysto ds phystokenc' 
 
-       ndex(1) = 0
-         i=itap 
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
-         CALL histwrite(physid,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
-c
-         i=itap
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
-         CALL histwrite(physid,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
 	
       ENDIF
 c
+      ndex2d = 0
+      ndex3d = 0
+      i=itap 
+      CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
+      CALL histwrite(physid,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex2d)
+c
+      i=itap
+      CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
+      CALL histwrite(physid,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex2d)
+
       iadvtr=iadvtr+1
 c
-c
-c   reinitialisation des champs cumules
       if (mod(iadvtr,istphy).eq.1.or.istphy.eq.1) then
 	print*,'reinitialisation des champs cumules 
@@ -140,6 +150,4 @@
          do k=1,klev
             do i=1,klon
-               frac_impa(i,k)=1.
-               frac_nucl(i,k)=1.
                mfu(i,k)=0.
                mfd(i,k)=0.
@@ -149,9 +157,10 @@
                de_d(i,k)=0.
                coefh(i,k)=0.
-		t(i,k)=0.
-            enddo
-         enddo
-         do i=1,klon
-            rain(i)=0.
+                t(i,k)=0.
+		fm_therm(i,k)=0.
+               entr_therm(i,k)=0.
+            enddo
+         enddo
+         do i=1,klon
             pyv1(i)=0.
             pyu1(i)=0.
@@ -169,6 +178,217 @@
       do k=1,klev
          do i=1,klon
-            frac_impa(i,k)=frac_impa(i,k)*pfrac_impa(i,k)
-            frac_nucl(i,k)=frac_nucl(i,k)*pfrac_nucl(i,k)
+            mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
+            mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
+            en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
+            de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
+            en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
+            de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
+            coefh(i,k)=coefh(i,k)+pcoefh(i,k)*pdtphys
+                t(i,k)=t(i,k)+pt(i,k)*pdtphys
+       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
+       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
+         enddo
+      enddo
+         do i=1,klon
+            pyv1(i)=pyv1(i)+yv1(i)*pdtphys
+            pyu1(i)=pyu1(i)+yu1(i)*pdtphys
+         end do
+         do k=1,nbsrf
+             do i=1,klon
+               pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
+               ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
+            enddo
+         enddo
+
+      dtcum=dtcum+pdtphys
+
+      IF(mod(iadvtr,istphy).eq.0) THEN 
+c
+c   normalisation par le temps cumule
+         do k=1,klev
+            do i=1,klon
+               mfu(i,k)=mfu(i,k)/dtcum
+               mfd(i,k)=mfd(i,k)/dtcum
+               en_u(i,k)=en_u(i,k)/dtcum
+               de_u(i,k)=de_u(i,k)/dtcum
+               en_d(i,k)=en_d(i,k)/dtcum
+               de_d(i,k)=de_d(i,k)/dtcum
+               coefh(i,k)=coefh(i,k)/dtcum
+c Unitel a enlever
+	      t(i,k)=t(i,k)/dtcum	
+               fm_therm(i,k)=fm_therm(i,k)/dtcum
+	       entr_therm(i,k)=entr_therm(i,k)/dtcum
+            enddo
+         enddo
+         do i=1,klon
+            pyv1(i)=pyv1(i)/dtcum
+            pyu1(i)=pyu1(i)/dtcum
+         end do
+         do k=1,nbsrf
+             do i=1,klon
+               pftsol(i,k)=pftsol(i,k)/dtcum
+               pftsol1(i) = pftsol(i,1)
+               pftsol2(i) = pftsol(i,2)
+               pftsol3(i) = pftsol(i,3)
+               pftsol4(i) = pftsol(i,4)
+
+               ppsrf(i,k)=ppsrf(i,k)/dtcum
+               ppsrf1(i) = ppsrf(i,1)
+               ppsrf2(i) = ppsrf(i,2)
+               ppsrf3(i) = ppsrf(i,3)
+               ppsrf4(i) = ppsrf(i,4)
+
+            enddo
+         enddo
+c
+c   ecriture des champs
+c
+         irec=irec+1
+
+ccccc
+         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
+         CALL histwrite(physid,"t",itap,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+
+         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfu, zx_tmp_3d)
+      CALL histwrite(physid,"mfu",itap,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+	CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfd, zx_tmp_3d)
+      CALL histwrite(physid,"mfd",itap,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_u, zx_tmp_3d)
+      CALL histwrite(physid,"en_u",itap,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_u, zx_tmp_3d)
+      CALL histwrite(physid,"de_u",itap,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_d, zx_tmp_3d)
+      CALL histwrite(physid,"en_d",itap,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_d, zx_tmp_3d)       
+      CALL histwrite(physid,"de_d",itap,zx_tmp_3d,    
+     .                                   iim*(jjm+1)*klev,ndex3d)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, coefh, zx_tmp_3d)         
+      CALL histwrite(physid,"coefh",itap,zx_tmp_3d,    
+     .                                   iim*(jjm+1)*klev,ndex3d)	
+
+c ajou...
+	do k=1,klev
+           do i=1,klon
+	 fm_therm1(i,k)=fm_therm(i,k)	
+	   enddo
+	enddo
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, fm_therm1, zx_tmp_3d)
+      CALL histwrite(physid,"fm_th",itap,zx_tmp_3d,
+     .                                 iim*(jjm+1)*klev,ndex3d)
+c
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, entr_therm, zx_tmp_3d)
+      CALL histwrite(physid,"en_th",itap,zx_tmp_3d,
+     .                                iim*(jjm+1)*klev,ndex3d)
+cccc
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_impa,zx_tmp_3d)
+        CALL histwrite(physid,"frac_impa",itap,zx_tmp_3d,
+     .  iim*(jjm+1)*klev,ndex3d)
+
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_nucl,zx_tmp_3d)
+        CALL histwrite(physid,"frac_nucl",itap,zx_tmp_3d,
+     .  iim*(jjm+1)*klev,ndex3d)
+ 
+        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyu1,zx_tmp_2d)
+      CALL histwrite(physid,"pyu1",itap,zx_tmp_2d,iim*(jjm+1),
+     .                                                ndex2d)
+	
+	CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyv1,zx_tmp_2d)
+      CALL histwrite(physid,"pyv1",itap,zx_tmp_2d,iim*(jjm+1)
+     .                                                ,ndex2d)
+	
+	CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol1, zx_tmp_2d)
+      CALL histwrite(physid,"ftsol1",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex2d)
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol2, zx_tmp_2d)
+      CALL histwrite(physid,"ftsol2",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex2d)
+          CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol3, zx_tmp_2d)
+      CALL histwrite(physid,"ftsol3",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex2d)
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol4, zx_tmp_2d)
+      CALL histwrite(physid,"ftsol4",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex2d)
+
+        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf1, zx_tmp_2d)
+      CALL histwrite(physid,"psrf1",itap,zx_tmp_2d,   
+     .                                   iim*(jjm+1),ndex2d)
+        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf2, zx_tmp_2d)
+      CALL histwrite(physid,"psrf2",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex2d)
+        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf3, zx_tmp_2d)
+      CALL histwrite(physid,"psrf3",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex2d)
+        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf4, zx_tmp_2d)
+      CALL histwrite(physid,"psrf4",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex2d)
+
+      if (ok_sync) call histsync(physid)
+c     if (ok_sync) call histsync
+	
+c
+cAA Test sur la valeur des coefficients de lessivage 
+c
+         zmin=1e33
+         zmax=-1e33
+         do k=1,klev
+            do i=1,klon
+                  zmax=max(zmax,frac_nucl(i,k))
+                  zmin=min(zmin,frac_nucl(i,k))
+            enddo
+         enddo
+         Print*,'------ coefs de lessivage (min et max) --------'
+         Print*,'facteur de nucleation ',zmin,zmax
+         zmin=1e33
+         zmax=-1e33
+         do k=1,klev
+            do i=1,klon
+                  zmax=max(zmax,frac_impa(i,k))
+                  zmin=min(zmin,frac_impa(i,k))
+            enddo
+         enddo
+         Print*,'facteur d impaction ',zmin,zmax
+
+      ENDIF 
+
+c   reinitialisation des champs cumules
+	go to 768
+      if (mod(iadvtr,istphy).eq.1) then
+         do k=1,klev
+            do i=1,klon
+               mfu(i,k)=0.
+               mfd(i,k)=0.
+               en_u(i,k)=0.
+               de_u(i,k)=0.
+               en_d(i,k)=0.
+               de_d(i,k)=0.
+               coefh(i,k)=0.
+	       t(i,k)=0.
+               fm_therm(i,k)=0.
+	       entr_therm(i,k)=0.
+            enddo
+         enddo
+         do i=1,klon
+            pyv1(i)=0.
+            pyu1(i)=0.
+         end do
+         do k=1,nbsrf
+             do i=1,klon
+               pftsol(i,k)=0.
+               ppsrf(i,k)=0.
+            enddo
+         enddo
+
+         dtcum=0.
+      endif
+
+      do k=1,klev
+         do i=1,klon
             mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
             mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
@@ -179,4 +399,6 @@
             coefh(i,k)=coefh(i,k)+pcoefh(i,k)*pdtphys
 		t(i,k)=t(i,k)+pt(i,k)*pdtphys
+       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
+       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
          enddo
       enddo
@@ -193,151 +415,5 @@
 
       dtcum=dtcum+pdtphys
-c
-      IF(mod(iadvtr,istphy).eq.0) THEN
-c
-c   normalisation par le temps cumule
-         do k=1,klev
-            do i=1,klon
-c              frac_impa=frac_impa : c'est la fraction cumulee qu'on stoke
-c              frac_nucl=frac_nucl : c'est la fraction cumulee qu'on stoke
-               mfu(i,k)=mfu(i,k)/dtcum
-               mfd(i,k)=mfd(i,k)/dtcum
-               en_u(i,k)=en_u(i,k)/dtcum
-               de_u(i,k)=de_u(i,k)/dtcum
-               en_d(i,k)=en_d(i,k)/dtcum
-               de_d(i,k)=de_d(i,k)/dtcum
-               coefh(i,k)=coefh(i,k)/dtcum
-		t(i,k)=t(i,k)/dtcum
-            enddo
-         enddo
-         do i=1,klon
-            rain(i)=rain(i)/dtcum
-            pyv1(i)=pyv1(i)/dtcum
-            pyu1(i)=pyu1(i)/dtcum
-         end do
-c modif abderr 23 11 00         do k=1,nbsrf
-             do i=1,klon
-	      do k=1,nbsrf
-               pftsol(i,k)=pftsol(i,k)/dtcum
-	       ppsrf(i,k)=ppsrf(i,k)/dtcum
-	      enddo
-               pftsol1(i) = pftsol(i,1)
-               pftsol2(i) = pftsol(i,2)
-               pftsol3(i) = pftsol(i,3)
-               pftsol4(i) = pftsol(i,4)
-
-c               ppsrf(i,k)=ppsrf(i,k)/dtcum
-               ppsrf1(i) = ppsrf(i,1)
-               ppsrf2(i) = ppsrf(i,2)
-               ppsrf3(i) = ppsrf(i,3)
-               ppsrf4(i) = ppsrf(i,4)
-
-            enddo
-c         enddo
-c
-c   ecriture des champs
-c
-         irec=irec+1
-
-ccccc
-      print*,'AVANT ECRITURE'
-         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
-         CALL histwrite(physid,"t",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-      print*,'APRES ECRITURE'
-
-         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfu, zx_tmp_3d)
-      CALL histwrite(physid,"mfu",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfd, zx_tmp_3d)
-      CALL histwrite(physid,"mfd",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_u, zx_tmp_3d)
-      CALL histwrite(physid,"en_u",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_u, zx_tmp_3d)
-      CALL histwrite(physid,"de_u",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_d, zx_tmp_3d)
-      CALL histwrite(physid,"en_d",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_d, zx_tmp_3d)
-      CALL histwrite(physid,"de_d",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, coefh, zx_tmp_3d)
-      CALL histwrite(physid,"coefh",itap,zx_tmp_3d,
-     .                                   iim*(jjm+1)*klev,ndex)
-cccc
-       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_impa,zx_tmp_3d)
-        CALL histwrite(physid,"frac_impa",itap,zx_tmp_3d,
-     .  iim*(jjm+1)*klev,ndex)
-
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_nucl,zx_tmp_3d)
-        CALL histwrite(physid,"frac_nucl",itap,zx_tmp_3d,
-     .  iim*(jjm+1)*klev,ndex)
-
-        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyu1,zx_tmp_2d)
-      CALL histwrite(physid,"pyu1",itap,zx_tmp_2d,iim*(jjm+1),ndex)
-
-        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyv1,zx_tmp_2d)
-      CALL histwrite(physid,"pyv1",itap,zx_tmp_2d,iim*(jjm+1),ndex)
-
-        CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol1, zx_tmp_2d)
-      CALL histwrite(physid,"ftsol1",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol2, zx_tmp_2d)
-      CALL histwrite(physid,"ftsol2",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-          CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol3, zx_tmp_2d)
-      CALL histwrite(physid,"ftsol3",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-
-c
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol4, zx_tmp_2d)
-      CALL histwrite(physid,"ftsol4",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-
-        CALL gr_fi_ecrit(1,klon,iim,jjm+1, rain, zx_tmp_2d)
-      CALL histwrite(physid,"rain",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-
-        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf1, zx_tmp_2d)
-      CALL histwrite(physid,"psrf1",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf2, zx_tmp_2d)
-      CALL histwrite(physid,"psrf2",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf3, zx_tmp_2d)
-      CALL histwrite(physid,"psrf3",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf4, zx_tmp_2d)
-      CALL histwrite(physid,"psrf4",itap,zx_tmp_2d,
-     .                                   iim*(jjm+1),ndex)
-
-c
-cAA Test sur la valeur des coefficients de lessivage
-c
-         zmin=1e33
-         zmax=-1e33
-         do k=1,klev
-            do i=1,klon
-                  zmax=max(zmax,frac_nucl(i,k))
-                  zmin=min(zmin,frac_nucl(i,k))
-            enddo
-         enddo
-         Print*,'------ coefs de lessivage (min et max) --------'
-         Print*,'facteur de nucleation ',zmin,zmax
-         zmin=1e33
-         zmax=-1e33
-         do k=1,klev
-            do i=1,klon
-                  zmax=max(zmax,frac_impa(i,k))
-                  zmin=min(zmin,frac_impa(i,k))
-            enddo
-         enddo
-         Print*,'facteur d impaction ',zmin,zmax
-
-      ENDIF
-
+768   continue
 
       RETURN
Index: /LMDZ4/trunk/libf/phylmd/phytrac.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/phytrac.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/phytrac.F	(revision 541)
@@ -23,4 +23,5 @@
      I                    pde_d,
      I                    coefh,
+     I                    fm_therm,entr_therm,
      I                    yu1,
      I                    yv1,
@@ -130,5 +131,6 @@
       logical debutphy       ! le flag de l'initialisation de la physique
       logical lafin          ! le flag de la fin de la physique
-
+c Olivia      
+      integer isplit,nsplit
       REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)   !--lessivage convection
       REAL prfl(klon,klev+1),   psfl(klon,klev+1)     !--lessivage large-scale
@@ -146,4 +148,12 @@
       REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
       REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
+
+c
+c   thermiques:
+c   -----------
+c
+      real fm_therm(klon,klev+1),entr_therm(klon,klev)
+	real fm_therm1(klon,klev)
+c
       REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
       REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
@@ -265,5 +275,8 @@
       REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs 
       REAL d_tr_cl(klon,klev) ! tendance de traceurs  couche limite
-      REAL d_tr_cv(klon,klev) ! tendance de traceurs  convection 
+      REAL d_tr_cli(klon,klev,nbtr) ! tendance de traceurs  CL pour chq traceur
+      REAL d_tr_cv(klon,klev) ! tendance de traceurs  convection
+      REAL d_tr_cvi(klon,klev,nbtr) ! tendance de traceurs  conv pour chq traceur
+      REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
       REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance 
 c                                   ! radioactive du rn - > pb 
@@ -280,4 +293,6 @@
       REAL flestottr(klon,klev,nbtr) ! flux de lessivage 
 c                                    ! dans chaque couche 
+      real zmasse(klon,klev)
+      real ztra_th(klon,klev)
 
 C
@@ -288,9 +303,14 @@
 c-------------
       logical first,couchelimite,convection,lessivage,sorties,
-     s        rnpb,inirnpb
-      save first,couchelimite,convection,lessivage,sorties,
-     s     inirnpb
-      data first,couchelimite,convection,lessivage,sorties
-     s     /.true.,.true.,.false.,.true.,.true./
+     s        rnpb,inirnpb,thermiques
+      save first,couchelimite,convection,lessivage,thermiques,
+     s        sorties,inirnpb
+c      data first,couchelimite,convection,lessivage,sorties
+c     s     /.true.,.true.,.false.,.true.,.true./
+c Olivia
+       data first,couchelimite,convection,lessivage,
+     s      thermiques,sorties
+     s     /.true.,.true.,.true.,.true.,.true.,.true./
+
 
 #ifdef INCA
@@ -522,4 +542,6 @@
       DO i = 1, klon
          tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
+	 d_tr_cvi(i,k,it)=d_tr_cv(i,k)
+c	 print*,'en k i d_tr_cv=',k,i,d_tr_cv(i,k)
       ENDDO
       ENDDO
@@ -541,4 +563,63 @@
 c      enddo
 
+
+c======================================================================
+c   Calcul de l'effet des thermiques
+c======================================================================
+
+      do k=1,klev
+         do i=1,klon
+            zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
+         enddo
+      enddo
+
+c      print*,'masse dans ph ',zmasse
+      do it=1,nqmax
+         do k=1,klev
+            do i=1,klon
+               d_tr_th(i,k,it)=0.
+               tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
+               tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
+            enddo
+         enddo
+      enddo
+
+      if (thermiques) then
+        print*,'calcul de leffet des thermiques'
+        nsplit=10
+        DO it=1, nqmax
+c        WRITE(itn,'(i1)') it
+c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
+c	     print*,'avant dqthermiquesretro'
+c             call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI      ')
+
+         do isplit=1,nsplit
+c  Abderr 25 11 02
+C Thermiques 
+c	print*,'Avant dans phytrac',avant
+            call dqthermcell(klon,klev,pdtphys/nsplit
+     .       ,fm_therm,entr_therm,zmasse
+     .       ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
+
+            do k=1,klev
+               do i=1,klon
+                  d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
+                  d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
+                  tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
+               enddo
+            enddo
+          enddo ! nsplit
+          print*,'apres thermiques'
+c             call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th      ')
+c            do k=1,klev
+c	print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
+c	   enddo
+
+c        WRITE(itn,'(i1)') it
+c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
+       ENDDO ! it
+       endif ! Thermiques
+c       print*,'ATTENTION: sdans thermniques'
+      
 c======================================================================
 c   Calcul de l'effet de la couche limite
@@ -576,4 +657,5 @@
             DO i = 1, klon
               tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
+	      d_tr_cli(i,k,it)=d_tr_cl(i,k)
             ENDDO
           ENDDO
@@ -608,4 +690,5 @@
                DO i = 1, klon
                   tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
+		  d_tr_cli(i,k,it)=d_tr_cl(i,k)
                ENDDO
             ENDDO
Index: /LMDZ4/trunk/libf/phylmd/read_pstoke.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/read_pstoke.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/read_pstoke.F	(revision 541)
@@ -7,5 +7,14 @@
      .   zrec,zklono,zklevo,airefi,phisfi,
      .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,
+     .   fm_therm,en_therm,
      .   frac_impa,frac_nucl,pyu1,pyv1,ftsol,psrf)
+
+C******************************************************************************
+C  Frederic HOURDIN, Abderrahmane IDELKADI
+C Lecture des parametres physique stockes online necessaires pour
+C recalculer offline le transport de traceurs sur une grille 2x plus fine que 
+C celle online 
+C A FAIRE : une seule routine au lieu de 2 (lectflux, redecoupe)!
+C******************************************************************************
 
 
@@ -26,28 +35,31 @@
 #include "dimphy.h"
 	
-	  integer*4 klono,klevo,imo,jmo
+	  integer 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)
+	  REAL phisfi(klono)
+          REAL phisfi2(imo,jmo+1),airefi2(imo,jmo+1)
+
+          REAL mfu(klono,klevo), mfd(klono,klevo)
+          REAL en_u(klono,klevo), de_u(klono,klevo)
+          REAL en_d(klono,klevo), de_d(klono,klevo)
+          REAL coefh(klono,klevo)
+           REAL fm_therm(klono,klevo),en_therm(klono,klevo)
+
+          REAL mfu2(imo,jmo+1,klevo), mfd2(imo,jmo+1,klevo)
+          REAL en_u2(imo,jmo+1,klevo), de_u2(imo,jmo+1,klevo)
+          REAL en_d2(imo,jmo+1,klevo), de_d2(imo,jmo+1,klevo)
+          REAL coefh2(imo,jmo+1,klevo)
+           REAL fm_therm2(imo,jmo+1,klevo)
+           REAL en_therm2(imo,jmo+1,klevo)
+
+          REAL pl(klevo)
           integer irec
-          integer*4 xid,yid,zid,tid
+          integer xid,yid,zid,tid
           real zrec,zklono,zklevo,zim,zjm
-          integer*4 ncrec,ncklono,ncklevo,ncim,ncjm
-
-          real*4 airefi(klono)
-          character namedim
+          integer ncrec,ncklono,ncklevo,ncim,ncjm
+
+          real airefi(klono)
+          character*20 namedim
 
 c  !! attention !!
@@ -56,20 +68,20 @@
 	  
 	 
-          REAL*4 frac_impa(klono,klevo), frac_nucl(klono,klevo)
-          REAL*4 frac_impa2(imo,jmo+1,klevo), 
+          REAL frac_impa(klono,klevo), frac_nucl(klono,klevo)
+          REAL 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),
+          REAL pyu1(klono), pyv1(klono)
+          REAL pyu12(imo,jmo+1), pyv12(imo,jmo+1)
+          REAL ftsol(klono,nbsrf)
+          REAL psrf(klono,nbsrf) 
+	  REAL ftsol1(klono),ftsol2(klono),ftsol3(klono),ftsol4(klono)
+          REAL psrf1(klono),psrf2(klono),psrf3(klono),psrf4(klono)
+          REAL 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),
+          REAL 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)	
+		REAL t(klono,klevo)
+		REAL t2(imo,jmo+1,klevo)	
 	  integer ncidp
           save ncidp
@@ -77,4 +89,5 @@
           integer varidmfu, varidmfd, varidps, varidenu, variddeu	
           integer varidend,varidded,varidch,varidfi,varidfn
+           integer varidfmth,varidenth
           integer varidyu1,varidyv1,varidpl,varidai,varididvt
           integer varidfts1,varidfts2,varidfts3,varidfts4
@@ -82,4 +95,5 @@
           save varidmfu, varidmfd, varidps, varidenu, variddeu
           save varidend,varidded,varidch,varidfi,varidfn
+           save varidfmth,varidenth
           save varidyu1,varidyv1,varidpl,varidai,varididvt
           save varidfts1,varidfts2,varidfts3,varidfts4
@@ -113,6 +127,8 @@
             print*,'ncidp,varidai',ncidp,varidai
 
+c A FAIRE: Es-il necessaire de stocke t?
 	        varidt=NCVID(ncidp,'t',rcode)
                 print*,'ncidp,varidt',ncidp,varidt
+
             varidmfu=NCVID(ncidp,'mfu',rcode)
             print*,'ncidp,varidmfu',ncidp,varidmfu
@@ -136,4 +152,11 @@
             print*,'ncidp,varidch',ncidp,varidch
 	
+c abder (pour thermiques)
+             varidfmth=NCVID(ncidp,'fm_th',rcode)
+             print*,'ncidp,varidfmth',ncidp,varidfmth
+
+             varidenth=NCVID(ncidp,'en_th',rcode)
+             print*,'ncidp,varidenth',ncidp,varidenth
+
 	    varidfi=NCVID(ncidp,'frac_impa',rcode)
             print*,'ncidp,varidfi',ncidp,varidfi
@@ -200,6 +223,9 @@
 
 c niveaux de pression
-
-            status=NF_GET_VARA_REAL(ncidp,varidpl,1,zklevo,pl)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpl,1,zklevo,pl)
+#else
+      status=NF_GET_VARA_REAL(ncidp,varidpl,1,zklevo,pl)
+#endif
 
 c lecture de aire et phis
@@ -216,15 +242,17 @@
 
 c phis
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidps,start,count,phisfi2)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,phisfi2,phisfi)
 
 c aire
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidai,start,count,airefi2)
+#else
       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')
+#endif
        call gr_ecrit_fi(1,klono,imo,jmo+1,airefi2,airefi)
       else
@@ -249,165 +277,190 @@
       count(4)=1
 
+
+C *** Lessivage******************************************************
 c frac_impa 
-
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfi,start,count,frac_impa2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_impa2,frac_impa)
 
 c frac_nucl 
-
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfn,start,count,frac_nucl2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_nucl2,frac_nucl)
 
+C*** Temperature ******************************************************
 c abder t
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidt,start,count,t2)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,t2,t)
 
+C*** Flux pour le calcul de la convection TIEDTK ***********************
 c mfu
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidmfu,start,count,mfu2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfu2,mfu)
 
 c mfd
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidmfd,start,count,mfd2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfd2,mfd)
 
 c en_u 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidenu,start,count,en_u2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_u2,en_u)
 
 c de_u 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,variddeu,start,count,de_u2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_u2,de_u)
 
 c en_d 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidend,start,count,en_d2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_d2,en_d)
 
 c de_d 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidded,start,count,de_d2)
+#else
       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')
+#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_d2,de_d)
 
+C **** Coeffecient du mellange turbulent**********************************
 c coefh 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidch,start,count,coefh2)
+#else
       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')
+#endif
        call gr_ecrit_fi(klevo,klono,imo,jmo+1,coefh2,coefh)
 
+C*** Flux ascendant et entrant pour les Thermiques************************
+cabder thermiques
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfmth,start,count,fm_therm2)
+#else
+      status=NF_GET_VARA_REAL(ncidp,varidfmth,start,count,fm_therm2)
+#endif
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,fm_therm2,fm_therm)
+
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidenth,start,count,en_therm2)
+#else
+      status=NF_GET_VARA_REAL(ncidp,varidenth,start,count,en_therm2)
+#endif
+      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_therm2,en_therm)
+
+C*** Vitesses aux sol ******************************************************
       start(3)=irec
       start(4)=0
       count(3)=1
       count(4)=0
-
 c pyu1
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidyu1,start,count,pyu12)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,pyu12,pyu1)
 
 c pyv1
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidyv1,start,count,pyv12)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,pyv12,pyv1)
 
+C*** Temperature au sol ********************************************
 c ftsol1
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts1,start,count,ftsol12)
+#else
       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')
+#endif
        call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol12,ftsol1)
 
 c ftsol2
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts2,start,count,ftsol22)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol22,ftsol2)
 
 c ftsol3
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts3,start,count,ftsol32)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol32,ftsol3)
 
 c ftsol4
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts4,start,count,ftsol42)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol42,ftsol4)
 
+C*** Nature du sol **************************************************
 c psrf1 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr1,start,count,psrf12)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,psrf12,psrf1)
 
 c psrf2 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr2,start,count,psrf22)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,psrf22,psrf2)
 
 c psrf3 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr3,start,count,psrf32)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,psrf32,psrf3)
 
 c psrf4 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr4,start,count,psrf42)
+#else
       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')
+#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,psrf42,psrf4)
 	
Index: /LMDZ4/trunk/libf/phylmd/read_pstoke0.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/read_pstoke0.F	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/read_pstoke0.F	(revision 541)
@@ -7,5 +7,13 @@
      .   zrec,zkon,zkev,airefi,phisfi,
      .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,
+     .   fm_therm,en_therm,
      .   frac_impa,frac_nucl,pyu1,pyv1,ftsol,psrf)
+
+C******************************************************************************
+C  Frederic HOURDIN, Abderrahmane IDELKADI
+C Lecture des parametres physique stockes online necessaires pour
+C recalculer offline le transport des traceurs sur la meme grille que online
+C A FAIRE : une seule routine au lieu de 2 (lectflux, redecoupe)!
+C******************************************************************************
 
 
@@ -26,29 +34,36 @@
 #include "dimphy.h"
 	  
-	  integer*4 kon,kev,zkon,zkev
+	  integer 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)
+	  REAL phisfi(kon)
+          REAL phisfi2(iim,jjm+1),airefi2(iim,jjm+1)
+
+          REAL mfu(kon,kev), mfd(kon,kev)
+          REAL en_u(kon,kev), de_u(kon,kev)
+          REAL en_d(kon,kev), de_d(kon,kev)
+          REAL coefh(kon,kev)
+
+c abd 25 11 02
+c Thermiques
+	 REAL fm_therm(kon,kev),en_therm(kon,kev)
+		REAL t(kon,kev)
+
+          REAL mfu2(iim,jjm+1,kev), mfd2(iim,jjm+1,kev)
+          REAL en_u2(iim,jjm+1,kev), de_u2(iim,jjm+1,kev)
+          REAL en_d2(iim,jjm+1,kev), de_d2(iim,jjm+1,kev)
+          REAL coefh2(iim,jjm+1,kev)
+		REAL t2(iim,jjm+1,kev)
+c Thermiques
+	 REAL fm_therm2(iim,jjm+1,kev)
+         REAL en_therm2(iim,jjm+1,kev)       
+
+          REAL 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
+          integer xid,yid,zid,tid
+          integer zrec,zim,zjm
+          integer ncrec,nckon,nckev,ncim,ncjm
+
+          real airefi(kon)
+          character*20 namedim
 
 c  !! attention !!
@@ -56,17 +71,17 @@
 c dim de phis??
 
-          REAL*4 frac_impa(kon,kev), frac_nucl(kon,kev)
-          REAL*4 frac_impa2(iim,jjm+1,kev), 
+          REAL frac_impa(kon,kev), frac_nucl(kon,kev)
+          REAL 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),
+          REAL pyu1(kon), pyv1(kon)
+          REAL pyu12(iim,jjm+1), pyv12(iim,jjm+1)
+          REAL ftsol(kon,nbsrf)
+          REAL psrf(kon,nbsrf) 
+	  REAL ftsol1(kon),ftsol2(kon),ftsol3(kon),ftsol4(kon)
+          REAL psrf1(kon),psrf2(kon),psrf3(kon),psrf4(kon)
+          REAL 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),
+          REAL psrf12(iim,jjm+1),psrf22(iim,jjm+1),psrf32(iim,jjm+1),
      .     psrf42(iim,jjm+1)
 	
@@ -76,4 +91,6 @@
 		integer varidt
           integer varidend,varidded,varidch,varidfi,varidfn
+c therm
+	  integer varidfmth,varidenth
           integer varidyu1,varidyv1,varidpl,varidai,varididvt
           integer varidfts1,varidfts2,varidfts3,varidfts4
@@ -82,4 +99,6 @@
 		save varidt
           save varidend,varidded,varidch,varidfi,varidfn
+c therm
+	   save varidfmth,varidenth
           save varidyu1,varidyv1,varidpl,varidai,varididvt
           save varidfts1,varidfts2,varidfts3,varidfts4
@@ -112,10 +131,10 @@
             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
 
-		varidt=NCVID(ncidp,'t',rcode)
-                print*,'ncidp,varidt',ncidp,varidt
-
             varidmfd=NCVID(ncidp,'mfd',rcode)
             print*,'ncidp,varidmfd',ncidp,varidmfd
@@ -135,4 +154,11 @@
             varidch=NCVID(ncidp,'coefh',rcode)
             print*,'ncidp,varidch',ncidp,varidch
+
+c Thermiques
+            varidfmth=NCVID(ncidp,'fm_th',rcode)
+            print*,'ncidp,varidfmth',ncidp,varidfmth
+
+            varidenth=NCVID(ncidp,'en_th',rcode)
+            print*,'ncidp,varidenth',ncidp,varidenth
 	
 	    varidfi=NCVID(ncidp,'frac_impa',rcode)
@@ -216,10 +242,20 @@
 
 c 
+C**** Geopotentiel au sol ***************************************
 c phis
-      status=NF_GET_VARA_REAL(ncidp,varidps,start,count,phisfi2)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidps,start,count,phisfi2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidps,start,count,phisfi2)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,phisfi2,phisfi)
 
+C**** Aires des mails aux sol ************************************
 c aire
-      status=NF_GET_VARA_REAL(ncidp,varidai,start,count,airefi2)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidai,start,count,airefi2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidai,start,count,airefi2)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,airefi2,airefi)
       else
@@ -244,46 +280,116 @@
       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**** Temperature ********************************************
+cA FAIRE : Es-ce necessaire ?
 
 c abder t
-      status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidt,start,count,t2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
+#endif
       call gr_ecrit_fi(kev,kon,iim,jjm+1,t2,t)
 
+C**** Flux pour la convection (Tiedtk) ********************************************
 c mfu
-      status=NF_GET_VARA_REAL(ncidp,varidmfu,start,count,mfu2)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidmfu,start,count,mfu2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidmfu,start,count,mfu2)
+#endif
       call gr_ecrit_fi(kev,kon,iim,jjm+1,mfu2,mfu)
 
 c mfd
-      status=NF_GET_VARA_REAL(ncidp,varidmfd,start,count,mfd2)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidmfd,start,count,mfd2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidmfd,start,count,mfd2)
+#endif
       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)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidenu,start,count,en_u2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidenu,start,count,en_u2)
+#endif
       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)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,variddeu,start,count,de_u2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,variddeu,start,count,de_u2)
+#endif
       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)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidend,start,count,en_d2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidend,start,count,en_d2)
+#endif
       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)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidded,start,count,de_d2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidded,start,count,de_d2)
+#endif
       call gr_ecrit_fi(kev,kon,iim,jjm+1,de_d2,de_d)
 
+C**** Coefficient de mellange turbulent *******************************************
 c coefh 
 	print*,'LECTURE de coefh a irec =',irec
-       status=NF_GET_VARA_REAL(ncidp,varidch,start,count,coefh2)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidch,start,count,coefh2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidch,start,count,coefh2)
+#endif
        call gr_ecrit_fi(kev,kon,iim,jjm+1,coefh2,coefh)
+c      call dump2d(iip1,jjp1,coefh2(1,2),'COEFH2READ   ')
+c      call dump2d(iim ,jjm ,coefh (2,2),'COEFH2READ   ')
+
+C**** Flux ascendants et entrant dans le thermique **********************************
+cThermiques
+       print*,'LECTURE de fm_therm a irec =',irec
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfmth,start,
+     .                         count,fm_therm2)
+#else
+       status=NF_GET_VARA_REAL(ncidp,varidfmth,start,
+     .                         count,fm_therm2)
+#endif
+       call gr_ecrit_fi(kev,kon,iim,jjm+1,fm_therm2,fm_therm)
+       print*,'LECTURE de en_therm a irec =',irec
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidenth,start,
+     .                          count,en_therm2)
+#else
+       status=NF_GET_VARA_REAL(ncidp,varidenth,start,
+     .                          count,en_therm2)
+#endif
+       call gr_ecrit_fi(kev,kon,iim,jjm+1,en_therm2,en_therm)
+
+C**** Coefficients de lessivage *******************************************
+c frac_impa
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfi,start,count,frac_impa2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidfi,start,count,frac_impa2)
+#endif
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,frac_impa2,frac_impa)
+
+c frac_nucl
+
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfn,start,count,frac_nucl2)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidfn,start,count,frac_nucl2)
+#endif
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,frac_nucl2,frac_nucl)
+
+C**** Vents aux sol ********************************************
 
       start(3)=irec
@@ -294,47 +400,89 @@
 c pyu1
 	print*,'LECTURE de yu1 a irec =',irec
-      status=NF_GET_VARA_REAL(ncidp,varidyu1,start,count,pyu12)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidyu1,start,count,pyu12)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidyu1,start,count,pyu12)
+#endif
       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)
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidyv1,start,count,pyv12)
+#else
+        status=NF_GET_VARA_REAL(ncidp,varidyv1,start,count,pyv12)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,pyv12,pyv1)
 
+C**** Temerature au sol ********************************************
 c ftsol1
         print*,'LECTURE de ftsol1 a irec =',irec
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts1,start,count,ftsol12)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidfts1,start,count,ftsol12)
+#endif
        call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol12,ftsol1)
 
 c ftsol2
         print*,'LECTURE de ftsol2 a irec =',irec
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts2,start,count,ftsol22)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidfts2,start,count,ftsol22)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol22,ftsol2)
 
 c ftsol3
 	 print*,'LECTURE de ftsol3 a irec =',irec
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts3,start,count,ftsol32)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidfts3,start,count,ftsol32)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol32,ftsol3)
 
 c ftsol4
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidfts4,start,count,ftsol42)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidfts4,start,count,ftsol42)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol42,ftsol4)
 
+C**** Nature sol ********************************************
 c psrf1 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr1,start,count,psrf12)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidpsr1,start,count,psrf12)
+#endif
 c      call dump2d(iip1-1,jjm+1,psrf12,'PSRF1NC')
       call gr_ecrit_fi(1,kon,iim,jjm+1,psrf12,psrf1)
 
 c psrf2 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr2,start,count,psrf22)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidpsr2,start,count,psrf22)
+#endif
 c      call dump2d(iip1-1,jjm+1,psrf22,'PSRF2NC')
       call gr_ecrit_fi(1,kon,iim,jjm+1,psrf22,psrf2)
 
 c psrf3 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr3,start,count,psrf32)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidpsr3,start,count,psrf32)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,psrf32,psrf3)
 
 c psrf4 
+#ifdef NC_DOUBLE
+      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr4,start,count,psrf42)
+#else
       status=NF_GET_VARA_REAL(ncidp,varidpsr4,start,count,psrf42)
+#endif
       call gr_ecrit_fi(1,kon,iim,jjm+1,psrf42,psrf4)
 	
@@ -344,4 +492,6 @@
         psrf(i,2) = psrf2(i)
         psrf(i,3) = psrf3(i)
+c test abderr
+c	print*,'Dans read_pstoke psrf3 =',psrf3(i),i
         psrf(i,4) = psrf4(i)
   
Index: /LMDZ4/trunk/libf/phylmd/ustarhb.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/ustarhb.F	(revision 541)
+++ /LMDZ4/trunk/libf/phylmd/ustarhb.F	(revision 541)
@@ -0,0 +1,53 @@
+!
+! $Header$
+!
+      SUBROUTINE ustarhb(knon,u,v,cd_m, ustar)
+      IMPLICIT none
+c======================================================================
+c Laurent Li (LMD/CNRS), le 30 septembre 1998
+c Couche limite non-locale. Adaptation du code du CCM3.
+c Code non teste, donc a ne pas utiliser.
+c======================================================================
+c Nonlocal scheme that determines eddy diffusivities based on a
+c diagnosed boundary layer height and a turbulent velocity scale.
+c Also countergradient effects for heat and moisture are included.
+c
+c For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
+c Local versus nonlocal boundary-layer diffusion in a global climate
+c model. J. of Climate, vol. 6, 1825-1842.
+c======================================================================
+#include "dimensions.h"
+#include "dimphy.h"
+#include "YOMCST.h"
+c
+c Arguments:
+c
+      INTEGER knon ! nombre de points a calculer
+      REAL u(klon,klev) ! vitesse U (m/s)
+      REAL v(klon,klev) ! vitesse V (m/s)
+      REAL cd_m(klon) ! coefficient de friction au sol pour vitesse
+      REAL ustar(klon)
+c
+      INTEGER i, k
+      REAL zxt, zxq, zxu, zxv, zxmod, taux, tauy
+      REAL zx_alf1, zx_alf2 ! parametres pour extrapolation
+      LOGICAL unssrf(klon)  ! unstb pbl w/lvls within srf pbl lyr
+      LOGICAL unsout(klon)  ! unstb pbl w/lvls in outer pbl lyr
+      LOGICAL check(klon)   ! True=>chk if Richardson no.>critcal
+c
+#include "YOETHF.h"
+#include "FCTTRE.h"
+      DO i = 1, knon
+        zx_alf1 = 1.0
+        zx_alf2 = 1.0 - zx_alf1
+        zxu = u(i,1)*zx_alf1+u(i,2)*zx_alf2
+        zxv = v(i,1)*zx_alf1+v(i,2)*zx_alf2
+        zxmod = 1.0+SQRT(zxu**2+zxv**2)
+        taux = zxu *zxmod*cd_m(i)
+        tauy = zxv *zxmod*cd_m(i)
+        ustar(i) = SQRT(taux**2+tauy**2)
+c       print*,'Ust ',zxu,zxmod,taux,ustar(i)
+      ENDDO
+c
+      return
+      end
Index: /LMDZ4/trunk/libf/phylmd/vdif_kcay.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/vdif_kcay.F	(revision 541)
+++ /LMDZ4/trunk/libf/phylmd/vdif_kcay.F	(revision 541)
@@ -0,0 +1,739 @@
+!
+! $Header$
+!
+      SUBROUTINE vdif_kcay(ngrid,dt,g,rconst,plev,temp
+     s   ,zlev,zlay,u,v,teta,cd,q2,q2diag,km,kn,ustar
+     s   ,l_mix)
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+#include "dimphy.h"
+c.......................................................................
+c
+c dt : pas de temps
+c g  : g
+c zlev : altitude a chaque niveau (interface inferieure de la couche
+c        de meme indice)
+c zlay : altitude au centre de chaque couche
+c u,v : vitesse au centre de chaque couche
+c       (en entree : la valeur au debut du pas de temps)
+c teta : temperature potentielle au centre de chaque couche
+c        (en entree : la valeur au debut du pas de temps)
+c cd : cdrag
+c      (en entree : la valeur au debut du pas de temps)
+c q2 : $q^2$ au bas de chaque couche
+c      (en entree : la valeur au debut du pas de temps)
+c      (en sortie : la valeur a la fin du pas de temps)
+c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
+c      couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c 
+c.......................................................................
+      REAL dt,g,rconst
+      real plev(klon,klev+1),temp(klon,klev)
+      real ustar(klon),snstable
+      REAL zlev(klon,klev+1)
+      REAL zlay(klon,klev)
+      REAL u(klon,klev)
+      REAL v(klon,klev)
+      REAL teta(klon,klev)
+      REAL cd(klon)
+      REAL q2(klon,klev+1),q2s(klon,klev+1)
+      REAL q2diag(klon,klev+1)
+      REAL km(klon,klev+1)
+      REAL kn(klon,klev+1)
+      real sq(klon),sqz(klon),zz(klon,klev+1),zq,long0(klon)
+
+      integer l_mix,iii
+c.......................................................................
+c
+c nlay : nombre de couches        
+c nlev : nombre de niveaux
+c ngrid : nombre de points de grille       
+c unsdz : 1 sur l'epaisseur de couche
+c unsdzdec : 1 sur la distance entre le centre de la couche et le
+c            centre de la couche inferieure
+c q : echelle de vitesse au bas de chaque couche
+c     (valeur a la fin du pas de temps)
+c
+c.......................................................................
+      INTEGER nlay,nlev,ngrid
+      REAL unsdz(klon,klev)
+      REAL unsdzdec(klon,klev+1)
+      REAL q(klon,klev+1)
+
+c.......................................................................
+c
+c kmpre : km au debut du pas de temps
+c qcstat : q : solution stationnaire du probleme couple
+c          (valeur a la fin du pas de temps)
+c q2cstat : q2 : solution stationnaire du probleme couple
+c           (valeur a la fin du pas de temps)
+c
+c.......................................................................
+      REAL kmpre(klon,klev+1)
+      REAL qcstat
+      REAL q2cstat
+      real sss,sssq
+c.......................................................................
+c
+c long : longueur de melange calculee selon Blackadar
+c
+c.......................................................................
+      REAL long(klon,klev+1)
+c.......................................................................
+c
+c kmq3 : terme en q^3 dans le developpement de km
+c        (valeur au debut du pas de temps)
+c kmcstat : valeur de km solution stationnaire du systeme {q2 ; du/dz}
+c           (valeur a la fin du pas de temps)
+c knq3 : terme en q^3 dans le developpement de kn
+c mcstat : valeur de m solution stationnaire du systeme {q2 ; du/dz}
+c          (valeur a la fin du pas de temps)
+c m2cstat : valeur de m2 solution stationnaire du systeme {q2 ; du/dz}
+c           (valeur a la fin du pas de temps)
+c m : valeur a la fin du pas de temps
+c mpre : valeur au debut du pas de temps
+c m2 : valeur a la fin du pas de temps
+c n2 : valeur a la fin du pas de temps
+c 
+c.......................................................................
+      REAL kmq3
+      REAL kmcstat
+      REAL knq3
+      REAL mcstat
+      REAL m2cstat
+      REAL m(klon,klev+1)
+      REAL mpre(klon,klev+1)
+      REAL m2(klon,klev+1)
+      REAL n2(klon,klev+1)
+c.......................................................................
+c
+c gn : intermediaire pour les coefficients de stabilite
+c gnmin : borne inferieure de gn (-0.23 ou -0.28)
+c gnmax : borne superieure de gn (0.0233)
+c gninf : vrai si gn est en dessous de sa borne inferieure
+c gnsup : vrai si gn est en dessus de sa borne superieure
+c gm : drole d'objet bien utile
+c ri : nombre de Richardson
+c sn : coefficient de stabilite pour n
+c snq2 : premier terme du developement limite de sn en q2
+c sm : coefficient de stabilite pour m
+c smq2 : premier terme du developement limite de sm en q2
+c
+c.......................................................................
+      REAL gn
+      REAL gnmin
+      REAL gnmax
+      LOGICAL gninf
+      LOGICAL gnsup
+      REAL gm
+c      REAL ri(klon,klev+1)
+      REAL sn(klon,klev+1)
+      REAL snq2(klon,klev+1)
+      REAL sm(klon,klev+1)
+      REAL smq2(klon,klev+1)
+c.......................................................................
+c
+c kappa : consatnte de Von Karman (0.4)
+c long00 : longueur de reference pour le calcul de long (160)
+c a1,a2,b1,b2,c1 : constantes d'origine pour les  coefficients
+c                  de stabilite (0.92/0.74/16.6/10.1/0.08)
+c cn1,cn2 : constantes pour sn
+c cm1,cm2,cm3,cm4 : constantes pour sm
+c
+c.......................................................................
+      REAL kappa
+      REAL long00
+      REAL a1,a2,b1,b2,c1
+      REAL cn1,cn2
+      REAL cm1,cm2,cm3,cm4
+c.......................................................................
+c
+c termq : termes en $q$ dans l'equation de q2
+c termq3 : termes en $q^3$ dans l'equation de q2
+c termqm2 : termes en $q*m^2$ dans l'equation de q2
+c termq3m2 : termes en $q^3*m^2$ dans l'equation de q2
+c
+c.......................................................................
+      REAL termq
+      REAL termq3
+      REAL termqm2
+      REAL termq3m2
+c.......................................................................
+c
+c q2min : borne inferieure de q2
+c q2max : borne superieure de q2
+c
+c.......................................................................
+      REAL q2min
+      REAL q2max
+c.......................................................................
+c knmin : borne inferieure de kn
+c kmmin : borne inferieure de km
+c.......................................................................
+      REAL knmin
+      REAL kmmin
+c.......................................................................
+      INTEGER ilay,ilev,igrid
+      REAL tmp1,tmp2
+c.......................................................................
+      PARAMETER (kappa=0.4E+0)
+      PARAMETER (long00=160.E+0)
+c     PARAMETER (gnmin=-10.E+0)
+      PARAMETER (gnmin=-0.28)
+      PARAMETER (gnmax=0.0233E+0)
+      PARAMETER (a1=0.92E+0)
+      PARAMETER (a2=0.74E+0)
+      PARAMETER (b1=16.6E+0)
+      PARAMETER (b2=10.1E+0)
+      PARAMETER (c1=0.08E+0)
+      PARAMETER (knmin=1.E-5)
+      PARAMETER (kmmin=1.E-5)
+      PARAMETER (q2min=1.e-5)
+      PARAMETER (q2max=1.E+2)
+      PARAMETER (nlay=klev)
+      PARAMETER (nlev=klev+1)
+c
+      PARAMETER (
+     &  cn1=a2*(1.E+0 -6.E+0 *a1/b1)
+     &          )
+      PARAMETER (
+     &  cn2=-3.E+0 *a2*(6.E+0 *a1+b2)
+     &          )
+      PARAMETER (
+     &  cm1=a1*(1.E+0 -3.E+0 *c1-6.E+0 *a1/b1)
+     &          )
+      PARAMETER (
+     &  cm2=a1*(-3.E+0 *a2*((b2-3.E+0 *a2)*(1.E+0 -6.E+0 *a1/b1)
+     &          -3.E+0 *c1*(b2+6.E+0 *a1)))
+     &          )
+      PARAMETER (
+     &  cm3=-3.E+0 *a2*(6.E+0 *a1+b2)
+     &          )
+      PARAMETER (
+     &  cm4=-9.E+0 *a1*a2
+     &          )
+
+      logical first
+      save first
+      data first/.true./
+c.......................................................................
+c  traitment des valeur de q2 en entree
+c.......................................................................
+c
+c   Initialisation de q2
+
+      call yamada(ngrid,dt,g,rconst,plev,temp
+     s   ,zlev,zlay,u,v,teta,cd,q2diag,km,kn,ustar
+     s   ,l_mix)
+      if (first.and.1.eq.1) then
+      first=.false.
+      q2=q2diag
+      endif
+
+      DO ilev=1,nlev
+                                                      DO igrid=1,ngrid 
+        q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
+        q(igrid,ilev)=sqrt(q2(igrid,ilev))
+                                                      ENDDO
+      ENDDO
+c
+                                                      DO igrid=1,ngrid 
+      tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
+      q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
+      q2(igrid,1)=amax1(q2(igrid,1),q2min)
+      q(igrid,1)=sqrt(q2(igrid,1))
+                                                      ENDDO
+c
+c.......................................................................
+c  les increments verticaux
+c.......................................................................
+c
+c!!!!! allerte !!!!!c
+c!!!!! zlev n'est pas declare a nlev !!!!!c
+c!!!!! ---->
+                                                      DO igrid=1,ngrid 
+            zlev(igrid,nlev)=zlay(igrid,nlay)
+     &             +( zlay(igrid,nlay) - zlev(igrid,nlev-1) )
+                                                      ENDDO            
+c!!!!! <----
+c!!!!! allerte !!!!!c
+c
+      DO ilay=1,nlay
+                                                      DO igrid=1,ngrid 
+        unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
+                                                      ENDDO
+      ENDDO
+                                                      DO igrid=1,ngrid 
+      unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
+                                                      ENDDO
+      DO ilay=2,nlay
+                                                      DO igrid=1,ngrid 
+        unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
+                                                      ENDDO
+      ENDDO
+                                                      DO igrid=1,ngrid 
+      unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
+                                                      ENDDO
+c
+c.......................................................................
+c  le cisaillement et le gradient de temperature
+c.......................................................................
+c
+                                                      DO igrid=1,ngrid 
+      m2(igrid,1)=(unsdzdec(igrid,1)
+     &                   *u(igrid,1))**2
+     &                 +(unsdzdec(igrid,1)
+     &                   *v(igrid,1))**2
+      m(igrid,1)=sqrt(m2(igrid,1))
+      mpre(igrid,1)=m(igrid,1)
+                                                      ENDDO
+c
+c-----------------------------------------------------------------------
+      DO ilev=2,nlev-1
+                                                      DO igrid=1,ngrid 
+c-----------------------------------------------------------------------
+c
+        n2(igrid,ilev)=g*unsdzdec(igrid,ilev)
+     &                   *(teta(igrid,ilev)-teta(igrid,ilev-1))
+     &                   /(teta(igrid,ilev)+teta(igrid,ilev-1)) *2.E+0
+c       n2(igrid,ilev)=0.
+c
+c --->
+c       on ne sais traiter que les cas stratifies. et l'ajustement
+c       convectif est cense faire en sorte que seul des configurations
+c       stratifiees soient rencontrees en entree de cette routine.
+c       mais, bon ... on sait jamais (meme on sait que n2 prends
+c       quelques valeurs negatives ... parfois) alors : 
+c<---
+c
+        IF (n2(igrid,ilev).lt.0.E+0) THEN
+          n2(igrid,ilev)=0.E+0
+        ENDIF
+c
+        m2(igrid,ilev)=(unsdzdec(igrid,ilev)
+     &                     *(u(igrid,ilev)-u(igrid,ilev-1)))**2
+     &                   +(unsdzdec(igrid,ilev)
+     &                     *(v(igrid,ilev)-v(igrid,ilev-1)))**2
+        m(igrid,ilev)=sqrt(m2(igrid,ilev))
+        mpre(igrid,ilev)=m(igrid,ilev)
+c
+c-----------------------------------------------------------------------
+                                                      ENDDO
+      ENDDO
+c-----------------------------------------------------------------------
+c
+                                                      DO igrid=1,ngrid 
+      m2(igrid,nlev)=m2(igrid,nlev-1)
+      m(igrid,nlev)=m(igrid,nlev-1)
+      mpre(igrid,nlev)=m(igrid,nlev)
+                                                      ENDDO
+c
+c.......................................................................
+c  calcul des fonctions de stabilite
+c.......................................................................
+c
+      if (l_mix.eq.4) then
+                                                      DO igrid=1,ngrid 
+         sqz(igrid)=1.e-10
+         sq(igrid)=1.e-10
+                                                      ENDDO
+         do ilev=2,nlev-1
+                                                      DO igrid=1,ngrid 
+           zq=sqrt(q2(igrid,ilev))
+           sqz(igrid)
+     .     =sqz(igrid)+zq*zlev(igrid,ilev)
+     .     *(zlay(igrid,ilev)-zlay(igrid,ilev-1))
+           sq(igrid)=sq(igrid)+zq*(zlay(igrid,ilev)-zlay(igrid,ilev-1))
+                                                      ENDDO
+         enddo
+                                                      DO igrid=1,ngrid 
+         long0(igrid)=0.2*sqz(igrid)/sq(igrid)
+                                                      ENDDO
+      else if (l_mix.eq.3) then
+         long0(igrid)=long00
+      endif
+
+c (abd 5 2)      print*,'LONG0=',long0
+
+c-----------------------------------------------------------------------
+      DO ilev=2,nlev-1
+                                                      DO igrid=1,ngrid 
+c-----------------------------------------------------------------------
+c
+        tmp1=kappa*(zlev(igrid,ilev)-zlev(igrid,1))
+        if (l_mix.ge.10) then
+            long(igrid,ilev)=l_mix
+        else
+           long(igrid,ilev)=tmp1/(1.E+0 + tmp1/long0(igrid))
+        endif
+        long(igrid,ilev)=max(min(long(igrid,ilev)
+     s    ,0.5*sqrt(q2(igrid,ilev))/sqrt(max(n2(igrid,ilev),1.e-10)))
+     s    ,5.)
+
+        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
+     &                                           * n2(igrid,ilev)
+        gm=long(igrid,ilev)**2 / q2(igrid,ilev)
+     &                                           * m2(igrid,ilev)
+c
+        gninf=.false.
+        gnsup=.false.
+        long(igrid,ilev)=long(igrid,ilev)
+        long(igrid,ilev)=long(igrid,ilev)
+c
+        IF (gn.lt.gnmin) THEN
+          gninf=.true.
+          gn=gnmin
+        ENDIF
+c
+        IF (gn.gt.gnmax) THEN
+          gnsup=.true.
+          gn=gnmax
+        ENDIF
+c
+        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
+        sm(igrid,ilev)=
+     &    (cm1+cm2*gn)
+     &   /( (1.E+0 +cm3*gn)
+     &     *(1.E+0 +cm4*gn) )
+c
+        IF ((gninf).or.(gnsup)) THEN
+          snq2(igrid,ilev)=0.E+0
+          smq2(igrid,ilev)=0.E+0
+        ELSE
+          snq2(igrid,ilev)=
+     &     -gn
+     &     *(-cn1*cn2/(1.E+0 +cn2*gn)**2 )
+          smq2(igrid,ilev)=
+     &     -gn
+     &     *( cm2*(1.E+0 +cm3*gn)
+     &           *(1.E+0 +cm4*gn)
+     &       -( cm3*(1.E+0 +cm4*gn)
+     &         +cm4*(1.E+0 +cm3*gn) )
+     &       *(cm1+cm2*gn)            )
+     &     /( (1.E+0 +cm3*gn)
+     &       *(1.E+0 +cm4*gn) )**2
+        ENDIF
+c
+c abd
+c        if(ilev.le.57.and.ilev.ge.37) then
+c            print*,'L=',ilev,'   GN=',gn,'  SM=',sm(igrid,ilev)
+c        endif
+c --->
+c       la decomposition de Taylor en q2 n'a de sens que
+c       dans les cas stratifies ou sn et sm sont quasi
+c       proportionnels a q2. ailleurs on laisse le meme
+c       algorithme car l'ajustement convectif fait le travail.
+c       mais c'est delirant quand sn et snq2 n'ont pas le meme
+c       signe : dans ces cas, on ne fait pas la decomposition.
+c<---
+c
+        IF (snq2(igrid,ilev)*sn(igrid,ilev).le.0.E+0)
+     &      snq2(igrid,ilev)=0.E+0
+        IF (smq2(igrid,ilev)*sm(igrid,ilev).le.0.E+0)
+     &      smq2(igrid,ilev)=0.E+0
+c
+C   Correction pour les couches stables.
+C   Schema repris de JHoltzlag Boville, lui meme venant de...
+
+        if (1.eq.1) then
+        snstable=1.-zlev(igrid,ilev)
+     s     /(700.*max(ustar(igrid),0.0001))
+        snstable=1.-zlev(igrid,ilev)/400.
+        snstable=max(snstable,0.)
+        snstable=snstable*snstable
+
+c abde       print*,'SN ',ilev,sn(1,ilev),snstable
+        if (sn(igrid,ilev).lt.snstable) then
+           sn(igrid,ilev)=snstable
+           snq2(igrid,ilev)=0.
+        endif
+
+        if (sm(igrid,ilev).lt.snstable) then
+           sm(igrid,ilev)=snstable
+           smq2(igrid,ilev)=0.
+        endif
+
+        endif
+
+c sn : coefficient de stabilite pour n
+c snq2 : premier terme du developement limite de sn en q2
+c-----------------------------------------------------------------------
+                                                      ENDDO
+      ENDDO
+c-----------------------------------------------------------------------
+c
+c.......................................................................
+c  calcul de km et kn au debut du pas de temps
+c.......................................................................
+c
+                                                      DO igrid=1,ngrid 
+      kn(igrid,1)=knmin
+      km(igrid,1)=kmmin
+      kmpre(igrid,1)=km(igrid,1)
+                                                      ENDDO
+c
+c-----------------------------------------------------------------------
+      DO ilev=2,nlev-1
+                                                      DO igrid=1,ngrid 
+c-----------------------------------------------------------------------
+c
+        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
+     &                                         *sn(igrid,ilev)
+        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
+     &                                         *sm(igrid,ilev)
+        kmpre(igrid,ilev)=km(igrid,ilev)
+c
+c-----------------------------------------------------------------------
+                                                      ENDDO
+      ENDDO
+c-----------------------------------------------------------------------
+c
+                                                      DO igrid=1,ngrid 
+      kn(igrid,nlev)=kn(igrid,nlev-1)
+      km(igrid,nlev)=km(igrid,nlev-1)
+      kmpre(igrid,nlev)=km(igrid,nlev)
+                                                      ENDDO
+c
+c.......................................................................
+c  boucle sur les niveaux 2 a nlev-1
+c.......................................................................
+c
+c---->
+      DO 10001 ilev=2,nlev-1
+c---->
+      DO 10002 igrid=1,ngrid 
+c
+c.......................................................................
+c
+c  calcul des termes sources et puits de l'equation de q2
+c  ------------------------------------------------------
+c
+        knq3=kn(igrid,ilev)*snq2(igrid,ilev)
+     &                                    /sn(igrid,ilev)
+        kmq3=km(igrid,ilev)*smq2(igrid,ilev)
+     &                                    /sm(igrid,ilev)
+c
+        termq=0.E+0
+        termq3=0.E+0
+        termqm2=0.E+0
+        termq3m2=0.E+0
+c
+        tmp1=dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
+        tmp2=dt*2.E+0 *kmq3*m2(igrid,ilev)
+        termqm2=termqm2
+     &    +dt*2.E+0 *km(igrid,ilev)*m2(igrid,ilev)
+     &    -dt*2.E+0 *kmq3*m2(igrid,ilev)
+        termq3m2=termq3m2
+     &    +dt*2.E+0 *kmq3*m2(igrid,ilev)
+c 
+        termq=termq
+     &    -dt*2.E+0 *kn(igrid,ilev)*n2(igrid,ilev)
+     &    +dt*2.E+0 *knq3*n2(igrid,ilev)
+        termq3=termq3
+     &    -dt*2.E+0 *knq3*n2(igrid,ilev)
+c
+        termq3=termq3
+     &    -dt*2.E+0 *q(igrid,ilev)**3 / (b1*long(igrid,ilev))
+c
+c.......................................................................
+c
+c  resolution stationnaire couplee avec le gradient de vitesse local
+c  -----------------------------------------------------------------
+c
+c  -----{on cherche le cisaillement qui annule l'equation de q^2
+c        supposee en q3}
+c
+        tmp1=termq+termq3
+        tmp2=termqm2+termq3m2
+        m2cstat=m2(igrid,ilev)
+     &      -(tmp1+tmp2)/(dt*2.E+0*km(igrid,ilev))
+        mcstat=sqrt(m2cstat)
+
+c  abde      print*,'M2 L=',ilev,mpre(igrid,ilev),mcstat
+c
+c  -----{puis on ecrit la valeur de q qui annule l'equation de m
+c        supposee en q3}
+c
+        IF (ilev.eq.2) THEN
+          kmcstat=1.E+0 / mcstat
+     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
+     &                        *mpre(igrid,ilev+1)
+     &      +unsdz(igrid,ilev-1)
+     &              *cd(igrid)
+     &              *( sqrt(u(igrid,3)**2+v(igrid,3)**2)
+     &                -mcstat/unsdzdec(igrid,ilev)
+     &                -mpre(igrid,ilev+1)/unsdzdec(igrid,ilev+1) )**2)
+     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
+        ELSE
+          kmcstat=1.E+0 / mcstat
+     &    *( unsdz(igrid,ilev)*kmpre(igrid,ilev+1)
+     &                        *mpre(igrid,ilev+1)
+     &      +unsdz(igrid,ilev-1)*kmpre(igrid,ilev-1)
+     &                          *mpre(igrid,ilev-1) )
+     &      /( unsdz(igrid,ilev)+unsdz(igrid,ilev-1) )
+        ENDIF
+        tmp2=kmcstat
+     &      /( sm(igrid,ilev)/q2(igrid,ilev) )
+     &      /long(igrid,ilev)
+        qcstat=tmp2**(1.E+0/3.E+0)
+        q2cstat=qcstat**2
+c
+c.......................................................................
+c
+c  choix de la solution finale
+c  ---------------------------
+c
+          q(igrid,ilev)=qcstat
+          q2(igrid,ilev)=q2cstat
+          m(igrid,ilev)=mcstat
+c abd       if(ilev.le.57.and.ilev.ge.37) then
+c           print*,'L=',ilev,'   M2=',m2(igrid,ilev),m2cstat,
+c     s     'N2=',n2(igrid,ilev)
+c abd       endif
+          m2(igrid,ilev)=m2cstat
+c
+c --->
+c       pour des raisons simples q2 est minore 
+c<---
+c
+        IF (q2(igrid,ilev).lt.q2min) THEN
+          q2(igrid,ilev)=q2min
+          q(igrid,ilev)=sqrt(q2min)
+        ENDIF
+c
+c.......................................................................
+c
+c  calcul final de kn et km
+c  ------------------------
+c
+        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
+     &                                           * n2(igrid,ilev)
+        IF (gn.lt.gnmin) gn=gnmin
+        IF (gn.gt.gnmax) gn=gnmax
+        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
+        sm(igrid,ilev)=
+     &    (cm1+cm2*gn)
+     &   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
+        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
+     &                 *sn(igrid,ilev)
+        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
+     &                 *sm(igrid,ilev)
+c abd
+c        if(ilev.le.57.and.ilev.ge.37) then
+c            print*,'L=',ilev,'   GN=',gn,'  SM=',sm(igrid,ilev)
+c        endif
+c
+c.......................................................................
+c
+10002 CONTINUE
+c
+10001 CONTINUE
+c
+c.......................................................................
+c
+c
+                                                      DO igrid=1,ngrid 
+      kn(igrid,1)=knmin
+      km(igrid,1)=kmmin
+c     kn(igrid,1)=cd(igrid)
+c     km(igrid,1)=cd(igrid)
+      q2(igrid,nlev)=q2(igrid,nlev-1)
+      q(igrid,nlev)=q(igrid,nlev-1)
+      kn(igrid,nlev)=kn(igrid,nlev-1)
+      km(igrid,nlev)=km(igrid,nlev-1)
+                                                      ENDDO
+c
+c  CALCUL DE LA DIFFUSION VERTICALE DE Q2
+      if (1.eq.1) then
+
+        do ilev=2,klev-1
+           sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
+           sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
+        enddo
+c        print*,'Q2moy avant',sssq/sss
+c       print*,'Q2q20 ',(q2(1,ilev),ilev=1,10)
+c       print*,'Q2km0 ',(km(1,ilev),ilev=1,10)
+c   ! C'est quoi ca qu'etait dans l'original???
+c       do igrid=1,ngrid
+c          q2(igrid,1)=10.
+c       enddo
+c        q2s=q2
+c       do iii=1,10
+c       call vdif_q2(dt,g,rconst,plev,temp,km,q2)
+c       do ilev=1,klev+1
+c          write(iii+49,*) q2(1,ilev),zlev(1,ilev)
+c       enddo
+c       enddo
+c       stop
+c       do ilev=1,klev
+c          print*,zlev(1,ilev),q2s(1,ilev),q2(1,ilev)
+c       enddo
+c        q2s=q2-q2s
+c       do ilev=1,klev
+c          print*,q2s(1,ilev),zlev(1,ilev)
+c       enddo
+        do ilev=2,klev-1
+           sss=sss+plev(1,ilev-1)-plev(1,ilev+1)
+           sssq=sssq+(plev(1,ilev-1)-plev(1,ilev+1))*q2(1,ilev)
+        enddo
+        print*,'Q2moy apres',sssq/sss
+c
+c
+        do ilev=1,nlev
+           do igrid=1,ngrid
+              q2(igrid,ilev)=max(q2(igrid,ilev),q2min)
+              q(igrid,ilev)=sqrt(q2(igrid,ilev))
+
+c.......................................................................
+c
+c  calcul final de kn et km
+c  ------------------------
+c
+        gn=-long(igrid,ilev)**2 / q2(igrid,ilev)
+     &                                           * n2(igrid,ilev)
+        IF (gn.lt.gnmin) gn=gnmin
+        IF (gn.gt.gnmax) gn=gnmax
+        sn(igrid,ilev)=cn1/(1.E+0 +cn2*gn)
+        sm(igrid,ilev)=
+     &    (cm1+cm2*gn)
+     &   /( (1.E+0 +cm3*gn)*(1.E+0 +cm4*gn) )
+C   Correction pour les couches stables.
+C   Schema repris de JHoltzlag Boville, lui meme venant de...
+
+        if (1.eq.1) then
+        snstable=1.-zlev(igrid,ilev)
+     s     /(700.*max(ustar(igrid),0.0001))
+        snstable=1.-zlev(igrid,ilev)/400.
+        snstable=max(snstable,0.)
+        snstable=snstable*snstable
+
+c abde      print*,'SN ',ilev,sn(1,ilev),snstable
+        if (sn(igrid,ilev).lt.snstable) then
+           sn(igrid,ilev)=snstable
+           snq2(igrid,ilev)=0.
+        endif
+
+        if (sm(igrid,ilev).lt.snstable) then
+           sm(igrid,ilev)=snstable
+           smq2(igrid,ilev)=0.
+        endif
+
+        endif
+
+c sn : coefficient de stabilite pour n
+        kn(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
+     &                 *sn(igrid,ilev)
+        km(igrid,ilev)=long(igrid,ilev)*q(igrid,ilev)
+c
+           enddo
+        enddo
+c       print*,'Q2km1 ',(km(1,ilev),ilev=1,10)
+
+      endif
+
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/phylmd/write_histrac.h
===================================================================
--- /LMDZ4/trunk/libf/phylmd/write_histrac.h	(revision 540)
+++ /LMDZ4/trunk/libf/phylmd/write_histrac.h	(revision 541)
@@ -234,4 +234,17 @@
      .                                   iim*(jjm+1)*klev,ndex3d)
       endif
+      
+c----Olivia
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,d_tr_th(1,1,it),zx_tmp_3d)
+       CALL histwrite(nid_tra,"d_tr_th_"//tnom(it+2),itau_w,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,d_tr_cvi(1,1,it),zx_tmp_3d)
+       CALL histwrite(nid_tra,"d_tr_cv_"//tnom(it+2),itau_w,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,d_tr_cli(1,1,it),zx_tmp_3d)
+       CALL histwrite(nid_tra,"d_tr_cl_"//tnom(it+2),itau_w,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex3d)
+c---fin Olivia      
+      
 #endif
       ENDDO
Index: /LMDZ4/trunk/libf/phylmd/yamada.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/yamada.F	(revision 541)
+++ /LMDZ4/trunk/libf/phylmd/yamada.F	(revision 541)
@@ -0,0 +1,169 @@
+!
+! $Header$
+!
+      SUBROUTINE yamada(ngrid,dt,g,rconst,plev,temp
+     s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,ustar
+     s   ,l_mix)
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+#include "dimphy.h"
+c.......................................................................
+c
+c dt : pas de temps
+c g  : g
+c zlev : altitude a chaque niveau (interface inferieure de la couche
+c        de meme indice)
+c zlay : altitude au centre de chaque couche
+c u,v : vitesse au centre de chaque couche
+c       (en entree : la valeur au debut du pas de temps)
+c teta : temperature potentielle au centre de chaque couche
+c        (en entree : la valeur au debut du pas de temps)
+c cd : cdrag
+c      (en entree : la valeur au debut du pas de temps)
+c q2 : $q^2$ au bas de chaque couche
+c      (en entree : la valeur au debut du pas de temps)
+c      (en sortie : la valeur a la fin du pas de temps)
+c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
+c      couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c
+c.......................................................................
+      REAL dt,g,rconst
+      real plev(klon,klev+1),temp(klon,klev)
+      real ustar(klon),snstable
+      REAL zlev(klon,klev+1)
+      REAL zlay(klon,klev)
+      REAL u(klon,klev)
+      REAL v(klon,klev)
+      REAL teta(klon,klev)
+      REAL cd(klon)
+      REAL q2(klon,klev+1)
+      REAL km(klon,klev+1)
+      REAL kn(klon,klev+1)
+      integer l_mix,ngrid
+
+
+      integer nlay,nlev
+      PARAMETER (nlay=klev)
+      PARAMETER (nlev=klev+1)
+
+      logical first
+      save first
+      data first/.true./
+
+
+      integer ig,k
+
+      real ri,zrif,zalpha,zsm
+      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
+
+      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
+      real l(klon,klev+1),l0(klon)
+
+      real sq(klon),sqz(klon),zz(klon,klev+1)
+      integer iter
+
+      real ric,rifc,b1,kap
+      save ric,rifc,b1,kap
+      data ric,rifc,b1,kap/0.195,0.191,16.6,0.3/
+
+      real frif,falpha,fsm
+
+      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
+      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
+      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
+
+      if (0.eq.1.and.first) then
+      do ig=1,1000
+         ri=(ig-800.)/500.
+         if (ri.lt.ric) then
+            zrif=frif(ri)
+         else
+            zrif=rifc
+         endif
+         if(zrif.lt.0.16) then
+            zalpha=falpha(zrif)
+            zsm=fsm(zrif)
+         else
+            zalpha=1.12
+            zsm=0.085
+         endif
+         print*,ri,rif,zalpha,zsm
+      enddo
+      first=.false.
+      endif
+
+c  Correction d'un bug sauvage a verifier.
+c      do k=2,nlev
+      do k=2,nlay
+                                                          do ig=1,ngrid
+         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
+         m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
+     s             /(dz(ig,k)*dz(ig,k))
+         n2(ig,k)=g*2.*(teta(ig,k)-teta(ig,k-1))
+     s            /(teta(ig,k-1)+teta(ig,k))  /dz(ig,k)
+         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
+         if (ri.lt.ric) then
+            rif(ig,k)=frif(ri)
+         else
+            rif(ig,k)=rifc
+         endif
+         if(rif(ig,k).lt.0.16) then
+            alpha(ig,k)=falpha(rif(ig,k))
+            sm(ig,k)=fsm(rif(ig,k))
+         else
+            alpha(ig,k)=1.12
+            sm(ig,k)=0.085
+         endif
+         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
+                                                          enddo
+      enddo
+
+c iterration pour determiner la longueur de melange
+
+                                                          do ig=1,ngrid
+      l0(ig)=100.
+                                                          enddo
+      do k=2,klev-1
+                                                          do ig=1,ngrid
+        l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
+                                                          enddo
+      enddo
+
+      do iter=1,10
+                                                          do ig=1,ngrid
+         sq(ig)=1.e-10
+         sqz(ig)=1.e-10
+                                                          enddo
+         do k=2,klev-1
+                                                          do ig=1,ngrid
+           q2(ig,k)=l(ig,k)**2*zz(ig,k)
+           l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
+     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
+           zq=sqrt(q2(ig,k))
+           sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
+           sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
+                                                          enddo
+         enddo
+                                                          do ig=1,ngrid
+         l0(ig)=0.2*sqz(ig)/sq(ig)
+                                                          enddo
+c(abd 3 5 2)         print*,'ITER=',iter,'  L0=',l0
+
+      enddo
+
+      do k=2,klev
+                                                          do ig=1,ngrid
+         l(ig,k)=min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
+     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10)))
+         q2(ig,k)=l(ig,k)**2*zz(ig,k)
+         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
+         kn(ig,k)=km(ig,k)*alpha(ig,k)
+                                                          enddo
+      enddo
+
+      return
+      end
Index: /LMDZ4/trunk/libf/phylmd/yamada4.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/yamada4.F	(revision 541)
+++ /LMDZ4/trunk/libf/phylmd/yamada4.F	(revision 541)
@@ -0,0 +1,462 @@
+!
+! $Header$
+!
+      SUBROUTINE yamada4(ngrid,dt,g,rconst,plev,temp
+     s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
+     s   ,iflag_pbl)
+      IMPLICIT NONE
+c.......................................................................
+#include "dimensions.h"
+#include "dimphy.h"
+c.......................................................................
+c
+c dt : pas de temps
+c g  : g
+c zlev : altitude a chaque niveau (interface inferieure de la couche
+c        de meme indice)
+c zlay : altitude au centre de chaque couche
+c u,v : vitesse au centre de chaque couche
+c       (en entree : la valeur au debut du pas de temps)
+c teta : temperature potentielle au centre de chaque couche
+c        (en entree : la valeur au debut du pas de temps)
+c cd : cdrag
+c      (en entree : la valeur au debut du pas de temps)
+c q2 : $q^2$ au bas de chaque couche
+c      (en entree : la valeur au debut du pas de temps)
+c      (en sortie : la valeur a la fin du pas de temps)
+c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
+c      couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
+c      (en sortie : la valeur a la fin du pas de temps)
+c
+c  iflag_pbl doit valoir entre 6 et 9
+c      l=6, on prend  systematiquement une longueur d'equilibre
+c    iflag_pbl=6 : MY 2.0
+c    iflag_pbl=7 : MY 2.0.Fournier
+c    iflag_pbl=8 : MY 2.5
+c    iflag_pbl=9 : un test ?
+
+c.......................................................................
+      REAL dt,g,rconst
+      real plev(klon,klev+1),temp(klon,klev)
+      real ustar(klon)
+      real kmin,qmin,pblhmin(klon),coriol(klon)
+      REAL zlev(klon,klev+1)
+      REAL zlay(klon,klev)
+      REAL u(klon,klev)
+      REAL v(klon,klev)
+      REAL teta(klon,klev)
+      REAL cd(klon)
+      REAL q2(klon,klev+1),qpre
+      REAL unsdz(klon,klev)
+      REAL unsdzdec(klon,klev+1)
+
+      REAL km(klon,klev+1)
+      REAL kmpre(klon,klev+1),tmp2
+      REAL mpre(klon,klev+1)
+      REAL kn(klon,klev+1)
+      REAL kq(klon,klev+1)
+      real ff(klon,klev+1),delta(klon,klev+1)
+      real aa(klon,klev+1),aa0,aa1
+      integer iflag_pbl,ngrid
+
+
+      integer nlay,nlev
+      PARAMETER (nlay=klev)
+      PARAMETER (nlev=klev+1)
+
+      logical first
+      integer ipas
+      save first,ipas
+      data first,ipas/.true.,0/
+
+
+      integer ig,k
+
+
+      real ri,zrif,zalpha,zsm,zsn
+      real rif(klon,klev+1),sm(klon,klev+1),alpha(klon,klev)
+
+      real m2(klon,klev+1),dz(klon,klev+1),zq,n2(klon,klev+1)
+      real dtetadz(klon,klev+1)
+      real m2cstat,mcstat,kmcstat
+      real l(klon,klev+1),l0(klon)
+      save l0
+
+      real sq(klon),sqz(klon),zz(klon,klev+1)
+      integer iter
+
+      real ric,rifc,b1,kap
+      save ric,rifc,b1,kap
+      data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
+
+      real frif,falpha,fsm
+      real fl,zzz,zl0,zq2,zn2
+
+      real rino(klon,klev+1),smyam(klon,klev),styam(klon,klev)
+     s  ,lyam(klon,klev),knyam(klon,klev)
+     s  ,w2yam(klon,klev),t2yam(klon,klev)
+      common/pbldiag/rino,smyam,styam,lyam,knyam,w2yam,t2yam
+
+      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
+      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
+      fsm(ri)=1.96*(0.1912-ri)*(0.2341-ri)/((1.-ri)*(0.2231-ri))
+      fl(zzz,zl0,zq2,zn2)=
+     s     max(min(l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
+     s     ,0.5*sqrt(q2(ig,k))/sqrt(max(n2(ig,k),1.e-10))) ,1.)
+
+      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.9)) then
+           stop'probleme de coherence dans appel a MY'
+      endif
+
+      ipas=ipas+1
+      if (0.eq.1.and.first) then
+      do ig=1,1000
+         ri=(ig-800.)/500.
+         if (ri.lt.ric) then
+            zrif=frif(ri)
+         else
+            zrif=rifc
+         endif
+         if(zrif.lt.0.16) then
+            zalpha=falpha(zrif)
+            zsm=fsm(zrif)
+         else
+            zalpha=1.12
+            zsm=0.085
+         endif
+c     print*,ri,rif,zalpha,zsm
+      enddo
+      endif
+
+c.......................................................................
+c  les increments verticaux
+c.......................................................................
+c
+c!!!!! allerte !!!!!c
+c!!!!! zlev n'est pas declare a nlev !!!!!c
+c!!!!! ---->
+                                                      DO ig=1,ngrid
+            zlev(ig,nlev)=zlay(ig,nlay)
+     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
+                                                      ENDDO
+c!!!!! <----
+c!!!!! allerte !!!!!c
+c
+      DO k=1,nlay
+                                                      DO ig=1,ngrid
+        unsdz(ig,k)=1.E+0/(zlev(ig,k+1)-zlev(ig,k))
+                                                      ENDDO
+      ENDDO
+                                                      DO ig=1,ngrid
+      unsdzdec(ig,1)=1.E+0/(zlay(ig,1)-zlev(ig,1))
+                                                      ENDDO
+      DO k=2,nlay
+                                                      DO ig=1,ngrid
+        unsdzdec(ig,k)=1.E+0/(zlay(ig,k)-zlay(ig,k-1))
+                                                     ENDDO
+      ENDDO
+                                                      DO ig=1,ngrid
+      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
+                                                     ENDDO
+c
+c.......................................................................
+
+      do k=2,klev
+                                                          do ig=1,ngrid
+         dz(ig,k)=zlay(ig,k)-zlay(ig,k-1)
+         m2(ig,k)=((u(ig,k)-u(ig,k-1))**2+(v(ig,k)-v(ig,k-1))**2)
+     s             /(dz(ig,k)*dz(ig,k))
+         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
+         n2(ig,k)=g*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
+c        n2(ig,k)=0.
+         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
+         if (ri.lt.ric) then
+            rif(ig,k)=frif(ri)
+         else
+            rif(ig,k)=rifc
+         endif
+         if(rif(ig,k).lt.0.16) then
+            alpha(ig,k)=falpha(rif(ig,k))
+            sm(ig,k)=fsm(rif(ig,k))
+         else
+            alpha(ig,k)=1.12
+            sm(ig,k)=0.085
+         endif
+         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
+c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
+
+
+                                                          enddo
+      enddo
+
+
+c====================================================================
+c   Au premier appel, on determine l et q2 de facon iterative.
+c iterration pour determiner la longueur de melange
+
+
+      if (first.or.iflag_pbl.eq.6) then
+                                                          do ig=1,ngrid
+      l0(ig)=10.
+                                                          enddo
+      do k=2,klev-1
+                                                          do ig=1,ngrid
+        l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
+                                                          enddo
+      enddo
+
+      do iter=1,10
+                                                          do ig=1,ngrid
+         sq(ig)=1.e-10
+         sqz(ig)=1.e-10
+                                                          enddo
+         do k=2,klev-1
+                                                          do ig=1,ngrid
+           q2(ig,k)=l(ig,k)**2*zz(ig,k)
+           l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
+           zq=sqrt(q2(ig,k))
+           sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
+           sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
+                                                          enddo
+         enddo
+                                                          do ig=1,ngrid
+         l0(ig)=0.2*sqz(ig)/sq(ig)
+c        l0(ig)=30.
+                                                          enddo
+c      print*,'ITER=',iter,'  L0=',l0
+
+      enddo
+
+c     print*,'Fin de l initialisation de q2 et l0'
+
+      endif ! first
+
+c====================================================================
+c  Calcul de la longueur de melange.
+c====================================================================
+
+c   Mise a jour de l0
+                                                          do ig=1,ngrid
+      sq(ig)=1.e-10
+      sqz(ig)=1.e-10
+                                                          enddo
+      do k=2,klev-1
+                                                          do ig=1,ngrid
+        zq=sqrt(q2(ig,k))
+        sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
+        sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
+                                                          enddo
+      enddo
+                                                          do ig=1,ngrid
+      l0(ig)=0.2*sqz(ig)/sq(ig)
+c        l0(ig)=30.
+                                                          enddo
+c      print*,'ITER=',iter,'  L0=',l0
+c   calcul de l(z)
+      do k=2,klev
+                                                          do ig=1,ngrid
+         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
+         if(first) then
+           q2(ig,k)=l(ig,k)**2*zz(ig,k)
+         endif
+                                                          enddo
+      enddo
+
+c====================================================================
+c   Yamada 2.0
+c====================================================================
+      if (iflag_pbl.eq.6) then
+
+      do k=2,klev
+                                                          do ig=1,ngrid
+         q2(ig,k)=l(ig,k)**2*zz(ig,k)
+                                                          enddo
+      enddo
+
+
+      else if (iflag_pbl.eq.7) then
+c====================================================================
+c   Yamada 2.Fournier
+c====================================================================
+
+c  Calcul de l,  km, au pas precedent
+      do k=2,klev
+                                                          do ig=1,ngrid
+c        print*,'SMML=',sm(ig,k),l(ig,k)
+         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
+         kmpre(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
+         mpre(ig,k)=sqrt(m2(ig,k))
+c        print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
+                                                          enddo
+      enddo
+
+      do k=2,klev-1
+                                                          do ig=1,ngrid
+        m2cstat=max(alpha(ig,k)*n2(ig,k)+delta(ig,k)/b1,1.e-12)
+        mcstat=sqrt(m2cstat)
+
+c        print*,'M2 L=',k,mpre(ig,k),mcstat
+c
+c  -----{puis on ecrit la valeur de q qui annule l'equation de m
+c        supposee en q3}
+c
+        IF (k.eq.2) THEN
+          kmcstat=1.E+0 / mcstat
+     &    *( unsdz(ig,k)*kmpre(ig,k+1)
+     &                        *mpre(ig,k+1)
+     &      +unsdz(ig,k-1)
+     &              *cd(ig)
+     &              *( sqrt(u(ig,3)**2+v(ig,3)**2)
+     &                -mcstat/unsdzdec(ig,k)
+     &                -mpre(ig,k+1)/unsdzdec(ig,k+1) )**2)
+     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
+        ELSE
+          kmcstat=1.E+0 / mcstat
+     &    *( unsdz(ig,k)*kmpre(ig,k+1)
+     &                        *mpre(ig,k+1)
+     &      +unsdz(ig,k-1)*kmpre(ig,k-1)
+     &                          *mpre(ig,k-1) )
+     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
+        ENDIF
+c       print*,'T2 L=',k,tmp2
+        tmp2=kmcstat
+     &      /( sm(ig,k)/q2(ig,k) )
+     &      /l(ig,k)
+        q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
+c       print*,'Q2 L=',k,q2(ig,k)
+c
+                                                          enddo
+      enddo
+
+      else if (iflag_pbl.ge.8) then
+c====================================================================
+c   Yamada 2.5 a la Didi
+c====================================================================
+
+
+c  Calcul de l,  km, au pas precedent
+      do k=2,klev
+                                                          do ig=1,ngrid
+c        print*,'SMML=',sm(ig,k),l(ig,k)
+         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
+         if (delta(ig,k).lt.1.e-20) then
+c     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
+            delta(ig,k)=1.e-20
+         endif
+         km(ig,k)=l(ig,k)*sqrt(q2(ig,k))*sm(ig,k)
+         aa0=
+     s   (m2(ig,k)-alpha(ig,k)*n2(ig,k)-delta(ig,k)/b1)
+         aa1=
+     s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
+c abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
+         aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
+c     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
+         qpre=sqrt(q2(ig,k))
+         if (iflag_pbl.eq.8 ) then
+            if (aa(ig,k).gt.0.) then
+               q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
+            else
+               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
+            endif
+         else ! iflag_pbl=9
+            if (aa(ig,k)*qpre.gt.0.9) then
+               q2(ig,k)=(qpre*10.)**2
+            else
+               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
+            endif
+         endif
+         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
+c     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
+                                                          enddo
+      enddo
+
+      endif ! Fin du cas 8
+
+c     print*,'OK8'
+
+c====================================================================
+c   Calcul des coefficients de mélange
+c====================================================================
+      do k=2,klev
+c     print*,'k=',k
+                                                          do ig=1,ngrid
+cabde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
+         zq=sqrt(q2(ig,k))
+         km(ig,k)=l(ig,k)*zq*sm(ig,k)
+         kn(ig,k)=km(ig,k)*alpha(ig,k)
+         kq(ig,k)=l(ig,k)*zq*0.2
+c     print*,'KML=',km(ig,k),kn(ig,k)
+                                                          enddo
+      enddo
+
+c     if (iflag_pbl.ge.7..and.0.eq.1) then
+c        q2(:,1)=q2(:,2)
+c        call vdif_q2(dt,g,rconst,plev,temp,kq,q2)
+c     endif
+
+c   Traitement des cas noctrunes avec l'introduction d'une longueur
+c   minilale.
+
+c====================================================================
+c   Traitement particulier pour les cas tres stables.
+c   D'apres Holtslag Boville.
+
+      print*,'YAMADA4 0'
+
+                                                          do ig=1,ngrid
+      coriol(ig)=1.e-4
+      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
+                                                          enddo
+
+       print*,'pblhmin ',pblhmin
+CTest a remettre 21 11 02
+c test abd 13 05 02      if(0.eq.1) then
+      if(1.eq.1) then
+      do k=2,klev
+         do ig=1,klon
+            if (teta(ig,2).gt.teta(ig,1)) then
+               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
+               kmin=kap*zlev(ig,k)*qmin
+            else
+               kmin=-1. ! kmin n'est utilise que pour les SL stables.
+            endif 
+            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
+c               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
+c     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
+               kn(ig,k)=kmin
+               km(ig,k)=kmin
+               kq(ig,k)=kmin
+c   la longueur de melange est suposee etre l= kap z
+c   K=l q Sm d'ou q2=(K/l Sm)**2
+               q2(ig,k)=(qmin/sm(ig,k))**2
+            endif
+         enddo
+      enddo
+      endif
+
+      print*,'YAMADA4 1'
+c   Diagnostique pour stokage
+
+      rino=rif
+      smyam(:,1:klev)=sm(:,1:klev)
+      styam=sm(:,1:klev)*alpha(:,1:klev)
+      lyam(1:klon,1:klev)=l(:,1:klev)
+      knyam(1:klon,1:klev)=kn(:,1:klev)
+
+c   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
+
+        if(1.eq.0)then
+      w2yam=q2(:,1:klev)*0.24
+     s    +lyam(:,1:klev)*5.17*kn(:,1:klev)*n2(:,1:klev)
+     s   /sqrt(q2(:,1:klev))
+
+      t2yam=9.1*kn(:,1:klev)*dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))
+     s  *lyam(:,1:klev)
+	endif
+
+c     print*,'OKFIN'
+      first=.false.
+      return
+      end
Index: /LMDZ4/trunk/makegcm
===================================================================
--- /LMDZ4/trunk/makegcm	(revision 540)
+++ /LMDZ4/trunk/makegcm	(revision 541)
@@ -779,4 +779,5 @@
      else
        set opt_link=" -C hopt -float0 $optdbl -P static -L$MODIPSLDIR $link_veget -lsxioipsl $NCDFLIB "
+     endif
    endif
    set mod_loc_dir="./"
Index: /LMDZ4/trunk/physiq.def
===================================================================
--- /LMDZ4/trunk/physiq.def	(revision 540)
+++ /LMDZ4/trunk/physiq.def	(revision 541)
@@ -80,3 +80,4 @@
 #cdhmax = 2.0E-3
 cdhmax = 0.002
-
+iflag_pbl = 1
+iflag_thermals = 0
