Index: /LMDZ5/branches/testing/libf/bibio/i1mach.F
===================================================================
--- /LMDZ5/branches/testing/libf/bibio/i1mach.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/bibio/i1mach.F	(revision 2220)
@@ -1,4 +1,5 @@
 *DECK I1MACH
       INTEGER FUNCTION I1MACH (I)
+      IMPLICIT NONE
 C***BEGIN PROLOGUE  I1MACH
 C***PURPOSE  Return integer machine dependent constants.
@@ -95,4 +96,5 @@
       SAVE IMACH
       EQUIVALENCE (IMACH(4),OUTPUT)
+      INTEGER I
 C***FIRST EXECUTABLE STATEMENT  I1MACH
       IMACH( 1) =         5
Index: /LMDZ5/branches/testing/libf/bibio/j4save.F
===================================================================
--- /LMDZ5/branches/testing/libf/bibio/j4save.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/bibio/j4save.F	(revision 2220)
@@ -1,4 +1,5 @@
 *DECK J4SAVE
       FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
+      IMPLICIT NONE
 C***BEGIN PROLOGUE  J4SAVE
 C***SUBSIDIARY
@@ -59,4 +60,5 @@
       DATA IPARAM(5)/1/
       DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
+      INTEGER J4SAVE,IWHICH,IVALUE
 C***FIRST EXECUTABLE STATEMENT  J4SAVE
       J4SAVE = IPARAM(IWHICH)
Index: /LMDZ5/branches/testing/libf/bibio/xercnt.F
===================================================================
--- /LMDZ5/branches/testing/libf/bibio/xercnt.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/bibio/xercnt.F	(revision 2220)
@@ -1,4 +1,5 @@
 *DECK XERCNT
       SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
+      IMPLICIT NONE
 C***BEGIN PROLOGUE  XERCNT
 C***SUBSIDIARY
@@ -56,4 +57,5 @@
 C***END PROLOGUE  XERCNT
       CHARACTER*(*) LIBRAR, SUBROU, MESSG
+      INTEGER NERR, LEVEL, KONTRL
 C***FIRST EXECUTABLE STATEMENT  XERCNT
       RETURN
Index: /LMDZ5/branches/testing/libf/bibio/xermsg.F
===================================================================
--- /LMDZ5/branches/testing/libf/bibio/xermsg.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/bibio/xermsg.F	(revision 2220)
@@ -1,4 +1,5 @@
 *DECK XERMSG
       SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
+      IMPLICIT NONE
 C***BEGIN PROLOGUE  XERMSG
 C***PURPOSE  Process error messages for SLATEC and other libraries.
@@ -189,4 +190,7 @@
       CHARACTER*72  TEMP
       CHARACTER*20  LFIRST
+      INTEGER NERR, LEVEL, LKNTRL
+      INTEGER J4SAVE, MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
+      INTEGER MKNTRL, LTEMP
 C***FIRST EXECUTABLE STATEMENT  XERMSG
       LKNTRL = J4SAVE (2, 0, .FALSE.)
Index: /LMDZ5/branches/testing/libf/bibio/xerprn.F
===================================================================
--- /LMDZ5/branches/testing/libf/bibio/xerprn.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/bibio/xerprn.F	(revision 2220)
@@ -1,4 +1,5 @@
 *DECK XERPRN
       SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
+      IMPLICIT NONE
 C***BEGIN PROLOGUE  XERPRN
 C***SUBSIDIARY
@@ -81,4 +82,6 @@
       CHARACTER*2 NEWLIN
       PARAMETER (NEWLIN = '$$')
+      INTEGER N, I1MACH, I, LPREF, LWRAP, LENMSG, NEXTC
+      INTEGER LPIECE, IDELTA
 C***FIRST EXECUTABLE STATEMENT  XERPRN
       CALL XGETUA(IU,NUNIT)
Index: /LMDZ5/branches/testing/libf/bibio/xersve.F
===================================================================
--- /LMDZ5/branches/testing/libf/bibio/xersve.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/bibio/xersve.F	(revision 2220)
@@ -2,4 +2,5 @@
       SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL,
      +   ICOUNT)
+      IMPLICIT NONE
 C***BEGIN PROLOGUE  XERSVE
 C***SUBSIDIARY
@@ -58,5 +59,5 @@
 C   920501  Reformatted the REFERENCES section.  (WRB)
 C***END PROLOGUE  XERSVE
-      PARAMETER (LENTAB=10)
+      INTEGER,PARAMETER :: LENTAB=10
       INTEGER LUN(5)
       CHARACTER*(*) LIBRAR, SUBROU, MESSG
@@ -66,4 +67,7 @@
       SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
       DATA KOUNTX/0/, NMSG/0/
+      INTEGER NERR,LEVEL,KONTRL
+      INTEGER NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
+      INTEGER KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I1MACH, I
 C***FIRST EXECUTABLE STATEMENT  XERSVE
 C
Index: /LMDZ5/branches/testing/libf/bibio/xgetua.F
===================================================================
--- /LMDZ5/branches/testing/libf/bibio/xgetua.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/bibio/xgetua.F	(revision 2220)
@@ -1,4 +1,5 @@
 *DECK XGETUA
       SUBROUTINE XGETUA (IUNITA, N)
+      IMPLICIT NONE
 C***BEGIN PROLOGUE  XGETUA
 C***PURPOSE  Return unit number(s) to which error messages are being
@@ -41,4 +42,5 @@
 C***END PROLOGUE  XGETUA
       DIMENSION IUNITA(5)
+      INTEGER IUNITA, N, J4SAVE, INDEX, I
 C***FIRST EXECUTABLE STATEMENT  XGETUA
       N = J4SAVE(5,0,.FALSE.)
Index: /LMDZ5/branches/testing/libf/dyn3d_common/grid_atob.F
===================================================================
--- /LMDZ5/branches/testing/libf/dyn3d_common/grid_atob.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/dyn3d_common/grid_atob.F	(revision 2220)
@@ -936,4 +936,5 @@
 c
       SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,im,jm,distance)
+      IMPLICIT NONE
 c
 c Auteur: Laurent Li (le 30 decembre 1996)
@@ -958,4 +959,5 @@
       REAL radius
       PARAMETER (radius=6371229.)
+      INTEGER i,j
 c
       pi = 4.0 * ATAN(1.0)
Index: /LMDZ5/branches/testing/libf/dyn3d_common/grid_noro.F
===================================================================
--- /LMDZ5/branches/testing/libf/dyn3d_common/grid_noro.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/dyn3d_common/grid_noro.F	(revision 2220)
@@ -54,8 +54,8 @@
 C=======================================================================
 
-      IMPLICIT INTEGER (I,J)
-      IMPLICIT REAL(X,Z) 
+      IMPLICIT NONE 
       
-	  parameter(iusn=2160,jusn=1080,iext=216, epsfra = 1.e-5)
+      INTEGER,PARAMETER :: iusn=2160,jusn=1080,iext=216
+      REAL,PARAMETER :: epsfra = 1.e-5
 #include "dimensions.h"
 	  REAL xusn(iusn+2*iext),yusn(jusn+2)	
@@ -89,4 +89,13 @@
       REAL a(2200),b(2200),c(1100),d(1100)
       logical masque_lu
+      INTEGER i,j,ii,jj
+      REAL xpi, rad, zdeltay, masque
+      REAL zdeltax, zlenx, zleny, xincr
+      REAL zbordnor, zbordsud, weighy, zbordest, zbordoue, weighx
+      REAL zllmmea, zllmstd, zllmsig, zllmgam, zllmpic, zllmval
+      REAL zllmthe, zminthe
+      REAL xk, xl, xm, xp, xq, xw
+      REAL zmeanor, zmeasud, zstdnor, zstdsud, zsignor, zsigsud
+      REAL zweinor, zweisud, zpicnor, zpicsud, zvalnor, zvalsud
 c
       print *,' parametres de l orographie a l echelle sous maille' 
@@ -455,9 +464,11 @@
 
       SUBROUTINE MVA9(X,IMAR,JMAR)
-
+      IMPLICIT NONE
 C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
-
+      INTEGER IMAR,JMAR
       REAL X(IMAR,JMAR),XF(IMAR,JMAR)
       real WEIGHTpb(-1:1,-1:1)
+      INTEGER I,J,IS,JS
+      REAL SUM
 
 
Index: /LMDZ5/branches/testing/libf/dyn3d_common/juldate.F
===================================================================
--- /LMDZ5/branches/testing/libf/dyn3d_common/juldate.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/dyn3d_common/juldate.F	(revision 2220)
@@ -7,5 +7,11 @@
 c	En entree:an,mois,jour,heure,min.,sec.
 c	En sortie:tjd
-	implicit real (a-h,o-z)
+        IMPLICIT NONE
+        INTEGER,INTENT(IN) :: ian,imoi,ijou,oh,om,os
+        REAL,INTENT(OUT) :: tjd,tjdsec
+        
+        REAL frac,year,rmon,cf,a,b
+        INTEGER ojou
+        
 	frac=((os/60.+om)/60.+oh)/24.
 	ojou=dble(ijou)+frac
Index: /LMDZ5/branches/testing/libf/dyn3d_common/massbar.F
===================================================================
--- /LMDZ5/branches/testing/libf/dyn3d_common/massbar.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/dyn3d_common/massbar.F	(revision 2220)
@@ -3,4 +3,5 @@
 !
       SUBROUTINE massbar(  masse, massebx, masseby )
+      IMPLICIT NONE
 c
 c **********************************************************************
@@ -24,4 +25,5 @@
       REAL    masse( ip1jmp1,llm ), massebx( ip1jmp1,llm )  ,
      *      masseby(   ip1jm,llm )
+      INTEGER ij,l
 c
 c
Index: /LMDZ5/branches/testing/libf/dyn3d_common/massbarxy.F
===================================================================
--- /LMDZ5/branches/testing/libf/dyn3d_common/massbarxy.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/dyn3d_common/massbarxy.F	(revision 2220)
@@ -3,4 +3,5 @@
 !
       SUBROUTINE massbarxy(  masse, massebxy )
+      IMPLICIT NONE
 c
 c **********************************************************************
@@ -23,4 +24,5 @@
 c
        REAL  masse( ip1jmp1,llm ), massebxy( ip1jm,llm )
+       INTEGER ij,l
 c
 
Index: /LMDZ5/branches/testing/libf/dyn3d_common/ppm3d.F
===================================================================
--- /LMDZ5/branches/testing/libf/dyn3d_common/ppm3d.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/dyn3d_common/ppm3d.F	(revision 2220)
@@ -66,5 +66,5 @@
      &                  JNP,j1,NLAY,AP,BP,PT,AE,fill,dum,Umax)
 
-c      implicit none
+      implicit none
 
 c     rajout de déclarations
@@ -270,5 +270,5 @@
 C     User modifiable parameters
 C
-      parameter (Jmax = 361, kmax = 150)
+      integer,parameter :: Jmax = 361, kmax = 150
 C
 C ****6***0*********0*********0*********0*********0*********0**********72
@@ -299,7 +299,13 @@
       data NDT0, NSTEP /0, 0/
       data cross /.true./
+      REAL DTDY, DTDY5, RCAP
+      INTEGER JS0, JN0, IML, JMR, IMJM
       SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML,
      &     DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK
 C
+      INTEGER NDT0, NSTEP, j2, k,j,i,ic,l,JS,JN,IMH
+      INTEGER IU,IIU,JT,iad,jad,krd
+      REAL r23,r3,PI,DL,DP,DT,CR1,MAXDT,ZTC,D5
+      REAL sum1,sum2,ru
             
       JMR = JNP -1
@@ -756,6 +762,8 @@
       subroutine FZPPM(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6,
      &                 flux,wk1,wk2,wz2,delp,KORD)
-      parameter ( kmax = 150 )
-      parameter ( R23 = 2./3., R3 = 1./3.)
+      implicit none
+      integer,parameter :: kmax = 150 
+      real,parameter :: R23 = 2./3., R3 = 1./3.
+      integer IMR,JNP,NLAY,J1,KORD
       real WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY),
      &     wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY),
@@ -764,4 +772,6 @@
       real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*),
      &     wz2(IMR,*)
+      integer JMR,IMJM,NLAYM1,LMT,K,I,J
+      real c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp
 C
       JMR = JNP - 1
@@ -922,7 +932,14 @@
       subroutine xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC,
      &               fx1,xmass,IORD)
+      implicit none
+      integer IMR,JNP,IML,j1,j2,JN,JS,IORD
+      real PU,DQ,Q,UC,fx1,xmass
+      real dc,qtmp
+      integer ISAVE(IMR)
       dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP)
      &    ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML)
-      dimension PU(IMR,JNP),Q(IMR,JNP),ISAVE(IMR)
+      dimension PU(IMR,JNP),Q(IMR,JNP)
+      integer jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp
+      real rut
 C
       IMP = IMR + 1
@@ -1031,8 +1048,11 @@
 C
       subroutine fxppm(IMR,IML,UT,P,DC,flux,IORD)
-      parameter ( R3 = 1./3., R23 = 2./3. )
+      implicit none
+      integer IMR,IML,IORD
+      real UT,P,DC,flux
+      real,parameter ::  R3 = 1./3., R23 = 2./3. 
       DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1)
-      DIMENSION AR(0:IMR),AL(0:IMR),A6(0:IMR)
-      integer LMT 
+      REAL :: AR(0:IMR),AL(0:IMR),A6(0:IMR)
+      integer LMT,IMP,JLVL,i
 c      logical first
 c      data first /.true./
@@ -1088,6 +1108,10 @@
 C
       subroutine xmist(IMR,IML,P,DC)
-      parameter( R24 = 1./24.)
-      dimension P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
+      implicit none
+      integer IMR,IML
+      real,parameter :: R24 = 1./24.
+      real :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
+      integer :: i
+      real :: tmp,pmax,pmin
 C
       do 10  i=1,IMR
@@ -1101,8 +1125,13 @@
       subroutine ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2
      &              ,ymass,fx,A6,AR,AL,JORD)
+      implicit none
+      integer :: IMR,JNP,j1,j2,JORD
+      real :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL
       dimension P(IMR,JNP),VC(IMR,JNP),ymass(IMR,JNP)
      &       ,DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP)
 C Work array
       DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
+      integer :: JMR,len,i,jt,j
+      real :: sum1,sum2
 C
       JMR = JNP - 1
@@ -1161,6 +1190,10 @@
 C
       subroutine  ymist(IMR,JNP,j1,P,DC,ID)
-      parameter ( R24 = 1./24. )
-      dimension P(IMR,JNP),DC(IMR,JNP)
+      implicit none
+      integer :: IMR,JNP,j1,ID
+      real,parameter :: R24 = 1./24.
+      real :: P(IMR,JNP),DC(IMR,JNP)
+      integer :: iimh,jmr,ijm3,imh,i
+      real :: pmax,pmin,tmp
 C
       IMH = IMR / 2
@@ -1239,9 +1272,12 @@
 C
       subroutine fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)
-      parameter ( R3 = 1./3., R23 = 2./3. )
+      implicit none
+      integer IMR,JNP,j1,j2,JORD
+      real,parameter :: R3 = 1./3., R23 = 2./3.
       real VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*)
 C Local work arrays.
       real AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
-      integer LMT
+      integer LMT,i
+      integer IMH,JMR,j11,IMJM1,len
 c      logical first
 C      data first /.true./
@@ -1315,6 +1351,10 @@
 C
 	subroutine yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)
+        implicit none
+        integer IMR,JNP,j1,j2,IAD
 	REAL p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP)
         REAL WK(IMR,-1:JNP+2)
+        INTEGER JMR,IMH,i,j,jp
+        REAL rv,a1,b1,sum1,sum2
 C
 	JMR = JNP-1
@@ -1401,5 +1441,9 @@
 C
 	subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)
+        implicit none
+        INTEGER IMR,JNP,j1,j2,JS,JN,IML,IAD
 	REAL p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP)
+        INTEGER JMR,j,i,ip,iu,iiu
+        REAL ru,a1,b1
 C
 	JMR = JNP-1
@@ -1489,4 +1533,5 @@
 C
       subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)
+      implicit none
 C
 C A6 =  CURVATURE OF THE TEST PARABOLA
@@ -1503,6 +1548,9 @@
 C LMT = 2: POSITIVE-DEFINITE CONSTRAINT
 C
-      parameter ( R12 = 1./12. )
-      dimension A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
+      real,parameter :: R12 = 1./12. 
+      real :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
+      integer :: IM,LMT
+      INTEGER i
+      REAL da1,da2,a6da,fmin
 C
       if(LMT.eq.0) then
@@ -1564,5 +1612,8 @@
 C
       subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
-      dimension U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*)
+      implicit none
+      integer IMR,JMR,j1,j2
+      real :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5
+      integer i,j
 C
       do 35 j=j1,j2
@@ -1579,5 +1630,9 @@
 C
       subroutine cosa(cosp,cose,JNP,PI,DP)
-      dimension cosp(*),cose(*)
+      implicit none
+      integer JNP
+      real :: cosp(*),cose(*),PI,DP
+      integer JMR,j,jeq
+      real ph5
       JMR = JNP-1
       do 55 j=2,JNP
@@ -1606,5 +1661,9 @@
 C
       subroutine cosc(cosp,cose,JNP,PI,DP)
-      dimension cosp(*),cose(*)
+      implicit none
+      integer JNP
+      real :: cosp(*),cose(*),PI,DP
+      real phi
+      integer j
 C
       phi = -0.5*PI
@@ -1628,7 +1687,10 @@
      &                   cross,IC,NSTEP)
 C
-      parameter( tiny = 1.E-60 )
-      DIMENSION Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)
+      real,parameter :: tiny = 1.E-60 
+      INTEGER :: IMR,JNP,NLAY,j1,j2,IC,NSTEP
+      REAL :: Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)
       logical cross
+      INTEGER :: NLAYM1,len,ip,L,icr,ipy,ipx,i
+      real :: qup,qly,dup,sum
 C
       NLAYM1 = NLAY-1
@@ -1730,5 +1792,9 @@
 C
       subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
-      dimension q(IMR,*),cosp(*),acosp(*)
+      implicit none
+      integer :: IMR,JNP,j1,j2,icr
+      real :: q(IMR,*),cosp(*),acosp(*),tiny
+      integer :: i,j
+      real :: dq,dn,d0,d1,ds,d2
       icr = 0
       do 65 j=j1+1,j2-1
@@ -1828,5 +1894,9 @@
 C
       subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
-      dimension q(IMR,*),cosp(*),acosp(*)
+      implicit none
+      integer :: IMR,JNP,j1,j2,ipy
+      real :: q(IMR,*),cosp(*),acosp(*),tiny
+      real :: DP,CAP1,dq,dn,d0,d1,ds,d2
+      INTEGER :: i,j
 c      logical first
 c      data first /.true./
@@ -1910,5 +1980,9 @@
 C
       subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)
-      dimension q(IMR,*),qtmp(JNP,IMR)
+      implicit none
+      integer :: IMR,JNP,j1,j2,ipx
+      real :: q(IMR,*),qtmp(JNP,IMR),tiny
+      integer :: i,j
+      real :: d0,d1,d2
 C
       ipx = 0
@@ -1983,8 +2057,11 @@
 C
       subroutine zflip(q,im,km,nc)
+      implicit none
 C This routine flip the array q (in the vertical).
+      integer :: im,km,nc
       real q(im,km,nc)
 C local dynamic array
       real qtmp(im,km)
+      integer IC,k,i
 C
       do 4000 IC = 1, nc
Index: /LMDZ5/branches/testing/libf/dyn3d_common/ran1.F
===================================================================
--- /LMDZ5/branches/testing/libf/dyn3d_common/ran1.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/dyn3d_common/ran1.F	(revision 2220)
@@ -3,11 +3,15 @@
 !
       FUNCTION RAN1(IDUM)
-      DIMENSION R(97)
-      save r
-      save iff,ix1,ix2,ix3
-      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
-      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
-      PARAMETER (M3=243000,IA3=4561,IC3=51349)
-      DATA IFF /0/
+      IMPLICIT NONE
+      REAL RAN1
+      REAL,SAVE :: R(97)
+      REAL,PARAMETER :: RM1=3.8580247E-6,RM2=7.4373773E-6
+      INTEGER,SAVE :: IFF=0
+      integer,save :: ix1,ix2,ix3
+      INTEGER,PARAMETER :: M1=259200,IA1=7141,IC1=54773
+      INTEGER,PARAMETER :: M2=134456,IA2=8121,IC2=28411
+      INTEGER,PARAMETER :: M3=243000,IA3=4561,IC3=51349
+      INTEGER :: IDUM,J
+
       IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
         IFF=1
Index: /LMDZ5/branches/testing/libf/filtrez/acc.F
===================================================================
--- /LMDZ5/branches/testing/libf/filtrez/acc.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/filtrez/acc.F	(revision 2220)
@@ -3,5 +3,10 @@
 !
         subroutine acc(vec,d,im)
-        dimension vec(im,im),d(im)
+        implicit none
+        integer :: im
+        real :: vec(im,im),d(im)
+        integer :: i,j
+        real ::sum
+        real,external :: ssum
         do j=1,im
           do i=1,im
Index: /LMDZ5/branches/testing/libf/filtrez/eigen.F
===================================================================
--- /LMDZ5/branches/testing/libf/filtrez/eigen.F	(revision 2219)
+++ /LMDZ5/branches/testing/libf/filtrez/eigen.F	(revision 2220)
@@ -3,7 +3,9 @@
 !
       SUBROUTINE eigen( e,d)
+      IMPLICIT NONE
 #include "dimensions.h"
-      dimension e( iim,iim ), d( iim )
-      dimension asm( iim )
+      real :: e( iim,iim ), d( iim )
+      real :: asm( iim )
+      integer :: im,i,j
       im=iim
 c
