Index: dynamico_lmdz/simple_physics/phyparam/param/coefdifv.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/coefdifv.F	(revision 4201)
+++ dynamico_lmdz/simple_physics/phyparam/param/coefdifv.F	(revision 4202)
@@ -3,5 +3,5 @@
      $               pcdzv,pcdzh)
       USE constlim, ONLY : ccdzh, cdzmin, cpgam, cdzconst, dgcdrag
-      USE vdif_mod, ONLY : lmixmin
+      USE turbulence, ONLY : lmixmin
       USE phys_const, ONLY : g
       IMPLICIT NONE
Index: dynamico_lmdz/simple_physics/phyparam/param/inifrict.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/inifrict.F	(revision 4201)
+++ 	(revision )
@@ -1,65 +1,0 @@
-      SUBROUTINE inifrict(timestep)
-      USE comgeomfi, ONLY : ngridmax,nlayermx
-      USE constlim
-      USE phys_const
-      IMPLICIT NONE
-c=======================================================================
-c
-c   Calcul des coefficients de ls diffusion verticale
-c
-c=======================================================================
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-      
-
-c   local:
-c   ------
-
-      REAL dgrad,cpr,rl2,rrr,timestep
-      INTEGER l,ij
-
-c-----------------------------------------------------------------------
-
-      call initconstlim
-
-      dtradia=timestep
-      PRINT*,'DTPHYS',dtradia
-      lmixmin=100.
-      ais1    = 1.
-      ais2    = ais1 - 1.
-      print*,'ais1',ais1,'ais2',ais2
-      cdzmin  = 1.e-6
-      OPEN(99,file='cdzmin',status='old',err=9999)
-      READ(99,*) cdzmin
-9999  CLOSE(99)
-      PRINT*,'cdzmin=',cdzmin
-
-      cpr     = cpp/ r
-      ccdzh   = 2.5*g
-c!!!  cpgam   = 5.e-3*cpp
-      cpgam=0.
-
-c-----------------------------------------------------------------------
-c   coefficient de diffusion dans l'atmosphere:
-c   -------------------------------------------
-
-      rl2=lmixmin**2
-      cdzconst(1)= 0.
-      DO 15 l=1,nlayermx - 1
-         cdzconst(l+1)= dtradia*g*g*cpr*rl2
-         print*,'cdzconst(',l+1,')  =  ',cdzconst(l+1)
-  15  CONTINUE
-
-c-----------------------------------------------------------------------
-c   couche limite de surface:
-c   -------------------------
-
-      cdrat   = 2.e-3
-      dgrad   = dtradia*g*cpp/r
-      DO 16 ij = 1, ngridmax
-         dgcdrag( ij ) = cdrat * dgrad
-  16  CONTINUE
-
-      RETURN
-      END
Index: dynamico_lmdz/simple_physics/phyparam/param/inifrict.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/inifrict.F90	(revision 4202)
+++ dynamico_lmdz/simple_physics/phyparam/param/inifrict.F90	(revision 4202)
@@ -0,0 +1,65 @@
+      SUBROUTINE inifrict(timestep) 
+      USE comgeomfi, ONLY : ngridmax,nlayermx 
+      USE constlim 
+      USE phys_const 
+      IMPLICIT NONE 
+!=======================================================================
+!                                                                       
+!   Calcul des coefficients de ls diffusion verticale                   
+!                                                                       
+!=======================================================================
+!-----------------------------------------------------------------------
+!   Declarations:                                                       
+!   -------------                                                       
+                                                                        
+                                                                        
+!   local:                                                              
+!   ------                                                              
+                                                                        
+      REAL dgrad,cpr,rl2,rrr,timestep 
+      INTEGER l,ij 
+                                                                        
+!-----------------------------------------------------------------------
+                                                                        
+      call initconstlim 
+                                                                        
+      dtradia=timestep 
+      PRINT*,'DTPHYS',dtradia 
+      lmixmin=100. 
+      ais1    = 1. 
+      ais2    = ais1 - 1. 
+      print*,'ais1',ais1,'ais2',ais2 
+      cdzmin  = 1.e-6 
+      OPEN(99,file='cdzmin',status='old',err=9999) 
+      READ(99,*) cdzmin 
+ 9999 CLOSE(99) 
+      PRINT*,'cdzmin=',cdzmin 
+                                                                        
+      cpr     = cpp/ r 
+      ccdzh   = 2.5*g 
+!!!!  cpgam   = 5.e-3*cpp                                               
+      cpgam=0. 
+                                                                        
+!-----------------------------------------------------------------------
+!   coefficient de diffusion dans l'atmosphere:                         
+!   -------------------------------------------                         
+                                                                        
+      rl2=lmixmin**2 
+      cdzconst(1)= 0. 
+      DO 15 l=1,nlayermx - 1 
+         cdzconst(l+1)= dtradia*g*g*cpr*rl2 
+         print*,'cdzconst(',l+1,')  =  ',cdzconst(l+1) 
+   15 END DO 
+                                                                        
+!-----------------------------------------------------------------------
+!   couche limite de surface:                                           
+!   -------------------------                                           
+                                                                        
+      cdrat   = 2.e-3 
+      dgrad   = dtradia*g*cpp/r 
+      DO 16 ij = 1, ngridmax 
+         dgcdrag( ij ) = cdrat * dgrad 
+   16 END DO 
+                                                                        
+      RETURN 
+      END                                           
Index: dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F	(revision 4201)
+++ dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F	(revision 4202)
@@ -14,5 +14,5 @@
       USE planet, ONLY : coefir, coefvis
       USE astronomy
-      USE vdif_mod, ONLY : lmixmin, emin_turb 
+      USE turbulence, ONLY : lmixmin, emin_turb 
       IMPLICIT NONE
 
