Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 388)
+++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F	(revision 389)
@@ -482,4 +482,5 @@
       DO j = 1, knon
          yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
+     $      + 0.000014 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))
          yrugm(j) = MAX(1.5e-05,yrugm(j))
       ENDDO
@@ -1100,5 +1101,7 @@
       REAL zcfm1(klon), zcfm2(klon)
       REAL zcfh1(klon), zcfh2(klon)
-      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
+c$$$      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn
+c$$$cPB differenciation coefficient de frottement drag et flux chaleur
+      REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn,zcdh, rugh
       REAL zscf, zucf, zcr
       REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs
@@ -1183,6 +1186,126 @@
 c Calculer le frottement au sol (Cdrag)
 c
-      CALL clcdrag(knon, nsrf, zxli, u, v, t, q, zgeop,
-     .             ts, qsurf, rugos, pcfm, pcfh, zcdn, zri)
+c$$$c PB
+c$$$c essais d'itération pour l'océan
+c
+      rugh = 1.3 e-4
+      IF (nsrf.EQ.is_oce) THEN 
+      DO k=1,10
+c$$$        WRITE(*,*) 'k',k
+c$$$        WRITE(*,*) rugos(100)
+      DO i = 1, knon
+         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
+         zdphi=zgeop(i,1)
+         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
+c$$$         ztsolv = ts(i) * (1.0+RETV*qsurf(i))
+         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
+     .       *(1.+RETV*q(i,1))
+         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
+         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
+c PB ajout drag neutre pour flux chaleur
+         zcdh = ckap**2/(log(1.+zgeop(i,1)/(RG*rugos(i)))
+     $       * log(1.+zgeop(i,1)/(RG*rugh)))
+         IF (zri(i) .ge. 0.) THEN ! situation stable
+           IF (.NOT.zxli) THEN
+           zscf=SQRT(1.+cd*ABS(zri(i)))
+           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
+!           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf)
+           zcfm1(i) = zcdn * FRIV
+           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
+!           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
+c PB avec drag neutre pour flux chaleur different
+c           zcfh1(i) = zcdn * FRIH
+           zcfh1(i) = zcdh * FRIH
+           pcfm(i,1) = zcfm1(i)
+           pcfh(i,1) = zcfh1(i)
+           ELSE
+           pcfm(i,1) = zcdn* fsta(zri(i))
+           pcfh(i,1) = zcdn* fsta(zri(i))
+           ENDIF
+         ELSE ! situation instable
+           IF (.NOT.zxli) THEN
+           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
+     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
+           zcfm2(i) = zcdn*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
+c PB ajout pour prendre drag neutre des flux chaleur different drag neutre vent
+           zucf=1./(1.+3.0*cb*cc*zcdh*SQRT(ABS(zri(i))
+     .              *(1.0+zgeop(i,1)/(RG*rugh))))
+c$$$           zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
+           zcfh2(i) = zcdh*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
+           pcfm(i,1) = zcfm2(i)
+           pcfh(i,1) = zcfh2(i)
+           ELSE
+           pcfm(i,1) = zcdn* fins(zri(i))
+           pcfh(i,1) = zcdn* fins(zri(i))
+           ENDIF
+           zcr = (0.0016/(zcdh*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
+           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdh*(1.0+zcr**1.25)**(1./1.25)
+         ENDIF
+C
+c$$$C PB test drag
+c$$$         pcfm(i,1)=zcdn
+c$$$         pcfh(i,1) = zcdn
+         rugos(i)= 0.018*pcfm(i,1) * zdu2/RG 
+     $       +0.11*0.000014/sqrt(pcfm(i,1) * zdu2)
+         rugh = 0.62*0.000014/sqrt(pcfm(i,1) * zdu2)+1.4e-5
+      END DO
+      END DO 
+
+      ELSE 
+      DO i = 1, knon
+         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
+         zdphi=zgeop(i,1)
+         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
+c         ztsolv = ts(i) * (1.0+RETV*qsurf(i))
+         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
+     .       *(1.+RETV*q(i,1))
+         zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd)
+         zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2
+c pb ajout pour avoir drage neutre flux chaleur differents
+         zcdh = ckap**2/(log(1.+zgeop(i,1)/(RG*rugos(i)))
+     $       * log(1.+zgeop(i,1)/(RG*1.3e-4)))
+         IF (zri(i) .ge. 0.) THEN ! situation stable
+           IF (.NOT.zxli) THEN
+           zscf=SQRT(1.+cd*ABS(zri(i)))
+           FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1)
+!           zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf)
+           zcfm1(i) = zcdn * FRIV
+           FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 )
+!           zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf)
+c$$$C Modif pour drag neutre flux chaleur differents
+c$$$           zcfh1(i) = zcdn * FRIH
+           zcfh1(i) = zcdh * FRIH
+           pcfm(i,1) = zcfm1(i)
+           pcfh(i,1) = zcfh1(i)
+           ELSE
+           pcfm(i,1) = zcdn* fsta(zri(i))
+           pcfh(i,1) = zcdn* fsta(zri(i))
+           ENDIF
+         ELSE ! situation instable
+           IF (.NOT.zxli) THEN
+           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
+     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
+           zcfm2(i) = zcdn*amax1((1.-2.0*cb*zri(i)*zucf),0.1)
+C PB ajout pour drage neutre flux chaleur differents
+           zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i))
+     .              *(1.0+zgeop(i,1)/(RG*rugos(i)))))
+c$$$           zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
+           zcfh2(i) = zcdh*amax1((1.-3.0*cb*zri(i)*zucf),0.1)
+           pcfm(i,1) = zcfm2(i)
+           pcfh(i,1) = zcfh2(i)
+           ELSE
+           pcfm(i,1) = zcdn* fins(zri(i))
+           pcfh(i,1) = zcdn* fins(zri(i))
+           ENDIF
+           zcr = (0.0016/(zcdn*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
+           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
+         ENDIF
+C
+c$$$C PB test drag
+c$$$         pcfm(i,1)=zcdn
+c$$$         pcfh(i,1) = zcdn
+      END DO
+c$$$CPB fin test iterations
+      ENDIF 
 c
 c Calculer les coefficients turbulents dans l'atmosphere