Index: /LMDZ5/branches/testing/libf/grid/dimension/makdim
===================================================================
--- /LMDZ5/branches/testing/libf/grid/dimension/makdim	(revision 2219)
+++ /LMDZ5/branches/testing/libf/grid/dimension/makdim	(revision 2220)
@@ -5,25 +5,32 @@
 if (( $# < 1 )) || (( $# > 3 ))
 then
- echo "Wrong number of parameters in $0 !!!"
- echo " Usage:"
- echo "  $0 [im] [jm] lm"
- echo " where im, jm and lm are the dimensions"
- exit
+    echo "Wrong number of parameters in $0 !!!"
+    echo " Usage:"
+    echo "  $0 [im] [jm] lm"
+    echo " where im, jm and lm are the dimensions"
+    exit 1
+fi
+
+if (($1 % 8 != 0)) && (( $# = 3 ))
+then
+    echo "The number of longitudes must be a multiple of 8."
+    echo "See the files dyn3d/groupe.F and dyn3dmem/groupe_loc.F."
+    exit 1
 fi
 
 # build "fichnom", the relevant 'dimensions.im.jm.lm' file name
 for i in $*
-  do
-  list=$list.$i
+do
+    list=$list.$i
 done
 fichdim=dimensions${list}
 
 if [ ! -f $fichdim ]
-    then
+then
 #    echo "$fichdim does not exist"
 
     # assign values of im, jm and lm
     if [ $# -ge 3 ]
-        then
+    then
         im=$1
         jm=$2
@@ -31,5 +38,5 @@
         ndm=1
     elif [ $# -ge 2 ]
-        then
+    then
         im=1
         jm=$1
@@ -37,5 +44,5 @@
         ndm=1
     elif [ $# -ge 1 ]
-        then
+    then
         im=1
         jm=1
@@ -45,5 +52,5 @@
 
 # since the file doesn't exist, we create it
-cat << EOF > $fichdim
+    cat << EOF > $fichdim
 !-----------------------------------------------------------------------
 !   INCLUDE 'dimensions.h'
@@ -64,5 +71,5 @@
 # remove 'old' dimensions.h file (if any) and replace it with new one
 if [ -f ../dimensions.h ] ; then
-\rm ../dimensions.h
+    \rm ../dimensions.h
 fi
 tar cf - $fichdim | ( cd .. ; tar xf - ; mv $fichdim dimensions.h )
Index: /LMDZ5/branches/testing/libf/phylmd/1DUTILS.h
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/1DUTILS.h	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/1DUTILS.h	(revision 2220)
@@ -99,4 +99,6 @@
 !             LS convergence imposed from  RICO (cst) 
 !         = 6 ==> forcing_amma = .true.
+!         = 10 ==> forcing_case = .true.
+!             initial profiles from case.nc file 
 !         = 40 ==> forcing_GCSSold = .true.
 !             initial profile from GCSS file
@@ -132,4 +134,84 @@
           CALL getin('turb_fcg',xTurb_fcg_gcssold)
         ENDIF
+
+!Paramètres de forçage
+!Config  Key  = tend_t
+!Config  Desc = forcage ou non par advection de T
+!Config  Def  = false
+!Config  Help = forcage ou non par advection de T
+       tend_t =0
+       CALL getin('tend_t',tend_t)
+
+!Config  Key  = tend_q
+!Config  Desc = forcage ou non par advection de q
+!Config  Def  = false
+!Config  Help = forcage ou non par advection de q
+       tend_q =0
+       CALL getin('tend_q',tend_q)
+
+!Config  Key  = tend_u
+!Config  Desc = forcage ou non par advection de u
+!Config  Def  = false
+!Config  Help = forcage ou non par advection de u
+       tend_u =0
+       CALL getin('tend_u',tend_u)
+
+!Config  Key  = tend_v
+!Config  Desc = forcage ou non par advection de v
+!Config  Def  = false
+!Config  Help = forcage ou non par advection de v
+       tend_v =0
+       CALL getin('tend_v',tend_v)
+
+!Config  Key  = tend_w
+!Config  Desc = forcage ou non par vitesse verticale
+!Config  Def  = false
+!Config  Help = forcage ou non par vitesse verticale
+       tend_w =0
+       CALL getin('tend_w',tend_w)
+
+!Config  Key  = tend_rayo
+!Config  Desc = forcage ou non par dtrad
+!Config  Def  = false
+!Config  Help = forcage ou non par dtrad
+       tend_rayo =0
+       CALL getin('tend_rayo',tend_rayo)
+
+
+!Config  Key  = nudge_t
+!Config  Desc = constante de nudging de T
+!Config  Def  = false
+!Config  Help = constante de nudging de T
+       nudge_t =0.
+       CALL getin('nudge_t',nudge_t)
+
+!Config  Key  = nudge_q
+!Config  Desc = constante de nudging de q
+!Config  Def  = false
+!Config  Help = constante de nudging de q
+       nudge_q =0.
+       CALL getin('nudge_q',nudge_q)
+
+!Config  Key  = nudge_u
+!Config  Desc = constante de nudging de u
+!Config  Def  = false
+!Config  Help = constante de nudging de u
+       nudge_u =0.
+       CALL getin('nudge_u',nudge_u)
+
+!Config  Key  = nudge_v
+!Config  Desc = constante de nudging de v
+!Config  Def  = false
+!Config  Help = constante de nudging de v
+       nudge_v =0.
+       CALL getin('nudge_v',nudge_v)
+
+!Config  Key  = nudge_w
+!Config  Desc = constante de nudging de w
+!Config  Def  = false
+!Config  Help = constante de nudging de w
+       nudge_w =0.
+       CALL getin('nudge_w',nudge_w)
+
 
 !Config  Key  = iflag_nudge
@@ -2431,4 +2513,162 @@
  
 !=====================================================================
+       SUBROUTINE interp_case_vertical(play,nlev_cas,plev_prof_cas            &
+     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
+     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
+     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
+     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
+     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
+     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
+ 
+       implicit none
+ 
+#include "dimensions.h"
+
+!-------------------------------------------------------------------------
+! Vertical interpolation of TOGA-COARE forcing data onto mod_casel levels
+!-------------------------------------------------------------------------
+ 
+       integer nlevmax
+       parameter (nlevmax=41)
+       integer nlev_cas,mxcalc
+!       real play(llm), plev_prof(nlevmax) 
+!       real t_prof(nlevmax),q_prof(nlevmax)
+!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
+!       real ht_prof(nlevmax),vt_prof(nlevmax)
+!       real hq_prof(nlevmax),vq_prof(nlevmax)
+ 
+       real play(llm), plev_prof_cas(nlev_cas) 
+       real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
+       real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
+       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas)
+       real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
+       real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
+       real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
+       real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
+ 
+       real t_mod_cas(llm),q_mod_cas(llm)
+       real u_mod_cas(llm),v_mod_cas(llm)
+       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm)
+       real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
+       real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
+       real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
+       real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
+ 
+       integer l,k,k1,k2
+       real frac,frac1,frac2,fact
+ 
+       do l = 1, llm
+
+        if (play(l).ge.plev_prof_cas(nlev_cas)) then
+ 
+        mxcalc=l
+         k1=0
+         k2=0
+
+         if (play(l).le.plev_prof_cas(1)) then
+
+         do k = 1, nlev_cas-1
+          if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
+            k1=k
+            k2=k+1
+          endif
+         enddo
+
+         if (k1.eq.0 .or. k2.eq.0) then
+          write(*,*) 'PB! k1, k2 = ',k1,k2
+          write(*,*) 'l,play(l) = ',l,play(l)/100
+         do k = 1, nlev_cas-1
+          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
+         enddo
+         endif
+
+         frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
+         t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
+         q_mod_cas(l)= q_prof_cas(k2) - frac*(q_prof_cas(k2)-q_prof_cas(k1))
+         u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
+         v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
+         ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
+         vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
+         w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
+         du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
+         hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
+         vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
+         dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
+         hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
+         vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
+         dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
+         ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
+         vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
+         dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
+         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
+         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
+     
+         else !play>plev_prof_cas(1)
+
+         k1=1
+         k2=2
+         frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
+         frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
+         t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
+         q_mod_cas(l)= frac1*q_prof_cas(k1) - frac2*q_prof_cas(k2)
+         u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
+         v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
+         ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
+         vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
+         w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
+         du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
+         hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
+         vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
+         dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
+         hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
+         vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
+         dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
+         ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
+         vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
+         dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
+         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
+         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
+
+         endif ! play.le.plev_prof_cas(1)
+
+        else ! above max altitude of forcing file
+ 
+!jyg
+         fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
+         fact = max(fact,0.)                                           !jyg
+         fact = exp(-fact)                                             !jyg
+         t_mod_cas(l)= t_prof_cas(nlev_cas)                                   !jyg
+         q_mod_cas(l)= q_prof_cas(nlev_cas)*fact                              !jyg
+         u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                              !jyg
+         v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                              !jyg
+         ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact                              !jyg
+         vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact                              !jyg
+         w_mod_cas(l)= 0.0                                                 !jyg
+         du_mod_cas(l)= du_prof_cas(nlev_cas)*fact 
+         hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                            !jyg
+         vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                            !jyg
+         dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact 
+         hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                            !jyg
+         vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                            !jyg
+         dt_mod_cas(l)= dt_prof_cas(nlev_cas) 
+         ht_mod_cas(l)= ht_prof_cas(nlev_cas)                                 !jyg
+         vt_mod_cas(l)= vt_prof_cas(nlev_cas)                                 !jyg
+         dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact 
+         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                            !jyg
+         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                            !jyg
+ 
+        endif ! play
+ 
+       enddo ! l
+
+!       do l = 1,llm
+!       print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
+!     $        l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
+!       enddo
+ 
+          return
+          end
+!***************************************************************************** 
+!=====================================================================
        SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof   &
      &         ,th_prof,qv_prof,u_prof,v_prof,o3_prof                     &
Index: /LMDZ5/branches/testing/libf/phylmd/1D_decl_cases.h
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/1D_decl_cases.h	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/1D_decl_cases.h	(revision 2220)
@@ -240,4 +240,15 @@
         real u_profa(nlev_astex),v_profa(nlev_astex),w_profa(nlev_astex)
         real tke_profa(nlev_astex)
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!Declarations specifiques au cas standard
+
+        real w_mod_cas(llm), t_mod_cas(llm),q_mod_cas(llm)
+        real ug_mod_cas(llm),vg_mod_cas(llm)
+        real u_mod_cas(llm),v_mod_cas(llm)
+        real ht_mod_cas(llm),vt_mod_cas(llm),dt_mod_cas(llm),dtrad_mod_cas(llm)
+        real hq_mod_cas(llm),vq_mod_cas(llm),dq_mod_cas(llm)
+        real hu_mod_cas(llm),vu_mod_cas(llm),du_mod_cas(llm)
+        real hv_mod_cas(llm),vv_mod_cas(llm),dv_mod_cas(llm)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
Index: /LMDZ5/branches/testing/libf/phylmd/1D_interp_cases.h
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/1D_interp_cases.h	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/1D_interp_cases.h	(revision 2220)
@@ -597,4 +597,161 @@
       enddo
       endif ! forcing_astex
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Interpolation forcing standard case
+!---------------------------------------------------------------------
+      if (forcing_case) then
+
+        print*,                                                             &
+     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
+     &    daytime,day1,(daytime-day1)*86400.,                               &
+     &    (daytime-day1)*86400/pdt_cas
+
+! time interpolation:
+        CALL interp_case_time(daytime,day1,annee_ref                        &
+     &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas       &
+     &       ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas,ug_cas,vg_cas    &
+     &       ,vitw_cas,du_cas,hu_cas,vu_cas           &
+     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas             &
+     &       ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas                 &
+     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas         &
+     &       ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
+     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
+     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas               &
+     &       ,sens_prof_cas)
+
+             ts_cur = ts_prof_cas 
+             psurf=plev_prof_cas(1)
+
+! vertical interpolation:
+      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
+     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
+     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
+     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
+     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
+     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
+     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
+
+
+!calcul de l'advection verticale a partir du omega
+!Calcul des gradients verticaux
+!initialisation
+      d_t_z(:)=0.
+      d_q_z(:)=0.
+      d_u_z(:)=0.
+      d_v_z(:)=0.
+      d_t_dyn_z(:)=0.
+      d_q_dyn_z(:)=0.
+      d_u_dyn_z(:)=0.
+      d_v_dyn_z(:)=0.
+      DO l=2,llm-1
+       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
+       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
+       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
+       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
+      ENDDO
+      d_t_z(1)=d_t_z(2)
+      d_q_z(1)=d_q_z(2)
+      d_u_z(1)=d_u_z(2)
+      d_v_z(1)=d_v_z(2)
+      d_t_z(llm)=d_t_z(llm-1)
+      d_q_z(llm)=d_q_z(llm-1)
+      d_u_z(llm)=d_u_z(llm-1)
+      d_v_z(llm)=d_v_z(llm-1)
+
+!Calcul de l advection verticale
+      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
+      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
+      d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
+      d_v_dyn_z(:)=w_mod_cas(:)*d_v_z(:)
+
+!wind nudging 
+      if (nudge_u.gt.0.) then
+        do l=1,llm
+           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
+        enddo
+      else
+        do l=1,llm
+        u(l) = u_mod_cas(l) 
+        enddo
+      endif
+
+      if (nudge_v.gt.0.) then
+        do l=1,llm
+           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
+        enddo
+      else
+        do l=1,llm
+        v(l) = v_mod_cas(l) 
+        enddo
+      endif
+
+      if (nudge_w.gt.0.) then
+        do l=1,llm
+           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
+        enddo
+      else
+        do l=1,llm
+        w(l) = w_mod_cas(l) 
+        enddo
+      endif
+
+!nudging of q and temp 
+      if (nudge_t.gt.0.) then
+        do l=1,llm
+           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
+        enddo
+      endif
+      if (nudge_q.gt.0.) then
+        do l=1,llm
+           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
+        enddo
+      endif
+
+      do l = 1, llm
+       omega(l) = w_mod_cas(l)
+       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
+       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
+
+!calcul advection
+        if ((tend_u.eq.1).and.(tend_w.eq.0)) then
+           d_u_adv(l)=du_mod_cas(l)
+        else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
+           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
+        endif
+
+        if ((tend_v.eq.1).and.(tend_w.eq.0)) then
+           d_v_adv(l)=dv_mod_cas(l)
+        else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
+           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
+        endif
+
+        if ((tend_t.eq.1).and.(tend_w.eq.0)) then
+!           d_th_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
+           d_th_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
+        else if ((tend_t.eq.1).and.(tend_w.eq.1)) then
+!           d_th_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
+           d_th_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
+        endif
+
+        if ((tend_q.eq.1).and.(tend_w.eq.0)) then
+!           d_q_adv(l,1)=dq_mod_cas(l)
+           d_q_adv(l,1)=-1*dq_mod_cas(l)
+        else if ((tend_q.eq.1).and.(tend_w.eq.1)) then
+!           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
+           d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
+        endif
+         
+        if (tend_rayo.eq.1) then
+           dt_cooling(l) = dtrad_mod_cas(l)
+        else
+           dt_cooling(l) = 0.0
+        endif
+      enddo
+
+      endif ! forcing_case
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
Index: /LMDZ5/branches/testing/libf/phylmd/1D_read_forc_cases.h
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/1D_read_forc_cases.h	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/1D_read_forc_cases.h	(revision 2220)
@@ -720,3 +720,63 @@
 
       endif ! forcing_astex
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Forcing from standard case :
+!---------------------------------------------------------------------
+
+      if (forcing_case) then
+
+         write(*,*),'avant call read_1D_cas'
+         call read_1D_cas
+         write(*,*) 'Forcing read'
+
+!Time interpolation for initial conditions using TOGA interpolation routine
+         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
+      CALL interp_case_time(day,day1,annee_ref                &
+     &         ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas       &
+     &         ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas               &
+     &         ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas           &
+     &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas             &
+     &         ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas                 &
+     &         ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas         &
+     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
+     &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas       &
+     &         ,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas)
+
+! vertical interpolation using TOGA interpolation routine:
+!      write(*,*)'avant interp vert', t_prof
+      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
+     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas    &
+     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
+     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
+     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas           &
+     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
+     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
+!       write(*,*) 'Profil initial forcing case interpole',t_mod
+
+! initial and boundary conditions :
+!      tsurf = ts_prof_cas
+      ts_cur = ts_prof_cas 
+      psurf=plev_prof_cas(1)
+      write(*,*) 'SST initiale: ',tsurf
+      do l = 1, llm
+       temp(l) = t_mod_cas(l)
+       q(l,1) = q_mod_cas(l)
+       q(l,2) = 0.0
+       u(l) = u_mod_cas(l)
+       v(l) = v_mod_cas(l)
+       omega(l) = w_mod_cas(l)
+       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
+
+       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
+!on applique le forcage total au premier pas de temps
+!attention: signe different de toga
+       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
+       d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
+       d_q_adv(l,2) = 0.0
+       d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
+       d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
+      enddo     
+       
+      endif !forcing_case
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Index: /LMDZ5/branches/testing/libf/phylmd/YOMCST2.h
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/YOMCST2.h	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/YOMCST2.h	(revision 2220)
@@ -3,4 +3,5 @@
       REAL  gammas, alphas, betas, Fmax, qqa1, qqa2, qqa3, scut
       REAL  Qcoef1max,Qcoef2max,Supcrit1,Supcrit2
+      REAL coef_clos_ls
 !
       COMMON/YOMCST2/gammas,    alphas, betas, Fmax, scut,              &
@@ -8,5 +9,5 @@
      &               Qcoef1max,Qcoef2max,                               &
      &               Supcrit1, Supcrit2,                                &
-     &               choice,iflag_mix
+     &               choice,iflag_mix,coef_clos_ls
 !$OMP THREADPRIVATE(/YOMCST2/)
 !    --------------------------------------------------------------------
Index: /LMDZ5/branches/testing/libf/phylmd/calcratqs.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/calcratqs.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/calcratqs.F90	(revision 2220)
@@ -1,4 +1,4 @@
 SUBROUTINE calcratqs(klon,klev,prt_level,lunout,       &
-           iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
+           iflag_ratqs,iflag_con,iflag_cldth,pdtphys, &
            ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,   &
            ptconv,ptconvth,clwcon0th, rnebcon0th,      &
@@ -19,5 +19,5 @@
 ! Input
 integer,intent(in) :: klon,klev,prt_level,lunout
-integer,intent(in) :: iflag_con,iflag_cldcon,iflag_ratqs
+integer,intent(in) :: iflag_con,iflag_cldth,iflag_ratqs
 real,intent(in) :: pdtphys,ratqsbas,ratqshaut,fact_cldcon,tau_ratqs
 real, dimension(klon,klev+1),intent(in) :: paprs
@@ -43,5 +43,5 @@
 !   ----------------
 !   on ecrase le tableau ratqsc calcule par clouds_gno
-      if (iflag_cldcon.eq.1) then
+      if (iflag_cldth.eq.1) then
          do k=1,klev
          do i=1,klon
@@ -58,5 +58,5 @@
 !  par nversion de la fonction log normale
 !-----------------------------------------------------------------------
-      else if (iflag_cldcon.eq.4) then
+      else if (iflag_cldth.eq.4) then
          ptconvth(:,:)=.false.
          ratqsc(:,:)=0.
@@ -136,5 +136,5 @@
 !  -----------
 
-      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2.or.iflag_cldcon.eq.4) then
+      if (iflag_cldth.eq.1 .or.iflag_cldth.eq.2.or.iflag_cldth.eq.4) then
 
 ! On ajoute une constante au ratqsc*2 pour tenir compte de 
@@ -165,5 +165,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
-      else if (iflag_cldcon<=6) then
+      else if (iflag_cldth<=6) then
 !   on ne prend que le ratqs stable pour fisrtilp
          ratqs(:,:)=ratqss(:,:)
@@ -174,5 +174,5 @@
              do i=1,klon
                 if (ratqsc(i,k).gt.1.e-10) then
-                   ratqs(i,k)=ratqs(i,k)*zfratqs2+(iflag_cldcon/100.)*ratqsc(i,k)*(1.-zfratqs2)
+                   ratqs(i,k)=ratqs(i,k)*zfratqs2+(iflag_cldth/100.)*ratqsc(i,k)*(1.-zfratqs2)
                 endif
                 ratqs(i,k)=min(ratqs(i,k)*zfratqs1+ratqss(i,k)*(1.-zfratqs1),0.5)
Index: /LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90	(revision 2220)
@@ -23,9 +23,9 @@
 
     USE dimphy 
-    USE surface_data, ONLY : type_ocean
+    USE surface_data, ONLY : type_ocean,version_ocean
     USE limit_read_mod
     USE pbl_surface_mod, ONLY : pbl_surface_newfrac
     USE cpl_mod, ONLY : cpl_receive_frac
-    USE ocean_slab_mod, ONLY : ocean_slab_frac
+    USE ocean_slab_mod, ONLY : fsic, ocean_slab_frac
     USE indice_sol_mod
 
@@ -146,4 +146,10 @@
           pctsrf(:,is_sic) = 0.
        END WHERE
+! Send fractions back to slab ocean if needed
+       IF (type_ocean == 'slab'.AND. version_ocean.NE.'sicINT') THEN
+           WHERE (1.-zmasq(:)>EPSFRA)
+               fsic(:)=pctsrf(:,is_sic)/(1.-zmasq(:))
+           END WHERE
+       END IF
 
 !****************************************************************************************
Index: /LMDZ5/branches/testing/libf/phylmd/clift.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/clift.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/clift.F90	(revision 2220)
@@ -3,4 +3,5 @@
 
 SUBROUTINE clift(p, t, rr, rs, plcl, dplcldt, dplcldq)
+IMPLICIT NONE
   ! ***************************************************************
   ! *                                                             *
@@ -41,4 +42,6 @@
 
   include "YOMCST.h"
+  real :: p,t,rr,rs,plcl,dplcldt,dplcldq,cpd,cpv,cl,cpvmcl,eps,alv0,a,b
+  real :: rh,chi,alv
 
   cpd = rcpd
Index: /LMDZ5/branches/testing/libf/phylmd/clouds_bigauss.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/clouds_bigauss.F90	(revision 2220)
+++ /LMDZ5/branches/testing/libf/phylmd/clouds_bigauss.F90	(revision 2220)
@@ -0,0 +1,119 @@
+
+! $Header$
+
+
+! ================================================================================
+
+SUBROUTINE clouds_bigauss(klon, nd, r, rs, qtc, sigt, ptconv, ratqsc, cldf)
+  IMPLICIT NONE
+
+  ! --------------------------------------------------------------------------------
+
+  ! Inputs:
+
+  ! ND----------: Number of vertical levels
+  ! R--------ND-: Domain-averaged mixing ratio of total water
+  ! RS-------ND-: Mean saturation humidity mixing ratio within the gridbox
+  ! QSUB-----ND-: Mixing ratio of condensed water within clouds associated
+  ! with SUBGRID-SCALE condensation processes (here, it is
+  ! predicted by the convection scheme)
+  ! Outputs:
+
+  ! PTCONV-----ND-: Point convectif = TRUE
+  ! RATQSC-----ND-: Largeur normalisee de la distribution
+  ! CLDF-----ND-: Fraction nuageuse
+
+  ! --------------------------------------------------------------------------------
+
+
+  INTEGER klon, nd
+  REAL r(klon, nd), rs(klon, nd), qtc(klon, nd), sigt(klon, nd)
+  LOGICAL ptconv(klon, nd)
+  REAL ratqsc(klon, nd)
+  REAL cldf(klon, nd)
+
+  ! -- parameters controlling the iteration:
+  ! --    nmax    : maximum nb of iterations (hopefully never reached)
+  ! --    epsilon : accuracy of the numerical resolution
+  ! --    vmax    : v-value above which we use an asymptotic expression for
+  ! ERF(v)
+
+  INTEGER nmax
+  PARAMETER (nmax=10)
+  REAL epsilon, vmax0, vmax(klon)
+  PARAMETER (epsilon=0.02, vmax0=2.0)
+
+  REAL min_mu, min_q
+  PARAMETER (min_mu=1.E-12, min_q=1.E-12)
+
+  INTEGER i, k, n, m
+  REAL mu, qsat, delta
+  REAL sigma1, sigma2, alpha, qconv
+  REAL xconv, xenv
+  REAL cconv, cenv
+  REAL pi, u, v
+  REAL erf
+  REAL sqrtpi, sqrt2
+  ! lconv = true si le calcul a converge (entre autre si qsub < min_q)
+  LOGICAL lconv(klon)
+
+
+  cldf(1:klon, 1:nd) = 0.0 ! cym
+  ratqsc(1:klon, 1:nd) = 0.0
+  ptconv(1:klon, 1:nd) = .FALSE.
+  ! cdir end arraycomb
+
+  pi = acos(-1.)
+  sqrtpi = sqrt(pi)
+  sqrt2 = sqrt(2.)
+
+
+  DO k = 1, nd
+
+  DO i = 1, klon ! vector
+
+      mu = r(i, k)
+      mu = max(mu, min_mu)
+      qsat = rs(i, k)
+      qsat = max(qsat, min_mu)
+      delta = log(mu/qsat)
+      qconv=qtc(i,k)
+      alpha=sigt(i,k)
+
+     IF (qconv<min_q) THEN
+        ptconv(i, k) = .FALSE.
+        ratqsc(i, k) = 0.
+
+        ! Rien on a deja initialise
+
+      ELSE
+    
+      sigma1=0.1*((qconv-mu)**2)**0.5+0.002*mu
+      sigma2=0.1*((qconv-mu)**2)**0.5+0.002*qconv 
+
+!      sigma2=0.09*((qconv-mu)**2)**0.5/(alpha+0.01)**0.5+0.002*qconv 
+!-----------------------------------------------------------------------------------------------------------------
+! Calcul de la couverture nuageuse et de ratqs
+!-----------------------------------------------------------------------------------------------------------------
+
+      xconv=(qsat-qconv)/(sqrt(2.)*sigma2)
+      xenv=(qsat-mu)/(sqrt(2.)*sigma1)
+
+      cconv=0.5*(1.-1.*erf(xconv))
+      cenv=0.5*(1.-1.*erf(xenv)) 
+      cldf(i,k)=alpha*cconv+(1.-1.*alpha)*cenv
+      ratqsc(i,k)= alpha*sigma1+(1.-1.*alpha)*sigma2
+      ptconv(i,k)= .TRUE.
+
+     END IF
+
+  END DO ! vector
+
+
+  END DO  ! K
+
+  RETURN
+  END SUBROUTINE clouds_bigauss
+
+
+
Index: /LMDZ5/branches/testing/libf/phylmd/compar1d.h
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/compar1d.h	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/compar1d.h	(revision 2220)
@@ -3,4 +3,6 @@
 !
       integer :: forcing_type
+      integer :: tend_u,tend_v,tend_w,tend_t,tend_q,tend_rayo
+      real :: nudge_u,nudge_v,nudge_w,nudge_t,nudge_q
       integer :: iflag_nudge
       real :: nat_surf
@@ -32,5 +34,6 @@
      & qsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,    &
      & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
-     & forcing_type,                                                    &
+     & forcing_type,tend_u,tend_v,tend_w,tend_t,tend_q,tend_rayo,       &
+     & nudge_u,nudge_v,nudge_w,nudge_t,nudge_q,                         &
      & iflag_nudge,                                                     &
      & restart,ok_old_disvert
Index: /LMDZ5/branches/testing/libf/phylmd/concvl.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/concvl.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/concvl.F90	(revision 2220)
@@ -15,5 +15,6 @@
                   dd_t, dd_q, lalim_conv, wght_th, &                 ! RomP
                   evap, ep, epmlmMm, eplaMm, &                       ! RomP
-                  wdtrainA, wdtrainM, wght)                          ! RomP+RL
+                  wdtrainA, wdtrainM, wght, qtc, sigt, &
+                  tau_cld_cv, coefw_cld_cv)                           ! RomP+RL, AJ
 !RomP <<<
 ! **************************************************************
@@ -29,4 +30,5 @@
   USE dimphy
   USE infotrac, ONLY: nbtr
+  USE phys_local_var_mod, ONLY: omega
   IMPLICIT NONE
 ! ======================================================================
@@ -135,4 +137,6 @@
   REAL dplcldt(klon), dplcldr(klon)
   REAL qcondc(klon, klev)
+  REAL qtc(klon, klev)
+  REAL sigt(klon, klev)
   REAL wd(klon)
   REAL plim1(klon), plim2(klon), asupmax(klon, klev)
@@ -141,4 +145,5 @@
   REAL sigd(klon)
   REAL zx_t, zdelta, zx_qs, zcor
+  REAL tau_cld_cv, coefw_cld_cv
 
 !   INTEGER iflag_mix
@@ -203,4 +208,5 @@
 !$OMP THREADPRIVATE(itap, igout)
 
+
   include "YOMCST.h"
   include "YOMCST2.h"
@@ -398,5 +404,5 @@
                     t, q, qs, t_wake, q_wake, qs_wake, s_wake, u, v, tra, &
                     em_p, em_ph, &
-                    Ale, Alp, &
+                    Ale, Alp, omega, &
                     em_sig1feed, em_sig2feed, em_wght, &
                     iflag, d_t, d_q, d_u, d_v, d_tra, rain, kbas, ktop, &
@@ -411,5 +417,6 @@
                     da, phi, mp, phi2, d1a, dam, sij, wght, &           ! RomP+RL
                     clw, elij, evap, ep, epmlmMm, eplaMm, &             ! RomP+RL
-                    wdtrainA, wdtrainM)                                 ! RomP
+                    wdtrainA, wdtrainM, qtc, sigt, &
+                    tau_cld_cv, coefw_cld_cv)                           ! RomP,AJ
 !AC!+!RomP+jyg
   END IF
Index: /LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90	(revision 2220)
@@ -15,5 +15,5 @@
        solarlong0,seuil_inversion, &
        fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
-       iflag_cldcon, &
+       iflag_cldth, &
        iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
        ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
@@ -81,5 +81,5 @@
     REAL                 :: bl95_b0, bl95_b1
     real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut,tau_ratqs
-    integer              :: iflag_cldcon
+    integer              :: iflag_cldth
     integer              :: iflag_ratqs
 
@@ -101,4 +101,8 @@
     REAL,SAVE           :: freq_COSP_omp
     real,SAVE           :: fact_cldcon_omp, facttemps_omp,ratqsbas_omp
+    real,SAVE           :: tau_cld_cv_omp, coefw_cld_cv_omp
+    integer,SAVE        :: iflag_cld_cv_omp
+
+
     real,SAVE           :: ratqshaut_omp
     real,SAVE           :: tau_ratqs_omp
@@ -107,5 +111,5 @@
     integer,SAVE        :: iflag_rrtm_omp
     integer,SAVE        :: NSW_omp
-    integer,SAVE        :: iflag_cldcon_omp, ip_ebil_phy_omp
+    integer,SAVE        :: iflag_cldth_omp, ip_ebil_phy_omp
     integer,SAVE        :: iflag_ratqs_omp
 
@@ -132,4 +136,5 @@
     integer,SAVE :: iflag_coupl_omp,iflag_clos_omp,iflag_wake_omp
     integer,SAVE :: iflag_cvl_sigd_omp
+    REAL, SAVE :: coef_clos_ls_omp
     REAL, SAVE :: supcrit1_omp, supcrit2_omp
     INTEGER, SAVE :: iflag_mix_omp
@@ -886,11 +891,45 @@
 
     !
-    !Config Key  = iflag_cldcon 
+    !Config Key  = iflag_cldth 
     !Config Desc =  
     !Config Def  = 1
     !Config Help = 
     !
-    iflag_cldcon_omp = 1
-    call getin('iflag_cldcon',iflag_cldcon_omp)
+    iflag_cldth_omp = 1
+! On lit deux fois avec l'ancien et le nouveau nom 
+! pour assurer une retrocompatiblite.
+! A abandonner un jour
+    call getin('iflag_cldcon',iflag_cldth_omp)
+    call getin('iflag_cldth',iflag_cldth_omp)
+
+    !
+    !Config Key  = iflag_cld_cv
+    !Config Desc =
+    !Config Def  = 1
+    !Config Help =
+    !
+    iflag_cld_cv_omp = 1
+    call getin('iflag_cld_cv',iflag_cld_cv_omp)
+
+    !
+    !Config Key  = tau_cld_cv
+    !Config Desc =
+    !Config Def  = 10.
+    !Config Help =
+    !
+    tau_cld_cv_omp = 10.
+    call getin('tau_cld_cv',tau_cld_cv_omp)
+
+    !
+    !Config Key  = coefw_cld_cv
+    !Config Desc =
+    !Config Def  = 0.1
+    !Config Help =
+    !
+    coefw_cld_cv_omp = 0.1
+    call getin('coefw_cld_cv',coefw_cld_cv_omp)
+
+
+
 
     !
@@ -1343,4 +1382,13 @@
     iflag_clos_omp = 1
     call getin('iflag_clos',iflag_clos_omp)
+    !
+    !Config Key  = coef_clos_ls
+    !Config Desc =  
+    !Config Def  = 0
+    !Config Help = 
+    !
+    coef_clos_ls_omp = 0.
+    call getin('coef_clos_ls',coef_clos_ls_omp)
+
     !
     !Config Key  = iflag_cvl_sigd
@@ -1925,5 +1973,8 @@
     iflag_rrtm = iflag_rrtm_omp
     NSW = NSW_omp
-    iflag_cldcon = iflag_cldcon_omp
+    iflag_cldth = iflag_cldth_omp
+    iflag_cld_cv = iflag_cld_cv_omp
+    tau_cld_cv = tau_cld_cv_omp
+    coefw_cld_cv = coefw_cld_cv_omp
     iflag_ratqs = iflag_ratqs_omp
     ip_ebil_phy = ip_ebil_phy_omp
@@ -1946,4 +1997,5 @@
     iflag_clos = iflag_clos_omp
     iflag_wake = iflag_wake_omp
+    coef_clos_ls = coef_clos_ls_omp
     alp_offset = alp_offset_omp
     iflag_cvl_sigd = iflag_cvl_sigd_omp
@@ -2076,5 +2128,8 @@
     write(lunout,*)' reevap_ice = ', reevap_ice
     write(lunout,*)' iflag_pdf = ', iflag_pdf
-    write(lunout,*)' iflag_cldcon = ', iflag_cldcon
+    write(lunout,*)' iflag_cldth = ', iflag_cldth
+    write(lunout,*)' iflag_cld_cv = ', iflag_cld_cv
+    write(lunout,*)' tau_cld_cv = ', tau_cld_cv
+    write(lunout,*)' coefw_cld_cv = ', coefw_cld_cv
     write(lunout,*)' iflag_radia = ', iflag_radia
     write(lunout,*)' iflag_rrtm = ', iflag_rrtm
@@ -2135,4 +2190,5 @@
     write(lunout,*)' iflag_thermals_closure = ', iflag_thermals_closure
     write(lunout,*)' iflag_clos = ', iflag_clos
+    write(lunout,*)' coef_clos_ls = ', coef_clos_ls
     write(lunout,*)' type_run = ',type_run 
     write(lunout,*)' ok_cosp = ',ok_cosp
Index: /LMDZ5/branches/testing/libf/phylmd/convect3.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/convect3.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/convect3.F90	(revision 2220)
@@ -17,5 +17,5 @@
   USE dimphy
   USE infotrac, ONLY: nbtr
-
+  IMPLICIT NONE
   include "dimensions.h"
   INTEGER na
@@ -73,5 +73,15 @@
 
 
-
+  REAL :: cpv,cl,cpvmcl,eps,alv0,rdcp,pbcrit,ptcrit,sigd,spfac
+  REAL :: tau,beta,alpha,dtcrit,dtovsh,ahm,rm,um,vm,dphinv
+  REAL :: a2,x,tvx,tvy,plcl,pden,dpbase,tvpbase,tvbase,tdif
+  REAL :: ath1,ath,delti,deltap,dcape,dlnp,sigold,dtmin,fac,w
+  REAL :: amu,rti,cpd,bf2,anum,denom,dei,altem,cwat,stemp,qp
+  REAL :: scrit,alt,smax,asij,wgh,sjmax,sjmin,smid,delp,delm
+  REAL :: asum,bsum,csum,wflux,tinv,wdtrain,awat,afac,afac1,afac2
+  REAL :: bfac,pr1,pr2,sigt,b6,c6,revap,tevap,delth,amfac,amp2
+  REAL :: xf,tf,af,bf,fac2,ur,sru,d,ampmax,dpinv,am,amde,cpinv
+  REAL :: amp1,ad,rat,ax,bx,cx,dx,ex,dsum
+  INTEGER :: nk,i,j,nopt,jn,k,im,jm,n
 
   REAL dnwd0(nd) !  precipitation driven unsaturated downdraft flux
Index: /LMDZ5/branches/testing/libf/phylmd/cv3_inicp.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cv3_inicp.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/cv3_inicp.F90	(revision 2220)
@@ -9,5 +9,5 @@
   ! modified by :                                               *
   ! **************************************************************
-
+  IMPLICIT NONE
   include "YOMCST2.h"
 
@@ -19,4 +19,5 @@
 
   REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f
+  REAL :: sumcoef,sigma,aire,pdf,mu,df,ff
 
   qcoef1(f) = tanh(f/gammas)
Index: /LMDZ5/branches/testing/libf/phylmd/cv3_inip.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cv3_inip.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/cv3_inip.F90	(revision 2220)
@@ -1,19 +1,47 @@
 SUBROUTINE cv3_inip()
-  ! **************************************************************
-  ! *
-  ! CV3_INIP Lecture des choix de lois de probabilité de mélange*
-  ! et calcul de leurs coefficients normalisés.        *
-  ! *
-  ! written by   : Jean-Yves Grandpeix, 06/06/2006, 19.39.27    *
-  ! modified by :                                               *
-  ! **************************************************************
+  ! *******************************************************************
+  ! *                                                                  *
+  ! CV3_INIP Input = choice of mixing probability laws                 *
+  !          Output = normalized coefficients of the probability laws. *
+  ! *                                                                  *
+  ! written by   : Jean-Yves Grandpeix, 06/06/2006, 19.39.27           *
+  ! modified by :                                                      *
+  ! *******************************************************************
+!
+!----------------------------------------------
+!  INPUT (from Common YOMCST2 in "YOMCST2.h") : 
+! iflag_mix
+! gammas    
+! alphas
+! betas
+! Fmax
+! scut
+!
+!----------------------------------------------
+!  INPUT/OUTPUT (from and to Common YOMCST2 in "YOMCST2.h") : 
+! qqa1
+! qqa2
+!
+!----------------------------------------------
+!  OUTPUT (to Common YOMCST2 in "YOMCST2.h") :
+! Qcoef1max
+! Qcoef2max
+!
+!----------------------------------------------
+
+  IMPLICIT NONE
 
   include "YOMCST2.h"
 
-  ! INTEGER iflag_mix
   include 'iniprint.h'
 
-  CHARACTER (LEN=20) :: modname = 'cv3_inip'
-  CHARACTER (LEN=80) :: abort_message
+!----------------------------------------------
+! Local variables :
+   CHARACTER (LEN=20) :: modname = 'cv3_inip'
+   CHARACTER (LEN=80) :: abort_message
+
+   REAL               :: sumcoef
+   REAL               :: sigma, aire, pdf, mu, df
+   REAL               :: ff
 
 
Index: /LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90	(revision 2220)
@@ -1,4 +1,5 @@
 
 ! $Id$
+
 
 
@@ -2763,5 +2764,5 @@
                      cbmf, upwd, dnwd, dnwd0, ma, mip, &
                      tls, tps, qcondc, wd, &
-                     ftd, fqd)
+                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
 
   IMPLICIT NONE
@@ -2800,4 +2801,5 @@
 !input/output:
       INTEGER iflag(nloc)
+      REAL,INTENT(in) :: tau_cld_cv, coefw_cld_cv
 !
 !outputs:
@@ -2811,4 +2813,5 @@
       REAL tls(nloc, nd), tps(nloc, nd)
       REAL qcondc(nloc, nd) ! cld
+      REAL qtc(nloc,nd), sigt(nloc,nd) ! cld
       REAL wd(nloc) ! gust
       REAL cbmf(nloc)
@@ -2830,5 +2833,6 @@
       REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
       REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
-
+      REAL sument(nloc), sigment(nloc,nd), qtment(nloc,nd) ! cld
+      REAL qnk(nloc)
       REAL sumdq !jyg
 ! 
@@ -2861,4 +2865,8 @@
       qcondc(il, i) = 0.0 ! cld
       qcond(il, i) = 0.0 ! cld
+      qtc(il, i) = 0.0 ! cld
+      qtment(il, i) = 0.0 ! cld
+      sigment(il, i) = 0.0 ! cld
+      sigt(il, i) = 0.0 ! cld
       nqcond(il, i) = 0.0 ! cld
     END DO
@@ -3234,5 +3242,6 @@
 ! (saturated updrafts resulting from mixing)                                   ! cld
           qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il))                ! cld