Index: dynamico_lmdz/simple_physics/phyparam/param/interpol.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/interpol.F	(revision 4201)
+++ 	(revision )
@@ -1,95 +1,0 @@
-      SUBROUTINE interpol(xinp,yinp,ninp,xout,yout,nout)
-      IMPLICIT NONE
-
-c=======================================================================
-c
-c  objet:
-c  interpolation lineaire
-c
-c  ATTENTION!:
-c  les tableaux x doivent etre croissants
-c
-c  Arguments:
-c  ----------
-c
-c  entree:
-c  -------
-c  ninp        INT  nombre de points des donnes
-c  xinp(ninp)  REAL abbcisses des donnees (doivent etre croissantes)
-c  yinp(ninp)  REAL valeurs des donnees
-c
-c  sorties:
-c  --------
-c  nout        INT  nombre de points interpolles
-c  xout(nout)  REAL abbcisses (doivent etre croissantes)
-c  yout(nout)  REAL valeurs interpollees
-c
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c    0. Declarations :
-c    -----------------
-
-      INTEGER ninp,nout
-      INTEGER iinp,iout,iout0,iout1
-      REAL xinp(ninp),yinp(ninp)
-      REAL xout(nout),yout(nout)
-
-c-----------------------------------------------------------------------
-c   traitement des valeurs telles que xout<xinp(1):
-c   -----------------------------------------------
-
-      DO 100 iout=1,nout
-         IF(xout(iout).lt.xinp(1)) THEN
-              yout(iout)=yinp(1)
-c   si on veut extrapoler aux bords:
-c    1        +(xout(iout)-xinp(1))*(yinp(1)-yinp(2))/(xinp(1)-xinp(2))
-          ELSE
-              iout0=iout
-              GO TO 200
-          ENDIF
-100   CONTINUE
-
-      RETURN
-
-200   CONTINUE
-
-c-----------------------------------------------------------------------
-c   traitement du centre:
-c   ---------------------
-
-      iinp=1
-      DO 300 iout=iout0,nout
-          IF(xout(iout).ge.xinp(ninp)) THEN 
-             iout1=iout
-             GO TO 400
-          ENDIF
-320       CONTINUE
-          IF(  (xout(iout).ge.xinp(iinp)  ) .AND.
-     1         (xout(iout).lt.xinp(iinp+1))        ) THEN
-               yout(iout)=yinp(iinp)+
-     1         (xout(iout)-xinp(iinp))*
-     2         (yinp(iinp+1)-yinp(iinp))/(xinp(iinp+1)-xinp(iinp))
-          ELSE
-               iinp=iinp+1
-               GO TO 320
-          ENDIF
-300   CONTINUE
-
-      RETURN
-
-400   CONTINUE
-
-c-----------------------------------------------------------------------
-c   traitement des points tels que xout>xinp(ninp):
-c   -----------------------------------------------
-
-      DO 500 iout=iout1,nout
-          yout(iout)=yinp(ninp)
-c   si on veut extrapoler aux bords:
-c    1    +(xout(iout)-xinp(ninp))*
-c    2    (yinp(ninp-1)-yinp(ninp))/(xinp(ninp-1)-xinp(ninp))
-500   CONTINUE
-
-      RETURN
-      END
Index: dynamico_lmdz/simple_physics/phyparam/param/mellor_yamada.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/mellor_yamada.F	(revision 4202)
+++ dynamico_lmdz/simple_physics/phyparam/param/mellor_yamada.F	(revision 4202)
@@ -0,0 +1,247 @@
+      subroutine my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
+
+***********************************************************************
+******* FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) *******
+** q2 au interfaces entre mailles.
+***********************************************************************
+
+
+************** DECLARATIONS *******************************************
+
+
+      integer imax,kmax
+
+      parameter (impmax=5)
+
+      real kappa,khmin,kmmin,kqmin,longc
+
+      parameter (kappa=0.4)
+      parameter (a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08,
+     a           e1=1.8, e2=1.33)
+      parameter (khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5,
+     a           q2min=0.001, q2lmin=0.001)
+      parameter (ghmax=0.023, ghmin=-0.28)
+
+      real cd(imax),teta(imax,kmax),u(imax,kmax),v(imax,kmax),
+     a     z(imax,kmax),zi(imax,kmax+1)
+      real kh(imax,kmax+1),km(imax,kmax+1),q2(imax,kmax+1),
+     a     long(imax,kmax+1)
+      real unsdz(imax,kmax),unsdzi(imax,kmax+1)
+      real kq(imax,kmax),
+     a     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)
+      real a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax),
+     a     alph(imax,kmax)
+      real ksdz2inf,ksdz2sup
+
+c      save q2,q2l
+
+      save imp 
+
+      data imp/0/
+
+***********************************************************************
+
+      imp=imp+1
+
+************** INCREMENTS VERTICAUX ***********************************
+
+      do 9 i=1,imax
+       zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax))
+ 9    continue
+      do 10 k=1,kmax
+       do 10 i=1,imax
+        unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k))
+ 10   continue
+      do 11 k=2,kmax
+       do 11 i=1,imax
+        unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1))
+ 11   continue
+      do 12 i=1,imax
+       unsdzi(i,1)=0.5/(z(i,1)-zi(i,1))
+       unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax))
+ 12   continue
+
+***********************************************************************
+
+
+************** DIFFUSIVITES KH, KM et KQ ******************************
+* Ci-dessous, une premiere estimation des diffusivites turbulentes km *
+* et kh est effectuee pour utilisation dans les taux de production    *
+* et de destruction de q2 et q2l. On calcule aussi kq.                *
+
+      do 100 k=2,kmax
+       do 100 i=1,imax
+        beta=2.0/(teta(i,k)+teta(i,k-1))
+        n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1))
+        n2(i,k)=amax1(0.0,n2(i,k))
+        du=unsdzi(i,k)*(u(i,k)-u(i,k-1))
+        dv=unsdzi(i,k)*(v(i,k)-v(i,k-1))
+        m2(i,k)=du*du+dv*dv
+        ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10)
+        ri(i,k)=amax1(-0.1,min(4.0,ri(i,k)))
+ 100  continue
+
+      do 110 k=2,kmax
+       do 110 i=1,imax
+        vt=kappa*(zi(i,k)-zi(i,1))
+        long(i,k)=vt/(1.0+vt/160.0)
+        if(n2(i,k).gt.0.0) then
+         long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
+        endif
+        gh=amax1(ghmin,  
+     a     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
+        sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*
+     a           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/
+     a        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
+        km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
+        sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
+        kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
+ 110  continue
+      do 111 i=1,imax
+       us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1)))
+       vt1=(b1**0.666667)*us*us
+       vt2=(b1**0.6666667)*kappa*kappa*
+     a     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))
+c       q2(i,1)=amax1(q2min,vt1,vt2)
+       q2(i,1)=amax1(q2min,vt1)
+       long(i,1)=0.0
+       long(i,kmax+1)=long(i,kmax)
+       sq=0.2
+       kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq)
+ 111  continue
+
+      do 120 k=2,kmax
+       do 120 i=1,imax
+        longc=0.5*(long(i,k)+long(i,k+1))
+        q2c=0.5*(q2(i,k)+q2(i,k+1))
+        sq=0.2
+        kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq)
+ 120  continue
+
+***********************************************************************
+
+
+************** CALCUL DE Q2 *******************************************
+
+      do 200 k=2,kmax
+       do 200 i=1,imax
+        prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k)))
+        dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k))
+     a           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))
+        if(k.lt.kmax) then
+         ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k)
+        else
+         ksdz2sup=0.0
+        endif
+        ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1)
+        b(i,k)=-ksdz2inf*dt
+        a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup)
+        c(i,k)=-ksdz2sup*dt
+        f(i,k)=q2(i,k)+dt*prod
+ 200  continue
+      do 201 i=1,imax
+       f(i,2)=f(i,2)
+     a       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)
+ 201  continue
+      
+      do 210 i=1,imax
+       alph(i,2)=a(i,2)
+ 210  continue
+      do 211 k=3,kmax
+       do 211 i=1,imax
+        bet=b(i,k)/alph(i,k-1)
+        alph(i,k)=a(i,k)-bet*c(i,k-1)
+        f(i,k)=f(i,k)-bet*f(i,k-1)
+ 211  continue
+      
+      do 220 i=1,imax
+       q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax))
+       q2(i,kmax+1)=q2(i,kmax)
+ 220  continue
+      do 221 k=kmax-1,2,-1
+       do 221 i=1,imax
+        q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k))
+ 221  continue
+      do 222 i=1,imax
+       q2(i,2)=amax1(q2(i,2),q2(i,1))
+ 222  continue
+
+***********************************************************************
+
+
+************** EVALUATION FINALE DE KH ET KM **************************
+
+      do 400 k=2,kmax
+       do 400 i=1,imax
+        if(n2(i,k).gt.0.0) then
+         long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))	
+        endif
+        gh=amax1(ghmin,  
+     a         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
+        sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*
+     a           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/
+     a        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
+        km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
+        sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
+        kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
+ 400  continue
+      do 401 i=1,imax
+       km(i,1)=kmmin
+       km(i,kmax+1)=km(i,kmax)
+       kh(i,1)=khmin
+       kh(i,kmax+1)=kh(i,kmax)
+ 401  continue
+ 
+***********************************************************************
+
+c      if(imp.eq.impmax) then
+       am=1.0/float(imax)
+       imp=0
+       do 1000 k=kmax,1,-1
+        au=0.0
+        ateta=0.0
+        az=0.0
+        adz=0.0
+        akq=0.0
+        acd=0.0
+        do 1001 i=1,imax
+         au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k))
+         ateta=ateta+am*teta(i,k)
+         az=az+am*z(i,k)
+         adz=adz+am*(1.0/unsdz(i,k))
+         akq=akq+am*kq(i,k)
+         acd=acd+am*cd(i)
+ 1001   continue
+c        write(*,2000) k,az,adz,au,ateta,akq,acd*1000.0
+ 2000   format(2x,i3,2x,6(f9.2,2x))
+ 1000  continue
+       
+       write(*,*)
+       write(*,*)
+
+       do 1002 k=kmax+1,1,-1
+        azi=0.0
+        adzi=0.0
+        aq2=0.0
+        al=0.0
+        akm=0.0
+        akh=0.0
+        am2=0.0
+        al0=0.0
+        do 1003 i=1,imax
+         azi=azi+am*zi(i,k)
+         adzi=adzi+am*(1.0/unsdzi(i,k))
+         aq2=aq2+am*q2(i,k)
+         al=al+am*long(i,k)
+         akm=akm+am*km(i,k)
+         akh=akh+am*kh(i,k)
+         am2=am2+am*m2(i,k)
+c         al0=al0+am*long0d(i)
+ 1003   continue
+c        write(*,2001) k,azi,aq2,al,akm,akh,am2*1.0e5
+ 2001   format(2x,i3,6(2x,f9.3))
+ 1002  continue
+c      endif
+
+      return
+      end
Index: dynamico_lmdz/simple_physics/phyparam/param/my_25.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/my_25.F	(revision 4201)
+++ 	(revision )
@@ -1,247 +1,0 @@
-      subroutine my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh)
-
-***********************************************************************
-******* FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) *******
-** q2 au interfaces entre mailles.
-***********************************************************************
-
-
-************** DECLARATIONS *******************************************
-
-
-      integer imax,kmax
-
-      parameter (impmax=5)
-
-      real kappa,khmin,kmmin,kqmin,longc
-
-      parameter (kappa=0.4)
-      parameter (a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08,
-     a           e1=1.8, e2=1.33)
-      parameter (khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5,
-     a           q2min=0.001, q2lmin=0.001)
-      parameter (ghmax=0.023, ghmin=-0.28)
-
-      real cd(imax),teta(imax,kmax),u(imax,kmax),v(imax,kmax),
-     a     z(imax,kmax),zi(imax,kmax+1)
-      real kh(imax,kmax+1),km(imax,kmax+1),q2(imax,kmax+1),
-     a     long(imax,kmax+1)
-      real unsdz(imax,kmax),unsdzi(imax,kmax+1)
-      real kq(imax,kmax),
-     a     m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1)
-      real a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax),
-     a     alph(imax,kmax)
-      real ksdz2inf,ksdz2sup
-
-c      save q2,q2l
-
-      save imp 
-
-      data imp/0/
-
-***********************************************************************
-
-      imp=imp+1
-
-************** INCREMENTS VERTICAUX ***********************************
-
-      do 9 i=1,imax
-       zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax))
- 9    continue
-      do 10 k=1,kmax
-       do 10 i=1,imax
-        unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k))
- 10   continue
-      do 11 k=2,kmax
-       do 11 i=1,imax
-        unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1))
- 11   continue
-      do 12 i=1,imax
-       unsdzi(i,1)=0.5/(z(i,1)-zi(i,1))
-       unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax))
- 12   continue
-
-***********************************************************************
-
-
-************** DIFFUSIVITES KH, KM et KQ ******************************
-* Ci-dessous, une premiere estimation des diffusivites turbulentes km *
-* et kh est effectuee pour utilisation dans les taux de production    *
-* et de destruction de q2 et q2l. On calcule aussi kq.                *
-
-      do 100 k=2,kmax
-       do 100 i=1,imax
-        beta=2.0/(teta(i,k)+teta(i,k-1))
-        n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1))
-        n2(i,k)=amax1(0.0,n2(i,k))
-        du=unsdzi(i,k)*(u(i,k)-u(i,k-1))
-        dv=unsdzi(i,k)*(v(i,k)-v(i,k-1))
-        m2(i,k)=du*du+dv*dv
-        ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10)
-        ri(i,k)=amax1(-0.1,min(4.0,ri(i,k)))
- 100  continue
-
-      do 110 k=2,kmax
-       do 110 i=1,imax
-        vt=kappa*(zi(i,k)-zi(i,1))
-        long(i,k)=vt/(1.0+vt/160.0)
-        if(n2(i,k).gt.0.0) then
-         long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))
-        endif
-        gh=amax1(ghmin,  
-     a     min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
-        sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*
-     a           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/
-     a        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
-        km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
-        sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
-        kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
- 110  continue
-      do 111 i=1,imax
-       us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1)))
-       vt1=(b1**0.666667)*us*us
-       vt2=(b1**0.6666667)*kappa*kappa*
-     a     m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1))
-c       q2(i,1)=amax1(q2min,vt1,vt2)
-       q2(i,1)=amax1(q2min,vt1)
-       long(i,1)=0.0
-       long(i,kmax+1)=long(i,kmax)
-       sq=0.2
-       kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq)
- 111  continue
-
-      do 120 k=2,kmax
-       do 120 i=1,imax
-        longc=0.5*(long(i,k)+long(i,k+1))
-        q2c=0.5*(q2(i,k)+q2(i,k+1))
-        sq=0.2
-        kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq)
- 120  continue
-
-***********************************************************************
-
-
-************** CALCUL DE Q2 *******************************************
-
-      do 200 k=2,kmax
-       do 200 i=1,imax
-        prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k)))
-        dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k))
-     a           +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k)))
-        if(k.lt.kmax) then
-         ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k)
-        else
-         ksdz2sup=0.0
-        endif
-        ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1)
-        b(i,k)=-ksdz2inf*dt
-        a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup)
-        c(i,k)=-ksdz2sup*dt
-        f(i,k)=q2(i,k)+dt*prod
- 200  continue
-      do 201 i=1,imax
-       f(i,2)=f(i,2)
-     a       +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1)
- 201  continue
-      
-      do 210 i=1,imax
-       alph(i,2)=a(i,2)
- 210  continue
-      do 211 k=3,kmax
-       do 211 i=1,imax
-        bet=b(i,k)/alph(i,k-1)
-        alph(i,k)=a(i,k)-bet*c(i,k-1)
-        f(i,k)=f(i,k)-bet*f(i,k-1)
- 211  continue
-      
-      do 220 i=1,imax
-       q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax))
-       q2(i,kmax+1)=q2(i,kmax)
- 220  continue
-      do 221 k=kmax-1,2,-1
-       do 221 i=1,imax
-        q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k))
- 221  continue
-      do 222 i=1,imax
-       q2(i,2)=amax1(q2(i,2),q2(i,1))
- 222  continue
-
-***********************************************************************
-
-
-************** EVALUATION FINALE DE KH ET KM **************************
-
-      do 400 k=2,kmax
-       do 400 i=1,imax
-        if(n2(i,k).gt.0.0) then
-         long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k))	
-        endif
-        gh=amax1(ghmin,  
-     a         min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k)))
-        sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh*
-     a           ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/
-     a        ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh))
-        km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm)
-        sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2))
-        kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh)
- 400  continue
-      do 401 i=1,imax
-       km(i,1)=kmmin
-       km(i,kmax+1)=km(i,kmax)
-       kh(i,1)=khmin
-       kh(i,kmax+1)=kh(i,kmax)
- 401  continue
- 
-***********************************************************************
-
-c      if(imp.eq.impmax) then
-       am=1.0/float(imax)
-       imp=0
-       do 1000 k=kmax,1,-1
-        au=0.0
-        ateta=0.0
-        az=0.0
-        adz=0.0
-        akq=0.0
-        acd=0.0
-        do 1001 i=1,imax
-         au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k))
-         ateta=ateta+am*teta(i,k)
-         az=az+am*z(i,k)
-         adz=adz+am*(1.0/unsdz(i,k))
-         akq=akq+am*kq(i,k)
-         acd=acd+am*cd(i)
- 1001   continue
-c        write(*,2000) k,az,adz,au,ateta,akq,acd*1000.0
- 2000   format(2x,i3,2x,6(f9.2,2x))
- 1000  continue
-       
-       write(*,*)
-       write(*,*)
-
-       do 1002 k=kmax+1,1,-1
-        azi=0.0
-        adzi=0.0
-        aq2=0.0
-        al=0.0
-        akm=0.0
-        akh=0.0
-        am2=0.0
-        al0=0.0
-        do 1003 i=1,imax
-         azi=azi+am*zi(i,k)
-         adzi=adzi+am*(1.0/unsdzi(i,k))
-         aq2=aq2+am*q2(i,k)
-         al=al+am*long(i,k)
-         akm=akm+am*km(i,k)
-         akh=akh+am*kh(i,k)
-         am2=am2+am*m2(i,k)
-c         al0=al0+am*long0d(i)
- 1003   continue
-c        write(*,2001) k,azi,aq2,al,akm,akh,am2*1.0e5
- 2001   format(2x,i3,6(2x,f9.3))
- 1002  continue
-c      endif
-
-      return
-      end
Index: dynamico_lmdz/simple_physics/phyparam/param/paramdef.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/paramdef.F	(revision 4201)
+++ dynamico_lmdz/simple_physics/phyparam/param/paramdef.F	(revision 4202)
@@ -1,5 +1,5 @@
       SUBROUTINE paramdef(ngrid,rnatur,albedo,inertie,emissiv,z0)
       USE comgeomfi, ONLY : lati,sinlat,coslat
