Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 332)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 336)
@@ -1176,8 +1176,177 @@
 c Calculer le frottement au sol (Cdrag)
 c
+      CALL clcdrag(knon, nsrf, is_oce, zxli, u, v, t, q, zgeop,
+     .             ts, qsurf, rugos, pcfm, pcfh, zcdn, zri)
+c
+c Calculer les coefficients turbulents dans l'atmosphere
+c
+      DO i = 1, knon
+         itop(i) = isommet
+      ENDDO
+
+      DO k = 2, isommet
+      DO i = 1, knon
+            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
+     .                     +(v(i,k)-v(i,k-1))**2)
+            zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
+            zdphi =zmgeom(i) / 2.0
+            zt = (t(i,k)+t(i,k-1)) * 0.5
+            zq = (q(i,k)+q(i,k-1)) * 0.5
+c
+c           calculer Qs et dQs/dT:
+c
+            IF (thermcep) THEN
+              zdelta = MAX(0.,SIGN(1.,RTT-zt))
+              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) 
+     .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
+              zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
+              zqs = MIN(0.5,zqs)
+              zcor = 1./(1.-RETV*zqs)
+              zqs = zqs*zcor
+              zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
+            ELSE
+              IF (zt .LT. t_coup) THEN
+                 zqs = qsats(zt) / pplay(i,k)
+                 zdqs = dqsats(zt,zqs)
+              ELSE
+                 zqs = qsatl(zt) / pplay(i,k)
+                 zdqs = dqsatl(zt,zqs)
+              ENDIF
+            ENDIF
+c
+c           calculer la fraction nuageuse (processus humide):
+c
+            zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
+            zfr = MAX(0.0,MIN(1.0,zfr))
+            IF (.NOT.richum) zfr = 0.0
+c
+c           calculer le nombre de Richardson:
+c
+            IF (tvirtu) THEN
+            ztvd =( t(i,k)
+     .             + zdphi/RCPD/(1.+RVTMP2*zq)
+     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
+     .            )*(1.+RETV*q(i,k))
+            ztvu =( t(i,k-1)
+     .             - zdphi/RCPD/(1.+RVTMP2*zq)
+     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
+     .            )*(1.+RETV*q(i,k-1))
+            zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
+            zri(i) = zri(i)
+     .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)
+     .               *(paprs(i,k)/101325.0)**RKAPPA
+     .               /(zdu2*0.5*(ztvd+ztvu))
+c
+            ELSE ! calcul de Ridchardson compatible LMD5
+c
+            zri(i) =(RCPD*(t(i,k)-t(i,k-1))
+     .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)
+     .               *(pplay(i,k)-pplay(i,k-1))
+     .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
+            zri(i) = zri(i) +
+     .             zmgeom(i)*zmgeom(i)*gamt(k)/RG
+cSB     .             /(paprs(i,k)/101325.0)**RKAPPA
+     .             *(paprs(i,k)/101325.0)**RKAPPA
+     .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))
+            ENDIF
+c
+c           finalement, les coefficients d'echange sont obtenus:
+c
+            zcdn=SQRT(zdu2) / zmgeom(i) * RG
+c
+          IF (opt_ec) THEN
+            z2geomf=zgeop(i,k-1)+zgeop(i,k)
+            zalm2=(0.5*ckap/RG*z2geomf
+     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
+            zalh2=(0.5*ckap/rg*z2geomf
+     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
+            IF (zri(i).LT.0.0) THEN  ! situation instable
+               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
+     .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
+               zscf = SQRT(-zri(i)*zscf)
+               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
+               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
+               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
+               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
+            ELSE ! situation stable
+               zscf=SQRT(1.+cd*zri(i))
+               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
+               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
+            ENDIF
+          ELSE
+            zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
+     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
+            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
+            pcfm(i,k)= zl2(i)* pcfm(i,k)
+            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
+          ENDIF
+      ENDDO
+      ENDDO
+c
+c Au-dela du sommet, pas de diffusion turbulente:
+c
+      DO i = 1, knon
+         IF (itop(i)+1 .LE. klev) THEN
+            DO k = itop(i)+1, klev
+               pcfh(i,k) = 0.0
+               pcfm(i,k) = 0.0
+            ENDDO
+         ENDIF
+      ENDDO
+c
+      RETURN
+      END
+
+      SUBROUTINE clcdrag(knon, nsrf, is_oce, zxli, 
+     .                   u, v, t, q, zgeop, 
+     .                   ts, qsurf, rugos,
+     .                   pcfm, pcfh, zcdn, zri)
+c ================================================================= c
+c Objet : calcul cdrags pour le moment et les flux chaleur sensible, latente (pcfm,pcfh)
+c         et du nombre de Richardson zri
+c ================================================================= c
+      IMPLICIT NONE
+#include "dimensions.h"
+#include "dimphy.h"
+#include "YOMCST.h"
+#include "indicesol.h"
+c
+      INTEGER knon, nsrf, is_oce
+      REAL ts(klon), qsurf(klon)
+      REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
+      REAL zgeop(klon,klev)
+      REAL rugos(klon), zri(klon)
+c
+      REAL zcdn
+      REAL pcfm(klon,klev), pcfh(klon,klev)
+c
+c Quelques constantes et options:
+c
+      REAL ckap, cb, cc, cd, cepdu2
+      PARAMETER (ckap=0.35)
+      PARAMETER (cb=5.0)
+      PARAMETER (cc=5.0)
+      PARAMETER (cd=5.0)
+      PARAMETER (cepdu2 =(0.1)**2)
+c
+c Variables locales
+      INTEGER i
+      REAL zdu2, zdphi, ztsolv, ztvd, zscf, zucf, zcr
+      REAL friv, frih
+      REAL zcfm1(klon), zcfm2(klon)
+      REAL zcfh1(klon), zcfh2(klon)
+c
+c Fonctions thermodynamiques et fonctions d'instabilite
+      REAL fsta, fins, x
+      LOGICAL zxli
+      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
+      fins(x) = SQRT(1.0-18.0*x)
+c
+c Calculer le frottement au sol (Cdrag)
+c
       DO i = 1, knon
          zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
          zdphi=zgeop(i,1)
-c         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
+c        ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
          ztsolv = ts(i) * (1.0+RETV*qsurf(i))
          ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
@@ -1215,125 +1384,8 @@
            IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
          ENDIF
-      ENDDO
-
-c
-c Calculer les coefficients turbulents dans l'atmosphere
-c
-      DO i = 1, knon
-         itop(i) = isommet
-      ENDDO
-
-      DO k = 2, isommet
-      DO i = 1, knon
-            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
-     .                     +(v(i,k)-v(i,k-1))**2)
-            zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
-            zdphi =zmgeom(i) / 2.0
-            zt = (t(i,k)+t(i,k-1)) * 0.5
-            zq = (q(i,k)+q(i,k-1)) * 0.5
-c
-c           calculer Qs et dQs/dT:
-c
-            IF (thermcep) THEN
-              zdelta = MAX(0.,SIGN(1.,RTT-zt))
-              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) 
-     .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
-              zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
-              zqs = MIN(0.5,zqs)
-              zcor = 1./(1.-RETV*zqs)
-              zqs = zqs*zcor
-              zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
-            ELSE
-              IF (zt .LT. t_coup) THEN
-                 zqs = qsats(zt) / pplay(i,k)
-                 zdqs = dqsats(zt,zqs)
-              ELSE
-                 zqs = qsatl(zt) / pplay(i,k)
-                 zdqs = dqsatl(zt,zqs)
-              ENDIF
-            ENDIF
-c
-c           calculer la fraction nuageuse (processus humide):
-c
-            zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
-            zfr = MAX(0.0,MIN(1.0,zfr))
-            IF (.NOT.richum) zfr = 0.0
-c
-c           calculer le nombre de Richardson:
-c
-            IF (tvirtu) THEN
-            ztvd =( t(i,k)
-     .             + zdphi/RCPD/(1.+RVTMP2*zq)
-     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
-     .            )*(1.+RETV*q(i,k))
-            ztvu =( t(i,k-1)
-     .             - zdphi/RCPD/(1.+RVTMP2*zq)
-     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
-     .            )*(1.+RETV*q(i,k-1))
-            zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
-            zri(i) = zri(i)
-     .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)
-     .               *(paprs(i,k)/101325.0)**RKAPPA
-     .               /(zdu2*0.5*(ztvd+ztvu))
-c
-            ELSE ! calcul de Ridchardson compatible LMD5
-c
-            zri(i) =(RCPD*(t(i,k)-t(i,k-1))
-     .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)
-     .               *(pplay(i,k)-pplay(i,k-1))
-     .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
-            zri(i) = zri(i) +
-     .             zmgeom(i)*zmgeom(i)*gamt(k)/RG
-cSB     .             /(paprs(i,k)/101325.0)**RKAPPA
-     .             *(paprs(i,k)/101325.0)**RKAPPA
-     .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))
-            ENDIF
-c
-c           finalement, les coefficients d'echange sont obtenus:
-c
-            zcdn=SQRT(zdu2) / zmgeom(i) * RG
-c
-          IF (opt_ec) THEN
-            z2geomf=zgeop(i,k-1)+zgeop(i,k)
-            zalm2=(0.5*ckap/RG*z2geomf
-     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
-            zalh2=(0.5*ckap/rg*z2geomf
-     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
-            IF (zri(i).LT.0.0) THEN  ! situation instable
-               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
-     .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
-               zscf = SQRT(-zri(i)*zscf)
-               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
-               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
-               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
-               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
-            ELSE ! situation stable
-               zscf=SQRT(1.+cd*zri(i))
-               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
-               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
-            ENDIF
-          ELSE
-            zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
-     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
-            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
-            pcfm(i,k)= zl2(i)* pcfm(i,k)
-            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
-          ENDIF
-      ENDDO
-      ENDDO
-c
-c Au-dela du sommet, pas de diffusion turbulente:
-c
-      DO i = 1, knon
-         IF (itop(i)+1 .LE. klev) THEN
-            DO k = itop(i)+1, klev
-               pcfh(i,k) = 0.0
-               pcfm(i,k) = 0.0
-            ENDDO
-         ENDIF
-      ENDDO
-c
+      END DO
       RETURN
       END
+
       SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
      .                  pcfm, pcfh)