-          nqcond(il, i) = nqcond(il, i) + 1. ! cld
+          qtment(il, i) = qtment(il, i) + qent(il,k,i)                         ! cld
+          nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
         END IF ! i
       END DO
@@ -3310,5 +3319,6 @@
 ! (saturated downdrafts resulting from mixing)                                 ! cld
           qcond(il, i) = qcond(il, i) + elij(il, k, i)                         ! cld
-          nqcond(il, i) = nqcond(il, i) + 1. ! cld
+          qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
+          nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
         END IF ! cld
       END DO ! cld
@@ -3319,4 +3329,5 @@
       IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
         qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
+        qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
         nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
       END IF                                                                   ! cld
@@ -3326,4 +3337,5 @@
       IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN            ! cld
         qcond(il, i) = qcond(il, i)/nqcond(il, i)                              ! cld
+        qtment(il, i) = qtment(il,i)/nqcond(il, i)                             ! cld
       END IF                                                                   ! cld
     END DO
@@ -3788,18 +3800,53 @@
       END IF                                                         ! cld
     END DO                                                           ! cld
-  END DO                                                             ! cld
-
-  DO i = 1, nl                                                       ! cld
+  END DO  
+                                                           ! cld
+  DO i = 1, nl  
+
+! 14/01/15 AJ je remets les parties manquantes cf JYG
+! Initialize sument to 0
+
+    DO il = 1,ncum
+     sument(il) = 0.
+    ENDDO
+
+! Sum mixed mass fluxes in sument
+
+    DO k = 1,nl
+      DO il = 1,ncum
+        IF  (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN   ! cld
+          sument(il) =sument(il) + abs(ment(il,k,i))
+        ENDIF
+      ENDDO     ! il
+    ENDDO       ! k
+
+! 14/01/15 AJ delta n'a rien à faire là...                                                 
     DO il = 1, ncum                                                  ! cld
       IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
-        siga(il, i) = mac(il, i)/wa(il, i) &                         ! cld
-        *rrd*tvp(il, i)/p(il, i)/100./delta                          ! cld
+        siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) &          ! cld
+        *rrd*tvp(il, i)/p(il, i)/100.                                ! cld
+
       siga(il, i) = min(siga(il,i), 1.0)                             ! cld
-! IM cf. FH                                                          
+
+! IM cf. FH 
+! 14/01/15 AJ ne correspond pas à ce qui a été codé par JYG et SB           
+                                                         
       IF (iflag_clw==0) THEN                                         ! cld
         qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) &       ! cld
           +(1.-siga(il,i))*qcond(il, i)                              ! cld
+
+
+        sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1))    ! cld
+        sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i))  ! cld
+        qtc(il, i) = (siga(il,i)*qnk(il)+sigment(il,i)*qtment(il,i)) & ! cld
+                     /(siga(il,i)+sigment(il,i))                     ! cld
+        sigt(il,i) = sigment(il, i) + siga(il, i)
+
+!        qtc(il, i) = siga(il,i)*qnk(il)+(1.-siga(il,i))*qtment(il,i) ! cld
+!     print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)  
+      		
       ELSE IF (iflag_clw==1) THEN                                    ! cld
         qcondc(il, i) = qcond(il, i)                                 ! cld
+        qtc(il,i) = qtment(il,i)                                     ! cld
       END IF                                                         ! cld
 
Index: /LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F90	(revision 2220)
@@ -4,8 +4,8 @@
     th1_wake, tra1, h1, lv1, lf1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
     h1_wake, lv1_wake, lf1_wake, cpn1_wake, tv1_wake, sig1, w01, ptop21, &
-    ale1, alp1, iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, hnk, unk, vnk, &
+    ale1, alp1, omega1, iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, hnk, unk, vnk, &
     wghti, pbase, buoybase, t, q, qs, t_wake, q_wake, qs_wake, s_wake, u, v, &
     gz, th, th_wake, tra, h, lv, lf, cpn, p, ph, tv, tp, tvp, clw, h_wake, &
-    lv_wake, lf_wake, cpn_wake, tv_wake, sig, w0, ptop2, ale, alp)
+    lv_wake, lf_wake, cpn_wake, tv_wake, sig, w0, ptop2, ale, alp, omega)
   ! **************************************************************
   ! *
@@ -40,4 +40,5 @@
   REAL sig1(len, nd), w01(len, nd), ptop21(len)
   REAL ale1(len), alp1(len)
+  REAL omega1(len,nd)
 
   ! outputs:
@@ -60,4 +61,5 @@
   REAL sig(len, nd), w0(len, nd), ptop2(len)
   REAL ale(len), alp(len)
+  REAL omega(len,nd)
 
   ! local variables:
@@ -102,4 +104,5 @@
         sig(nn, k) = sig1(i, k)
         w0(nn, k) = w01(i, k)
+        omega(nn, k) = omega1(i, k)
       END IF
     END DO
Index: /LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F90	(revision 2220)
@@ -6,4 +6,5 @@
     , clw, elij, evap, ep, epmlmmm, eplamm & ! RomP
     , wdtraina, wdtrainm &         ! RomP
+    , qtc, sigt          &
 
     , iflag1, kbas1, ktop1, precip1, cbmf1, plcl1, plfc1, wbeff1, sig1, w01, &
@@ -13,5 +14,6 @@
     , da1, phi1, mp1, phi21, d1a1, dam1, sigij1 & ! RomP+AC+jyg
     , clw1, elij1, evap1, ep1, epmlmmm1, eplamm1 & ! RomP
-    , wdtraina1, wdtrainm1) ! RomP
+    , wdtraina1, wdtrainm1 & ! RomP
+    , qtc1, sigt1)
 
   ! **************************************************************
@@ -56,4 +58,5 @@
   REAL evap(nloc, nd), ep(nloc, nd) !RomP
   REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd) !RomP+jyg
+  REAL qtc(nloc, nd), sigt(nloc, nd) !RomP
   REAL wdtraina(nloc, nd), wdtrainm(nloc, nd) !RomP
 
@@ -84,4 +87,5 @@
   REAL evap1(len, nd), ep1(len, nd) !RomP
   REAL epmlmmm1(len, nd, nd), eplamm1(len, nd) !RomP+jyg
+  REAL qtc1(len, nd), sigt1(len, nd) !RomP
   REAL wdtraina1(len, nd), wdtrainm1(len, nd) !RomP
 
@@ -141,4 +145,6 @@
       wdtraina1(idcum(i), k) = wdtraina(i, k) !RomP
       wdtrainm1(idcum(i), k) = wdtrainm(i, k) !RomP
+      qtc1(idcum(i), k) = qtc(i, k)
+      sigt1(idcum(i), k) = sigt(i, k)
 
     END DO
Index: /LMDZ5/branches/testing/libf/phylmd/cv3p1_closure.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cv3p1_closure.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/cv3p1_closure.F90	(revision 2220)
@@ -3,5 +3,5 @@
 
 SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
-    tvp, buoy, supmax, ok_inhib, ale, alp, sig, w0, ptop2, cape, cin, m, &
+    tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
     iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmf, plfc, &
     wbeff)
@@ -37,4 +37,5 @@
   LOGICAL ok_inhib ! enable convection inhibition by dryness
   REAL ale(nloc), alp(nloc)
+  REAL omega(nloc,nd)
 
   ! input/output:
@@ -52,5 +53,5 @@
 
   ! local variables:
-  INTEGER il, i, j, k, icbmax, i0(nloc)
+  INTEGER il, i, j, k, icbmax, i0(nloc), klfc
   REAL deltap, fac, w, amu
   REAL rhodp
@@ -523,4 +524,13 @@
   END DO
 
+!CR:Compute k at plfc
+  DO k=1,nl
+     DO il=1,ncum
+        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
+           klfc=k
+        endif
+     ENDDO
+  ENDDO
+!RC
 
   DO il = 1, ncum
@@ -528,4 +538,10 @@
     ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
     cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
+!CR: Add large-scale component to the mass-flux
+!encore connu sous le nom "Experience du tube de dentifrice"
+    if (coef_clos_ls.gt.0.) then 
+       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc))
+    endif
+!RC
     IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
       WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &
Index: /LMDZ5/branches/testing/libf/phylmd/cva_driver.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/cva_driver.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/cva_driver.F90	(revision 2220)
@@ -7,5 +7,5 @@
                       u1, v1, tra1, &
                       p1, ph1, &
-                      Ale1, Alp1, &
+                      Ale1, Alp1, omega1, &
                       sig1feed1, sig2feed1, wght1, &
                       iflag1, ft1, fq1, fu1, fv1, ftra1, &
@@ -24,5 +24,6 @@
                       da1, phi1, mp1, phi21, d1a1, dam1, sigij1, wghti1, & ! RomP, RL
                       clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &        ! RomP, RL
-                      wdtrainA1, wdtrainM1)                                ! RomP
+                      wdtrainA1, wdtrainM1, qtc1, sigt1, tau_cld_cv, &
+                      coefw_cld_cv)                                        ! RomP, AJ
 ! **************************************************************
 ! *
@@ -55,4 +56,6 @@
 ! iflag_ice_thermo Integer        Input        accounting for ice thermodynamics (0/1)
 ! iflag_clos    Integer        Input        version of closure (0/1)
+! tau_cld_cv    Real           Input        characteristic time of dissipation of mixing fluxes
+! coefw_cld_cv  Real           Input        coefficient for updraft velocity in convection
 ! ok_conserv_q  Logical        Input        when true corrections for water conservation are swtiched on
 ! delt          Real           Input        time step
@@ -119,5 +122,6 @@
 ! phi1          Real           Output     used in tracer transport (cvltr)
 ! mp1           Real           Output     used in tracer transport (cvltr)
-                                         
+! qtc1          Real           Output     specific humidity in convection
+! sigt1         Real           Output     surface fraction in adiabatic updrafts                                         
 ! phi21         Real           Output     used in tracer transport (cvltr)
                                          
@@ -163,4 +167,6 @@
   INTEGER iflag_clos
   LOGICAL ok_conserv_q
+  REAL tau_cld_cv
+  REAL coefw_cld_cv
   REAL delt
   REAL t1(len, nd)
@@ -178,4 +184,5 @@
   REAL Ale1(len)
   REAL Alp1(len)
+  REAL omega1(len,nd)
   REAL sig1feed1 ! pressure at lower bound of feeding layer
   REAL sig2feed1 ! pressure at upper bound of feeding layer
@@ -225,4 +232,7 @@
   REAL asupmaxmin1(len)
   INTEGER lalim_conv(len)
+  REAL qtc1(len, nd)         ! cld
+  REAL sigt1(len, nd)        ! cld
+
 ! RomP >>>
   REAL wdtrainA1(len, nd), wdtrainM1(len, nd)
@@ -455,4 +465,5 @@
   REAL supmax(nloc, klev)
   REAL Ale(nloc), Alp(nloc), coef_clos(nloc)
+  REAL omega(nloc,klev)
   REAL sigd(nloc)
 ! real mp(nloc,klev), qp(nloc,klev), up(nloc,klev), vp(nloc,klev)
@@ -487,4 +498,7 @@
   REAL hnk(nloc), unk(nloc), vnk(nloc)
 
+  REAL qtc(nloc, klev)         ! cld
+  REAL sigt(nloc, klev)        ! cld
+ 
 ! RomP >>>
   REAL wdtrainA(nloc, klev), wdtrainM(nloc, klev)