-      USE vdif_mod, ONLY : lmixmin, emin_turb
+      USE turbulence, ONLY : lmixmin, emin_turb
       USE planet, ONLY : coefir, coefvis
       USE astronomy
Index: dynamico_lmdz/simple_physics/phyparam/param/phyparam.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/phyparam.F	(revision 4201)
+++ dynamico_lmdz/simple_physics/phyparam/param/phyparam.F	(revision 4202)
@@ -13,5 +13,5 @@
       USE planet
       USE astronomy
-      USE vdif_mod, ONLY : vdif
+      USE turbulence, ONLY : vdif
       USE solar, ONLY : solang, zenang, mucorr
       USE radiative_sw, ONLY : sw
Index: dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90	(revision 4202)
+++ dynamico_lmdz/simple_physics/phyparam/physics/turbulence.F90	(revision 4202)
@@ -0,0 +1,446 @@
+MODULE turbulence
+  
+#include "use_logging.h"
+
+  IMPLICIT NONE
+  SAVE
+  PRIVATE
+
+  REAL, PARAMETER :: karman=0.4
+  REAL :: lmixmin=100., emin_turb=1e-8  
+
+  PUBLIC :: vdif, lmixmin, emin_turb
+  
+CONTAINS
+  
+  PURE SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
+    !=======================================================================
+    !
+    !   Subject: computation of the surface drag coefficient using the
+    !   -------  approch developed by Loui for ECMWF.
+    !
+    !   Author: Frederic Hourdin  15 /10 /93
+    !   -------
+    !
+    !   Arguments:
+    !   ----------
+    !
+    !   inputs:
+    !   ------
+    !     ngrid            size of the horizontal grid
+    !     pz0(ngrid)       roughness length ?
+    !     pg               gravity (m s -2)
+    !     pz(ngrid)        height of the first atmospheric layer
+    !     pu(ngrid)        u component of the wind in that layer
+    !     pv(ngrid)        v component of the wind in that layer
+    !     pts(ngrid)       surfacte temperature
+    !     ph(ngrid)        potential temperature T*(p/ps)^kappa
+    !
+    !   outputs:
+    !   --------
+    !     pcdv(ngrid)      Cd for the wind
+    !     pcdh(ngrid)      Cd for potential temperature
+    !
+    !=======================================================================
+    !
+    !-----------------------------------------------------------------------
+    !   Declarations:
+    !   -------------
+    
+    !   Arguments:
+    !   ----------
+    
+    INTEGER, INTENT(IN) ::  ngrid
+    REAL, INTENT(IN)  ::  pg, pz0(ngrid), pz(ngrid), &
+         pu(ngrid),pv(ngrid), pts(ngrid),ph(ngrid)
+    REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid)
+    
+    !   Local:
+    !   ------
+    
+    REAL, PARAMETER :: b=5., c=5., d=5., umin2=1e-12, &
+         c2b=2.*b, c3bc=3.*b*c, c3b=3.*b
+    INTEGER ig
+    REAL zu2,z1,zri,zcd0,zz    
+    
+    !-----------------------------------------------------------------------
+    !   couche de surface:
+    !   ------------------
+    
+    !     DO ig=1,ngrid
+    !        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
+    !        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
+    !        pcdh(ig)=pcdv(ig)
+    !     ENDDO
+    !     RETURN
+    
+    
+!!!! WARNING, verifier la formule originale de Louis!
+    DO ig=1,ngrid
+       zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
+       zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
+       z1=1.+pz(ig)/pz0(ig)
+       zcd0=karman/log(z1)
+       zcd0=zcd0*zcd0*sqrt(zu2)
+       IF(zri.LT.0.) THEN
+          z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
+          pcdv(ig)=zcd0*(1.-2.*z1)
+          pcdh(ig)=zcd0*(1.-3.*z1)
+       ELSE
+          zz=sqrt(1.+d*zri)
+          pcdv(ig)=zcd0/(1.+c2b*zri/zz)
+          pcdh(ig)=zcd0/(1.+c3b*zri*zz)
+       ENDIF
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    
+  END SUBROUTINE vdif_cd
+  
+  PURE SUBROUTINE vdif_k(ngrid,nlay,   &
+       ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pkv,pkh)
+    ! FIXME : pkh := pkv
+    USE planet
+    INTEGER, INTENT(IN) :: ngrid,nlay
+    REAL, INTENT(IN) ::  ptimestep, pg, &
+         pzlay(ngrid,nlay),pzlev(ngrid,nlay+1), pz0(ngrid), &
+         pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+    REAL, INTENT(OUT) :: pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
+    
+    INTEGER ig,il
+    REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
+    
+    DO ig=1,ngrid
+       pkv(ig,1)=0.
+       pkh(ig,1)=0.
+       pkv(ig,nlay+1)=0.
+       pkh(ig,nlay+1)=0.
+    ENDDO
+    
+    DO il=2,nlay
+       DO ig=1,ngrid
+          z1=pzlev(ig,il)+pz0(ig)
+          lmix=karman*z1/(1.+karman*z1/lmixmin)
+          !           lmix=lmixmin
+          ! WARNING test lmix=lmixmin
+          zdu=pu(ig,il)-pu(ig,il-1)
+          zdv=pv(ig,il)-pv(ig,il-1)
+          zdz=pzlay(ig,il)-pzlay(ig,il-1)
+          zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)
+          IF(zdvodz2.LT.1.e-5) THEN
+             pkv(ig,il)=lmix*sqrt(emin_turb)
+          ELSE
+             zri=2.*pg*(ph(ig,il)-ph(ig,il-1)) &
+                  / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
+             pkv(ig,il)= &
+                  lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
+          ENDIF
+          pkh(ig,il)=pkv(ig,il)
+          !           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
+          !    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
+          !    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
+          !    s      ph(ig,il),ph(ig,il-1)
+       ENDDO
+    ENDDO
+    
+  END SUBROUTINE vdif_k
+
+  
+  SUBROUTINE multipl(n,x1,x2,y)
+!=======================================================================
+!   multiplication de deux vecteurs
+!=======================================================================
+    INTEGER n,i
+    REAL x1(n),x2(n),y(n)
+
+    DO i=1,n
+       y(i)=x1(i)*x2(i)
+    END DO
+  END SUBROUTINE multipl
+
+  SUBROUTINE vdif(ngrid,nlay,ptime, &
+       ptimestep,pcapcal,pz0, &
+       pplay,pplev,pzlay,pzlev, &
+       pu,pv,ph,ptsrf,pemis, &
+       pdufi,pdvfi,pdhfi,pfluxsrf, &
+       pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l, &
+       lwrite)
+    USE phys_const
+    
+!=======================================================================
+!
+!   Diffusion verticale
+!   Shema implicite
+!   On commence par rajouter au variables x la tendance physique
+!   et on resoult en fait:
+!      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
+!
+!   arguments:
+!   ----------
+!
+!   entree:
+!   -------
+!
+!
+!=======================================================================
+
+!
+!   arguments:
+!   ----------
+
+    INTEGER ngrid,nlay
+    REAL ptime,ptimestep
+    REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+    REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
+    REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
+    REAL ptsrf(ngrid),pemis(ngrid)
+    REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
+    REAL pfluxsrf(ngrid)
+    REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
+    REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
+    REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
+    LOGICAL lwrite
+    !
+    !   local:
+    !   ------
+    
+    INTEGER ilev,ig,ilay,nlev
+    INTEGER unit,ierr,it1,it2
+    INTEGER cluvdb,putdat,putvdim,setname,setvdim
+    REAL z4st,zdplanck(ngrid),zu2
+    REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
+    REAL zcdv(ngrid),zcdh(ngrid)
+    REAL zu(ngrid,nlay),zv(ngrid,nlay)
+    REAL zh(ngrid,nlay)
+    REAL ztsrf2(ngrid)
+    REAL z1(ngrid),z2(ngrid)
+    REAL za(ngrid,nlay),zb(ngrid,nlay)
+    REAL zb0(ngrid,nlay)
+    REAL zc(ngrid,nlay),zd(ngrid,nlay)
+    REAL zcst1
+    
+    !
+    !-----------------------------------------------------------------------
+    !   initialisations:
+    !   ----------------
+    
+    nlev=nlay+1
+    
+    !   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
+    !   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
+    !   ---------------------------------
+    
+    DO ilay=1,nlay
+       DO ig=1,ngrid
+          za(ig,ilay) = (pplev(ig,ilay)-pplev(ig,ilay+1))/g
+       ENDDO
+    ENDDO
+    
+    zcst1=4.*g*ptimestep/(r*r)
+    DO ilev=2,nlev-1
+       DO ig=1,ngrid
+          zb0(ig,ilev)=pplev(ig,ilev) &
+               *(pplev(ig,1)/pplev(ig,ilev))**rcp &
+               /(ph(ig,ilev-1)+ph(ig,ilev))
+          zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev) &
+               / (pplay(ig,ilev-1)-pplay(ig,ilev))
+       ENDDO
+    ENDDO
+    DO ig=1,ngrid
+       zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
+    ENDDO
+    IF(lwrite) THEN
+       ig=ngrid/2+1
+       WRITELOG(*,*) 'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
+       DO ilay=1,nlay
+          WRITELOG(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay), &
+               pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
+       ENDDO
+       WRITELOG(*,*) 'Pression (mbar) ,altitude (km),zb'
+       DO ilev=1,nlay
+          WRITELOG(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev), &
+               zb0(ig,ilev)
+       ENDDO
+       LOG_DBG('vdif')
+    ENDIF
+    
+    !-----------------------------------------------------------------------
+    !   2. ajout des tendances physiques:
+    !   ------------------------------
+    
+    DO ilev=1,nlay
+       DO ig=1,ngrid
+          zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
+          zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
+          zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   3. calcul de  cd :
+    !   ----------------
+    !
+    CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
+    CALL vdif_k(ngrid,nlay,ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zkv,zkh)
+    
+    DO ig=1,ngrid
+       zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
+       zcdv(ig)=zcdv(ig)*sqrt(zu2)
+       zcdh(ig)=zcdh(ig)*sqrt(zu2)
+    ENDDO
+    
+    IF(lwrite) THEN
+       
+       WRITELOG(*,*) 'Diagnostique diffusion verticale'
+       WRITELOG(*,*) 'LMIXMIN',lmixmin
+       PRINT*,'coefficients Cd pour v et h'
+       PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
+       PRINT*,'coefficients K pour v et h'
+       DO ilev=1,nlay
+          PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
+       ENDDO
+       LOG_DBG('vdif')
+    ENDIF
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale pour u:
+    !   -----------------------------
+    
+    CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
+    CALL multipl(ngrid,zcdv,zb0,zb)
+    DO ig=1,ngrid
+       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+       zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
+       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+    ENDDO
+    
+    DO ilay=nlay-1,1,-1
+       DO ig=1,ngrid
+          z1(ig) = 1./(za(ig,ilay)+zb(ig,ilay) &
+               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+          zc(ig,ilay) = (za(ig,ilay)*zu(ig,ilay) &
+               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+          zd(ig,ilay) = zb(ig,ilay)*z1(ig)
+       ENDDO
+    ENDDO
+    
+    DO ig=1,ngrid
+       zu(ig,1)=zc(ig,1)
+    ENDDO
+    DO ilay=2,nlay
+       DO ig=1,ngrid
+          zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale pour v:
+    !   -----------------------------
+    !
+    DO ig=1,ngrid
+       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+       zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
+       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+    ENDDO
+    
+    DO ilay=nlay-1,1,-1
+       DO ig=1,ngrid
+          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
+               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+          zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay) &
+               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+       ENDDO
+    ENDDO
+    
+    DO ig=1,ngrid
+       zv(ig,1)=zc(ig,1)
+    ENDDO
+    DO ilay=2,nlay
+       DO ig=1,ngrid
+          zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale pour h:
+    !   -----------------------------
+    !
+    CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
+    CALL multipl(ngrid,zcdh,zb0,zb)
+    DO ig=1,ngrid
+       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
+       zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
+       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
+    ENDDO
+    
+    DO ilay=nlay-1,1,-1
+       DO ig=1,ngrid
+          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
+               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
+          zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay) &
+               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
+          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   rajout eventuel de planck dans le shema implicite:
+    !   --------------------------------------------------
+    
+    z4st=4.*5.67e-8*ptimestep
+    DO ig=1,ngrid
+       zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   calcul le l'evolution de la temperature du sol:
+    !   -----------------------------------------------
+    
+    DO ig=1,ngrid
+       z1(ig) = pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) &
+            + zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
+       z2(ig) = pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
+       ztsrf2(ig) = z1(ig)/z2(ig)
+       zh(ig,1) = zc(ig,1)+zd(ig,1)*ztsrf2(ig)
+       pdtsrf(ig) = (ztsrf2(ig)-ptsrf(ig))/ptimestep
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   integration verticale finale:
+    !   -----------------------------
+    
+    DO ilay=2,nlay
+       DO ig=1,ngrid
+          zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
+       ENDDO
+    ENDDO
+    
+    !-----------------------------------------------------------------------
+    !   calcul final des tendances de la diffusion verticale:
+    !   -----------------------------------------------------
+    
+    DO ilev = 1, nlay
+       DO ig=1,ngrid
+          pdudif(ig,ilev) = ( zu(ig,ilev) &
+               - (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
+          pdvdif(ig,ilev) = ( zv(ig,ilev) &
+               - (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
+          pdhdif(ig,ilev) = ( zh(ig,ilev) &
+               - (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep
+       ENDDO
+    ENDDO
+    
+    IF(lwrite) THEN
+       PRINT*
+       PRINT*,'Diagnostique de la diffusion verticale'
+       PRINT*,'h avant et apres diffusion verticale'
+       PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
+       DO  ilev=1,nlay
+          PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
+       ENDDO
+    END IF
+
+  END SUBROUTINE vdif
+  
+END MODULE turbulence
Index: dynamico_lmdz/simple_physics/phyparam/physics/vdif_mod.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/vdif_mod.F90	(revision 4201)
+++ 	(revision )
@@ -1,446 +1,0 @@
-MODULE vdif_mod
-  
-#include "use_logging.h"
-
-  IMPLICIT NONE
-  SAVE
-  PRIVATE
-
-  REAL, PARAMETER :: karman=0.4
-  REAL :: lmixmin=100., emin_turb=1e-8  
-
-  PUBLIC :: vdif, lmixmin, emin_turb
-  
-CONTAINS
-  
-  PURE SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
-    !=======================================================================
-    !
-    !   Subject: computation of the surface drag coefficient using the
-    !   -------  approch developed by Loui for ECMWF.
-    !
-    !   Author: Frederic Hourdin  15 /10 /93
-    !   -------
-    !
-    !   Arguments:
-    !   ----------
-    !
-    !   inputs:
-    !   ------
-    !     ngrid            size of the horizontal grid
-    !     pz0(ngrid)       roughness length ?
-    !     pg               gravity (m s -2)
-    !     pz(ngrid)        height of the first atmospheric layer
-    !     pu(ngrid)        u component of the wind in that layer
-    !     pv(ngrid)        v component of the wind in that layer
-    !     pts(ngrid)       surfacte temperature
-    !     ph(ngrid)        potential temperature T*(p/ps)^kappa
-    !
-    !   outputs:
-    !   --------
-    !     pcdv(ngrid)      Cd for the wind
-    !     pcdh(ngrid)      Cd for potential temperature
-    !
-    !=======================================================================
-    !
-    !-----------------------------------------------------------------------
-    !   Declarations:
-    !   -------------
-    
-    !   Arguments:
-    !   ----------
-    
-    INTEGER, INTENT(IN) ::  ngrid
-    REAL, INTENT(IN)  ::  pg, pz0(ngrid), pz(ngrid), &
-         pu(ngrid),pv(ngrid), pts(ngrid),ph(ngrid)
-    REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid)
-    
-    !   Local:
-    !   ------
-    
-    REAL, PARAMETER :: b=5., c=5., d=5., umin2=1e-12, &
-         c2b=2.*b, c3bc=3.*b*c, c3b=3.*b
-    INTEGER ig
-    REAL zu2,z1,zri,zcd0,zz    
-    
-    !-----------------------------------------------------------------------
-    !   couche de surface:
-    !   ------------------
-    
-    !     DO ig=1,ngrid
-    !        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
-    !        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
-    !        pcdh(ig)=pcdv(ig)
-    !     ENDDO
-    !     RETURN
-    
-    
-!!!! WARNING, verifier la formule originale de Louis!
-    DO ig=1,ngrid
-       zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
-       zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
-       z1=1.+pz(ig)/pz0(ig)
-       zcd0=karman/log(z1)
-       zcd0=zcd0*zcd0*sqrt(zu2)
-       IF(zri.LT.0.) THEN
-          z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
-          pcdv(ig)=zcd0*(1.-2.*z1)
-          pcdh(ig)=zcd0*(1.-3.*z1)
-       ELSE
-          zz=sqrt(1.+d*zri)
-          pcdv(ig)=zcd0/(1.+c2b*zri/zz)
-          pcdh(ig)=zcd0/(1.+c3b*zri*zz)
-       ENDIF
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    
-  END SUBROUTINE vdif_cd
-  
-  PURE SUBROUTINE vdif_k(ngrid,nlay,   &
-       ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pkv,pkh)
-    ! FIXME : pkh := pkv
-    USE planet
-    INTEGER, INTENT(IN) :: ngrid,nlay
-    REAL, INTENT(IN) ::  ptimestep, pg, &
-         pzlay(ngrid,nlay),pzlev(ngrid,nlay+1), pz0(ngrid), &
-         pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
-    REAL, INTENT(OUT) :: pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
-    
-    INTEGER ig,il
-    REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
-    
-    DO ig=1,ngrid
-       pkv(ig,1)=0.
-       pkh(ig,1)=0.
-       pkv(ig,nlay+1)=0.
-       pkh(ig,nlay+1)=0.
-    ENDDO
-    
-    DO il=2,nlay
-       DO ig=1,ngrid
-          z1=pzlev(ig,il)+pz0(ig)
-          lmix=karman*z1/(1.+karman*z1/lmixmin)
-          !           lmix=lmixmin
-          ! WARNING test lmix=lmixmin
-          zdu=pu(ig,il)-pu(ig,il-1)
-          zdv=pv(ig,il)-pv(ig,il-1)
-          zdz=pzlay(ig,il)-pzlay(ig,il-1)
-          zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)
-          IF(zdvodz2.LT.1.e-5) THEN
-             pkv(ig,il)=lmix*sqrt(emin_turb)
-          ELSE
-             zri=2.*pg*(ph(ig,il)-ph(ig,il-1)) &
-                  / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
-             pkv(ig,il)= &
-                  lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
-          ENDIF
-          pkh(ig,il)=pkv(ig,il)
-          !           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
-          !    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
-          !    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
-          !    s      ph(ig,il),ph(ig,il-1)
-       ENDDO
-    ENDDO
-    
-  END SUBROUTINE vdif_k
-
-  
-  SUBROUTINE multipl(n,x1,x2,y)
-!=======================================================================
-!   multiplication de deux vecteurs
-!=======================================================================
-    INTEGER n,i
-    REAL x1(n),x2(n),y(n)
-
-    DO i=1,n
-       y(i)=x1(i)*x2(i)
-    END DO
-  END SUBROUTINE multipl
-
-  SUBROUTINE vdif(ngrid,nlay,ptime, &
-       ptimestep,pcapcal,pz0, &
-       pplay,pplev,pzlay,pzlev, &
-       pu,pv,ph,ptsrf,pemis, &
-       pdufi,pdvfi,pdhfi,pfluxsrf, &
-       pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l, &
-       lwrite)
-    USE phys_const
-    
-!=======================================================================
-!
-!   Diffusion verticale
-!   Shema implicite
-!   On commence par rajouter au variables x la tendance physique
-!   et on resoult en fait:
-!      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
-!
-!   arguments:
-!   ----------
-!
-!   entree:
-!   -------
-!
-!
-!=======================================================================
-
-!
-!   arguments:
-!   ----------
-
-    INTEGER ngrid,nlay
-    REAL ptime,ptimestep
-    REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
-    REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
-    REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
-    REAL ptsrf(ngrid),pemis(ngrid)
-    REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
-    REAL pfluxsrf(ngrid)
-    REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
-    REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
-    REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
-    LOGICAL lwrite
-    !
-    !   local:
-    !   ------
-    
-    INTEGER ilev,ig,ilay,nlev
-    INTEGER unit,ierr,it1,it2
-    INTEGER cluvdb,putdat,putvdim,setname,setvdim
-    REAL z4st,zdplanck(ngrid),zu2
-    REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
-    REAL zcdv(ngrid),zcdh(ngrid)
-    REAL zu(ngrid,nlay),zv(ngrid,nlay)
-    REAL zh(ngrid,nlay)
-    REAL ztsrf2(ngrid)
-    REAL z1(ngrid),z2(ngrid)
-    REAL za(ngrid,nlay),zb(ngrid,nlay)
-    REAL zb0(ngrid,nlay)
-    REAL zc(ngrid,nlay),zd(ngrid,nlay)
-    REAL zcst1
-    
-    !
-    !-----------------------------------------------------------------------
-    !   initialisations:
-    !   ----------------
-    
-    nlev=nlay+1
-    
-    !   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
-    !   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
-    !   ---------------------------------
-    
-    DO ilay=1,nlay
-       DO ig=1,ngrid
-          za(ig,ilay) = (pplev(ig,ilay)-pplev(ig,ilay+1))/g
-       ENDDO
-    ENDDO
-    
-    zcst1=4.*g*ptimestep/(r*r)
-    DO ilev=2,nlev-1
-       DO ig=1,ngrid
-          zb0(ig,ilev)=pplev(ig,ilev) &
-               *(pplev(ig,1)/pplev(ig,ilev))**rcp &
-               /(ph(ig,ilev-1)+ph(ig,ilev))
-          zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev) &
-               / (pplay(ig,ilev-1)-pplay(ig,ilev))
-       ENDDO
-    ENDDO
-    DO ig=1,ngrid
-       zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
-    ENDDO
-    IF(lwrite) THEN
-       ig=ngrid/2+1
-       WRITELOG(*,*) 'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
-       DO ilay=1,nlay
-          WRITELOG(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay), &
-               pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
-       ENDDO
-       WRITELOG(*,*) 'Pression (mbar) ,altitude (km),zb'
-       DO ilev=1,nlay
-          WRITELOG(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev), &
-               zb0(ig,ilev)
-       ENDDO
-       LOG_DBG('vdif')
-    ENDIF
-    
-    !-----------------------------------------------------------------------
-    !   2. ajout des tendances physiques:
-    !   ------------------------------
-    
-    DO ilev=1,nlay
-       DO ig=1,ngrid
-          zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
-          zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
-          zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
-       ENDDO
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   3. calcul de  cd :
-    !   ----------------
-    !
-    CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
-    CALL vdif_k(ngrid,nlay,ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zkv,zkh)
-    
-    DO ig=1,ngrid
-       zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
-       zcdv(ig)=zcdv(ig)*sqrt(zu2)
-       zcdh(ig)=zcdh(ig)*sqrt(zu2)
-    ENDDO
-    
-    IF(lwrite) THEN
-       
-       WRITELOG(*,*) 'Diagnostique diffusion verticale'
-       WRITELOG(*,*) 'LMIXMIN',lmixmin
-       PRINT*,'coefficients Cd pour v et h'
-       PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
-       PRINT*,'coefficients K pour v et h'
-       DO ilev=1,nlay
-          PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
-       ENDDO
-       LOG_DBG('vdif')
-    ENDIF
-    
-    !-----------------------------------------------------------------------
-    !   integration verticale pour u:
-    !   -----------------------------
-    
-    CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
-    CALL multipl(ngrid,zcdv,zb0,zb)
-    DO ig=1,ngrid
-       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-       zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
-       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-    ENDDO
-    
-    DO ilay=nlay-1,1,-1
-       DO ig=1,ngrid
-          z1(ig) = 1./(za(ig,ilay)+zb(ig,ilay) &
-               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-          zc(ig,ilay) = (za(ig,ilay)*zu(ig,ilay) &
-               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-          zd(ig,ilay) = zb(ig,ilay)*z1(ig)
-       ENDDO
-    ENDDO
-    
-    DO ig=1,ngrid
-       zu(ig,1)=zc(ig,1)
-    ENDDO
-    DO ilay=2,nlay
-       DO ig=1,ngrid
-          zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
-       ENDDO
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   integration verticale pour v:
-    !   -----------------------------
-    !
-    DO ig=1,ngrid
-       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-       zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
-       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-    ENDDO
-    
-    DO ilay=nlay-1,1,-1
-       DO ig=1,ngrid
-          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
-               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-          zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay) &
-               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
-       ENDDO
-    ENDDO
-    
-    DO ig=1,ngrid
-       zv(ig,1)=zc(ig,1)
-    ENDDO
-    DO ilay=2,nlay
-       DO ig=1,ngrid
-          zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
-       ENDDO
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   integration verticale pour h:
-    !   -----------------------------
-    !
-    CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
-    CALL multipl(ngrid,zcdh,zb0,zb)
-    DO ig=1,ngrid
-       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
-       zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
-       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
-    ENDDO
-    
-    DO ilay=nlay-1,1,-1
-       DO ig=1,ngrid
-          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
-               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
-          zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay) &
-               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
-          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
-       ENDDO
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   rajout eventuel de planck dans le shema implicite:
-    !   --------------------------------------------------
-    
-    z4st=4.*5.67e-8*ptimestep
-    DO ig=1,ngrid
-       zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   calcul le l'evolution de la temperature du sol:
-    !   -----------------------------------------------
-    
-    DO ig=1,ngrid
-       z1(ig) = pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) &
-            + zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
-       z2(ig) = pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
-       ztsrf2(ig) = z1(ig)/z2(ig)
-       zh(ig,1) = zc(ig,1)+zd(ig,1)*ztsrf2(ig)
-       pdtsrf(ig) = (ztsrf2(ig)-ptsrf(ig))/ptimestep
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   integration verticale finale:
-    !   -----------------------------
-    
-    DO ilay=2,nlay
-       DO ig=1,ngrid
-          zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
-       ENDDO
-    ENDDO
-    
-    !-----------------------------------------------------------------------
-    !   calcul final des tendances de la diffusion verticale:
-    !   -----------------------------------------------------
-    
-    DO ilev = 1, nlay
-       DO ig=1,ngrid
-          pdudif(ig,ilev) = ( zu(ig,ilev) &
-               - (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
-          pdvdif(ig,ilev) = ( zv(ig,ilev) &
-               - (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
-          pdhdif(ig,ilev) = ( zh(ig,ilev) &
-               - (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep
-       ENDDO
-    ENDDO
-    
-    IF(lwrite) THEN
-       PRINT*
-       PRINT*,'Diagnostique de la diffusion verticale'
-       PRINT*,'h avant et apres diffusion verticale'
-       PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
-       DO  ilev=1,nlay
-          PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
-       ENDDO
-    END IF
-
-  END SUBROUTINE vdif
-  
-END MODULE vdif_mod