@@ -593,4 +607,6 @@
 
 ! RomP >>>
+  sigt1(:, :) = 0.
+  qtc1(:, :) = 0.
   wdtrainA1(:, :) = 0.
   wdtrainM1(:, :) = 0.
@@ -776,5 +792,5 @@
                          h1_wake, lv1_wake, lf1_wake, cpn1_wake, tv1_wake, &
                          sig1, w01, ptop21, &
-                         Ale1, Alp1, &
+                         Ale1, Alp1, omega1, &
                          iflag, nk, icb, icbs, &
                          plcl, tnk, qnk, gznk, hnk, unk, vnk, &
@@ -786,5 +802,5 @@
                          h_wake, lv_wake, lf_wake, cpn_wake, tv_wake, &
                          sig, w0, ptop2, &
-                         Ale, Alp)
+                         Ale, Alp, omega)
 
 ! print*,'tv ',tv
@@ -877,5 +893,5 @@
         CALL cv3p1_closure(nloc, ncum, nd, icb, inb, &         ! na->nd
                            pbase, plcl, p, ph, tv, tvp, buoy, &
-                           supmax, ok_inhib, Ale, Alp, &
+                           supmax, ok_inhib, Ale, Alp, omega, &
                            sig, w0, ptop2, cape, cin, m, iflag, coef_clos, &
                            Plim1, plim2, asupmax, supmax0, &
@@ -978,5 +994,5 @@
                      cbmf, upwd, dnwd, dnwd0, ma, mip, &
                      tls, tps, qcondc, wd, &
-                     ftd, fqd)
+                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
     END IF
 
@@ -1032,4 +1048,5 @@
                            clw, elij, evap, ep, epmlmMm, eplaMm, &       ! RomP
                            wdtrainA, wdtrainM, &                         ! RomP
+                           qtc, sigt, &
                            iflag1, kbas1, ktop1, &
                            precip1, cbmf1, plcl1, plfc1, wbeff1, sig1, w01, ptop21, &
@@ -1043,5 +1060,6 @@
                            da1, phi1, mp1, phi21, d1a1, dam1, sigij1,  & ! RomP
                            clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, & ! RomP
-                           wdtrainA1, wdtrainM1)                         ! RomP
+                           wdtrainA1, wdtrainM1,                       & ! RomP
+                           qtc1, sigt1)
     END IF
 
Index: /LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90	(revision 2220)
@@ -8,5 +8,5 @@
      frac_impa, frac_nucl, beta,                        &
      prfl, psfl, rhcl, zqta, fraca,                     &
-     ztv, zpspsk, ztla, zthl, iflag_cldcon,             &
+     ztv, zpspsk, ztla, zthl, iflag_cldth,             &
      iflag_ice_thermo)
 
@@ -82,5 +82,5 @@
   INTEGER ninter ! sous-intervals pour la precipitation
   INTEGER ncoreczq  
-  INTEGER iflag_cldcon
+  INTEGER iflag_cldth
   INTEGER iflag_ice_thermo
   PARAMETER (ninter=5)
@@ -545,5 +545,5 @@
            enddo
 
-           if (iflag_cldcon>=5) then
+           if (iflag_cldth>=5) then
 
               call cloudth(klon,klev,k,ztv, &
@@ -559,11 +559,11 @@
            endif
 
-           if (iflag_cldcon <= 4) then
+           if (iflag_cldth <= 4) then
               lognormale = .true.
-           elseif (iflag_cldcon >= 6) then
+           elseif (iflag_cldth >= 6) then
               ! lognormale en l'absence des thermiques
               lognormale = fraca(:,k) < 1e-10
            else
-              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
+              ! Dans le cas iflag_cldth=5, on prend systématiquement la
               ! bi-gaussienne
               lognormale = .false.
Index: /LMDZ5/branches/testing/libf/phylmd/hines_gwd.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/hines_gwd.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/hines_gwd.F90	(revision 2220)
@@ -625,5 +625,4 @@
     mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, il2, &
     nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod)
-
   ! Smooth cutoff wavenumbers and total rms velocity in the vertical
   ! direction NSMAX times, using FLUX_U as temporary work array.
@@ -715,5 +714,5 @@
     i_alpha, mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, &
     il2, nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod)
-
+  IMPLICIT NONE
   ! This routine calculates the cutoff vertical wavenumber and velocity
   ! variances on a longitude by altitude grid for the Hines' Doppler
@@ -766,5 +765,5 @@
 
   INTEGER naz, levbot, levtop, il1, il2, nlons, nlevs, nazmth
-  REAL slope, kstar(nlons), f1, f2, f3
+  REAL slope, kstar(nlons), f1, f2, f3, f2mfac
   REAL m_alpha(nlons, nlevs, nazmth)
   REAL sigma_alpha(nlons, nlevs, nazmth)
@@ -938,5 +937,5 @@
 SUBROUTINE hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, &
     nlons, nlevs, nazmth)
-
+    IMPLICIT NONE
   ! This routine calculates the azimuthal horizontal background wind
   ! components
@@ -1034,5 +1033,5 @@
     m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, nlons, &
     nlevs, nazmth, lorms)
-
+    IMPLICIT NONE
   ! Calculate zonal and meridional components of the vertical flux
   ! of horizontal momentum and corresponding wave drag (force per unit mass)
@@ -1089,5 +1088,5 @@
   ! Internal variables.
 
-  INTEGER i, l, lev1p, lev2m
+  INTEGER i, l, lev1p, lev2m, lev2p
   REAL cos45, prod2, prod4, prod6, prod8, dendz, dendz2
   DATA cos45/0.7071068/
@@ -1234,5 +1233,5 @@
     bvfreq, density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, &
     naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth)
-
+  IMPLICIT NONE
   ! This routine calculates the gravity wave induced heating and
   ! diffusion coefficient on a longitude by altitude grid for
@@ -1355,5 +1354,5 @@
 SUBROUTINE hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, lev, il1, &
     il2, nlons, nlevs, nazmth)
-
+  IMPLICIT NONE
   ! This routine calculates the total rms and azimuthal rms horizontal
   ! velocities at a given level on a longitude by altitude grid for
@@ -1450,5 +1449,5 @@
 SUBROUTINE hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, lev, &
     il1, il2, nlons, nlevs, nazmth, lorms)
-
+  IMPLICIT NONE
   ! This routine calculates the vertical wavenumber integral
   ! for a single vertical level at each azimuth on a longitude grid
@@ -1643,5 +1642,5 @@
     alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, nlons, &
     nazmth, coslat)
-
+  IMPLICIT NONE
   ! This routine specifies various parameters needed for the
   ! the Hines' Doppler spread gravity wave drag parameterization scheme.
@@ -1772,5 +1771,5 @@
     sigma_alpha, v_alpha, m_alpha, iu_print, iv_print, nmessg, ilprt1, &
     ilprt2, levprt1, levprt2, naz, nlons, nlevs, nazmth)
-
+  IMPLICIT NONE
   ! Print out altitude profiles of various quantities from
   ! Hines' Doppler spread gravity wave drag parameterization scheme.
@@ -1864,5 +1863,5 @@
 SUBROUTINE hines_exp(data, data_zmax, alt, alt_exp, iorder, il1, il2, lev1, &
     lev2, nlons, nlevs)
-
+  IMPLICIT NONE
   ! This routine exponentially damps a longitude by altitude array
   ! of data above a specified altitude.
@@ -1941,5 +1940,5 @@
 SUBROUTINE vert_smooth(data, work, coeff, nsmooth, il1, il2, lev1, lev2, &
     nlons, nlevs)
-
+  IMPLICIT NONE
   ! Smooth a longitude by altitude array in the vertical over a
   ! specified number of levels using a three point smoother.
Index: /LMDZ5/branches/testing/libf/phylmd/ini_wake.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/ini_wake.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/ini_wake.F90	(revision 2220)
@@ -4,4 +4,5 @@
 SUBROUTINE ini_wake(wape, fip, it_wape_prescr, wape_prescr, fip_prescr, &
     alp_bl_prescr, ale_bl_prescr)
+  IMPLICIT NONE
   ! **************************************************************
   ! *
@@ -39,7 +40,10 @@
   include 'iniprint.h'
   ! declarations
+  REAL wape, fip, wape_prescr, fip_prescr
+  INTEGER it_wape_prescr
   REAL ale_bl_prescr
   REAL alp_bl_prescr
   REAL it
+  REAL w,f,alebl,alpbl
 
   ! FH A mettre si besoin dans physiq.def
Index: /LMDZ5/branches/testing/libf/phylmd/limit_read_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/limit_read_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/limit_read_mod.F90	(revision 2220)
@@ -93,5 +93,6 @@
 !****************************************************************************************
 
-    IF (type_ocean == 'couple') THEN
+IF (type_ocean == 'couple'.OR. &
+         (type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
        ! limit.nc has not yet been read. Do it now!
        CALL limit_read_tot(itime, dtime, jour, is_modified)
Index: /LMDZ5/branches/testing/libf/phylmd/limit_slab.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/limit_slab.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/limit_slab.F90	(revision 2220)
@@ -1,5 +1,5 @@
 ! $Header$
 
-SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, diff_sst)
+SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
 
   USE dimphy
@@ -20,10 +20,11 @@
   INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
   REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
-  REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst
+  REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst, diff_siv
 
 ! Locals variables with attribute SAVE
 !****************************************************************************************
   REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, diff_sst_save
-!$OMP THREADPRIVATE(bils_save, diff_sst_save)
+  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: diff_siv_save
+!$OMP THREADPRIVATE(bils_save, diff_sst_save, diff_siv_save)
 
 ! Locals variables
@@ -33,6 +34,7 @@
   INTEGER, DIMENSION(2)    :: start, epais 
   REAL, DIMENSION(klon_glo):: bils_glo, sst_l_glo, sst_lp1_glo, diff_sst_glo
+  REAL, DIMENSION(klon_glo):: siv_l_glo, siv_lp1_glo, diff_siv_glo
   CHARACTER (len = 20)     :: modname = 'limit_slab'
-  LOGICAL                  :: read_bils,read_sst
+  LOGICAL                  :: read_bils,read_sst,read_siv
 
 ! End declaration
@@ -49,4 +51,5 @@
         read_bils=.TRUE.
         read_sst=.TRUE.
+        read_siv=.TRUE.
         
         ierr = NF90_OPEN ('limit_slab.nc', NF90_NOWRITE, nid)
@@ -54,4 +57,5 @@
             read_bils=.FALSE.
             read_sst=.FALSE.
+            read_siv=.FALSE.
         ELSE ! read file
         
@@ -63,5 +67,5 @@
 
 !****************************************************************************************
-! 2) Read bils and SST tendency
+! 2) Read bils and SST/ ice volume tendency
 !
 !****************************************************************************************
@@ -89,4 +93,19 @@
         END IF
 
+! Read siv_glo for this day
+        ierr = NF90_INQ_VARID(nid, 'SICV', nvarid)
+        IF (ierr /= NF90_NOERR)  THEN
+            read_siv=.FALSE.
+        ELSE
+            start(2) = jour
+            ierr = NF90_GET_VAR(nid,nvarid,siv_l_glo,start,epais)
+            IF (ierr /= NF90_NOERR) read_siv=.FALSE.
+! Read siv_glo for one day ahead
+            start(2) = jour + 1
+            IF (start(2) > 360) start(2)=1
+            ierr = NF90_GET_VAR(nid,nvarid,siv_lp1_glo,start,epais)
+            IF (ierr /= NF90_NOERR) read_siv=.FALSE.
+        END IF
+
 !****************************************************************************************
 ! 5) Close file and distribute variables to all processus
@@ -98,12 +117,14 @@
         IF (read_sst) THEN
 ! Calculate difference in temperature between this day and one ahead
-!            DO i=1, klon_glo-1
-!               diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
-!            END DO
-!            diff_sst_glo(klon_glo) = sst_lp1_glo(klon_glo) - sst_l_glo(1)
             DO i=1, klon_glo
                diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
             END DO
         END IF !read_sst
+        IF (read_siv) THEN
+! Calculate difference in temperature between this day and one ahead
+            DO i=1, klon_glo
+               diff_siv_glo(i) = siv_lp1_glo(i) - siv_l_glo(i)
+            END DO
+        END IF !read_siv
      ENDIF ! is_mpi_root
 
@@ -111,9 +132,9 @@
        
      IF (.NOT. ALLOCATED(bils_save)) THEN
-        ALLOCATE(bils_save(klon), diff_sst_save(klon), stat=ierr)
+        ALLOCATE(bils_save(klon), diff_sst_save(klon), diff_siv_save(klon), stat=ierr)
         IF (ierr /= 0) CALL abort_gcm('limit_slab', 'pb in allocation',1)
      END IF
 
-! Giveddefault values if needed
+! Give default values if needed
      IF (read_bils) THEN
          CALL Scatter(bils_glo, bils_save)
@@ -126,4 +147,9 @@
          diff_sst_save(:)=0.
      END IF
+     IF (read_siv) THEN
+         CALL Scatter(diff_siv_glo, diff_siv_save)
+     ELSE
+         diff_siv_save(:)=0.
+     END IF
      
   ENDIF ! time to read
@@ -131,4 +157,6 @@
   lmt_bils(:) = bils_save(:)
   diff_sst(:) = diff_sst_save(:)
+  diff_siv(:) = diff_siv_save(:)
+
   
 END SUBROUTINE limit_slab
Index: /LMDZ5/branches/testing/libf/phylmd/lmdz1d.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/lmdz1d.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/lmdz1d.F90	(revision 2220)
@@ -22,4 +22,5 @@
       USE phyaqua_mod
       USE mod_1D_cases_read
+      USE mod_1D_amma_read
 
       implicit none
@@ -118,4 +119,5 @@
         logical :: forcing_astex   = .false.
         logical :: forcing_fire    = .false.
+        logical :: forcing_case    = .false.
         integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
 !                                                            (cf read_tsurf1d.F)
@@ -171,4 +173,5 @@
       real :: dt_cooling(llm),d_th_adv(llm),d_t_nudge(llm)
       real :: d_u_nudge(llm),d_v_nudge(llm)
+      real :: du_adv(llm),dv_adv(llm)
       real :: alpha
       real :: ttt
@@ -285,4 +288,7 @@
 !             Different stages: soil model alone, atm. model alone
 !             then both models coupled
+!forcing_type = 10 ==> forcing_case = .true.
+!             initial profiles and large scale forcings in cas.nc
+!             LS convergence, omega and SST imposed from CINDY-DYNAMO files 
 !forcing_type = 40 ==> forcing_GCSSold = .true.
 !             initial profile from GCSS file
@@ -321,4 +327,6 @@
       elseif (forcing_type .eq.7) THEN
        forcing_dice = .true.
+      elseif (forcing_type .eq.10) THEN
+       forcing_case = .true.
       elseif (forcing_type .eq.40) THEN
        forcing_GCSSold = .true.
@@ -428,4 +436,11 @@
      & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice             & 
      & ,day_ju_ini_dice)
+      ELSEIF (forcing_type .eq.10) THEN
+! Convert the initial date to Julian day
+      print*,'time cindy',year_ini_cas,mth_ini_cas,day_ini_cas
+      call ymds2ju                                                         &
+     & (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas &
+     & ,day_ju_ini_cas)
+      print*,'time cindy 2',day_ju_ini_cas
       ELSEIF (forcing_type .eq.59) THEN
 ! Convert the initial date of Sandu case to Julian day
@@ -943,9 +958,9 @@
         u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
      &              du_phys(1:mxcalc)                                       &
-     &             +du_age(1:mxcalc)                                        &
+     &             +du_age(1:mxcalc)+du_adv(1:mxcalc)                       &
      &             +d_u_nudge(1:mxcalc) )            
         v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
      &              dv_phys(1:mxcalc)                                       &
-     &             +dv_age(1:mxcalc)                                        &
+     &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
      &             +d_v_nudge(1:mxcalc) )
         q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
Index: /LMDZ5/branches/testing/libf/phylmd/mod_1D_amma_read.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/mod_1D_amma_read.F90	(revision 2220)
+++ /LMDZ5/branches/testing/libf/phylmd/mod_1D_amma_read.F90	(revision 2220)
@@ -0,0 +1,500 @@
+MODULE mod_1D_amma_read
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!Declarations specifiques au cas AMMA
+        character*80 :: fich_amma
+! Option du cas AMMA ou on impose la discretisation verticale (Ap,Bp)
+        integer nlev_amma, nt_amma
+
+
+        integer year_ini_amma, day_ini_amma, mth_ini_amma
+        real heure_ini_amma
+        real day_ju_ini_amma   ! Julian day of amma first day
+        parameter (year_ini_amma=2006)
+        parameter (mth_ini_amma=7)
+        parameter (day_ini_amma=10)  ! 10 = 10Juil2006
+        parameter (heure_ini_amma=0.) !0h en secondes
+        real dt_amma
+        parameter (dt_amma=1800.)
+
+!profils initiaux:
+        real, allocatable::  plev_amma(:)
+
+        real, allocatable::  z_amma(:)
+        real, allocatable::  th_amma(:),q_amma(:)
+        real, allocatable::  u_amma(:)
+        real, allocatable::  v_amma(:)
+
+        real, allocatable::  th_ammai(:),q_ammai(:)
+        real, allocatable::  u_ammai(:)
+        real, allocatable::  v_ammai(:)
+        real, allocatable::  vitw_ammai(:)
+        real, allocatable::  ht_ammai(:)
+        real, allocatable::  hq_ammai(:)
+        real, allocatable::  vt_ammai(:)
+        real, allocatable::  vq_ammai(:)
+
+!forcings
+        real, allocatable::  ht_amma(:,:)
+        real, allocatable::  hq_amma(:,:)
+        real, allocatable::  vitw_amma(:,:)
+        real, allocatable::  lat_amma(:),sens_amma(:)
+
+!champs interpoles
+        real, allocatable::  vitw_profamma(:)
+        real, allocatable::  ht_profamma(:)
+        real, allocatable::  hq_profamma(:)
+        real lat_profamma,sens_profamma
+        real, allocatable::  vt_profamma(:)
+        real, allocatable::  vq_profamma(:)
+        real, allocatable::  th_profamma(:)
+        real, allocatable::  q_profamma(:)
+        real, allocatable::  u_profamma(:)
+        real, allocatable::  v_profamma(:)
+
+
+CONTAINS
+
+SUBROUTINE read_1D_cases
+      implicit none
+
+#include "netcdf.inc"
+
+      INTEGER nid,rid,ierr
+
+      fich_amma='amma.nc'
+      print*,'fich_amma ',fich_amma
+      ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
+      print*,'fich_amma,NF_NOWRITE,nid ',fich_amma,NF_NOWRITE,nid
+      if (ierr.NE.NF_NOERR) then
+         write(*,*) 'ERROR: GROS Pb opening forcings nc file '
+         write(*,*) NF_STRERROR(ierr)
+         stop ""
+      endif
+!.......................................................................
+      ierr=NF_INQ_DIMID(nid,'lev',rid)
+      IF (ierr.NE.NF_NOERR) THEN
+         print*, 'Oh probleme lecture dimension zz'
+      ENDIF
+      ierr=NF_INQ_DIMLEN(nid,rid,nlev_amma)
+      print*,'OK nid,rid,nlev_amma',nid,rid,nlev_amma
+!.......................................................................
+      ierr=NF_INQ_DIMID(nid,'time',rid)
+      print*,'nid,rid',nid,rid
+      nt_amma=0
+      IF (ierr.NE.NF_NOERR) THEN
+        stop 'probleme lecture dimension sens'
+      ENDIF
+      ierr=NF_INQ_DIMLEN(nid,rid,nt_amma)
+      print*,'nid,rid,nlev_amma',nid,rid,nt_amma
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!profils initiaux:
+        allocate(plev_amma(nlev_amma))
+        
+        allocate(z_amma(nlev_amma))
+        allocate(th_amma(nlev_amma),q_amma(nlev_amma))
+        allocate(u_amma(nlev_amma))
+        allocate(v_amma(nlev_amma))
+
+!forcings
+        allocate(ht_amma(nlev_amma,nt_amma))
+        allocate(hq_amma(nlev_amma,nt_amma))
+        allocate(vitw_amma(nlev_amma,nt_amma))
+        allocate(lat_amma(nt_amma),sens_amma(nt_amma))
+
+!profils initiaux:
+        allocate(th_ammai(nlev_amma),q_ammai(nlev_amma))
+        allocate(u_ammai(nlev_amma))
+        allocate(v_ammai(nlev_amma))
+        allocate(vitw_ammai(nlev_amma) )
+        allocate(ht_ammai(nlev_amma))
+        allocate(hq_ammai(nlev_amma))
+        allocate(vt_ammai(nlev_amma))
+        allocate(vq_ammai(nlev_amma))
+
+!champs interpoles
+        allocate(vitw_profamma(nlev_amma))
+        allocate(ht_profamma(nlev_amma))
+        allocate(hq_profamma(nlev_amma))
+        allocate(vt_profamma(nlev_amma))
+        allocate(vq_profamma(nlev_amma))
+        allocate(th_profamma(nlev_amma))
+        allocate(q_profamma(nlev_amma))
+        allocate(u_profamma(nlev_amma))
+        allocate(v_profamma(nlev_amma))
+
+        print*,'Allocations OK'
+        call read_amma(nid,nlev_amma,nt_amma                                  &
+     &     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma         &
+     &     ,ht_amma,hq_amma,sens_amma,lat_amma)
+
+END SUBROUTINE read_1D_cases
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+SUBROUTINE deallocate_1D_cases
+!profils initiaux:
+        deallocate(plev_amma)
+        
+        deallocate(z_amma)
+        deallocate(th_amma,q_amma)
+        deallocate(u_amma)
+        deallocate(v_amma)
+
+        deallocate(th_ammai,q_ammai)
+        deallocate(u_ammai)
+        deallocate(v_ammai)
+        deallocate(vitw_ammai )
+        deallocate(ht_ammai)
+        deallocate(hq_ammai)
+        deallocate(vt_ammai)
+        deallocate(vq_ammai)
+        
+!forcings
+        deallocate(ht_amma)
+        deallocate(hq_amma)
+        deallocate(vitw_amma)
+        deallocate(lat_amma,sens_amma)
+
+!champs interpoles
+        deallocate(vitw_profamma)
+        deallocate(ht_profamma)
+        deallocate(hq_profamma)
+        deallocate(vt_profamma)
+        deallocate(vq_profamma)
+        deallocate(th_profamma)
+        deallocate(q_profamma)
+        deallocate(u_profamma)
+        deallocate(v_profamma)
+END SUBROUTINE deallocate_1D_cases
+
+
+END MODULE mod_1D_amma_read
+!=====================================================================
+      subroutine read_amma(nid,nlevel,ntime                          &
+     &     ,zz,pp,temp,qv,u,v,dw                   &
+     &     ,dt,dq,sens,flat)
+
+!program reading forcings of the AMMA case study
+      implicit none
+#include "netcdf.inc"
+
+      integer ntime,nlevel
+
+      real zz(nlevel)
+      real temp(nlevel),pp(nlevel)
+      real qv(nlevel),u(nlevel)
+      real v(nlevel)
+      real dw(nlevel,ntime)
+      real dt(nlevel,ntime)
+      real dq(nlevel,ntime)
+      real flat(ntime),sens(ntime)
+
+
+      integer nid, ierr,rid
+      integer nbvar3d
+      parameter(nbvar3d=30)
+      integer var3didin(nbvar3d)
+
+       ierr=NF_INQ_VARID(nid,"zz",var3didin(1)) 
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'lev'
+         endif
+
+
+      ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'temp'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'qv'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'u'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'v'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'dw'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'dt'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'dq'
+         endif
+      
+      ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'sens'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'flat'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+      endif
+
+!dimensions lecture
+!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture z ok',zz
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture th ok',temp
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture qv ok',qv
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture u ok',u
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture v ok',v
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture w ok',dw
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dt ok',dt
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dq ok',dq
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture sens ok',sens
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture flat ok',flat
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture pp ok',pp
+
+         return 
+         end subroutine read_amma
+!======================================================================
+        SUBROUTINE interp_amma_time(day,day1,annee_ref                     &
+     &         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma       &
+     &         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma               &
+     &         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
+        implicit none
+
+!---------------------------------------------------------------------------------------
+! Time interpolation of a 2D field to the timestep corresponding to day
+!
+! day: current julian day (e.g. 717538.2)
+! day1: first day of the simulation
+! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
+! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
+!---------------------------------------------------------------------------------------
+
+#include "compar1d.h"
+
+! inputs:
+        integer annee_ref
+        integer nt_amma,nlev_amma
+        integer year_ini_amma
+        real day, day1,day_ini_amma,dt_amma
+        real vitw_amma(nlev_amma,nt_amma)
+        real ht_amma(nlev_amma,nt_amma)
+        real hq_amma(nlev_amma,nt_amma)
+        real lat_amma(nt_amma)
+        real sens_amma(nt_amma)
+! outputs:
+        real vitw_prof(nlev_amma)
+        real ht_prof(nlev_amma)
+        real hq_prof(nlev_amma)
+        real lat_prof,sens_prof
+! local:
+        integer it_amma1, it_amma2,k
+        real timeit,time_amma1,time_amma2,frac
+
+
+        if (forcing_type.eq.6) then
+! Check that initial day of the simulation consistent with AMMA case:
+       if (annee_ref.ne.2006) then
+        print*,'Pour AMMA, annee_ref doit etre 2006'
+        print*,'Changer annee_ref dans run.def'
+        stop
+       endif
+       if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
+        print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
+        print*,'Changer dayref dans run.def'
+        stop
+       endif
+       if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
+        print*,'AMMA a fini le 11 juillet'
+        print*,'Changer dayref ou nday dans run.def'
+        stop
+       endif
+       endif
+
+! Determine timestep relative to the 1st day of AMMA:
+!       timeit=(day-day1)*86400.
+!       if (annee_ref.eq.1992) then
+!        timeit=(day-day_ini_toga)*86400.
+!       else
+!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
+!       endif
+      timeit=(day-day_ini_amma)*86400
+
+! Determine the closest observation times:
+!       it_amma1=INT(timeit/dt_amma)+1
+!       it_amma2=it_amma1 + 1
+!       time_amma1=(it_amma1-1)*dt_amma
+!       time_amma2=(it_amma2-1)*dt_amma
+
+       it_amma1=INT(timeit/dt_amma)+1
+       IF (it_amma1 .EQ. nt_amma) THEN
+       it_amma2=it_amma1 
+       ELSE
+       it_amma2=it_amma1 + 1
+       ENDIF
+       time_amma1=(it_amma1-1)*dt_amma
+       time_amma2=(it_amma2-1)*dt_amma
+
+       if (it_amma1 .gt. nt_amma) then
+        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
+     &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
+        stop
+       endif
+
+! time interpolation:
+       frac=(time_amma2-timeit)/(time_amma2-time_amma1)
+       frac=max(frac,0.0)
+
+       lat_prof = lat_amma(it_amma2)                                       &
+     &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1)) 
+       sens_prof = sens_amma(it_amma2)                                     &
+     &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
+
+       do k=1,nlev_amma
+        vitw_prof(k) = vitw_amma(k,it_amma2)                               &
+     &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
+        ht_prof(k) = ht_amma(k,it_amma2)                                   &
+     &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
+        hq_prof(k) = hq_amma(k,it_amma2)                                   &
+     &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
+        enddo
+
+        return
+        END
+
Index: /LMDZ5/branches/testing/libf/phylmd/mod_1D_cases_read.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/mod_1D_cases_read.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/mod_1D_cases_read.F90	(revision 2220)
@@ -2,59 +2,79 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!Declarations specifiques au cas AMMA
-        character*80 :: fich_amma
-! Option du cas AMMA ou on impose la discretisation verticale (Ap,Bp)
-        integer nlev_amma, nt_amma
-
-
-        integer year_ini_amma, day_ini_amma, mth_ini_amma
-        real heure_ini_amma
-        real day_ju_ini_amma   ! Julian day of amma first day
-        parameter (year_ini_amma=2006)
-        parameter (mth_ini_amma=7)
-        parameter (day_ini_amma=10)  ! 10 = 10Juil2006
-        parameter (heure_ini_amma=0.) !0h en secondes
-        real dt_amma
-        parameter (dt_amma=1800.)
-
-!profils initiaux:
-        real, allocatable::  plev_amma(:)
-
-        real, allocatable::  z_amma(:)
-        real, allocatable::  th_amma(:),q_amma(:)
-        real, allocatable::  u_amma(:)
-        real, allocatable::  v_amma(:)
-
-        real, allocatable::  th_ammai(:),q_ammai(:)
-        real, allocatable::  u_ammai(:)
-        real, allocatable::  v_ammai(:)
-        real, allocatable::  vitw_ammai(:)
-        real, allocatable::  ht_ammai(:)
-        real, allocatable::  hq_ammai(:)
-        real, allocatable::  vt_ammai(:)
-        real, allocatable::  vq_ammai(:)
-
-!forcings
-        real, allocatable::  ht_amma(:,:)
-        real, allocatable::  hq_amma(:,:)
-        real, allocatable::  vitw_amma(:,:)
-        real, allocatable::  lat_amma(:),sens_amma(:)
+!Declarations specifiques au cas standard
+        character*80 :: fich_cas
+! Discrétisation 
+        integer nlev_cas, nt_cas
+
+
+        integer year_ini_cas, day_ini_cas, mth_ini_cas
+        real heure_ini_cas
+        real day_ju_ini_cas   ! Julian day of case first day
+        parameter (year_ini_cas=2011)
+        parameter (mth_ini_cas=10)
+        parameter (day_ini_cas=1)  ! 10 = 10Juil2006
+        parameter (heure_ini_cas=0.) !0h en secondes
+        real pdt_cas
+        parameter (pdt_cas=3.*3600)
+
+!CR ATTENTION TEST AMMA
+!        parameter (year_ini_cas=2006)
+!        parameter (mth_ini_cas=7)
+!        parameter (day_ini_cas=10)  ! 10 = 10Juil2006
+!        parameter (heure_ini_cas=0.) !0h en secondes
+!        parameter (pdt_cas=1800.)
+
+!profils environnementaux
+        real, allocatable::  plev_cas(:,:)
+
+        real, allocatable::  z_cas(:,:)
+        real, allocatable::  t_cas(:,:),q_cas(:,:),rh_cas(:,:)
+        real, allocatable::  th_cas(:,:),rv_cas(:,:)
+        real, allocatable::  u_cas(:,:)
+        real, allocatable::  v_cas(:,:)
+
+!forcing
+        real, allocatable::  ht_cas(:,:),vt_cas(:,:),dt_cas(:,:),dtrad_cas(:,:)
+        real, allocatable::  hth_cas(:,:),vth_cas(:,:),dth_cas(:,:)
+        real, allocatable::  hq_cas(:,:),vq_cas(:,:),dq_cas(:,:)
+        real, allocatable::  hr_cas(:,:),vr_cas(:,:),dr_cas(:,:)
+        real, allocatable::  hu_cas(:,:),vu_cas(:,:),du_cas(:,:)
+        real, allocatable::  hv_cas(:,:),vv_cas(:,:),dv_cas(:,:)
+        real, allocatable::  vitw_cas(:,:)
+        real, allocatable::  ug_cas(:,:),vg_cas(:,:)
+        real, allocatable::  lat_cas(:),sens_cas(:),ts_cas(:)
 
 !champs interpoles
-        real, allocatable::  vitw_profamma(:)
-        real, allocatable::  ht_profamma(:)
-        real, allocatable::  hq_profamma(:)
-        real lat_profamma,sens_profamma
-        real, allocatable::  vt_profamma(:)
-        real, allocatable::  vq_profamma(:)
-        real, allocatable::  th_profamma(:)
-        real, allocatable::  q_profamma(:)
-        real, allocatable::  u_profamma(:)
-        real, allocatable::  v_profamma(:)
+        real, allocatable::  plev_prof_cas(:)
+        real, allocatable::  t_prof_cas(:)
+        real, allocatable::  q_prof_cas(:)
+        real, allocatable::  u_prof_cas(:)
+        real, allocatable::  v_prof_cas(:)        
+
+        real, allocatable::  vitw_prof_cas(:)
+        real, allocatable::  ug_prof_cas(:)
+        real, allocatable::  vg_prof_cas(:)
+        real, allocatable::  ht_prof_cas(:)
+        real, allocatable::  hq_prof_cas(:)
+        real, allocatable::  vt_prof_cas(:)
+        real, allocatable::  vq_prof_cas(:)
+        real, allocatable::  dt_prof_cas(:)
+        real, allocatable::  dtrad_prof_cas(:)
+        real, allocatable::  dq_prof_cas(:)
+        real, allocatable::  hu_prof_cas(:)
+        real, allocatable::  hv_prof_cas(:)
+        real, allocatable::  vu_prof_cas(:)
+        real, allocatable::  vv_prof_cas(:)
+        real, allocatable::  du_prof_cas(:)
+        real, allocatable::  dv_prof_cas(:)
+
+
+        real lat_prof_cas,sens_prof_cas,ts_prof_cas
+      
 
 
 CONTAINS
 
-SUBROUTINE read_1D_cases
+SUBROUTINE read_1D_cas
       implicit none
 
@@ -62,9 +82,10 @@
 
       INTEGER nid,rid,ierr
-
-      fich_amma='amma.nc'
-      print*,'fich_amma ',fich_amma
-      ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
-      print*,'fich_amma,NF_NOWRITE,nid ',fich_amma,NF_NOWRITE,nid
+      INTEGER ii,jj
+
+      fich_cas='setup/cas.nc'
+      print*,'fich_cas ',fich_cas
+      ierr = NF_OPEN(fich_cas,NF_NOWRITE,nid)
+      print*,'fich_cas,NF_NOWRITE,nid ',fich_cas,NF_NOWRITE,nid
       if (ierr.NE.NF_NOERR) then
          write(*,*) 'ERROR: GROS Pb opening forcings nc file '
@@ -73,62 +94,89 @@
       endif
 !.......................................................................
+      ierr=NF_INQ_DIMID(nid,'lat',rid)
+      IF (ierr.NE.NF_NOERR) THEN
+         print*, 'Oh probleme lecture dimension lat'
+      ENDIF
+      ierr=NF_INQ_DIMLEN(nid,rid,ii)
+      print*,'OK nid,rid,lat',nid,rid,ii
+!.......................................................................
+      ierr=NF_INQ_DIMID(nid,'lon',rid)
+      IF (ierr.NE.NF_NOERR) THEN
+         print*, 'Oh probleme lecture dimension lon'
+      ENDIF
+      ierr=NF_INQ_DIMLEN(nid,rid,jj)
+      print*,'OK nid,rid,lat',nid,rid,jj
+!.......................................................................
       ierr=NF_INQ_DIMID(nid,'lev',rid)
       IF (ierr.NE.NF_NOERR) THEN
          print*, 'Oh probleme lecture dimension zz'
       ENDIF
-      ierr=NF_INQ_DIMLEN(nid,rid,nlev_amma)
-      print*,'OK nid,rid,nlev_amma',nid,rid,nlev_amma
+      ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
+      print*,'OK nid,rid,nlev_cas',nid,rid,nlev_cas
 !.......................................................................
       ierr=NF_INQ_DIMID(nid,'time',rid)
       print*,'nid,rid',nid,rid
-      nt_amma=0
+      nt_cas=0
       IF (ierr.NE.NF_NOERR) THEN
         stop 'probleme lecture dimension sens'
       ENDIF
-      ierr=NF_INQ_DIMLEN(nid,rid,nt_amma)
-      print*,'nid,rid,nlev_amma',nid,rid,nt_amma
+      ierr=NF_INQ_DIMLEN(nid,rid,nt_cas)
+      print*,'nid,rid,nlev_cas',nid,rid,nt_cas
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!profils initiaux:
-        allocate(plev_amma(nlev_amma))
-        
-        allocate(z_amma(nlev_amma))
-        allocate(th_amma(nlev_amma),q_amma(nlev_amma))
-        allocate(u_amma(nlev_amma))
-        allocate(v_amma(nlev_amma))
-
-!forcings
-        allocate(ht_amma(nlev_amma,nt_amma))
-        allocate(hq_amma(nlev_amma,nt_amma))
-        allocate(vitw_amma(nlev_amma,nt_amma))
-        allocate(lat_amma(nt_amma),sens_amma(nt_amma))
-
-!profils initiaux:
-        allocate(th_ammai(nlev_amma),q_ammai(nlev_amma))
-        allocate(u_ammai(nlev_amma))
-        allocate(v_ammai(nlev_amma))
-        allocate(vitw_ammai(nlev_amma) )
-        allocate(ht_ammai(nlev_amma))
-        allocate(hq_ammai(nlev_amma))
-        allocate(vt_ammai(nlev_amma))
-        allocate(vq_ammai(nlev_amma))
+!profils moyens:
+        allocate(plev_cas(nlev_cas,nt_cas))        
+        allocate(z_cas(nlev_cas,nt_cas))
+        allocate(t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
+        allocate(th_cas(nlev_cas,nt_cas),rv_cas(nlev_cas,nt_cas))
+        allocate(u_cas(nlev_cas,nt_cas))
+        allocate(v_cas(nlev_cas,nt_cas))
+
+!forcing
+        allocate(ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas),dt_cas(nlev_cas,nt_cas),dtrad_cas(nlev_cas,nt_cas))
+        allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
+        allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
+        allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
+        allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
+        allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
+        allocate(vitw_cas(nlev_cas,nt_cas))
+        allocate(ug_cas(nlev_cas,nt_cas))
+        allocate(vg_cas(nlev_cas,nt_cas))
+        allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas))
+
 
 !champs interpoles
-        allocate(vitw_profamma(nlev_amma))
-        allocate(ht_profamma(nlev_amma))
-        allocate(hq_profamma(nlev_amma))
-        allocate(vt_profamma(nlev_amma))
-        allocate(vq_profamma(nlev_amma))
-        allocate(th_profamma(nlev_amma))
-        allocate(q_profamma(nlev_amma))
-        allocate(u_profamma(nlev_amma))
-        allocate(v_profamma(nlev_amma))
+        allocate(plev_prof_cas(nlev_cas))
+        allocate(t_prof_cas(nlev_cas))
+        allocate(q_prof_cas(nlev_cas))
+        allocate(u_prof_cas(nlev_cas))
+        allocate(v_prof_cas(nlev_cas))
+
+        allocate(vitw_prof_cas(nlev_cas))
+        allocate(ug_prof_cas(nlev_cas))
+        allocate(vg_prof_cas(nlev_cas))
+        allocate(ht_prof_cas(nlev_cas))
+        allocate(hq_prof_cas(nlev_cas))
+        allocate(hu_prof_cas(nlev_cas))
+        allocate(hv_prof_cas(nlev_cas))
+        allocate(vt_prof_cas(nlev_cas))
+        allocate(vq_prof_cas(nlev_cas))
+        allocate(vu_prof_cas(nlev_cas))
+        allocate(vv_prof_cas(nlev_cas))
+        allocate(dt_prof_cas(nlev_cas))
+        allocate(dtrad_prof_cas(nlev_cas))
+        allocate(dq_prof_cas(nlev_cas))
+        allocate(du_prof_cas(nlev_cas))
+        allocate(dv_prof_cas(nlev_cas))
 
         print*,'Allocations OK'
-        call read_amma(nid,nlev_amma,nt_amma                                  &
-     &     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma         &
-     &     ,ht_amma,hq_amma,sens_amma,lat_amma)
-
-END SUBROUTINE read_1D_cases
+        call read_cas(nid,nlev_cas,nt_cas                                  &
+     &     ,z_cas,plev_cas,t_cas,q_cas,rh_cas,th_cas,rv_cas,u_cas,v_cas     &
+     &     ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas,dv_cas,hv_cas,vv_cas &
+     &     ,dt_cas,dtrad_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas             &
+     &     ,dth_cas,hth_cas,vth_cas,dr_cas,hr_cas,vr_cas,sens_cas,lat_cas,ts_cas)
+
+
+END SUBROUTINE read_1D_cas
 
 
@@ -136,37 +184,53 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 SUBROUTINE deallocate_1D_cases
-!profils initiaux:
-        deallocate(plev_amma)
+!profils environnementaux:
+        deallocate(plev_cas)
         
-        deallocate(z_amma)
-        deallocate(th_amma,q_amma)
-        deallocate(u_amma)
-        deallocate(v_amma)
-
-        deallocate(th_ammai,q_ammai)
-        deallocate(u_ammai)
-        deallocate(v_ammai)
-        deallocate(vitw_ammai )
-        deallocate(ht_ammai)
-        deallocate(hq_ammai)
-        deallocate(vt_ammai)
-        deallocate(vq_ammai)
+        deallocate(z_cas)
+        deallocate(t_cas,q_cas,rh_cas)
+        deallocate(th_cas,rv_cas)
+        deallocate(u_cas)
+        deallocate(v_cas)
         
-!forcings
-        deallocate(ht_amma)
-        deallocate(hq_amma)
-        deallocate(vitw_amma)
-        deallocate(lat_amma,sens_amma)
+!forcing
+        deallocate(ht_cas,vt_cas,dt_cas,dtrad_cas)
+        deallocate(hq_cas,vq_cas,dq_cas)
+        deallocate(hth_cas,vth_cas,dth_cas)
+        deallocate(hr_cas,vr_cas,dr_cas)
+        deallocate(hu_cas,vu_cas,du_cas)
+        deallocate(hv_cas,vv_cas,dv_cas)
+        deallocate(vitw_cas)
+        deallocate(ug_cas)
+        deallocate(vg_cas)
+        deallocate(lat_cas,sens_cas,ts_cas)
 
 !champs interpoles
-        deallocate(vitw_profamma)
-        deallocate(ht_profamma)
-        deallocate(hq_profamma)
-        deallocate(vt_profamma)
-        deallocate(vq_profamma)
-        deallocate(th_profamma)
-        deallocate(q_profamma)
-        deallocate(u_profamma)
-        deallocate(v_profamma)
+        deallocate(plev_prof_cas)
+        deallocate(t_prof_cas)
+        deallocate(q_prof_cas)
+        deallocate(u_prof_cas)
+        deallocate(v_prof_cas)
+
+        deallocate(vitw_prof_cas)
+        deallocate(ug_prof_cas)
+        deallocate(vg_prof_cas)
+        deallocate(ht_prof_cas)
+        deallocate(hq_prof_cas)
+        deallocate(hu_prof_cas)
+        deallocate(hv_prof_cas)
+        deallocate(vt_prof_cas)
+        deallocate(vq_prof_cas)
+        deallocate(vu_prof_cas)
+        deallocate(vv_prof_cas)
+        deallocate(dt_prof_cas)
+        deallocate(dtrad_prof_cas)
+        deallocate(dq_prof_cas)
+        deallocate(du_prof_cas)
+        deallocate(dv_prof_cas)
+        deallocate(t_prof_cas)
+        deallocate(q_prof_cas)
+        deallocate(u_prof_cas)
+        deallocate(v_prof_cas)
+
 END SUBROUTINE deallocate_1D_cases
 
@@ -174,9 +238,10 @@
 END MODULE mod_1D_cases_read
 !=====================================================================
-      subroutine read_amma(nid,nlevel,ntime                          &
-     &     ,zz,pp,temp,qv,u,v,dw                   &
-     &     ,dt,dq,sens,flat)
-
-!program reading forcings of the AMMA case study
+      subroutine read_cas(nid,nlevel,ntime                          &
+     &     ,zz,pp,temp,qv,rh,theta,rv,u,v,ug,vg,w,                   &
+     &     du,hu,vu,dv,hv,vv,dt,dtrad,ht,vt,dq,hq,vq,                     &
+     &     dth,hth,vth,dr,hr,vr,sens,flat,ts)
+
+!program reading forcing of the case study
       implicit none
 #include "netcdf.inc"
@@ -184,17 +249,26 @@
       integer ntime,nlevel
 
-      real zz(nlevel)
-      real temp(nlevel),pp(nlevel)
-      real qv(nlevel),u(nlevel)
-      real v(nlevel)
-      real dw(nlevel,ntime)
-      real dt(nlevel,ntime)
-      real dq(nlevel,ntime)
-      real flat(ntime),sens(ntime)
+      real zz(nlevel,ntime)
+      real pp(nlevel,ntime)
+      real temp(nlevel,ntime),qv(nlevel,ntime),rh(nlevel,ntime)
+      real theta(nlevel,ntime),rv(nlevel,ntime)
+      real u(nlevel,ntime)
+      real v(nlevel,ntime)
+      real ug(nlevel,ntime)
+      real vg(nlevel,ntime)
+      real w(nlevel,ntime)
+      real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
+      real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
+      real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
+      real dtrad(nlevel,ntime)
+      real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
+      real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime)
+      real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
+      real flat(ntime),sens(ntime),ts(ntime)
 
 
       integer nid, ierr,rid
       integer nbvar3d
-      parameter(nbvar3d=30)
+      parameter(nbvar3d=34)
       integer var3didin(nbvar3d)
 
@@ -204,7 +278,13 @@
            stop 'lev'
          endif
-
-
-      ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
+      
+      ierr=NF_INQ_VARID(nid,"pp",var3didin(2)) 
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'plev'
+         endif
+
+
+      ierr=NF_INQ_VARID(nid,"temp",var3didin(3))
          if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
@@ -212,5 +292,5 @@
          endif
 
-      ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
+      ierr=NF_INQ_VARID(nid,"qv",var3didin(4))
          if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
@@ -218,5 +298,24 @@
          endif
 
-      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
+      ierr=NF_INQ_VARID(nid,"rh",var3didin(5))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'rh'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'theta'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"rv",var3didin(7))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'rv'
+         endif
+
+
+      ierr=NF_INQ_VARID(nid,"u",var3didin(8))
          if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
@@ -224,5 +323,5 @@
          endif
 
-      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
+      ierr=NF_INQ_VARID(nid,"v",var3didin(9))
          if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
@@ -230,23 +329,137 @@
          endif
 
-      ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'dw'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'dt'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'dq'
+       ierr=NF_INQ_VARID(nid,"ug",var3didin(10))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'ug'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"vg",var3didin(11))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vg'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"w",var3didin(12))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'w'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"advu",var3didin(13))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'advu'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"hu",var3didin(14))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hu'
+         endif
+
+       ierr=NF_INQ_VARID(nid,"vu",var3didin(15))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vu'
+         endif
+
+       ierr=NF_INQ_VARID(nid,"advv",var3didin(16))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'advv'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"hv",var3didin(17))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hv'
+         endif
+
+       ierr=NF_INQ_VARID(nid,"vv",var3didin(18))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vv'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"advT",var3didin(19))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'advT'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"hT",var3didin(20))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hT'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"vT",var3didin(21))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vT'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"advq",var3didin(22))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'advq'
          endif
       
-      ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
+      ierr=NF_INQ_VARID(nid,"hq",var3didin(23))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hq'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"vq",var3didin(24))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vq'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"advth",var3didin(25))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'advth'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"hth",var3didin(26))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hth'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"vth",var3didin(27))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vth'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"advr",var3didin(28))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'advr'
+         endif
+      
+      ierr=NF_INQ_VARID(nid,"hr",var3didin(29))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hr'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"vr",var3didin(30))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vr'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"radT",var3didin(31))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'radT'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"sens",var3didin(32))
          if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
@@ -254,5 +467,5 @@
          endif
 
-      ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
+      ierr=NF_INQ_VARID(nid,"flat",var3didin(33))
          if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
@@ -260,11 +473,9 @@
          endif
 
-      ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-      endif
-
-!dimensions lecture
-!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
+      ierr=NF_INQ_VARID(nid,"ts",var3didin(34))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'ts'
+         endif
  
 #ifdef NC_DOUBLE
@@ -280,18 +491,30 @@
 
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture th ok',temp
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),pp)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(2),pp)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture pp ok',pp
+
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),temp)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(3),temp)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture T ok',temp
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),qv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(4),qv)
 #endif
          if(ierr/=NF_NOERR) then
@@ -302,7 +525,40 @@
  
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),rh)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(5),rh)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture rh ok',rh
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),theta)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(6),theta)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture theta ok',theta
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),rv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(7),rv)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture rv ok',rv
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),u)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(8),u)
 #endif
          if(ierr/=NF_NOERR) then
@@ -313,7 +569,7 @@
 
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),v)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(9),v)
 #endif
          if(ierr/=NF_NOERR) then
@@ -324,18 +580,106 @@
 
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture w ok',dw
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),ug)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(10),ug)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture ug ok',ug
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),vg)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(11),vg)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vg ok',vg
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),w)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(12),w)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture w ok',w
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),du)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(13),du)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture du ok',du
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),hu)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(14),hu)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hu ok',hu
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),vu)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(15),vu)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vu ok',vu
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),dv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(16),dv)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dv ok',dv
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(17),hv)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hv ok',hv
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),vv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(18),vv)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vv ok',vv
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),dt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(19),dt)
 #endif
          if(ierr/=NF_NOERR) then
@@ -346,7 +690,29 @@
 
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),ht)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(20),ht)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture ht ok',ht
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),vt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(21),vt)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vt ok',vt
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),dq)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(22),dq)
 #endif
          if(ierr/=NF_NOERR) then
@@ -357,7 +723,106 @@
 
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(23),hq)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(23),hq)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hq ok',hq
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(24),vq)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(24),vq)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vq ok',vq
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(25),dth)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(25),dth)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dth ok',dth
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(26),hth)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(26),hth)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hth ok',hth
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(27),vth)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(27),vth)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vth ok',vth
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(28),dr)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(28),dr)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dr ok',dr
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(29),hr)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(29),hr)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hr ok',hr
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(30),vr)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(30),vr)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vr ok',vr
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(31),dtrad)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(31),dtrad)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dtrad ok',dtrad
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(32),sens)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(32),sens)
 #endif
          if(ierr/=NF_NOERR) then
@@ -368,7 +833,7 @@
 
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(33),flat)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(33),flat)
 #endif
          if(ierr/=NF_NOERR) then
@@ -379,21 +844,31 @@
 
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture pp ok',pp
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(34),ts)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(34),ts)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture ts ok',ts
 
          return 
-         end subroutine read_amma
+         end subroutine read_cas
 !======================================================================
-        SUBROUTINE interp_amma_time(day,day1,annee_ref                     &
-     &         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma       &
-     &         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma               &
-     &         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
+        SUBROUTINE interp_case_time(day,day1,annee_ref                &
+     &         ,year_ini_cas,day_ini_cas,nt_cas,pdt_cas,nlev_cas      &
+     &         ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas               &
+     &         ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas           &
+     &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas   &
+     &         ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas          &
+     &         ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas       &
+     &         ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas         &
+     &         ,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
+     &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas       &
+     &         ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas    &
+     &         ,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas)
+          
+
         implicit none
 
@@ -403,6 +878,6 @@
 ! day: current julian day (e.g. 717538.2)
 ! day1: first day of the simulation
-! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
-! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
+! nt_cas: total nb of data in the forcing 
+! pdt_cas: total time interval (in sec) between 2 forcing data
 !---------------------------------------------------------------------------------------
 
@@ -411,23 +886,61 @@
 ! inputs:
         integer annee_ref
-        integer nt_amma,nlev_amma
-        integer year_ini_amma
-        real day, day1,day_ini_amma,dt_amma
-        real vitw_amma(nlev_amma,nt_amma)
-        real ht_amma(nlev_amma,nt_amma)
-        real hq_amma(nlev_amma,nt_amma)
-        real lat_amma(nt_amma)
-        real sens_amma(nt_amma)
+        integer nt_cas,nlev_cas
+        integer year_ini_cas
+        real day_ini_cas
+        real day, day1,pdt_cas
+        real ts_cas(nt_cas)
+        real plev_cas(nlev_cas,nt_cas)
+        real t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas)
+        real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
+        real ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
+        real vitw_cas(nlev_cas,nt_cas)
+        real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
+        real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
+        real dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
+        real dtrad_cas(nlev_cas,nt_cas)
+        real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
+        real lat_cas(nt_cas)
+        real sens_cas(nt_cas)
+
 ! outputs:
-        real vitw_prof(nlev_amma)
-        real ht_prof(nlev_amma)
-        real hq_prof(nlev_amma)
-        real lat_prof,sens_prof
+        real plev_prof_cas(nlev_cas)
+        real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
+        real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
+        real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas)
+        real vitw_prof_cas(nlev_cas)
+        real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
+        real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
+        real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
+        real dtrad_prof_cas(nlev_cas)
+        real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
+        real lat_prof_cas,sens_prof_cas,ts_prof_cas
 ! local:
-        integer it_amma1, it_amma2,k
-        real timeit,time_amma1,time_amma2,frac
-
-
-        if (forcing_type.eq.6) then
+        integer it_cas1, it_cas2,k
+        real timeit,time_cas1,time_cas2,frac
+
+
+        print*,'Check time',day1,day_ini_cas,day_ini_cas+90
+
+        if ((forcing_type.eq.10).and.(1.eq.1)) then
+! Check that initial day of the simulation consistent with the case:
+       if (annee_ref.ne.2011) then
+        print*,'Pour CINDY, annee_ref doit etre 2011'
+        print*,'Changer annee_ref dans run.def'
+        stop
+       endif
+       if (annee_ref.eq.2011 .and. day1.lt.day_ini_cas) then
+        print*,'CINDY a débuté le 1 octobre 2011',day1,day_ini_cas
+        print*,'Changer dayref dans run.def'
+        stop
+       endif
+       if (annee_ref.eq.2011 .and. day1.gt.day_ini_cas+90) then
+        print*,'CINDY a fini le 31 decembre'
+        print*,'Changer dayref ou nday dans run.def',day1,day_ini_cas+90
+        stop
+       endif
+       endif
+!CR:ATTENTION TEST AMMA
+      if ((forcing_type.eq.10).and.(1.eq.0)) then
 ! Check that initial day of the simulation consistent with AMMA case:
        if (annee_ref.ne.2006) then
@@ -436,10 +949,10 @@
         stop
        endif
-       if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
-        print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
+       if (annee_ref.eq.2006 .and. day1.lt.day_ini_cas) then
+        print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_cas
         print*,'Changer dayref dans run.def'
         stop
        endif
-       if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
+       if (annee_ref.eq.2006 .and. day1.gt.day_ini_cas+1) then
         print*,'AMMA a fini le 11 juillet'
         print*,'Changer dayref ou nday dans run.def'
@@ -448,50 +961,88 @@
        endif
 
-! Determine timestep relative to the 1st day of AMMA:
+! Determine timestep relative to the 1st day:
 !       timeit=(day-day1)*86400.
 !       if (annee_ref.eq.1992) then
-!        timeit=(day-day_ini_toga)*86400.
+!        timeit=(day-day_ini_cas)*86400.
 !       else
 !        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
 !       endif
-      timeit=(day-day_ini_amma)*86400
+      timeit=(day-day_ini_cas)*86400
 
 ! Determine the closest observation times:
-!       it_amma1=INT(timeit/dt_amma)+1
-!       it_amma2=it_amma1 + 1
-!       time_amma1=(it_amma1-1)*dt_amma
-!       time_amma2=(it_amma2-1)*dt_amma
-
-       it_amma1=INT(timeit/dt_amma)+1
-       IF (it_amma1 .EQ. nt_amma) THEN
-       it_amma2=it_amma1 
+!       it_cas1=INT(timeit/pdt_cas)+1
+!       it_cas2=it_cas1 + 1
+!       time_cas1=(it_cas1-1)*pdt_cas
+!       time_cas2=(it_cas2-1)*pdt_cas
+
+       it_cas1=INT(timeit/pdt_cas)+1
+       IF (it_cas1 .EQ. nt_cas) THEN
+       it_cas2=it_cas1 
        ELSE
-       it_amma2=it_amma1 + 1
+       it_cas2=it_cas1 + 1
        ENDIF
-       time_amma1=(it_amma1-1)*dt_amma
-       time_amma2=(it_amma2-1)*dt_amma
-
-       if (it_amma1 .gt. nt_amma) then
-        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
-     &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
+       time_cas1=(it_cas1-1)*pdt_cas
+       time_cas2=(it_cas2-1)*pdt_cas
+
+       if (it_cas1 .gt. nt_cas) then
+        write(*,*) 'PB-stop: day, it_cas1, it_cas2, timeit: '            &
+     &        ,day,day_ini_cas,it_cas1,it_cas2,timeit/86400.
         stop
        endif
 
 ! time interpolation:
-       frac=(time_amma2-timeit)/(time_amma2-time_amma1)
+       frac=(time_cas2-timeit)/(time_cas2-time_cas1)
        frac=max(frac,0.0)
 
-       lat_prof = lat_amma(it_amma2)                                       &
-     &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1)) 
-       sens_prof = sens_amma(it_amma2)                                     &
-     &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
-
-       do k=1,nlev_amma
-        vitw_prof(k) = vitw_amma(k,it_amma2)                               &
-     &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
-        ht_prof(k) = ht_amma(k,it_amma2)                                   &
-     &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
-        hq_prof(k) = hq_amma(k,it_amma2)                                   &
-     &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
+       lat_prof_cas = lat_cas(it_cas2)                                       &
+     &          -frac*(lat_cas(it_cas2)-lat_cas(it_cas1)) 
+       sens_prof_cas = sens_cas(it_cas2)                                     &
+     &          -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
+       ts_prof_cas = ts_cas(it_cas2)                                     &
+     &          -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
+
+       do k=1,nlev_cas
+        plev_prof_cas(k) = plev_cas(k,it_cas2)                               &
+     &          -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
+        t_prof_cas(k) = t_cas(k,it_cas2)                               &
+     &          -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
+        q_prof_cas(k) = q_cas(k,it_cas2)                               &
+     &          -frac*(q_cas(k,it_cas2)-q_cas(k,it_cas1))
+        u_prof_cas(k) = u_cas(k,it_cas2)                               &
+     &          -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
+        v_prof_cas(k) = v_cas(k,it_cas2)                               &
+     &          -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
+        ug_prof_cas(k) = ug_cas(k,it_cas2)                               &
+     &          -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
+        vg_prof_cas(k) = vg_cas(k,it_cas2)                               &
+     &          -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
+        vitw_prof_cas(k) = vitw_cas(k,it_cas2)                               &
+     &          -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
+        du_prof_cas(k) = du_cas(k,it_cas2)                                   &
+     &          -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
+        hu_prof_cas(k) = hu_cas(k,it_cas2)                                   &
+     &          -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
+        vu_prof_cas(k) = vu_cas(k,it_cas2)                                   &
+     &          -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
+        dv_prof_cas(k) = dv_cas(k,it_cas2)                                   &
+     &          -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
+        hv_prof_cas(k) = hv_cas(k,it_cas2)                                   &
+     &          -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
+        vv_prof_cas(k) = vv_cas(k,it_cas2)                                   &
+     &          -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
+        dt_prof_cas(k) = dt_cas(k,it_cas2)                                   &
+     &          -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
+        ht_prof_cas(k) = ht_cas(k,it_cas2)                                   &
+     &          -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
+        vt_prof_cas(k) = vt_cas(k,it_cas2)                                   &
+     &          -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
+        dtrad_prof_cas(k) = dtrad_cas(k,it_cas2)                                   &
+     &          -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
+        dq_prof_cas(k) = dq_cas(k,it_cas2)                                   &
+     &          -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
+        hq_prof_cas(k) = hq_cas(k,it_cas2)                                   &
+     &          -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
+        vq_prof_cas(k) = vq_cas(k,it_cas2)                                   &
+     &          -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
         enddo
 
@@ -499,2 +1050,3 @@
         END
 
+!**********************************************************************************************
Index: /LMDZ5/branches/testing/libf/phylmd/nuage.h
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/nuage.h	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/nuage.h	(revision 2220)
@@ -5,9 +5,11 @@
       REAL exposant_glace
       REAL rei_min,rei_max
+      REAL tau_cld_cv,coefw_cld_cv
 
-      INTEGER iflag_t_glace
+      INTEGER iflag_t_glace,iflag_cld_cv
 
       common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max,     &
      &                  t_glace_min,exposant_glace,rei_min,rei_max,     &
-     &                  iflag_t_glace
+     &                  tau_cld_cv,coefw_cld_cv,                        &
+     &                  iflag_t_glace,iflag_cld_cv
 !$OMP THREADPRIVATE(/nuagecom/)
Index: /LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90	(revision 2220)
@@ -8,8 +8,13 @@
   USE dimphy
   USE indice_sol_mod
+  USE surface_data
 
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice!, ocean_slab_ice
+  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
+
+!****************************************************************************************
+! Global saved variables
+!****************************************************************************************
 
   INTEGER, PRIVATE, SAVE                           :: cpl_pas
@@ -21,10 +26,56 @@
   REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
   !$OMP THREADPRIVATE(tslab)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE  :: pctsrf
-  !$OMP THREADPRIVATE(pctsrf)
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
+  !$OMP THREADPRIVATE(fsic)
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
+  !$OMP THREADPRIVATE(tice)
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
+  !$OMP THREADPRIVATE(seaice)
   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
   !$OMP THREADPRIVATE(slab_bils)
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
   !$OMP THREADPRIVATE(bils_cum)
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
+  !$OMP THREADPRIVATE(slab_bilg)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
+  !$OMP THREADPRIVATE(bilg_cum)
+
+!****************************************************************************************
+! Parameters (could be read in def file: move to slab_init)
+!****************************************************************************************
+! snow and ice physical characteristics:
+    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
+    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
+    REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
+    REAL, PARAMETER :: ice_den=917.
+    REAL, PARAMETER :: sea_den=1025.
+    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
+    REAL, PARAMETER :: sno_cond=0.31*sno_den
+    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
+    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
+
+! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
+    REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
+    REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
+    REAL, PARAMETER :: ice_frac_min=0.001
+    REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction 
+    REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness 
+    REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness 
+    ! below ice_thin, priority is melt lateral / grow height
+    ! ice_thin is also height of new ice
+    REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness 
+    ! above ice_thick, priority is melt height / grow lateral
+    REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice 
+    REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height 
+
+! albedo  and radiation parameters
+    REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
+    REAL, PARAMETER :: alb_sno_del=0.3  !max snow albedo = min + del
+    REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
+    REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
+    REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
+    REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
+
+!****************************************************************************************
 
 CONTAINS
@@ -56,14 +107,21 @@
 ! Allocate surface fraction read from restart file
 !****************************************************************************************
-    ALLOCATE(pctsrf(klon,nbsrf), stat = error)
+    ALLOCATE(fsic(klon), stat = error)
     IF (error /= 0) THEN
        abort_message='Pb allocation tmp_pctsrf_slab'
        CALL abort_gcm(modname,abort_message,1)
     ENDIF
-    pctsrf(:,:) = pctsrf_rst(:,:)
-
-!****************************************************************************************
-! Allocate local variables
-!****************************************************************************************
+    fsic(:)=0.
+    WHERE (1.-zmasq(:)>EPSFRA)
+        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
+    END WHERE
+
+!****************************************************************************************
+! Allocate saved variables
+!****************************************************************************************
+    ALLOCATE(tslab(klon,nslay), stat=error)
+       IF (error /= 0) CALL abort_gcm &
+         (modname,'pb allocation tslab', 1)
+
     ALLOCATE(slab_bils(klon), stat = error)
     IF (error /= 0) THEN
@@ -79,4 +137,32 @@
     bils_cum(:) = 0.0   
 
+    IF (version_ocean=='sicINT') THEN
+        ALLOCATE(slab_bilg(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_bilg'
+           CALL abort_gcm(modname,abort_message,1)
+        ENDIF
+        slab_bilg(:) = 0.0   
+        ALLOCATE(bilg_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_bilg_cum'
+           CALL abort_gcm(modname,abort_message,1)
+        ENDIF
+        bilg_cum(:) = 0.0   
+        ALLOCATE(tice(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_tice'
+           CALL abort_gcm(modname,abort_message,1)
+        ENDIF
+        ALLOCATE(seaice(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_seaice'
+           CALL abort_gcm(modname,abort_message,1)
+        ENDIF
+    END IF
+
+!****************************************************************************************
+! Define some parameters
+!****************************************************************************************
 ! Layer thickness
     ALLOCATE(slabh(nslay), stat = error)
@@ -94,4 +180,5 @@
     CALL getin('cpl_pas',cpl_pas)
     print *,'cpl_pas',cpl_pas
+
  END SUBROUTINE ocean_slab_init
 !
@@ -101,5 +188,4 @@
 
     USE limit_read_mod
-    USE surface_data
 
 !    INCLUDE "clesphys.h"
@@ -119,5 +205,6 @@
        CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
     ELSE
-       pctsrf_chg(:,:)=pctsrf(:,:)
+       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
+       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
        is_modified=.TRUE.
     END IF
@@ -133,10 +220,9 @@
        AcoefU, AcoefV, BcoefU, BcoefV, &
        ps, u1, v1, tsurf_in, &
-       radsol, snow, agesno, &
+       radsol, snow, &
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
        tsurf_new, dflux_s, dflux_l, qflux)
     
     USE calcul_fluxs_mod
-    USE surface_data
 
     INCLUDE "iniprint.h"
@@ -158,10 +244,9 @@
     REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
     REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
+    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
 
 ! In/Output arguments
 !****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
     REAL, DIMENSION(klon), INTENT(INOUT) :: snow
-    REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
     
 ! Output arguments
@@ -177,8 +262,13 @@
 !****************************************************************************************
     INTEGER               :: i,ki
+    ! surface fluxes
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd
-    REAL, DIMENSION(klon) :: diff_sst, lmt_bils
+    ! flux correction
+    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
+    ! surface wind stress
     REAL, DIMENSION(klon) :: u0, v0
     REAL, DIMENSION(klon) :: u1_lay, v1_lay
+    ! ice creation
+    REAL                  :: e_freeze, h_new, dfsic
 
 !****************************************************************************************
@@ -189,5 +279,4 @@
     beta(:)     = 1.
     dif_grnd(:) = 0.
-    agesno(:)   = 0.
     
 ! Suppose zero surface speed
@@ -215,5 +304,6 @@
     DO i=1,knon
         ki=knindex(i)
-        slab_bils(ki)=(fluxlat(i)+fluxsens(i)+radsol(i))*pctsrf(ki,is_oce)/(1.-zmasq(ki))
+        slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
+                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
         bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
 ! Also taux, tauy, saved vars...
@@ -225,7 +315,7 @@
 !****************************************************************************************
     lmt_bils(:)=0.
-    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst)  ! global pour un processus
-    ! lmt_bils and diff_sst saved by limit_slab
-    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.
+    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
+    ! lmt_bils and diff_sst,siv saved by limit_slab
+    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
     ! qflux = total QFlux correction (in W/m2)
 
@@ -248,28 +338,51 @@
                 tsurf_new(i)=tslab(ki,1)
             END DO
-        CASE('sicOBS') ! check for sea ice or tsurf below freezing
+        CASE('sicOBS') ! check for sea ice or tslab below freezing
             DO i=1,knon
                 ki=knindex(i)
-                IF ((tslab(ki,1).LT.t_freeze).OR.(pctsrf(ki,is_sic).GT.epsfra)) THEN
-                    tsurf_new(i)=t_freeze
+                IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
                     tslab(ki,1)=t_freeze
-                ELSE
-                    tsurf_new(i)=tslab(ki,1)
                 END IF
+                tsurf_new(i)=tslab(ki,1)
             END DO
         CASE('sicINT')
             DO i=1,knon
                 ki=knindex(i)
-                IF (pctsrf(ki,is_sic).LT.epsfra) THEN ! Free of ice
-                    IF (tslab(ki,1).GT.t_freeze) THEN 
+                IF (fsic(ki).LT.epsfra) THEN ! Free of ice
+                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
+                        ! quantity of new ice formed
+                        e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
+                        ! new ice
+                        tice(ki)=t_freeze
+                        fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
+                        IF (fsic(ki).GT.ice_frac_min) THEN
+                            seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
+                            tslab(ki,1)=t_freeze
+                        ELSE
+                            fsic(ki)=0.
+                        END IF
+                        tsurf_new(i)=t_freeze
+                    ELSE
                         tsurf_new(i)=tslab(ki,1)
-                    ELSE
-                        tsurf_new(i)=t_freeze
-                        ! Call new ice routine
-                        tslab(ki,1)=t_freeze
-                    END IF
-                ELSE ! ice present, tslab update completed in slab_ice
+                    END IF 
+                ELSE ! ice present
                     tsurf_new(i)=t_freeze
-                END IF !ice free
+                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
+                        ! quantity of new ice formed over open ocean
+                        e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
+                                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
+                        ! new ice height and fraction
+                        h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
+                        dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
+                        h_new=MIN(e_freeze/dfsic,h_ice_max)
+                        ! update tslab to freezing over open ocean only
+                        tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
+                        ! update sea ice
+                        seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
+                                   /(dfsic+fsic(ki))
+                        fsic(ki)=fsic(ki)+dfsic
+                        ! update snow?
+                    END IF !freezing
+                END IF ! sea ice present
             END DO
         END SELECT
@@ -279,54 +392,308 @@
 !
 !****************************************************************************************
-!
-!  SUBROUTINE ocean_slab_ice(   &
-!       itime, dtime, jour, knon, knindex, &
-!       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
-!       AcoefH, AcoefQ, BcoefH, BcoefQ, &
-!       AcoefU, AcoefV, BcoefU, BcoefV, &
-!       ps, u1, v1, &
-!       radsol, snow, qsurf, qsol, agesno, tsoil, &
-!       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-!       tsurf_new, dflux_s, dflux_l)
-!
+
+  SUBROUTINE ocean_slab_ice(   &
+       itime, dtime, jour, knon, knindex, &
+       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+       AcoefH, AcoefQ, BcoefH, BcoefQ, &
+       AcoefU, AcoefV, BcoefU, BcoefV, &
+       ps, u1, v1, &
+       radsol, snow, qsurf, qsol, agesno, &
+       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
+       tsurf_new, dflux_s, dflux_l, swnet)
+
+   USE calcul_fluxs_mod
+
+   INCLUDE "YOMCST.h"
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                  :: itime, jour, knon
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
+    REAL, INTENT(IN)                     :: dtime
+    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
+    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
+    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
+    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
+    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
+    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
+    REAL, DIMENSION(klon), INTENT(IN)    :: ps
+    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
+    REAL, DIMENSION(klon), INTENT(IN)    :: swnet
+
+! In/Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
+    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
+    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
+    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
+    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
+    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l
+
+! Local variables
+!****************************************************************************************
+    INTEGER               :: i,ki
+    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
+    REAL, DIMENSION(klon) :: u0, v0
+    REAL, DIMENSION(klon) :: u1_lay, v1_lay
+    ! intermediate heat fluxes:
+    REAL                  :: f_cond, f_swpen
+    ! for snow/ice albedo:
+    REAL                  :: alb_snow, alb_ice, alb_pond
+    REAL                  :: frac_snow, frac_ice, frac_pond
+    ! for ice melt / freeze
+    REAL                  :: e_melt, snow_evap, h_test
+    ! dhsic, dfsic change in ice mass, fraction.
+    REAL                  :: dhsic, dfsic, frac_mf
+
 !****************************************************************************************
 ! 1) Flux calculation
 !****************************************************************************************
-! set beta, cal etc. depends snow / ice surf ?
+! Suppose zero surface speed
+    u0(:)=0.0
+    v0(:)=0.0
+    u1_lay(:) = u1(:) - u0(:)
+    v1_lay(:) = v1(:) - v0(:)
+
+! set beta, cal, compute conduction fluxes inside ice/snow 
+    slab_bilg(:)=0.
+    dif_grnd(:)=0.
+    beta(:) = 1.
+    DO i=1,knon
+    ki=knindex(i)
+        IF (snow(i).GT.snow_min) THEN
+            ! snow-layer heat capacity
+            cal(i)=2.*RCPD/(snow(i)*ice_cap)
+            ! snow conductive flux
+            f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
+            ! all shortwave flux absorbed
+            f_swpen=0.
+            ! bottom flux (ice conduction)
+            slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
+            ! update ice temperature
+            tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
+                     *(slab_bilg(ki)+f_cond)*dtime
+       ELSE ! bare ice
+            ! ice-layer heat capacity
+            cal(i)=2.*RCPD/(seaice(ki)*ice_cap)
+            ! conductive flux
+            f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
+            ! penetrative shortwave flux...
+            f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den)
+            slab_bilg(ki)=f_swpen-f_cond
+        END IF
+        radsol(i)=radsol(i)+f_cond-f_swpen
+    END DO
+    ! weight fluxes to ocean by sea ice fraction
+    slab_bilg(:)=slab_bilg(:)*fsic(:)
+
 ! calcul_fluxs (sens, lat etc)
+    CALL calcul_fluxs(knon, is_sic, dtime, &
+        tsurf_in, p1lay, cal, beta, cdragh, ps, &
+        precip_rain, precip_snow, snow, qsurf,  &
+        radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
+        AcoefH, AcoefQ, BcoefH, BcoefQ, &
+        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
+    DO i=1,knon 
+        IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
+    END DO
+
 ! calcul_flux_wind
-
-!****************************************************************************************
-! 2) Update surface
-!****************************************************************************************
-! neige, fonte
-! flux glace-ocean
-! update temperature
-! neige precip, evap
-! Melt snow & ice from above
+    CALL calcul_flux_wind(knon, dtime, &
+         u0, v0, u1, v1, cdragm, &
+         AcoefU, AcoefV, BcoefU, BcoefV, &
+         p1lay, temp_air, &
+         flux_u1, flux_v1)
+
+!****************************************************************************************
+! 2) Update snow and ice surface
+!****************************************************************************************
+! snow precip
+    DO i=1,knon
+        ki=knindex(i)
+        IF (precip_snow(i) > 0.) THEN
+            snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
+        END IF
+! snow and ice sublimation
+        IF (evap(i) > 0.) THEN 
+           snow_evap = MIN (snow(i) / dtime, evap(i))
+           snow(i) = snow(i) - snow_evap * dtime
+           snow(i) = MAX(0.0, snow(i))
+           seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
+        ENDIF
+! Melt / Freeze from above if Tsurf>0
+        IF (tsurf_new(i).GT.t_melt) THEN
+            ! energy available for melting snow (in kg/m2 of snow)
+            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
+               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
+            ! remove snow
+            IF (snow(i).GT.e_melt) THEN
+                snow(i)=snow(i)-e_melt
+                tsurf_new(i)=t_melt
+            ELSE ! all snow is melted
+                ! add remaining heat flux to ice
+                e_melt=e_melt-snow(i)
+                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
+                tsurf_new(i)=tice(ki)
+            END IF
+        END IF
+! melt ice from above if Tice>0
+        IF (tice(ki).GT.t_melt) THEN
+            ! quantity of ice melted (kg/m2)
+            e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. & 
+             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
+            ! melt from above, height only
+            dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
+            e_melt=e_melt-dhsic
+            IF (e_melt.GT.0) THEN
+            ! lateral melt if ice too thin
+            dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
+            ! if all melted add remaining heat to ocean
+            e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min)
+            slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
+            ! update height and fraction
+            fsic(ki)=fsic(ki)-dfsic
+            END IF
+            seaice(ki)=seaice(ki)-dhsic
+            ! surface temperature at melting point
+            tice(ki)=t_melt
+            tsurf_new(i)=t_melt
+        END IF
+        ! convert snow to ice if below floating line
+        h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den
+        IF (h_test.GT.0.) THEN !snow under water
+            ! extra snow converted to ice (with added frozen sea water)
+            dhsic=h_test/(sea_den-ice_den+sno_den)
+            seaice(ki)=seaice(ki)+dhsic
+            snow(i)=snow(i)-dhsic*sno_den/ice_den
+            ! available energy (freeze sea water + bring to tice)
+            e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
+                   ice_cap/2.*(t_freeze-tice(ki)))
+            ! update ice temperature
+            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
+        END IF
+    END DO
+
 ! New albedo
-
-!****************************************************************************************
-! 3) Recalculate new ocean temperature
+    DO i=1,knon
+        ki=knindex(i)
+       ! snow albedo: update snow age
+        IF (snow(i).GT.0.0001) THEN
+             agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
+                         * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.)
+        ELSE
+            agesno(i)=0.0
+        END IF
+        ! snow albedo
+        alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
+        ! ice albedo (varies with ice tkickness and temp)
+        alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1)
+        IF (tice(ki).GT.t_freeze-0.01) THEN
+            alb_ice=MIN(alb_ice,alb_ice_wet)
+        ELSE
+            alb_ice=MIN(alb_ice,alb_ice_dry)
+        END IF
+        ! pond albedo
+        alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
+        ! pond fraction
+        frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
+        ! snow fraction
+        frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
+        ! ice fraction
+        frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
+        ! total albedo
+        alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond
+    END DO
+    alb2_new(:) = alb1_new(:)
+
+!****************************************************************************************
+! 3) Recalculate new ocean temperature (add fluxes below ice)
 !    Melt / freeze from below
 !***********************************************o*****************************************
-
-
-!  END SUBROUTINE ocean_slab_ice
+    !cumul fluxes
+    bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
+    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
+        ! Add cumulated surface fluxes
+        tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime
+        DO i=1,knon
+            ki=knindex(i)
+            ! split lateral/top melt-freeze
+            frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
+            IF (tslab(ki,1).LE.t_freeze) THEN 
+               ! ****** Form new ice from below *******
+               ! quantity of new ice
+                e_melt=(t_freeze-tslab(ki,1))/cyang & 
+                       /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
+               ! first increase height to h_thin
+               dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
+               seaice(ki)=dhsic+seaice(ki)
+               e_melt=e_melt-fsic(ki)*dhsic
+               IF (e_melt.GT.0.) THEN
+               ! frac_mf fraction used for lateral increase
+               dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
+               fsic(ki)=fsic(ki)+dfsic
+               e_melt=e_melt-dfsic*seaice(ki)
+               ! rest used to increase height
+               seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
+               END IF
+               tslab(ki,1)=t_freeze
+           ELSE ! slab temperature above freezing
+               ! ****** melt ice from below *******
+               ! quantity of melted ice
+               e_melt=(tslab(ki,1)-t_freeze)/cyang & 
+                       /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
+               ! first decrease height to h_thick
+               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
+               seaice(ki)=seaice(ki)-dhsic
+               e_melt=e_melt-fsic(ki)*dhsic
+               IF (e_melt.GT.0) THEN
+               ! frac_mf fraction used for height decrease
+               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
+               seaice(ki)=seaice(ki)-dhsic
+               e_melt=e_melt-fsic(ki)*dhsic
+               ! rest used to decrease fraction (up to 0!)
+               dfsic=MIN(fsic(ki),e_melt/seaice(ki))
+               ! keep remaining in ocean
+               e_melt=e_melt-dfsic*seaice(ki)
+               END IF
+               tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang
+               fsic(ki)=fsic(ki)-dfsic
+           END IF
+        END DO
+        bilg_cum(:)=0.
+    END IF ! coupling time
+    
+    !tests fraction glace
+    WHERE (fsic.LT.ice_frac_min)
+        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
+        tice=t_melt
+        fsic=0.
+        seaice=0.
+    END WHERE
+
+  END SUBROUTINE ocean_slab_ice
 !
 !****************************************************************************************
 !
   SUBROUTINE ocean_slab_final
-  !, seaice_rst etc
-
-    ! For ok_xxx vars (Ekman...)
-    INCLUDE "clesphys.h"
 
 !****************************************************************************************
 ! Deallocate module variables
-!
-!****************************************************************************************
-    IF (ALLOCATED(pctsrf)) DEALLOCATE(pctsrf)
+!****************************************************************************************
     IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
+    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
+    IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
+    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
+    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
+    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
+    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
 
   END SUBROUTINE ocean_slab_final
Index: /LMDZ5/branches/testing/libf/phylmd/ozonecm_m.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/ozonecm_m.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/ozonecm_m.F90	(revision 2220)
@@ -91,5 +91,5 @@
     ozonecm = max(ozonecm, 1e-12)
 
-    print*,'ozonecm Version2'
+!    print*,'ozonecm Version2'
 
   END function ozonecm
Index: /LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2220)
@@ -33,5 +33,5 @@
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
   !$OMP THREADPRIVATE(fder)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: snow   ! snow at surface
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: snow   ! snow at surface
   !$OMP THREADPRIVATE(snow)
   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
@@ -172,5 +172,5 @@
        debut,     lafin,                              &
        rlon,      rlat,      rugoro,   rmu0,          &
-       zsig,      sollwd_m,  pphi,     cldt,          &
+       zsig,      lwdown_m,  pphi,     cldt,          &
        rain_f,    snow_f,    solsw_m,  sollw_m,       &
        t,         q,         u,        v,             &
@@ -182,5 +182,5 @@
        pplay,     paprs,     pctsrf,                  &
        ts,        alb1, alb2,ustar, u10m, v10m,wstar, &
-       lwdown_m,  cdragh,    cdragm,   zu1,    zv1,   &
+       cdragh,    cdragm,   zu1,    zv1,              &
        alb1_m,    alb2_m,    zxsens,   zxevap,        &
        alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
@@ -327,5 +327,5 @@
 ! Martin
     REAL, DIMENSION(klon),        INTENT(IN)        :: zsig    ! slope
-    REAL, DIMENSION(klon),        INTENT(IN)        :: sollwd_m ! net longwave radiation at mean s    
+    REAL, DIMENSION(klon),        INTENT(IN)        :: lwdown_m ! downward longwave radiation at mean s    
     REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud fraction
     REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pphi    ! geopotential (m2/s2)
@@ -367,5 +367,4 @@
 ! Output variables
 !****************************************************************************************
-    REAL, DIMENSION(klon),        INTENT(OUT)       :: lwdown_m   ! Downcoming longwave radiation
     REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
     REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
@@ -780,5 +779,4 @@
     ! Martin
     REAL, DIMENSION(klon, nbsrf)       :: sollwd ! net longwave radiation at surface
-    REAL, DIMENSION(klon)              :: ysollwd
     REAL, DIMENSION(klon)              :: ytoice
     REAL, DIMENSION(klon)              :: ysnowhgt, yqsnow, ysissnow, yrunoff
@@ -855,5 +853,4 @@
 ! 2a) Initialization of all argument variables with INTENT(OUT)
 !****************************************************************************************
- lwdown_m(:)=0.
  cdragh(:)=0. ; cdragm(:)=0.
  zu1(:)=0. ; zv1(:)=0.
@@ -938,5 +935,5 @@
     ! Martin
     ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
-    yalb3_new = 0.0  ; ysissnow = 0.0  ; ysollwd = 0.0
+    yalb3_new = 0.0  ; ysissnow = 0.0 
     ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
     ! Martin
@@ -1109,17 +1106,16 @@
        DO i = 1, klon
           sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
-          ! Martin
-          sollwd(i,nsrf)= sollwd_m(i)
-          ! Martin
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!         ! Martin
+! Apparently introduced for sisvat but not used
+!         sollwd(i,nsrf)= sollwd_m(i)
+!         ! Martin
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
           solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
        ENDDO
     ENDDO
 
-
-! Downwelling longwave radiation at mean surface
-    lwdown_m(:) = 0.0
-    DO i = 1, klon
-       lwdown_m(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
-    ENDDO
 
 !****************************************************************************************
@@ -1180,4 +1176,5 @@
           yagesno(j) = agesno(i,nsrf)
           yfder(j)   = fder(i)
+          ylwdown(j) = lwdown_m(i)
           ysolsw(j)  = solsw(i,nsrf)
           ysollw(j)  = sollw(i,nsrf)
@@ -1703,9 +1700,5 @@
      
        CASE(is_ter)
-          ! ylwdown : to be removed, calculation is now done at land surface in surf_land
-          ylwdown(:)=0.0
-          DO i=1,knon
-             ylwdown(i)=lwdown_m(ni(i))
-          END DO
+!          print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
           CALL surf_land(itap, dtime, date0, jour, knon, ni,&
                rlon, rlat, &
@@ -1746,5 +1739,5 @@
           CALL surf_landice(itap, dtime, knon, ni, &
                rlon, rlat, debut, lafin, &
-               yrmu0, ysollwd, yalb, ypphi(:,1), &
+               yrmu0, ylwdown, yalb, ypphi(:,1), &
                ysolsw, ysollw, yts, ypplay(:,1), &
                ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
Index: /LMDZ5/branches/testing/libf/phylmd/phyaqua_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phyaqua_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phyaqua_mod.F90	(revision 2220)
@@ -238,5 +238,5 @@
     rugos = rugos_omp
     WRITE (*, *) 'iniaqua: rugos=', rugos
-    zmasq(:) = pctsrf(:, is_oce)
+    zmasq(:) = pctsrf(:, is_ter)
 
     ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
@@ -538,34 +538,73 @@
       dims(2) = ntim
 
-      ! cc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
+#else
       ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
+#endif
       ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
-      ! cc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
+
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
+#else
       ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
+#endif
       ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
         'Nature du sol (0,1,2,3)')
-      ! cc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
+
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
+#else
       ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
+#endif
       ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
         'Temperature superficielle de la mer')
-      ! cc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
+
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
+#else
       ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
+#endif
       ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
         'Reference flux de chaleur au sol')
-      ! cc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
+
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
+#else
       ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
+#endif
       ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
-      ! cc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
+
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
+#else
       ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
+#endif
       ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
 
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
+#else
       ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
-      ierr = nf_put_att_text(nid, id_fter, 'title', 8, 'Frac. Terre')
+#endif
+      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
+#else
       ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
-      ierr = nf_put_att_text(nid, id_foce, 'title', 8, 'Frac. Terre')
+#endif
+      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
+#else
       ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
-      ierr = nf_put_att_text(nid, id_fsic, 'title', 8, 'Frac. Terre')
+#endif
+      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
+#ifdef NC_DOUBLE
+      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
+#else
       ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
-      ierr = nf_put_att_text(nid, id_flic, 'title', 8, 'Frac. Terre')
+#endif
+      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
 
       ierr = nf_enddef(nid)
Index: /LMDZ5/branches/testing/libf/phylmd/phyetat0.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phyetat0.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phyetat0.F90	(revision 2220)
@@ -8,9 +8,9 @@
   USE fonte_neige_mod,  ONLY : fonte_neige_init
   USE pbl_surface_mod,  ONLY : pbl_surface_init
-  USE surface_data,     ONLY : type_ocean
+  USE surface_data,     ONLY : type_ocean, version_ocean
   USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
        du_gwd_rando, dv_gwd_rando, entr_therm, f0, falb1, falb2, fm_therm, &
        ftsol, pbl_tke, pctsrf, q_ancien, radpas, radsol, rain_fall, ratqs, &
-       rlat, rlon, rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, &
+       rlat, rlon, rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
        solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
@@ -23,5 +23,5 @@
   USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
   USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
-  USE ocean_slab_mod, ONLY: tslab, ocean_slab_init
+  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
 
   IMPLICIT none
@@ -37,4 +37,5 @@
   include "thermcell.h"
   include "compbl.h"
+  include "YOMCST.h"
   !======================================================================
   CHARACTER*(*) fichnom
@@ -53,4 +54,5 @@
   REAL fractint(klon)
   REAL trs(klon, nbtr)
+  REAL zts(klon)
 
   CHARACTER*6 ocean_in
@@ -513,4 +515,25 @@
   PRINT*, 'Rayonnement IF au sol sollw:', xmin, xmax
 
+  CALL get_field("sollwdown", sollwdown, found)
+  IF (.NOT. found) THEN
+     PRINT*, 'phyetat0: Le champ <sollwdown> est absent'
+     PRINT*, 'mis a zero'
+     sollwdown = 0.
+     zts=0.
+     do nsrf=1,nbsrf
+        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
+     enddo
+     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
+  ENDIF
+!  print*,'TS SOLL',zts(klon/2),sollw(klon/2),sollwdown(klon/2)
+  xmin = 1.0E+20
+  xmax = -1.0E+20
+  DO i = 1, klon
+     xmin = MIN(sollwdown(i), xmin)
+     xmax = MAX(sollwdown(i), xmax)
+  ENDDO
+  PRINT*, 'Rayonnement IF au sol sollwdown:', xmin, xmax
+
+
   ! Lecture derive des flux:
 
@@ -1137,7 +1160,7 @@
   ! Initialize Slab variables
   IF ( type_ocean == 'slab' ) THEN
-      ALLOCATE(tslab(klon,nslay), stat=ierr)
-      IF (ierr /= 0) CALL abort_gcm &
-         ('phyetat0', 'pb allocation tslab', 1)
+      print*, "calling slab_init"
+      CALL ocean_slab_init(dtime, pctsrf)
+      ! tslab
       CALL get_field("tslab", tslab, found)
       IF (.NOT. found) THEN 
@@ -1145,9 +1168,25 @@
           PRINT*, "Initialisation a tsol_oce"
           DO i=1,nslay
-              tslab(:,i)=ftsol(:,is_oce)
+              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
           END DO
       END IF 
-      print*, "calling slab_init"
-      CALL ocean_slab_init(dtime, pctsrf)
+      ! Sea ice variables
+      IF (version_ocean == 'sicINT') THEN
+          CALL get_field("slab_tice", tice, found)
+          IF (.NOT. found) THEN 
+              PRINT*, "phyetat0: Le champ <tice> est absent"
+              PRINT*, "Initialisation a tsol_sic"
+                  tice(:)=ftsol(:,is_sic)
+          END IF 
+          CALL get_field("seaice", seaice, found)
+          IF (.NOT. found) THEN
+              PRINT*, "phyetat0: Le champ <seaice> est absent"
+              PRINT*, "Initialisation a 0/1m suivant fraction glace"
+              seaice(:)=0.
+              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
+                  seaice=917.
+              END WHERE
+          END IF
+      END IF !sea ice INT
   END IF ! Slab        
 
Index: /LMDZ5/branches/testing/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phyredem.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phyredem.F90	(revision 2220)
@@ -16,5 +16,5 @@
   USE indice_sol_mod
   USE surface_data
-  USE ocean_slab_mod, ONLY : tslab
+  USE ocean_slab_mod, ONLY : tslab, seaice, tice, fsic
 
   IMPLICIT none
@@ -107,4 +107,12 @@
   ! BP ajout des fraction de chaque sous-surface
 
+  ! Get last fractions from slab ocean
+  IF (type_ocean == 'slab' .AND. version_ocean == "sicINT") THEN
+      WHERE (1.-zmasq(:).GT.EPSFRA)
+          pctsrf(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
+          pctsrf(:,is_sic)=fsic(:)*(1.-zmasq(:))
+      END WHERE
+  END IF
+
   ! 1. fraction de terre 
 
@@ -209,4 +217,6 @@
 
   CALL put_field("sollw", "Rayonnement IF a la surface", sollw)
+
+  CALL put_field("sollwdown", "Rayonnement down IF a la surface", sollw)
 
   CALL put_field("fder", "Derive de flux", fder)
@@ -350,4 +360,8 @@
   IF (type_ocean == 'slab') THEN
       CALL put_field("tslab", "Slab ocean temperature", tslab)
+      IF (version_ocean == 'sicINT') THEN
+          CALL put_field("seaice", "Slab seaice (kg/m2)", seaice)
+          CALL put_field("slab_tice", "Slab sea ice temperature", tice)
+      END IF
   END IF
 
Index: /LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90	(revision 2220)
@@ -35,4 +35,6 @@
       REAL, SAVE, ALLOCATABLE :: d_t_lsc(:,:),d_q_lsc(:,:),d_ql_lsc(:,:),d_qi_lsc(:,:)
       !$OMP THREADPRIVATE(d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc)
+      REAL, SAVE, ALLOCATABLE :: d_t_lwr(:,:),d_t_lw0(:,:),d_t_swr(:,:),d_t_sw0(:,:)
+      !$OMP THREADPRIVATE(d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0)
       REAL, SAVE, ALLOCATABLE :: d_t_ajsb(:,:), d_q_ajsb(:,:)
       !$OMP THREADPRIVATE(d_t_ajsb, d_q_ajsb)
@@ -393,4 +395,6 @@
       allocate(d_t_wake(klon,klev),d_q_wake(klon,klev))
       allocate(d_t_lsc(klon,klev),d_q_lsc(klon,klev))
+      allocate(d_t_lwr(klon,klev),d_t_lw0(klon,klev))
+      allocate(d_t_swr(klon,klev),d_t_sw0(klon,klev))
       allocate(d_ql_lsc(klon,klev),d_qi_lsc(klon,klev))
       allocate(d_t_ajsb(klon,klev),d_q_ajsb(klon,klev))
@@ -599,4 +603,6 @@
       deallocate(d_t_wake,d_q_wake)
       deallocate(d_t_lsc,d_q_lsc)
+      deallocate(d_t_lwr,d_t_lw0)
+      deallocate(d_t_swr,d_t_sw0)
       deallocate(d_ql_lsc,d_qi_lsc)
       deallocate(d_t_ajsb,d_q_ajsb)
Index: /LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2220)
@@ -465,8 +465,14 @@
   TYPE(ctrl_out), SAVE :: o_slab_bils = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     'slab_bils', 'flux atmos - slab ponderes foce', 'W/m2', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_bilg = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'slab_bilg', 'flux glace - slab ponderes fsic', 'W/m2', (/ ('', i=1, 9) /))
   TYPE(ctrl_out), SAVE :: o_slab_qflux = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     'slab_qflux', 'Correction flux slab', 'W/m2', (/ ('', i=1, 9) /))
   TYPE(ctrl_out), SAVE :: o_tslab = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     'tslab', 'Temperature ocean slab', 'K', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_tice = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'slab_tice', 'Temperature banquise slab', 'K', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_slab_sic = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
+    'seaice', 'Epaisseur banquise slab', 'kg/m2', (/ ('', i=1, 9) /))
   TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), &
     'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /))
@@ -1008,4 +1014,10 @@
       ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'ages_oce',"Snow age", "day", (/ ('', i=1, 9) /)), &
       ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'ages_sic',"Snow age", "day", (/ ('', i=1, 9) /)) /)
+
+  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_snow_srf     = (/ &
+      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_ter', "Snow", "kg/m2", (/ ('', i=1, 9) /)), &
+      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_lic', "Snow", "kg/m2", (/ ('', i=1, 9) /)), &
+      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_oce',"Snow", "kg/m2", (/ ('', i=1, 9) /)), &
+      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_sic',"Snow", "kg/m2", (/ ('', i=1, 9) /)) /)
 
   TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_rugs_srf     = (/ &
Index: /LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90	(revision 2220)
@@ -80,4 +80,5 @@
          o_alp_bl_conv, o_alp_bl_stat, &
          o_slab_qflux, o_tslab, o_slab_bils, &
+         o_slab_bilg, o_slab_sic, o_slab_tice, &
          o_weakinv, o_dthmin, o_cldtau, &
          o_cldemi, o_pr_con_l, o_pr_con_i, &
@@ -113,5 +114,5 @@
          o_rnebls, o_rhum, o_ozone, o_ozone_light, &
          o_dtphy, o_dqphy, o_albe_srf, o_rugs_srf, &
-         o_ages_srf, o_alb1, o_alb2, o_tke, &
+         o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, &
          o_tke_max, o_kz, o_kz_max, o_clwcon, &
          o_dtdyn, o_dqdyn, o_dudyn, o_dvdyn, &
@@ -166,6 +167,6 @@
          wake_deltaq, ftd, fqd, ale_bl_trig, albsol1, &
          rnebcon, wo, falb1, albsol2, coefh, clwcon0, &
-         ratqs, entr_therm, zqasc, detr_therm, f0, heat, &
-         heat0, cool, cool0, lwup, lwdn, lwup0, coefm, &
+         ratqs, entr_therm, zqasc, detr_therm, f0, &
+         lwup, lwdn, lwup0, coefm, &
          swupp, lwupp, swup0p, lwup0p, swdnp, lwdnp, &
          swdn0p, lwdn0p, tnondef, O3sumSTD, uvsumSTD, &
@@ -215,4 +216,5 @@
          d_u_ajs, d_v_ajs, &
          d_u_con, d_v_con, d_q_con, d_q_ajs, d_t_lsc, &
+         d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
          d_t_eva, d_q_lsc, beta_prec, d_t_lscth, &
          d_t_lscst, d_q_lscth, d_q_lscst, plul_th, &
@@ -226,9 +228,10 @@
          bils_ec,bils_ech, bils_tke, bils_kinetic, bils_latent, bils_enthalp, &
          itau_con, nfiles, clef_files, nid_files, zvstr_gwd_rando
-    USE ocean_slab_mod, only: tslab, slab_bils
+    USE ocean_slab_mod, only: tslab, slab_bils, slab_bilg, tice, seaice
+    USE pbl_surface_mod, only: snow
     USE indice_sol_mod, only: nbsrf
     USE infotrac, only: nqtot, nqo, type_trac
     USE comgeomphy, only: airephy
-    USE surface_data, only: type_ocean, ok_veget, ok_snow
+    USE surface_data, only: type_ocean, version_ocean, ok_veget, ok_snow
 !    USE aero_mod, only: naero_spc
     USE aero_mod, only: naero_tot, id_STRAT_phy
@@ -399,5 +402,5 @@
        CALL histwrite_phy(o_pluc, zx_tmp_fi2d)
        CALL histwrite_phy(o_snow, snow_fall)
-       CALL histwrite_phy(o_msnow, snow_o)
+       CALL histwrite_phy(o_msnow, zxsnow)
        CALL histwrite_phy(o_fsnow, zfra_o)
        CALL histwrite_phy(o_evap, evap)
@@ -513,5 +516,5 @@
 
        IF (ok_snow) THEN
-          CALL histwrite_phy(o_snowsrf, zxsnow)
+          CALL histwrite_phy(o_snowsrf, snow_o)
           CALL histwrite_phy(o_qsnow, qsnow)
           CALL histwrite_phy(o_snowhgt,snowhgt)
@@ -754,4 +757,9 @@
               CALL histwrite_phy(o_tslab, tslab)
           END IF
+          IF (version_ocean=='sicINT') THEN
+              CALL histwrite_phy(o_slab_bilg, slab_bilg)
+              CALL histwrite_phy(o_slab_tice, tice)
+              CALL histwrite_phy(o_slab_sic, seaice)
+          END IF
        ENDIF !type_ocean == force/slab
        CALL histwrite_phy(o_weakinv, weak_inversion)
@@ -969,4 +977,6 @@
           IF (vars_defined) zx_tmp_fi2d(1 : klon) = agesno( 1 : klon, nsrf)
           CALL histwrite_phy(o_ages_srf(nsrf), zx_tmp_fi2d)
+          IF (vars_defined) zx_tmp_fi2d(1 : klon) = snow( 1 : klon, nsrf)
+          CALL histwrite_phy(o_snow_srf(nsrf), zx_tmp_fi2d)
        ENDDO !nsrf=1, nbsrf
        CALL histwrite_phy(o_alb1, albsol1)
@@ -1127,11 +1137,11 @@
        IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dqajs, zx_tmp_fi3d)
-       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY
+       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dtswr, zx_tmp_fi3d)
-       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=heat0(1:klon,1:klev)/RDAY
+       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_sw0(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dtsw0, zx_tmp_fi3d)
-       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=-1.*cool(1:klon,1:klev)/RDAY
+       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lwr(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dtlwr, zx_tmp_fi3d)
-       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=-1.*cool0(1:klon,1:klev)/RDAY
+       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lw0(1:klon,1:klev)/pdtphys
        CALL histwrite_phy(o_dtlw0, zx_tmp_fi3d)
        IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_ec(1:klon,1:klev)/pdtphys
@@ -1193,6 +1203,6 @@
        CALL histwrite_phy(o_tnt, zx_tmp_fi3d)
        IF(vars_defined) THEN
-          zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY - &
-               cool(1:klon,1:klev)/RDAY
+          zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys + &
+               d_t_lwr(1:klon,1:klev)/pdtphys
        ENDIF
        CALL histwrite_phy(o_tntr, zx_tmp_fi3d)
Index: /LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90	(revision 2220)
@@ -60,4 +60,6 @@
       REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:)
 !$OMP THREADPRIVATE(clwcon,rnebcon)
+      REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:)
+!$OMP THREADPRIVATE(qtc_cv,sigt_cv)
       REAL, ALLOCATABLE, SAVE :: ratqs(:,:)
 !$OMP THREADPRIVATE(ratqs)
@@ -416,4 +418,5 @@
 !!! Rom P <<<
       ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
+      ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev))
       ALLOCATE(ratqs(klon,klev))
       ALLOCATE(pbl_tke(klon,klev+1,nbsrf+1))
@@ -566,4 +569,5 @@
       deallocate(zthe, zpic, zval)
       deallocate(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
+      deallocate(qtc_cv,sigt_cv)
       deallocate(        u_ancien, v_ancien                 )
       deallocate(        tr_ancien)                           !RomP
Index: /LMDZ5/branches/testing/libf/phylmd/physiq.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/physiq.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/physiq.F90	(revision 2220)
@@ -636,7 +636,7 @@
   !$OMP THREADPRIVATE(fact_cldcon,facttemps)
 
-  integer iflag_cldcon
-  save iflag_cldcon
-  !$OMP THREADPRIVATE(iflag_cldcon)
+  integer iflag_cldth
+  save iflag_cldth
+  !$OMP THREADPRIVATE(iflag_cldth)
   logical ptconv(klon,klev)
   !IM cf. AM 081204 BEG
@@ -913,5 +913,5 @@
           solarlong0,seuil_inversion, &
           fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
-          iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
+          iflag_cldth,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
           ok_ade, ok_aie, ok_cdnc, aerosol_couple,  &
           flag_aerosol, flag_aerosol_strat, new_aod, &
@@ -1014,6 +1014,6 @@
      print*,'CYCLE_DIURNE', cycle_diurne
      !
-     IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
-        abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
+     IF (iflag_con.EQ.2.AND.iflag_cldth.GT.-1) THEN
+        abort_message = 'Tiedtke needs iflag_cldth=-2 or -1'
         CALL abort_gcm (modname,abort_message,1)
      ENDIF
@@ -1130,5 +1130,5 @@
                 ,alp_bl_prescr, ale_bl_prescr)
            ! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
-           !        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
+           !        print*,'apres ini_wake iflag_cldth=', iflag_cldth
         endif
 
@@ -1308,4 +1308,6 @@
              annee_ref, &
              day_ref,  &
+             day_ini, &
+             start_time, &
              itau_phy, &
              io_lon, &
@@ -1811,7 +1813,7 @@
 !>nrlmd+jyg
           pplay,     paprs,     pctsrf,             &
-          ftsol,falb1,falb2,ustar,u10m,v10m,wstar, &
-          sollwdown, cdragh,    cdragm,  u1,    v1, &
-          albsol1,   albsol2,   sens,    evap,   &
+          ftsol,falb1,falb2,ustar,u10m,v10m,wstar,  &
+          cdragh,    cdragm,  u1,    v1,            &
+          albsol1,   albsol2,   sens,    evap,      &
           albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
           zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
@@ -2168,5 +2170,6 @@
              ftd,fqd,lalim_conv,wght_th, &
              ev, ep,epmlmMm,eplaMm, &
-             wdtrainA,wdtrainM,wght_cvfd)
+             wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
+             tau_cld_cv,coefw_cld_cv)
         ! RomP <<<
 
@@ -2218,6 +2221,12 @@
      !   calcul des proprietes des nuages convectifs
      clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
+     IF (iflag_cld_cv <= 1) THEN
      call clouds_gno &
           (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
+     ELSE
+     call clouds_bigauss &
+          (klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
+     ENDIF
+
 
      ! =================================================================== c
@@ -2452,5 +2461,5 @@
   END IF
 
-  !      print*,'apres callwake iflag_cldcon=', iflag_cldcon
+  !      print*,'apres callwake iflag_cldth=', iflag_cldth
   !
   !===================================================================
@@ -2622,5 +2631,5 @@
               enddo
 
-           ELSE IF (iflag_trig_bl.eq.2) then
+           ELSE IF (iflag_trig_bl.ge.2) then
 
               do i=1,klon
@@ -2773,5 +2782,5 @@
   ! water distribution
   CALL  calcratqs(klon,klev,prt_level,lunout,        &
-       iflag_ratqs,iflag_con,iflag_cldcon,pdtphys,  &
+       iflag_ratqs,iflag_con,iflag_cldth,pdtphys,  &
        ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,   &
        ptconv,ptconvth,clwcon0th, rnebcon0th,     &
@@ -2795,5 +2804,5 @@
        frac_impa, frac_nucl, beta_prec_fisrt, &
        prfl, psfl, rhcl,  &
-       zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon, &
+       zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldth, &
        iflag_ice_thermo)
   !
@@ -2851,10 +2860,10 @@
   !
   !IM cf FH
-  !     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
-  IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
+  !     IF (iflag_cldth.eq.-1) THEN ! seulement pour Tiedtke
+  IF (iflag_cldth.le.-1) THEN ! seulement pour Tiedtke
      snow_tiedtke=0.
      !     print*,'avant calcul de la pseudo precip '
-     !     print*,'iflag_cldcon',iflag_cldcon
-     if (iflag_cldcon.eq.-1) then
+     !     print*,'iflag_cldth',iflag_cldth
+     if (iflag_cldth.eq.-1) then
         rain_tiedtke=rain_con
      else
@@ -2889,5 +2898,5 @@
      ENDDO
 
-  ELSE IF (iflag_cldcon.ge.3) THEN
+  ELSE IF (iflag_cldth.ge.3) THEN
      !  On prend pour les nuages convectifs le max du calcul de la
      !  convection et du calcul du pas de temps precedent diminue d'un facteur
@@ -2932,4 +2941,6 @@
              tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
              tausum_aero, tau3d_aero)
+
+             CALL aeropt_lw_rrtm
 #else
 
@@ -2976,5 +2987,5 @@
      !   On prend la somme des fractions nuageuses et des contenus en eau
 
-     if (iflag_cldcon>=5) then
+     if (iflag_cldth>=5) then
 
         do k=1,klev
@@ -3132,5 +3143,5 @@
      calday = REAL(days_elapsed + 1) + jH_cur
 
-     call chemtime(itap+itau_phy-1, date0, dtime)
+     call chemtime(itap+itau_phy-1, date0, dtime, itap)
      IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
         CALL AEROSOL_METEO_CALC( &
@@ -3455,10 +3466,11 @@
   ! Ajouter la tendance des rayonnements (tous les pas)
   !
-  DO k = 1, klev
-     DO i = 1, klon
-        t_seri(i,k) = t_seri(i,k) &
-             + (heat(i,k)-cool(i,k)) * dtime/RDAY
-     ENDDO
-  ENDDO
+  d_t_swr(:,:)=heat(:,:)*dtime/RDAY
+  d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
+  d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
+  d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
+  CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW')
+  CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW')
+
   !
   if (mydebug) then
Index: /LMDZ5/branches/testing/libf/phylmd/phytrac_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/phytrac_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/phytrac_mod.F90	(revision 2220)
@@ -310,6 +310,11 @@
     !RomP >>>
     INTEGER,SAVE  :: iflag_lscav_omp,iflag_lscav
+    REAL, SAVE ::   ccntrAA_in,ccntrAA_omp
+    REAL, SAVE ::   ccntrENV_in,ccntrENV_omp
+    REAL, SAVE ::   coefcoli_in,coefcoli_omp
+
     LOGICAL,SAVE  :: convscav_omp,convscav
 !$OMP THREADPRIVATE(iflag_lscav)
+!$OMP THREADPRIVATE(ccntrAA_in,ccntrENV_in,coefcoli_in)
 !$OMP THREADPRIVATE(convscav)
     !RomP <<<
@@ -412,7 +417,16 @@
        iflag_lscav_omp=1
        call getin('iflag_lscav', iflag_lscav_omp)
+       ccntrAA_omp=1
+       ccntrENV_omp=1.
+       coefcoli_omp=0.001
+       call getin('ccntrAA', ccntrAA_omp)
+       call getin('ccntrENV', ccntrENV_omp)
+       call getin('coefcoli', coefcoli_omp)
 !$OMP END MASTER
 !$OMP BARRIER
        iflag_lscav=iflag_lscav_omp
+       ccntrAA_in=ccntrAA_omp
+       ccntrENV_in=ccntrENV_omp
+       coefcoli_in=coefcoli_omp
        !
        SELECT CASE(iflag_lscav)
@@ -463,7 +477,7 @@
                 IF (convscav.and.aerosol(it)) THEN
                    flag_cvltr(it)=.true.
-                   ccntrAA(it) =1.0         !--a modifier par JYG a lire depuis fichier 
-                   ccntrENV(it)=1.0 
-                   coefcoli(it)=0.001
+                   ccntrAA(it) =ccntrAA_in    !--a modifier par JYG a lire depuis fichier 
+                   ccntrENV(it)=ccntrENV_in
+                   coefcoli(it)=coefcoli_in
                 ELSE 
                    flag_cvltr(it)=.false.
@@ -613,4 +627,6 @@
                 !--with the full array tr_seri even if only item it is processed
 
+                print*,'CV SCAV ',it,ccntrAA(it),ccntrENV(it)
+
                 CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep,         &
                      sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,              &     
@@ -747,4 +763,6 @@
           !
           DO it = 1, nbtr
+
+             IF (aerosol(it)) THEN
              !  incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
              ! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
@@ -763,4 +781,6 @@
              ENDDO
              CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//solsym(it))
+             ENDIF
+
           END DO  !tr
 
Index: /LMDZ5/branches/testing/libf/phylmd/radlwsw_m.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/radlwsw_m.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/radlwsw_m.F90	(revision 2220)
@@ -765,5 +765,5 @@
 ! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
          RII0=solaire/zdist/zdist
-	print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
+!print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
 !  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ! Ancien appel a RECMWF (celui du cy25)
@@ -819,5 +819,5 @@
          ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat) ! flags aerosols
             
-         print *,'RADLWSW: apres RECMWF'
+!        print *,'RADLWSW: apres RECMWF'
       if(lldebug) then
         CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
@@ -954,5 +954,5 @@
          zsollwdown(i)= ZFLDN(i,1)
       ENDDO
-      print*,'OK2'
+!     print*,'OK2'
 
 ! extrait de SW_AR4
Index: /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_lw_rrtm.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_lw_rrtm.F90	(revision 2220)
+++ /LMDZ5/branches/testing/libf/phylmd/rrtm/aeropt_lw_rrtm.F90	(revision 2220)
@@ -0,0 +1,13 @@
+!
+! aeropt_lw_rrtm.F90 2014-05-13 C. Kleinschmitt
+!
+
+SUBROUTINE AEROPT_LW_RRTM
+
+  USE phys_state_var_mod, ONLY : tau_aero_lw_rrtm
+
+  IMPLICIT NONE
+
+  tau_aero_lw_rrtm(:,:,:,:)=0.0
+
+END SUBROUTINE AEROPT_LW_RRTM
Index: /LMDZ5/branches/testing/libf/phylmd/rrtm/radlsw.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/rrtm/radlsw.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/rrtm/radlsw.F90	(revision 2220)
@@ -1062,5 +1062,5 @@
 !                 ------------------------------------
 
-print *,'RADLSW: LPHYLIN, LRRTM',LPHYLIN, LRRTM
+!print *,'RADLSW: LPHYLIN, LRRTM',LPHYLIN, LRRTM
 IF (.NOT.LPHYLIN) THEN
   IF ( .NOT. LRRTM) THEN
@@ -1074,5 +1074,5 @@
      & ZEMIT , PFLUX , PFLUC &
      & )  
-    print *,'RADLSW: apres CALL LW'
+!   print *,'RADLSW: apres CALL LW'
     IF(LLDEBUG) THEN
     call writefield_phy('radlsw_flux1',PFLUX(:,1,:),klev+1)
@@ -1101,5 +1101,5 @@
     ENDDO
 
-    print *,'RADLSW: avant CALL RRTM_RRTM_140GP,PAP=',PAP(1,:)
+!   print *,'RADLSW: avant CALL RRTM_RRTM_140GP,PAP=',PAP(1,:)
     CALL RRTM_RRTM_140GP &
      & ( KIDIA , KFDIA , KLON  , KLEV,&
@@ -1111,5 +1111,5 @@
      & PTAU_LW,&
      & ZEMIT , PFLUX , PFLUC , ZTCLEAR )
-    print *,'RADLSW: apres CALL RRTM_RRTM_140GP'
+!   print *,'RADLSW: apres CALL RRTM_RRTM_140GP'
 
   ENDIF
@@ -1118,5 +1118,5 @@
   PFLUX(:,:,:)= 0.0_JPRB
   PFLUC(:,:,:)= 0.0_JPRB
-  print *,'RADLSW: ZEMIT,PFLUX et PFLUC = 0'
+! print *,'RADLSW: ZEMIT,PFLUX et PFLUC = 0'
 ENDIF
 
Index: /LMDZ5/branches/testing/libf/phylmd/rrtm/sw1s.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/rrtm/sw1s.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/rrtm/sw1s.F90	(revision 2220)
@@ -304,5 +304,5 @@
 
 ELSEIF (NSW == 6) THEN
-print *,'... dans SW1S: NSW=',NSW
+!print *,'... dans SW1S: NSW=',NSW
 
 !*         3.2   SIX SPECTRAL INTERVALS
Index: /LMDZ5/branches/testing/libf/phylmd/surf_land_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/surf_land_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/surf_land_mod.F90	(revision 2220)
@@ -85,5 +85,4 @@
     REAL, DIMENSION(klon) :: pref_tmp
     REAL, DIMENSION(klon) :: swdown     ! downwelling shortwave radiation at land surface
-    REAL, DIMENSION(klon) :: lwdown     ! downwelling longwave radiation at land surface
     REAL, DIMENSION(klon) :: epot_air           ! potential air temperature
     REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
@@ -106,11 +105,9 @@
        pref_tmp(1:knon)  = pref(1:knon)/100.
 ! 
-!* Calculate incoming flux for SW and LW interval: swdown, lwdown
+!* Calculate incoming flux for SW and LW interval: swdown
 !
        swdown(:) = 0.0
-       lwdown(:) = 0.0
        DO i = 1, knon
           swdown(i) = swnet(i)/(1-albedo(i))
-          lwdown(i) = lwnet(i) + RSIGMA*tsurf(i)**4
        END DO
 !
Index: /LMDZ5/branches/testing/libf/phylmd/surf_ocean_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/surf_ocean_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/surf_ocean_mod.F90	(revision 2220)
@@ -123,5 +123,5 @@
             AcoefU, AcoefV, BcoefU, BcoefV, &
             ps, u1, v1, tsurf_in, &
-            radsol, snow, agesno, &
+            radsol, snow, &
             qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
             tsurf_new, dflux_s, dflux_l, lmt_bils)
Index: /LMDZ5/branches/testing/libf/phylmd/surf_seaice_mod.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/surf_seaice_mod.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/surf_seaice_mod.F90	(revision 2220)
@@ -25,4 +25,5 @@
   USE ocean_forced_mod, ONLY : ocean_forced_ice
   USE ocean_cpl_mod, ONLY    : ocean_cpl_ice
+  USE ocean_slab_mod, ONLY   : ocean_slab_ice
   USE indice_sol_mod
 
@@ -108,15 +109,15 @@
             tsurf_new, dflux_s, dflux_l)
        
-!    ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
-!       CALL ocean_slab_ice( & 
-!          itime, dtime, jour, knon, knindex, &
-!          debut, &
-!          tsurf, p1lay, cdragh, precip_rain, precip_snow, temp_air, spechum,&
-!          AcoefH, AcoefQ, BcoefH, BcoefQ, &
-!          ps, u1, v1, &
-!          radsol, snow, qsurf, qsol, agesno, tsoil, &
-!          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
-!          tsurf_new, dflux_s, dflux_l)
-!
+    ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
+       CALL ocean_slab_ice( & 
+          itime, dtime, jour, knon, knindex, &
+          tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum,&
+          AcoefH, AcoefQ, BcoefH, BcoefQ, &
+            AcoefU, AcoefV, BcoefU, BcoefV, &
+          ps, u1, v1, &
+          radsol, snow, qsurf, qsol, agesno, &
+          alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
+          tsurf_new, dflux_s, dflux_l, swnet)
+
       ELSE ! type_ocean=force or slab +sicOBS or sicNO
        CALL ocean_forced_ice( &
Index: /LMDZ5/branches/testing/libf/phylmd/surface_data.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/surface_data.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/surface_data.F90	(revision 2220)
@@ -7,6 +7,4 @@
   REAL, PARAMETER        :: tau_gl=86400.*5.
   REAL, PARAMETER        :: calsno=1./(2.3867e+06*.15)
-  REAL, PARAMETER        :: t_freeze=271.35
-  REAL, PARAMETER        :: t_melt=273.15
   
   LOGICAL, SAVE          :: ok_veget      ! true for use of vegetation model ORCHIDEE
Index: /LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90	(revision 2220)
@@ -976,4 +976,8 @@
       enddo
 
+! Ale sec (max de wmax/2 sous la zone d'inhibition) dans
+! le cas iflag_trig_bl=3
+      IF (iflag_trig_bl==3) Ale_bl(:)=0.5*wmax_sec(:)**2
+
 !test:calcul de la ponderation des couches pour KE
 !initialisations
Index: /LMDZ5/branches/testing/libf/phylmd/thermcell_plume.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/thermcell_plume.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/thermcell_plume.F90	(revision 2220)
@@ -359,6 +359,4 @@
                    zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
     &                    +(1.-fact_shell)*zbuoy(ig,l)
-         print*,'on est pass?? par l??1',l,lt,ztv1,ztv2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
-   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
                 elseif (zlmelup.ge.zinv) then
                  ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
@@ -370,12 +368,8 @@
     &            ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
 
-         print*,'on est pass?? par l??2',l,lt,ztv_est1,ztv_est2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
-   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
                 else 
                    ztv_est(ig,l)=atv1*zlmel+btv1
                    zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & 
     &                           +(1.-fact_shell)*zbuoy(ig,l)
-                 print*,'on est pass?? par l??3',l,lt,ztv1,ztv2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
-   &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l) 
                 endif
 
@@ -392,6 +386,4 @@
                 endif
 
-                 print*,'on est pass?? par l??4',l,lt,ztv1,ztv2,ztv(ig,lt),ztv(ig,l),ztva_est(ig,l), &
-   &                                   zlmelup,zbuoy(ig,l),zbuoyjam(ig,l) 
 !          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & 
 !    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
@@ -417,6 +409,4 @@
     &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
     &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
-         print*,'on est pass?? par l??',l,lt,ztv(ig,lt),ztva_est(ig,l),ztv(ig,l), &
-   &                                   zbuoy(ig,l),zbuoyjam(ig,l)
         endif !   if (iflag_thermals_ed.lt.8) then
 
Index: /LMDZ5/branches/testing/libf/phylmd/tilft43.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/tilft43.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/tilft43.F90	(revision 2220)
@@ -3,7 +3,11 @@
 
 SUBROUTINE tlift43(p, t, q, qs, gz, icb, nk, tvp, tpk, clw, nd, nl, kk)
+  IMPLICIT NONE
   REAL gz(nd), tpk(nd), clw(nd), p(nd)
   REAL t(nd), q(nd), qs(nd), tvp(nd), lv0
-
+  INTEGER icb, nk, nd, nl, kk
+  REAL cpd, cpv,  cl, g, rowl, gravity, cpvmcl, eps, epsi
+  REAL ah0, cpp, cpinv, tg, qg, alv, s, ahg, tc, denom, es
+  INTEGER i, nst, nsb, j 
   ! ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
 
Index: /LMDZ5/branches/testing/libf/phylmd/tlift.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/tlift.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/tlift.F90	(revision 2220)
@@ -4,10 +4,10 @@
 SUBROUTINE tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tpk, clw, nd, nl, &
     dtvpdt1, dtvpdq1)
-
+  IMPLICIT NONE
   ! Argument NK ajoute (jyg) = Niveau de depart de la
   ! convection
-
-  PARAMETER (na=60)
-  REAL gz(nd), tpk(nd), clw(nd)
+  INTEGER icb, nk, nd, nl
+  INTEGER,PARAMETER :: na=60
+  REAL gz(nd), tpk(nd), clw(nd), plcl
   REAL t(nd), rr(nd), rs(nd), tvp(nd), p(nd)
   REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual
@@ -17,5 +17,8 @@
   REAL dtpdt1(na), dtpdq1(na) ! Derivatives of parcel temperature
   ! wrt T1 and Q1
-
+  REAL gravity, cpd, cpv, cl, ci, cpvmcl, clmci, eps, alv0, alf0
+  REAL cpp, cpinv, ah0, alf, tg, s, ahg, tc, denom, alv, es, esi
+  REAL qsat_new, snew
+  INTEGER icbl, i, imin, j, icb1
 
   LOGICAL ice_conv
Index: /LMDZ5/branches/testing/libf/phylmd/wake.F90
===================================================================
--- /LMDZ5/branches/testing/libf/phylmd/wake.F90	(revision 2219)
+++ /LMDZ5/branches/testing/libf/phylmd/wake.F90	(revision 2220)
@@ -1756,5 +1756,5 @@
   ! a une humidite positive dans la zone (x) et dans la zone (w).
   ! ------------------------------------------------------
-
+  IMPLICIT NONE
 
   ! Input
@@ -1772,4 +1772,5 @@
   REAL epsilon
   ! DATA epsilon/1.e-15/
+  INTEGER i,k
 
   DO k = 1, nl
Index: /LMDZ5/branches/testing/makelmdz
===================================================================
--- /LMDZ5/branches/testing/makelmdz	(revision 2219)
+++ /LMDZ5/branches/testing/makelmdz	(revision 2220)
@@ -311,4 +311,10 @@
 else
   FLAG_PARA="$paramem"
+  if [[ $paramem == par ]]
+  then
+      echo "The version of the dynamics in dyn3dpar is no longer updated."
+      echo "You should use option \"-mem\"."
+      exit 1
+  fi     
 fi
 
@@ -470,4 +476,9 @@
 cd $LIBFGCM/grid/dimension
 ./makdim $dim
+if (($? != 0))
+then
+    exit 1
+fi
+
 cat $LIBFGCM/grid/dimensions.h
 cd $LMDGCM
Index: /LMDZ5/branches/testing/makelmdz_fcm
===================================================================
--- /LMDZ5/branches/testing/makelmdz_fcm	(revision 2219)
+++ /LMDZ5/branches/testing/makelmdz_fcm	(revision 2220)
@@ -12,4 +12,5 @@
 #
 ##set -x
+set -e
 ########################################################################
 # options par defaut pour la commande make
@@ -461,5 +462,9 @@
   if [[ "$paramem" == "mem" ]]
   then
-   SUFF_NAME=${SUFF_NAME}_${paramem}
+      SUFF_NAME=${SUFF_NAME}_${paramem}
+  else
+      echo "The version of the dynamics in dyn3dpar is no longer updated."
+      echo "You should use option \"-mem\"."
+      exit 1
   fi
 else
