Index: /LMDZ6/trunk/libf/phylmd/cdrag_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4721)
+++ /LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4722)
@@ -15,7 +15,7 @@
  SUBROUTINE cdrag(knon,  nsrf,   &
      speed, t1,    q1,    zgeop1, &
-     psol,  tsurf, qsurf, z0m, z0h,  &
+     psol, pblh,  tsurf, qsurf, z0m, z0h,  &
      ri_in, iri_in, &
-     cdm,  cdh,  zri,   pref)
+     cdm,  cdh,  zri,   pref, prain, tsol , pat1)
 
   USE dimphy
@@ -81,14 +81,20 @@
   
   INTEGER, INTENT(IN)                      :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
-  REAL, DIMENSION(klon), INTENT(IN)        :: speed ! module du vent au 1er niveau du modele
-  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1! geopotentiel au 1er niveau du modele
-  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf ! Surface temperature (K)
-  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf ! Surface humidity (Kg/Kg)
-  REAL, DIMENSION(klon), INTENT(IN)        :: z0m, z0h ! Rugosity at surface (m) 
+  REAL, DIMENSION(klon), INTENT(IN)        :: speed    ! module du vent au 1er niveau du modele
+  REAL, DIMENSION(klon), INTENT(IN)        :: zgeop1   ! geopotentiel au 1er niveau du modele
+  REAL, DIMENSION(klon), INTENT(IN)        :: tsurf    ! Surface temperature (K)
+  REAL, DIMENSION(klon), INTENT(IN)        :: qsurf    ! Surface humidity (Kg/Kg)
+  REAL, DIMENSION(klon), INTENT(INOUT)     :: z0m, z0h ! Rugosity at surface (m) 
   REAL, DIMENSION(klon), INTENT(IN)        :: ri_in ! Input Richardson 1st layer for first guess calculations of screen var. 
   INTEGER, INTENT(IN)                      :: iri_in! iflag to activate cdrag calculation using ri1
-  REAL, DIMENSION(klon), INTENT(IN)        :: t1  ! Temperature au premier niveau (K) 
-  REAL, DIMENSION(klon), INTENT(IN)        :: q1  ! humidite specifique au premier niveau (kg/kg)
-  REAL, DIMENSION(klon), INTENT(IN)        :: psol ! pression au sol
+  REAL, DIMENSION(klon), INTENT(IN)        :: t1       ! Temperature au premier niveau (K) 
+  REAL, DIMENSION(klon), INTENT(IN)        :: q1       ! humidite specifique au premier niveau (kg/kg)
+  REAL, DIMENSION(klon), INTENT(IN)        :: psol     ! pression au sol
+
+!------------------ Rajout pour COARE (OT2018) --------------------
+  REAL, DIMENSION(klon), INTENT(INOUT)     :: pblh  !hauteur de CL
+  REAL, DIMENSION(klon), INTENT(IN)        :: prain !rapport au precipitation
+  REAL, DIMENSION(klon), INTENT(IN)        :: tsol  !SST imposé sur la surface oceanique
+  REAL, DIMENSION(klon), INTENT(IN)        :: pat1  !pression premier lev
 
 
@@ -99,5 +105,5 @@
   REAL, DIMENSION(klon), INTENT(OUT)       :: cdh   ! Drag coefficient for heat flux
   REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
-  REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
+  REAL, DIMENSION(klon), INTENT(INOUT)       :: pref  ! Pression au niveau zgeop/RG
 
 ! Variables Locales
@@ -116,12 +122,10 @@
   REAL                   MU, CM, CH, B, CMstar, CHstar
   REAL                   PM, PH, BPRIME
-  REAL                   C
   INTEGER                ng_q1                      ! Number of grids that q1 < 0.0
   INTEGER                ng_qsurf                   ! Number of grids that qsurf < 0.0
-  INTEGER                i
+  INTEGER                i, k
   REAL                   zdu2, ztsolv
   REAL                   ztvd, zscf
   REAL                   zucf, zcr
-  REAL                   friv, frih
   REAL, DIMENSION(klon) :: FM, FH                   ! stability functions
   REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
@@ -129,4 +133,38 @@
   REAL, DIMENSION(klon) :: sm, prandtl              ! Stability function from atke turbulence scheme
   REAL                  ri0, ri1, cn                ! to have dimensionless term in sm and prandtl
+
+!------------------ Rajout (OT2018) --------------------
+!------------------ Rajout pour les appelles BULK (OT) --------------------
+  REAL, DIMENSION(klon,2) :: rugos_itm 
+  REAL, DIMENSION(klon,2) :: rugos_ith 
+  REAL, PARAMETER         :: tol_it_rugos=1.e-4
+  REAL, DIMENSION(3)      :: coeffs
+  LOGICAL                 :: mixte
+  REAL                    :: z_0m
+  REAL                    :: z_0h
+  REAL, DIMENSION(klon)   :: cdmm
+  REAL, DIMENSION(klon)   :: cdhh
+
+!------------------RAJOUT POUR ECUME -------------------
+
+  REAL, DIMENSION(klon) :: PQSAT 
+  REAL, DIMENSION(klon) :: PSFTH
+  REAL, DIMENSION(klon) :: PFSTQ 
+  REAL, DIMENSION(klon) :: PUSTAR
+  REAL, DIMENSION(klon) :: PCD      ! Drag coefficient for momentum
+  REAL, DIMENSION(klon) :: PCDN     ! Drag coefficient for momentum
+  REAL, DIMENSION(klon) :: PCH      ! Drag coefficient for momentum
+  REAL, DIMENSION(klon) :: PCE      ! Drag coefficient for momentum
+  REAL, DIMENSION(klon) :: PRI 
+  REAL, DIMENSION(klon) :: PRESA
+  REAL, DIMENSION(klon) :: PSSS
+
+  LOGICAL               :: OPRECIP
+  LOGICAL               :: OPWEBB
+  LOGICAL               :: OPERTFLUX
+  LOGICAL               :: LPRECIP
+  LOGICAL               :: LPWG
+
+
 
   LOGICAL, SAVE :: firstcall = .TRUE.
@@ -136,4 +174,6 @@
   INTEGER, SAVE :: iflag_corr_insta
   !$OMP THREADPRIVATE(iflag_corr_insta)
+  LOGICAL, SAVE :: ok_cdrag_iter
+  !$OMP THREADPRIVATE(ok_cdrag_iter)
 
 !===================================================================c
@@ -170,7 +210,7 @@
 
 ! Consistent with atke scheme
-cn=(1./sqrt(cepsilon))**(2/3)
-ri0=2./rpi*(cinf - cn)*ric/cn
-ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
+  cn=(1./sqrt(cepsilon))**(2./3.)
+  ri0=2./rpi*(cinf - cn)*ric/cn
+  ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
 
 
@@ -209,13 +249,24 @@
 ! On choisit les fonctions de stabilite utilisees au premier appel
 !**************************************************************************
-  IF (firstcall) THEN
+ IF (firstcall) THEN
    iflag_corr_sta=2
    iflag_corr_insta=2
+   ok_cdrag_iter = .FALSE.
  
    CALL getin_p('iflag_corr_sta',iflag_corr_sta)
    CALL getin_p('iflag_corr_insta',iflag_corr_insta)
+   CALL getin_p('ok_cdrag_iter',ok_cdrag_iter)
 
    firstcall = .FALSE.
  ENDIF
+
+!------------------ Rajout (OT2018) --------------------
+!---------    Rajout pour itération sur rugosité     ----------------
+  rugos_itm(:,2) = z0m
+  rugos_itm(:,1) = 3.*tol_it_rugos*z0m
+
+  rugos_ith(:,2) = z0h !cp nouveau rugos_it
+  rugos_ith(:,1) = 3.*tol_it_rugos*z0h
+!--------------------------------------------------------------------
 
 !xxxxxxxxxxxxxxxxxxxxxxx
@@ -227,211 +278,348 @@
 !***********************
      
-
      zdu2 = MAX(CEPDU2, speed(i)**2)
      pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
-                 (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero 
+          (1.+ RETV * max(q1(i),0.0))))           ! negative q1 set to zero
      ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0))    ! negative qsurf set to zero
-     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) &
+     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
           *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
-     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
+     
+!------------------ Rajout (OT2018) --------------------
+!--------------   ON DUPLIQUE POUR METTRE ITERATION SUR OCEAN    ------------------------     
      IF (iri_in.EQ.1) THEN
        zri(i) = ri_in(i)
      ENDIF
 
+     IF (nsrf == is_oce) THEN
+        
+!------------------  Pour Core 2 choix Core Pur ou Core Mixte  --------------------
+       IF ((choix_bulk > 1 .AND. choix_bulk < 4) .AND. (nsrf .eq. is_oce)) THEN
+         IF(choix_bulk .eq. 2) THEN
+           mixte = .false.
+         ELSE
+            mixte = .true.
+         ENDIF
+         call clc_core_cp ( sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),t1(i),q1(i),&
+             zgeop1(i)/RG, zgeop1(i)/RG,  zgeop1(i)/RG,&
+             psol(i),nit_bulk,mixte,&
+             coeffs,z_0m,z_0h)
+         cdmm(i) = coeffs(1)
+         cdhh(i) = coeffs(2)
+         cdm(i)=cdmm(i)
+         cdh(i)=cdhh(i)
+         write(*,*) "clc_core cd ch",cdmm(i),cdhh(i)
+
+!------------------                 Pour ECUME                 --------------------
+       ELSE IF ((choix_bulk .eq. 4) .and. (nsrf .eq. is_oce)) THEN
+         OPRECIP = .false.
+         OPWEBB  = .false.
+         OPERTFLUX = .false.
+         IF (nsrf .eq. is_oce) THEN
+           PSSS = 0.0
+         ENDIF
+         call ini_csts
+         call ecumev6_flux( z_0m,t1(i),tsurf(i),&
+             q1(i),qsurf(i),sqrt(zdu2),zgeop1(i)/RG,PSSS,zgeop1(i)/RG,&
+             psol(i),pat1(i), OPRECIP, OPWEBB,&
+             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,&
+             PCH,PCE,PRI,PRESA,prain,&
+             z_0h,OPERTFLUX,coeffs)
+
+         cdmm(i) = coeffs(1)
+         cdhh(i) = coeffs(2)
+         cdm(i)=cdmm(i)
+         cdh(i)=cdhh(i)
+   
+         write(*,*) "ecume cd ch",cdmm(i),cdhh(i)
+
+!------------------                 Pour COARE CNRM                 --------------------
+       ELSE IF ((choix_bulk .eq. 5) .and. (nsrf .eq. is_oce)) THEN
+         LPRECIP = .false.
+         LPWG    = .false.
+         call ini_csts
+         call coare30_flux_cnrm(z_0m,t1(i),tsurf(i), q1(i),  &
+             sqrt(zdu2),zgeop1(i)/RG,zgeop1(i)/RG,psol(i),qsurf(i),PQSAT, &
+             PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, &
+             PRESA,prain,pat1(i),z_0h, LPRECIP, LPWG, coeffs) 
+         cdmm(i) = coeffs(1)
+         cdhh(i) = coeffs(2)
+         cdm(i)=cdmm(i)
+         cdh(i)=cdhh(i)
+         write(*,*) "coare CNRM cd ch",cdmm(i),cdhh(i)
+
+!------------------                 Pour COARE Maison                 --------------------
+       ELSE IF ((choix_bulk .eq. 1) .and. (nsrf .eq. is_oce)) THEN
+         IF ( pblh(i) .eq. 0. ) THEN
+           pblh(i) = 1500.
+         ENDIF
+         write(*,*) "debug size",size(coeffs)
+         call coare_cp(sqrt(zdu2),t1(i)-tsurf(i),q1(i)-qsurf(i),&
+               t1(i),q1(i),&
+               zgeop1(i)/RG,zgeop1(i)/RG,zgeop1(i)/RG,&
+               psol(i), pblh(i),&
+               nit_bulk,&
+               coeffs,z_0m,z_0h)
+         cdmm(i) = coeffs(1)
+         cdhh(i) = coeffs(3)
+         cdm(i)=cdmm(i)
+         cdh(i)=cdhh(i)
+         write(*,*) "coare cd ch",cdmm(i),cdhh(i)
+       ELSE
+!------------------                 Pour La param LMDZ (ocean)              --------------------
+         zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
+         IF (iri_in.EQ.1) THEN
+           zri(i) = ri_in(i)
+         ENDIF
+       ENDIF
+     
+
+!----------------------- Rajout des itérations --------------
+       DO  k=1,nit_bulk
 
 ! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
 !********************************************************************
-
-     zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
-     cdmn(i) = zzzcd*zzzcd
-     cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
+         zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*rugos_itm(i,2)))
+         cdmn(i) = zzzcd*zzzcd
+         cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*rugos_ith(i,2))))
+
+! Calcul des fonctions de stabilite FMs, FHs, FMi, FHi :
+!*******************************************************
+         IF (zri(i) .LT. 0.) THEN    
+           SELECT CASE (iflag_corr_insta)
+             CASE (1) ! Louis 1979 + Mascart 1995
+                  MU=LOG(MAX(z0m(i)/z0h(i),0.01))
+                  CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
+                  PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
+                  CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
+                  PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
+                  CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
+                     & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
+                     & * ((zgeop1(i)/(RG*z0h(i)))**PH)
+                  CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
+                     & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
+                     & * ((zgeop1(i)/(RG*z0m(i)))**PM)
+                  FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
+                  FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
+             CASE (2) ! Louis 1982
+                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
+                        *(1.0+zgeop1(i)/(RG*z0m(i)))))
+                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
+                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
+             CASE (3) ! Laurent Li
+                  FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
+                  FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
+             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
+                  sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
+                  prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
+                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
+                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
+             CASE default ! Louis 1982
+                  zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
+                         *(1.0+zgeop1(i)/(RG*z0m(i)))))
+                  FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
+                  FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
+           END SELECT
+! Calcul des drags
+           cdmm(i)=cdmn(i)*FM(i)
+           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
+! Traitement particulier des cas oceaniques
+! on applique Miller et al 1992 en l'absence de gustiness 
+           IF (nsrf == is_oce) THEN
+!            cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
+             IF (iflag_gusts==0) THEN
+               zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
+               cdhh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
+             ENDIF
+             cdmm(i)=MIN(cdmm(i),cdmmax)
+             cdhh(i)=MIN(cdhh(i),cdhmax)
+!             write(*,*) "LMDZ cd ch",cdmm(i),cdhh(i)
+           END IF
+         ELSE                          
+
+!'''''''''''''''
+! Cas stables :
+!'''''''''''''''
+           zri(i) = MIN(20.,zri(i))
+           SELECT CASE (iflag_corr_sta)
+             CASE (1) ! Louis 1979 + Mascart 1995
+                  FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
+                  FH(i)=FM(i)
+             CASE (2) ! Louis 1982
+                  zscf = SQRT(1.+CD*ABS(zri(i)))
+                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
+                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
+             CASE (3) ! Laurent Li
+                  FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
+                  FH(i)=FM(i)
+             CASE (4)  ! King 2001
+                  IF (zri(i) .LT. C2/2.) THEN
+                    FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
+                    FH(i)=  FM(i)
+                  ELSE
+                    FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
+                    FH(i)= FM(i)
+                  ENDIF
+             CASE (5) ! MO
+                  if (zri(i) .LT. 1./alpha) then
+                    FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
+                    FH(i)=FM(i)
+                  else 
+                    FM(i)=MAX(1E-7,f_ri_cd_min)
+                    FH(i)=FM(i)
+                  endif
+             CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
+                  sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
+                  prandtl(i) = pr_neut + zri(i) * pr_slope
+                  FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
+                  FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
+             CASE default ! Louis 1982
+                  zscf = SQRT(1.+CD*ABS(zri(i)))
+                  FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
+                  FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
+           END SELECT
+
+! Calcul des drags
+           
+           cdmm(i)=cdmn(i)*FM(i)
+           cdhh(i)=f_cdrag_ter*cdhn(i)*FH(i)
+           IF (choix_bulk == 0) THEN
+             cdm(i)=cdmn(i)*FM(i)
+             cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
+           ENDIF
+
+           IF (nsrf.EQ.is_oce) THEN
+             cdhh(i)=f_cdrag_oce*cdhn(i)*FH(i)
+             cdmm(i)=MIN(cdm(i),cdmmax)
+             cdhh(i)=MIN(cdh(i),cdhmax)
+           ENDIF
+           IF (ok_cdrag_iter) THEN
+             rugos_itm(i,1) = rugos_itm(i,2)
+             rugos_ith(i,1) = rugos_ith(i,2)
+             rugos_itm(i,2) =  0.018*cdmm(i) * (speed(i))/RG  &
+                              +  0.11*14e-6 / SQRT(cdmm(i) * zdu2)
+
+!---------- Version SEPARATION DES Z0 ----------------------
+             IF (iflag_z0_oce==0) THEN
+               rugos_ith(i,2) = rugos_itm(i,2)
+             ELSE IF (iflag_z0_oce==1) THEN
+               rugos_ith(i,2) = 0.40*14e-6 / SQRT(cdmm(i) * zdu2)
+             ENDIF
+           ENDIF
+         ENDIF
+         IF (ok_cdrag_iter) THEN
+           rugos_itm(i,2) = MAX(1.5e-05,rugos_itm(i,2))
+           rugos_ith(i,2) = MAX(1.5e-05,rugos_ith(i,2))
+         ENDIF
+       ENDDO
+       IF (nsrf.EQ.is_oce) THEN
+         cdm(i)=MIN(cdmm(i),cdmmax)
+         cdh(i)=MIN(cdhh(i),cdhmax)
+       ENDIF
+       z0m = rugos_itm(:,2)
+       z0h = rugos_ith(:,2)
+     ELSE ! (nsrf == is_oce)
+       zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
+       IF (iri_in.EQ.1) THEN
+         zri(i) = ri_in(i)
+       ENDIF
+
+! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
+!********************************************************************
+       zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i)))
+       cdmn(i) = zzzcd*zzzcd
+       cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
 
 
 ! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
 !*******************************************************
-
-
-
 !''''''''''''''
 ! Cas instables
 !''''''''''''''
-
- IF (zri(i) .LT. 0.) THEN    
-
-
-        SELECT CASE (iflag_corr_insta)
-    
-        CASE (1) ! Louis 1979 + Mascart 1995
-
-           MU=LOG(MAX(z0m(i)/z0h(i),0.01))
-           CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
-           PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
-           CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
-           PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
-           CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
-            & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
-            & * ((zgeop1(i)/(RG*z0h(i)))**PH)
-           CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
-            & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
-            & * ((zgeop1(i)/(RG*z0m(i)))**PM)
-         
-
-
-
-           FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
-           FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
-
-        CASE (2) ! Louis 1982
-
-           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
-                *(1.0+zgeop1(i)/(RG*z0m(i)))))
-           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
-           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
-
-
-        CASE (3) ! Laurent Li
-
-           
-           FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
-           FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
-
-         CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
-         
-           sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
-           prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
-           FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
-           FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
-
-         CASE default ! Louis 1982
-
-           zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
-                *(1.0+zgeop1(i)/(RG*z0m(i)))))
-           FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
-           FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
-
-
+       IF (zri(i) .LT. 0.) THEN    
+         SELECT CASE (iflag_corr_insta)
+           CASE (1) ! Louis 1979 + Mascart 1995
+                MU=LOG(MAX(z0m(i)/z0h(i),0.01))
+                CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3)
+                PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3)
+                CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3)
+                PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3)
+                CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
+                   & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i)))       &
+                   & * ((zgeop1(i)/(RG*z0h(i)))**PH)
+                CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) &
+                   & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i)))         &
+                   & * ((zgeop1(i)/(RG*z0m(i)))**PM)
+                FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i))))
+                FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i))))
+           CASE (2) ! Louis 1982
+                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
+                       *(1.0+zgeop1(i)/(RG*z0m(i)))))
+                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
+                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
+           CASE (3) ! Laurent Li
+                FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
+                FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min)
+           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
+                 sm(i) = 2./rpi * (cinf - cn) * atan(-zri(i)/ri0) + cn
+                 prandtl(i) = -2./rpi * (pr_asym - pr_neut) * atan(zri(i)/ri1) + pr_neut
+                 FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
+                 FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
+           CASE default ! Louis 1982
+                zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) &
+                      *(1.0+zgeop1(i)/(RG*z0m(i)))))
+                FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
+                FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
          END SELECT
-
-
-
 ! Calcul des drags
-
-
-       cdm(i)=cdmn(i)*FM(i)
-       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
-
-
-! Traitement particulier des cas oceaniques
-! on applique Miller et al 1992 en l'absence de gustiness 
-
-  IF (nsrf == is_oce) THEN
-!        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
-
-        IF(iflag_gusts==0) THEN
-           zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
-           cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25)
-        ENDIF
-
-
-        cdm(i)=MIN(cdm(i),cdmmax)
-        cdh(i)=MIN(cdh(i),cdhmax)
-
-  END IF
-
-
-
- ELSE                          
-
+         cdm(i)=cdmn(i)*FM(i)
+         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
+       ELSE                          
 !'''''''''''''''
 ! Cas stables :
 !'''''''''''''''
-        zri(i) = MIN(20.,zri(i))
-
-       SELECT CASE (iflag_corr_sta)
-    
-        CASE (1) ! Louis 1979 + Mascart 1995
-   
-           FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
-           FH(i)=FM(i)
-
-
-        CASE (2) ! Louis 1982
-           
-           zscf = SQRT(1.+CD*ABS(zri(i)))
-           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
-           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
-
-
-        CASE (3) ! Laurent Li
-   
-           FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
-           FH(i)=FM(i)
-
-
-        CASE (4)  ! King 2001
-          
-           if (zri(i) .LT. C2/2.) then
-           FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
-           FH(i)=  FM(i)
-
-
-           else
-           FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
-           FH(i)= FM(i)
-           endif
-
-
-        CASE (5) ! MO
-   
-          if (zri(i) .LT. 1./alpha) then
-
-           FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
-           FH(i)=FM(i)
-
-           else 
-
-
-           FM(i)=MAX(1E-7,f_ri_cd_min)
-           FH(i)=FM(i)
-
-          endif
-
-        CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
-          sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
-          prandtl(i) = pr_neut + zri(i) * pr_slope
-          FM(i) = MAX((sm(i)**(3/2) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1/2)),f_ri_cd_min)
-          FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
-
-        CASE default ! Louis 1982
-
-           zscf = SQRT(1.+CD*ABS(zri(i)))
-           FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
-           FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
-
-
-
-   END SELECT
-
+         zri(i) = MIN(20.,zri(i))
+         SELECT CASE (iflag_corr_sta)
+           CASE (1) ! Louis 1979 + Mascart 1995
+                FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min)
+                FH(i)=FM(i)
+           CASE (2) ! Louis 1982
+                zscf = SQRT(1.+CD*ABS(zri(i)))
+                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
+                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
+           CASE (3) ! Laurent Li
+                FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min)
+                FH(i)=FM(i)
+           CASE (4)  ! King 2001
+                if (zri(i) .LT. C2/2.) then
+                  FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min)
+                  FH(i)=  FM(i)
+                else
+                  FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min)
+                  FH(i)= FM(i)
+                endif
+           CASE (5) ! MO
+                if (zri(i) .LT. 1./alpha) then
+                  FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min)
+                  FH(i)=FM(i)
+                else 
+                  FM(i)=MAX(1E-7,f_ri_cd_min)
+                  FH(i)=FM(i)
+                endif
+           CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
+                sm(i) = MAX(0.,cn*(1.-zri(i)/ric))
+                prandtl(i) = pr_neut + zri(i) * pr_slope
+                FM(i) = MAX((sm(i)**(3./2.) * sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1./2.)),f_ri_cd_min)
+                FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
+           CASE default ! Louis 1982
+                zscf = SQRT(1.+CD*ABS(zri(i)))
+                FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min)
+                FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min )
+         END SELECT
 ! Calcul des drags
-
-
-       cdm(i)=cdmn(i)*FM(i)
-       cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
-
-       IF(nsrf.EQ.is_oce) THEN
-
-        cdh(i)=f_cdrag_oce*cdhn(i)*FH(i)
-        cdm(i)=MIN(cdm(i),cdmmax)
-        cdh(i)=MIN(cdh(i),cdhmax)
-
-      ENDIF
-
-
- ENDIF
-
-!xxxxxxxxxxx
+         cdm(i)=cdmn(i)*FM(i)
+         cdh(i)=f_cdrag_ter*cdhn(i)*FH(i)
+       ENDIF
+     ENDIF ! fin du if (nsrf == is_oce)
   END DO   !  Fin de la boucle sur l'horizontal
-!xxxxxxxxxxx
-! ================================================================= c
      
 END SUBROUTINE cdrag
Index: /LMDZ6/trunk/libf/phylmd/clc_core_cp.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/clc_core_cp.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/clc_core_cp.F90	(revision 4722)
@@ -0,0 +1,232 @@
+!
+! $Id$
+!
+      SUBROUTINE clc_core_cp(du, dt, dq, t1, q1, &
+     				zu, zt, zq, &
+     					P,&
+     				n_it, mixte, &
+    			 coeffs, rugosm, rugosh)
+
+      IMPLICIT NONE
+
+
+  !  INTEGER     :: i,j,k
+  REAL, INTENT(IN) :: du,dt,dq,t1,q1
+  REAL, INTENT(IN) :: zu,zt,zq, P
+  INTEGER, INTENT(IN):: n_it
+  LOGICAL, INTENT(IN) :: mixte
+  REAL, DIMENSION(3), INTENT(OUT) :: coeffs
+  real, intent(out) :: rugosm
+  real, intent(out) :: rugosh
+
+
+!  REAL :: du,dt,t1,dq,q1
+  REAL :: dt_u, dq_u
+
+
+
+  REAL :: cd, ch, cd_rt, cd_n10, ce_n10, ce,&
+       u_N, X,&
+       rho,cpa,le
+
+  REAL :: usr,tsr,qsr
+  REAL, DIMENSION(3) :: zeta                           ! dans l'ordre: it courant, it precedente, var
+  REAL, DIMENSION(2,3) :: psi
+  REAL, PARAMETER :: g = 9.81, k=.41, sqrt_k = .640312 !avec g= gravité, k = cste von karman, sqrt_k = racine(k)
+  REAL :: logzu_10,logzt_10,logzq_10,logzt_zu,logzq_zu,&
+       u,cd_n10_rt,ch_n10_rt,ch_n10,chi,z_0m,z_0h,tv
+  INTEGER :: i,j,m1,m2,m3,m4,m5
+
+
+  REAL :: phih,phid,phie
+
+  
+  REAL, parameter :: alpha=2.7e-3,&           !alpha = a1 ds formulation smith 
+       beta = 1.42e-4,&                       !beta = a2 ds la formulation smith
+       gamma = 7.64e-5,&        !Large and yaeger 2006              !gamma = a3 ds formulation smith
+!       gamma = 13.09e-3,&         !Large and yaeger 2004
+       q0 = 1.64474                           ! q0 dzfinis mais pas reutilisée ds la suite ??
+
+  REAL, dimension(3) :: z
+  integer, dimension(2) :: tmp_shape
+  z = [zu, zt, zq]
+ 
+
+
+  logzu_10 = LOG(zu/10.)
+  logzt_10 = LOG(zt/10.)
+  logzq_10 = LOG(zq/10.)
+  logzt_zu = logzt_10 - logzu_10
+  logzq_zu = logzq_10 - logzu_10
+
+  cpa = 1004.67                               !!! pas reutilise ds nvelle formulation...
+
+  tv = t1 * (1+.608*q1)                       !!!!! tv pas reutiliser ds nvelle formulation...
+  rho = p / (287.1 * t1 * (1.+.61*q1))        !!!! rho pas reutiliser ds nvelle formulation....
+
+
+  le = (2.501-.00237*(t1-273.15-dt))*1.e6     !!!!!! le pas reutiliser ds nvelle formulation....
+  
+  u = max(du,.5)                              !!!! u = vitesse du vent 1 er niveau
+  u_N = u                                     !!!!!!! u_N utilisée pour le calcul des premier calcul...
+  
+
+  !!!!! initialisation des parametres et premier calcul de z_0 cd et ch
+
+
+  if (mixte) then 
+     z_0m = 1.e-4                             !!!!! pour l'initialisation on fixe z_0 a 10-4
+     cd_n10 = k**2 / (log(10/z_0m))**2        !!!!! on calcul ensuite cd a 10 m d'après smith 
+     cd_n10_rt = sqrt(cd_n10)   
+  else
+! Large and yaeger 2006
+     cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur
+! Large and yaeger 2004
+!     cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) 
+     cd_n10_rt = sqrt(cd_n10)                  !!!!!!
+     z_0m = 10. * exp(-k/cd_n10_rt)            !!!!!!! on calcul le premier z_0
+  endif
+
+
+  if ( dt .gt. 0. ) then
+     ch_n10 = 0.018 * cd_n10_rt                !!!!!!! si regime stable calcul du ch
+!     z_0h = 10. * exp(-k/ch_n10)
+  else
+     ch_n10 = 0.0327 * cd_n10_rt               !!!!!!! si regime instable calcul du ch
+!     z_0h = 10. * exp(-k/ch_n10)
+  endif
+
+  ce_n10 = 3.46e-2 *cd_n10_rt                  !!!!!! calcul ensuite de ce
+
+
+  dt_u = dt                                    !!! def
+  dq_u = dq                                    !!! def
+
+  cd = cd_n10                                  !!! def
+  ch = ch_n10                                  !!! def
+  ce = ce_n10                                  !!! def
+
+  cd_rt = sqrt(cd)                             !!! def
+
+  !  open(7,file="err_rel.dat",position="append")
+
+
+!!!!!!!!! début des itérations....
+
+  do i = 1,n_it
+
+     usr = cd_rt * u 
+     tsr = ch / cd_rt *  dt_u
+     qsr = ce / cd_rt * dq_u
+
+
+     zeta = k*g / (usr**2) * tsr / t1&
+          * z
+
+     zeta = sign( min(abs(zeta),10.0), zeta )
+
+
+!!!! calcul du zeta pour connaitre le signe et donc en deduire le regime...
+
+
+
+
+     IF ( zeta(1) .GT. 0 ) THEN !!! si regime stable alors valeurs pr chi et psi 
+
+        chi = 0.018 
+
+
+        psi(1,1) = -5.*zeta(1)
+        psi(2,1) = psi(1,1)
+
+     ELSE !!!!! si regime instable alors valeurs pr psi et chi !!!
+
+        chi = 0.0327
+
+        X = sqrt( sqrt( 1-16.*zeta(1) ))
+        psi(1,1) =  3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X)
+        psi(2,1) =  2.*LOG( (1 + X**2) / 2 )
+
+     END IF
+
+     do j =2,3
+        IF ( zeta(j) .GT. 0 ) THEN
+           psi(1,j) = -5.*zeta(j)
+           psi(2,j) = psi(1,1)
+
+        ELSE
+           X = sqrt( sqrt( 1-16.*zeta(j) ))
+           psi(1,j) =  3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X)
+           psi(2,j) =  2.*LOG( (1 + X**2) / 2 )
+
+        END IF
+     enddo
+
+
+     dt_u = dt - tsr / k * &
+          ( logzt_zu + psi(2,1) - psi(2,2) )
+
+     dq_u = dq - qsr / k * &
+          ( logzq_zu + psi(2,1) - psi(2,3) )
+
+
+
+     u_N = u / ( 1 + cd_n10_rt / k * ( logzu_10 - psi(1,1)) )
+
+     if (mixte) then
+        z_0m=0.018*usr**2/9.81 + 0.11*14e-6 / (usr)
+!        z_0h=(0.11*14e-6 / (usr)) + 1.4e-5 
+        cd_n10 = k**2 / (log(10/z_0m))**2
+        cd_n10_rt = sqrt(cd_n10)
+!        ch_n10 = k**2 / ((log(10/z_0m))*(log(10/z_0h)))
+!        ch_n10_rt = sqrt(ch_n10)
+        ch_n10 = 1.e-3
+     else
+! Large and yaeger 2006
+        cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur
+! Large and yaeger 2004
+!        cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) 
+        !        cd_n10_dun = -alpha / u_n**2 + gamma
+        cd_n10_rt = sqrt(cd_n10)
+        z_0m = 10. * exp(-k/cd_n10_rt) !!!!! c'est sur cette ligne la qu'il y a un bug... ou sur la ligne du u_N
+!        ch_n10 = chi * cd_n10_rt
+!        z_0h = 10. * exp(-k/ch_n10)
+!        ch_n10_rt = sqrt(ch_n10)
+        ch_n10 = chi * cd_n10_rt
+     endif
+
+
+     ce_n10 = 3.46e-2*cd_n10_rt
+
+
+
+     phid = 1+cd_n10_rt/k * (logzu_10 - psi(1,1))
+     phih = 1+chi/k * (logzu_10 - psi(2,1))
+     phie = 1+3.46e-2/k * (logzu_10 - psi(2,1))
+
+
+     cd_rt = cd_n10_rt / abs( phid )
+     ch = ch_n10 / ( phih * phid )
+     ce = ce_n10 / ( phie * phid )
+
+     cd=cd_rt**2
+
+  END DO
+
+  usr = cd_rt * u
+  tsr = ch / cd_rt *  dt_u
+  qsr = ce / cd_rt * dq_u
+
+
+  ! coeffs(:,m) = [ &
+  !      rho * usr**2,&
+  !      rho * cpa * usr *tsr,&
+  !      rho * le * usr *qsr]
+
+  coeffs = [ cd_rt**2, ch , ce ]
+
+  rugosm =  z_0m
+  rugosh =  z_0h
+
+      RETURN
+      END subroutine clc_core_cp
Index: /LMDZ6/trunk/libf/phylmd/clesphys.h
===================================================================
--- /LMDZ6/trunk/libf/phylmd/clesphys.h	(revision 4721)
+++ /LMDZ6/trunk/libf/phylmd/clesphys.h	(revision 4722)
@@ -43,4 +43,6 @@
 !IM seuils cdrm, cdrh
        REAL cdmmax, cdhmax
+!IM pour les params différentes Olivier Torres
+       INTEGER choix_bulk, nit_bulk, kz0
 !IM param. stabilite s/ terres et en dehors
        REAL ksta, ksta_ter, f_ri_cd_min
@@ -111,4 +113,5 @@
        COMMON/clesphys/                                                 &
 ! REAL FIRST
+! rajout choix_bulk et nit_bulk kz0 par Olivier Torres
      &       co2_ppm, solaire                                           &
      &     , RCO2, RCH4, RN2O, RCFC11, RCFC12                           &
@@ -138,4 +141,5 @@
      &     , ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad                &
      &     , iflag_con, nbapp_cv, nbapp_wk                              &
+     &     , choix_bulk, nit_bulk, kz0                                  &
      &     , iflag_ener_conserv                                         &
      &     , ok_suntime_rrtm                                            & 
Index: /LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90	(revision 4722)
@@ -0,0 +1,631 @@
+!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
+!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!SFX_LIC for details. version 1.
+!     #########
+
+!------------------------ Modif Olivier Torres pour rajout fonction et pas appel ----------------
+
+
+
+real function PSIFCTT(zet)
+
+  real, intent(in) :: zet
+
+  if(zet<0) then
+     x=(1.-(15*zet))**.5 
+     psik=2*log((1+x)/2) 
+     x=(1.-(34.15*zet))**.3333 
+     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 
+     f=zet*zet/(1+zet*zet) 
+     PSIFCTT=(1-f)*psik+f*psic   
+
+  else
+     c=min(50.,.35*zet) 
+     PSIFCTT=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
+  endif
+end FUNCTION PSIFCTT
+
+real function PSIFCTU(zet)
+
+  real, intent(in) :: zet
+
+  if (zet<0) then 
+
+     x=(1.-15.*zet)**.25 
+     psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) 
+     x=(1.-10.15*zet)**.3333 
+     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 
+     f=zet*zet/(1+zet*zet) 
+     PSIFCTU=(1-f)*psik+f*psic                                                
+  else
+     c=min(50.,.35*zet) 
+     PSIFCTU=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
+  endif
+END FUNCTION PSIFCTU
+
+!----------------------------------- Fin Modif ---------------------------------------------------
+
+
+SUBROUTINE COARE30_FLUX_CNRM(PZ0SEA,PTA,PSST,PQA,  &
+            PVMOD,PZREF,PUREF,PPS,PQSATA,PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI,&
+            PRESA,PRAIN,PPA,PZ0HSEA,LPRECIP, LPWG,coeffs )  
+!     #######################################################################
+!
+!
+!!****  *COARE25_FLUX*  
+!!
+!!    PURPOSE
+!!    -------
+!      Calculate the surface fluxes of heat, moisture, and momentum over
+!      sea surface with bulk algorithm COARE3.0. 
+!     
+!!**  METHOD
+!!    ------
+!      transfer coefficients were obtained using a dataset which combined COARE
+!      data with those from three other ETL field experiments, and reanalysis of
+!      the HEXMAX data (DeCosmos et al. 1996). 
+!      ITERMAX=3 
+!      Take account of the surface gravity waves on the velocity roughness and 
+!      hence the momentum transfer coefficient
+!        NGRVWAVES=0 no gravity waves action (Charnock) !default value
+!        NGRVWAVES=1 wave age parameterization of Oost et al. 2002
+!        NGRVWAVES=2 model of Taylor and Yelland 2001
+!
+!!    EXTERNAL
+!!    --------
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------ 
+!!      
+!!    REFERENCE
+!!    ---------
+!!      Fairall et al (2003), J. of Climate, vol. 16, 571-591
+!!      Fairall et al (1996), JGR, 3747-3764
+!!      Gosnell et al (1995), JGR, 437-442
+!!      Fairall et al (1996), JGR, 1295-1308
+!!      
+!!    AUTHOR
+!!    ------
+!!     C. Lebeaupin  *Météo-France* (adapted from C. Fairall's code)
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original     1/06/2006
+!!      B. Decharme    06/2009 limitation of Ri
+!!      B. Decharme    09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90
+!!      B. Decharme    06/2013 bug in z0 (output) computation 
+!!      J.Escobar      06/2013  for REAL4/8 add EPSILON management
+!!      C. Lebeaupin   03/2014 bug if PTA=PSST and PEXNA=PEXNS: set a minimum value
+!!                             add abort if no convergence
+!-------------------------------------------------------------------------------
+!
+!*       0.     DECLARATIONS
+!               ------------
+!
+!
+!USE MODD_SEAFLUX_n, ONLY : SEAFLUX_t
+!
+
+!----------Rajout Olive ---------
+USE dimphy
+USE indice_sol_mod
+
+!--------------------------------
+
+USE MODD_CSTS,       ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
+                            XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, &
+                            XP00
+
+!USE MODD_SURF_ATM,   ONLY : XVZ0CM
+!
+!USE MODD_SURF_PAR,   ONLY : XUNDEF, XSURF_EPSILON
+!USE MODD_WATER_PAR
+!
+!USE MODI_SURFACE_RI
+!USE MODI_WIND_THRESHOLD
+!USE MODE_COARE30_PSI
+!
+!USE MODE_THERMOS
+!
+!
+!USE MODI_ABOR1_SFX
+!
+!
+!USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
+!USE PARKIND1  ,ONLY : JPRB
+!
+IMPLICIT NONE
+!
+!*      0.1    declarations of arguments
+!
+!
+!
+!TYPE(SEAFLUX_t), INTENT(INOUT) :: S
+!
+REAL, DIMENSION(klon), INTENT(IN)       :: PTA   ! air temperature at atm. level (K)
+REAL, DIMENSION(klon), INTENT(IN)       :: PQA   ! air humidity at atm. level (kg/kg)
+!REAL, DIMENSION(:), INTENT(IN)       :: PEXNA ! Exner function at atm. level
+!REAL, DIMENSION(:), INTENT(IN)       :: PRHOA ! air density at atm. level
+REAL, DIMENSION(klon), INTENT(IN)       :: PVMOD ! module of wind at atm. wind level (m/s)
+REAL, DIMENSION(klon), INTENT(IN)       :: PZREF ! atm. level for temp. and humidity (m)
+REAL, DIMENSION(klon), INTENT(IN)       :: PUREF ! atm. level for wind (m)
+REAL, DIMENSION(klon), INTENT(IN)       :: PSST  ! Sea Surface Temperature (K)
+!REAL, DIMENSION(:), INTENT(IN)       :: PEXNS ! Exner function at sea surface
+REAL, DIMENSION(klon), INTENT(IN)       :: PPS   ! air pressure at sea surface (Pa)
+REAL, DIMENSION(klon), INTENT(IN)       :: PRAIN !precipitation rate (kg/s/m2)
+REAL, DIMENSION(klon), INTENT(IN)       :: PPA        ! air pressure at atm level (Pa)
+REAL, DIMENSION(klon), INTENT(IN)       :: PQSATA        ! air pressure at atm level (Pa)
+!
+REAL, DIMENSION(klon), INTENT(INOUT)    :: PZ0SEA! roughness length over the ocean
+!                                                                                 
+!  surface fluxes : latent heat, sensible heat, friction fluxes
+REAL, DIMENSION(klon), INTENT(OUT)      :: PSFTH ! heat flux (W/m2)
+REAL, DIMENSION(klon), INTENT(OUT)      :: PSFTQ ! water flux (kg/m2/s)
+REAL, DIMENSION(klon), INTENT(OUT)      :: PUSTAR! friction velocity (m/s)
+!
+! diagnostics
+REAL, DIMENSION(klon), INTENT(OUT)      :: PQSAT ! humidity at saturation
+REAL, DIMENSION(klon), INTENT(OUT)      :: PCD   ! heat drag coefficient
+REAL, DIMENSION(klon), INTENT(OUT)      :: PCDN  ! momentum drag coefficient
+REAL, DIMENSION(klon), INTENT(OUT)      :: PCH   ! neutral momentum drag coefficient
+REAL, DIMENSION(klon), INTENT(OUT)      :: PCE  !transfer coef. for latent heat flux
+REAL, DIMENSION(klon), INTENT(OUT)      :: PRI   ! Richardson number
+REAL, DIMENSION(klon), INTENT(OUT)      :: PRESA ! aerodynamical resistance
+REAL, DIMENSION(klon), INTENT(OUT)      :: PZ0HSEA ! heat roughness length
+
+LOGICAL,            INTENT(IN)    :: LPRECIP   !
+LOGICAL,            INTENT(IN)    :: LPWG    ! 
+real, dimension(3), intent(inout)      :: coeffs
+!
+!
+
+!INCLUDE "YOMCST.h"
+!INCLUDE "clesphys.h"
+
+!*      0.2    declarations of local variables
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZVMOD    ! wind intensity
+REAL, DIMENSION(SIZE(PTA))      :: ZPA      ! Pressure at atm. level
+REAL, DIMENSION(SIZE(PTA))      :: ZTA      ! Temperature at atm. level
+REAL, DIMENSION(SIZE(PTA))      :: ZQASAT   ! specific humidity at saturation  at atm. level (kg/kg)
+!
+!rajout
+REAL, DIMENSION(SIZE(PTA))       :: PEXNA      ! Exner function at atm level
+REAL, DIMENSION(SIZE(PTA))       :: PEXNS      ! Exner function at atm level
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZO       ! rougness length ref 
+REAL, DIMENSION(SIZE(PTA))      :: ZWG      ! gustiness factor (m/s)
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZDU,ZDT,ZDQ,ZDUWG !differences
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZUSR        !velocity scaling parameter "ustar" (m/s) = friction velocity
+REAL, DIMENSION(SIZE(PTA))      :: ZTSR        !temperature sacling parameter "tstar" (degC)
+REAL, DIMENSION(SIZE(PTA))      :: ZQSR        !humidity scaling parameter "qstar" (kg/kg)
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZU10,ZT10   !vertical profils (10-m height) 
+REAL, DIMENSION(SIZE(PTA))      :: ZVISA       !kinematic viscosity of dry air
+REAL, DIMENSION(SIZE(PTA))      :: ZO10,ZOT10  !roughness length at 10m
+REAL, DIMENSION(SIZE(PTA))      :: ZCD,ZCT,ZCC
+REAL, DIMENSION(SIZE(PTA))      :: ZCD10,ZCT10 !transfer coef. at 10m
+REAL, DIMENSION(SIZE(PTA))      :: ZRIBU,ZRIBCU
+REAL, DIMENSION(SIZE(PTA))      :: ZETU,ZL10
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZCHARN                      !Charnock number depends on wind module
+REAL, DIMENSION(SIZE(PTA))      :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZZL,ZZTL!,ZZQL    !Obukhovs stability 
+                                                     !param. z/l for u,T,q
+REAL, DIMENSION(SIZE(PTA))      :: ZRR
+REAL, DIMENSION(SIZE(PTA))      :: ZOT,ZOQ           !rougness length ref
+REAL, DIMENSION(SIZE(PTA))      :: ZPUZ,ZPTZ,ZPQZ    !PHI funct. for u,T,q 
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZBF               !constants to compute gustiness factor
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZTAU       !momentum flux (W/m2)
+REAL, DIMENSION(SIZE(PTA))      :: ZHF        !sensible heat flux (W/m2)
+REAL, DIMENSION(SIZE(PTA))      :: ZEF        !latent heat flux (W/m2)
+REAL, DIMENSION(SIZE(PTA))      :: ZWBAR      !diag for webb correction but not used here after
+REAL, DIMENSION(SIZE(PTA))      :: ZTAUR      !momentum flux due to rain (W/m2)
+REAL, DIMENSION(SIZE(PTA))      :: ZRF        !sensible heat flux due to rain (W/m2)
+REAL, DIMENSION(SIZE(PTA))      :: ZCHN,ZCEN  !neutral coef. for heat and vapor
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZLV      !latent heat constant
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact
+REAL, DIMENSION(SIZE(PTA))      :: ZXLR                           ! vaporisation  heat  at a given temperature
+REAL, DIMENSION(SIZE(PTA))      :: ZCPLW                          ! specific heat for water at a given temperature 
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZUSTAR2  ! square of friction velocity
+!
+REAL, DIMENSION(SIZE(PTA))      :: ZDIRCOSZW! orography slope cosine (=1 on water!)
+REAL, DIMENSION(SIZE(PTA))      :: ZAC      ! Aerodynamical conductance
+!
+!
+INTEGER, DIMENSION(SIZE(PTA))   :: ITERMAX             ! maximum number of iterations
+!
+REAL    :: ZRVSRDM1,ZRDSRV,ZR2 ! thermodynamic constants
+REAL    :: ZBETAGUST           !gustiness factor
+REAL    :: ZZBL                !atm. boundary layer depth (m)
+REAL    :: ZVISW               !m2/s kinematic viscosity of water
+REAL    :: ZS                  !height of rougness length ref
+REAL    :: ZCH10               !transfer coef. at 10m
+
+REAL    :: QSAT_SEAWATER
+REAL    :: QSATSEAW_1D
+REAL, EXTERNAL :: PSIFCTU, PSIFCTT
+!
+INTEGER :: J, JLOOP    !loop indice
+!REAL(KIND=JPRB) :: ZHOOK_HANDLE
+
+!--------- Modif Olive -----------------
+REAL, DIMENSION(SIZE(PTA))        :: PRHOA
+REAL, PARAMETER                   :: XUNDEF = 1.E+20
+
+REAL       :: XVCHRNK = 0.021
+REAL       :: XVZ0CM = 1.0E-5 
+!REAL       :: XRIMAX
+
+
+INTEGER :: PREF             ! reference pressure for exner function
+INTEGER :: NGRVWAVES        ! Pour le choix du z0
+
+INCLUDE "YOMCST.h"
+INCLUDE "clesphys.h"
+
+!--------------------------------------
+
+
+PRHOA(:) = PPS(:) / (287.1 * PTA(:) * (1.+.61*PQA(:)))
+
+PREF = 100000.                     ! = 1000 hPa
+NGRVWAVES = 1
+
+PEXNA = (PPA/PREF)**(RD/RCPD)
+PEXNS = (PPS/PREF)**(RD/RCPD)
+
+!
+!-------------------------------------------------------------------------------
+!
+!       1.     Initializations
+!              ---------------
+!
+!       1.1   Constants and parameters
+!
+!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',0,ZHOOK_HANDLE)
+!
+ZRVSRDM1  = XRV/XRD-1. ! 0.607766
+ZRDSRV    = XRD/XRV    ! 0.62198
+ZR2       = 1.-ZRDSRV  ! pas utilisé dans cette routine
+ZBETAGUST = 1.2        ! value based on TOGA-COARE experiment
+ZZBL      = 600.       ! Set a default value for boundary layer depth
+ZS        = 10.        ! Standard heigth =10m
+ZCH10     = 0.00115
+!
+ZVISW     = 1.E-6
+!
+!       1.2   Array initialization by undefined values
+!
+PSFTH (:)=XUNDEF
+PSFTQ (:)=XUNDEF
+PUSTAR(:)=XUNDEF
+!
+PCD(:) = XUNDEF
+PCDN(:) = XUNDEF
+PCH(:) = XUNDEF
+PCE(:) =XUNDEF
+PRI(:) = XUNDEF
+!
+PRESA(:)=XUNDEF
+!
+!-------------------------------------------------------------------------------
+!       2. INITIAL GUESS FOR THE ITERATIVE METHOD 
+!          -------------------------------------
+!
+!       2.0     Temperature 
+!
+! Set a non-zero value for the temperature gradient
+!
+WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.) 
+      ZTA(:)=PTA(:)-1E-3
+ELSEWHERE
+      ZTA(:)=PTA(:)      
+ENDWHERE
+
+!       2.1     Wind and humidity 
+!
+! Sea surface specific humidity 
+!
+!PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:))  
+PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:)) 
+
+
+       
+!              
+! Set a minimum value to wind 
+!
+!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))
+
+
+ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) )    !set a minimum value to wind
+!ZVMOD = PVMOD    !set a minimum value to wind
+
+!
+! Specific humidity at saturation at the atm. level 
+!
+ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD))
+!ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:)) 
+ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:)) 
+
+
+!
+!
+ZO(:)  = 0.0001
+ZWG(:) = 0.
+IF (LPWG) ZWG(:) = 0.5
+!
+ZCHARN(:) = 0.011  
+!
+DO J=1,SIZE(PTA)
+  !
+  !      2.2       initial guess
+  !    
+  ZDU(J) = ZVMOD(J)   !wind speed difference with surface current(=0) (m/s)
+                      !initial guess for gustiness factor
+  ZDT(J) = -(ZTA(J)/PEXNA(J)) + (PSST(J)/PEXNS(J)) !potential temperature difference
+  ZDQ(J) = PQSAT(J)-PQA(J)                         !specific humidity difference
+  !
+  ZDUWG(J) = SQRT(ZDU(J)**2+ZWG(J)**2)     !wind speed difference including gustiness ZWG
+  !
+  !      2.3   initialization of neutral coefficients
+  !
+  ZU10(J)  = ZDUWG(J)*LOG(ZS/ZO(J))/LOG(PUREF(J)/ZO(J))
+  ZUSR(J)  = 0.035*ZU10(J)
+  ZVISA(J) = 1.326E-5*(1.+6.542E-3*(ZTA(J)-XTT)+&
+             8.301E-6*(ZTA(J)-XTT)**2-4.84E-9*(ZTA(J)-XTT)**3) !Andrea (1989) CRREL Rep. 89-11
+  ! 
+  ZO10(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J)
+  ZCD(J)  = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2  !drag coefficient
+  ZCD10(J)= (XKARMAN/LOG(ZS/ZO10(J)))**2
+  ZCT10(J)= ZCH10/SQRT(ZCD10(J))
+  ZOT10(J)= ZS/EXP(XKARMAN/ZCT10(J))
+  !
+  !-------------------------------------------------------------------------------
+  !             Grachev and Fairall (JAM, 1997)
+  ZCT(J) = XKARMAN/LOG(PZREF(J)/ZOT10(J))      !temperature transfer coefficient
+  ZCC(J) = XKARMAN*ZCT(J)/ZCD(J)               !z/L vs Rib linear coef.
+  !
+  ZRIBCU(J) = -PUREF(J)/(ZZBL*0.004*ZBETAGUST**3) !saturation or plateau Rib
+  !ZRIBU(J) =-XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*(ZTA(J)-XTT)*ZDQ)/&
+  !     &((ZTA(J)-XTT)*ZDUWG(J)**2)
+  ZRIBU(J)  = -XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*ZTA(J)*ZDQ(J))/&
+               (ZTA(J)*ZDUWG(J)**2)  
+  !
+  IF (ZRIBU(J)<0.) THEN
+    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+ZRIBU(J)/ZRIBCU(J))    !Unstable G and F
+  ELSE
+    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable 
+  ENDIF
+  !
+  ZL10(J) = PUREF(J)/ZETU(J) !MO length
+  !
+ENDDO
+!
+!  First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL)
+ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF/ZL10))
+ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF/ZL10))
+ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF/ZL10))
+!
+ZZL(:) = 0.0
+!
+DO J=1,SIZE(PTA)
+  !
+  IF (ZETU(J)>50.) THEN
+    ITERMAX(J) = 1
+  ELSE
+    ITERMAX(J) = 3 !number of iterations
+  ENDIF
+  !
+  !then modify Charnork for high wind speeds Chris Fairall's data
+  IF (ZDUWG(J)>10.) ZCHARN(J) = 0.011 + (0.018-0.011)*(ZDUWG(J)-10.)/(18.-10.)
+  IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018
+  !
+  !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR 
+  !                -------------------------------------------
+  !
+  ZHWAVE(J) = 0.018*ZVMOD(J)*ZVMOD(J)*(1.+0.015*ZVMOD(J))
+  ZTWAVE(J) = 0.729*ZVMOD(J)
+  ZCWAVE(J) = XG*ZTWAVE(J)/(2.*XPI)
+  ZLWAVE(J) = ZTWAVE(J)*ZCWAVE(J)
+  !
+ENDDO
+!
+   
+!
+DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop
+  !
+  DO J=1,SIZE(PTA)
+    !
+    IF (JLOOP.GT.ITERMAX(J)) CYCLE
+    !
+    IF (NGRVWAVES==0) THEN
+      ZO(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG + 0.11*ZVISA(J)/ZUSR(J) !Smith 1988
+    ELSE IF (NGRVWAVES==1) THEN
+      ZO(J) = (50./(2.*XPI))*ZLWAVE(J)*(ZUSR(J)/ZCWAVE(J))**4.5 &
+              + 0.11*ZVISA(J)/ZUSR(J)                       !Oost et al. 2002  
+    ELSE IF (NGRVWAVES==2) THEN
+      ZO(J) = 1200.*ZHWAVE(J)*(ZHWAVE(J)/ZLWAVE(J))**4.5 &
+              + 0.11*ZVISA(J)/ZUSR(J)                       !Taulor and Yelland 2001  
+    ENDIF
+    !
+    ZRR(J) = ZO(J)*ZUSR(J)/ZVISA(J)
+    ZOQ(J) = MIN(1.15E-4 , 5.5E-5/ZRR(J)**0.6)
+    ZOT(J) = ZOQ(J)
+    !
+    ZZL(J) = XKARMAN * XG * PUREF(J) * &
+              ( ZTSR(J)*(1.+ZRVSRDM1*PQA(J)) + ZRVSRDM1*ZTA(J)*ZQSR(J) ) / &
+              ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) )  
+    ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J)  ! for T 
+!    ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J)  ! for Q
+  ENDDO
+  !
+  ZPUZ(:) = PSIFCTU(ZZL(:))     
+  ZPTZ(:) = PSIFCTT(ZZTL(:))
+  !
+  DO J=1,SIZE(PTA)
+    !
+    ! ZPQZ(J)=PSIFCTT(ZZQL(J))    
+    ZPQZ(J) = ZPTZ(J)
+    !
+    !             3.1 scale parameters
+    !
+    ZUSR(J) = ZDUWG(J)*XKARMAN/(LOG(PUREF(J)/ZO(J)) -ZPUZ(J))
+    ZTSR(J) = -ZDT(J)  *XKARMAN/(LOG(PZREF(J)/ZOT(J))-ZPTZ(J))
+    ZQSR(J) = -ZDQ(J)  *XKARMAN/(LOG(PZREF(J)/ZOQ(J))-ZPQZ(J))
+    !
+    !             3.2 Gustiness factor (ZWG)
+    !
+    IF(LPWG) THEN
+      ZBF(J) = -XG/ZTA(J)*ZUSR(J)*(ZTSR(J)+ZRVSRDM1*ZTA(J)*ZQSR(J))
+      IF (ZBF(J)>0.) THEN
+        ZWG(J) = ZBETAGUST*(ZBF(J)*ZZBL)**(1./3.)
+      ELSE
+        ZWG(J) = 0.2
+      ENDIF
+    ENDIF  
+    ZDUWG(J) = SQRT(ZVMOD(J)**2 + ZWG(J)**2)
+    !
+  ENDDO
+  !
+ENDDO
+!-------------------------------------------------------------------------------
+!
+!            4.  COMPUTE transfer coefficients PCD, PCH, ZCE and SURFACE FLUXES
+!                --------------------------------------------------------------
+!
+ZTAU(:) = XUNDEF
+ZHF(:)  = XUNDEF
+ZEF(:)  = XUNDEF
+!
+ZWBAR(:) = 0.
+ZTAUR(:) = 0.
+ZRF(:)   = 0.
+!
+DO J=1,SIZE(PTA)
+  !
+  !
+  !            4. transfert coefficients PCD, PCH and PCE 
+  !                 and neutral PCDN, ZCHN, ZCEN 
+  !
+  PCD(J) = (ZUSR(J)/ZDUWG(J))**2.
+  PCH(J) = ZUSR(J)*ZTSR(J)/(ZDUWG(J)*(ZTA(J)*PEXNS(J)/PEXNA(J)-PSST(J)))
+  PCE(J) = ZUSR(J)*ZQSR(J)/(ZDUWG(J)*(PQA(J)-PQSAT(J)))
+  !
+  PCDN(J) = (XKARMAN/LOG(ZS/ZO(J)))**2.
+  ZCHN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOT(J)))
+  ZCEN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOQ(J)))
+  !
+  ZLV(J) = XLVTT + (XCPV-XCL)*(PSST(J)-XTT)
+  !
+  !            4. 2 surface fluxes 
+  !
+  IF (ABS(PCDN(J))>1.E-2) THEN   !!!! secure COARE3.0 CODE 
+    write(*,*) 'pb PCDN in COARE30: ',PCDN(J)
+    write(*,*) 'point: ',J,"/",SIZE(PTA)
+    write(*,*) 'roughness: ', ZO(J)
+    write(*,*) 'ustar: ',ZUSR(J)
+    write(*,*) 'wind: ',ZDUWG(J)
+    CALL abort_physic('COARE30',': PCDN too large -> no convergence',1)
+  ELSE
+    ZTSR(J) = -ZTSR(J)
+    ZQSR(J) = -ZQSR(J)
+    ZTAU(J) = -PRHOA(J)*ZUSR(J)*ZUSR(J)*ZVMOD(J)/ZDUWG(J)
+    ZHF(J)  =  PRHOA(J)*XCPD*ZUSR(J)*ZTSR(J)
+    ZEF(J)  =  PRHOA(J)*ZLV(J)*ZUSR(J)*ZQSR(J)
+    !    
+    !           4.3 Contributions to surface  fluxes due to rainfall
+    !
+    ! SB: a priori, le facteur ZRDSRV=XRD/XRV est introduit pour
+    !     adapter la formule de Clausius-Clapeyron (pour l'air
+    !     sec) au cas humide.
+    IF (LPRECIP) THEN
+      ! 
+      ! heat surface  fluxes
+      !
+      ZTAC(J)  = ZTA(J)-XTT
+      !
+      ZXLR(J)  = XLVTT + (XCPV-XCL)* ZTAC(J)                            ! latent heat of rain vaporization
+      ZDQSDT(J)= ZQASAT(J) * ZXLR(J) / (XRD*ZTA(J)**2)                  ! Clausius-Clapeyron relation
+      ZDTMP(J) = (1.0 + 3.309e-3*ZTAC(J) -1.44e-6*ZTAC(J)*ZTAC(J)) * &  !heat diffusivity
+                  0.02411 / (PRHOA(J)*XCPD)
+      !
+      ZDWAT(J) = 2.11e-5 * (XP00/ZPA(J)) * (ZTA(J)/XTT)**1.94           ! water vapour diffusivity from eq (13.3)
+      !                                                                 ! of Pruppacher and Klett (1978)      
+      ZALFAC(J)= 1.0 / (1.0 + &                                         ! Eq.11 in GoF95
+                   ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD))   ! ZALFAC=wet-bulb factor (sans dim)     
+      ZCPLW(J) = 4224.8482 + ZTAC(J) * &
+                              ( -4.707 + ZTAC(J) * &
+                                (0.08499 + ZTAC(J) * &
+                                  (1.2826e-3 + ZTAC(J) * &
+                                    (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat  
+      !       
+      ZRF(J)   = PRAIN(J) * ZCPLW(J) * ZALFAC(J) * &                    !Eq.12 in GoF95 !SIGNE?
+                   (PSST(J) - ZTA(J) + (PQSAT(J)-PQA(J))*ZXLR(J)/XCPD )
+      !
+      ! Momentum flux due to rainfall  
+      !
+      ZTAUR(J)=-0.85*(PRAIN(J) *ZVMOD(J)) !pp3752 in FBR96
+      !
+    ENDIF
+    !
+    !             4.4   Webb correction to latent heat flux
+    ! 
+    ZWBAR(J)=- (1./ZRDSRV)*ZUSR(J)*ZQSR(J) / (1.0+(1./ZRDSRV)*PQA(J)) &
+               - ZUSR(J)*ZTSR(J)/ZTA(J)                        ! Eq.21*rhoa in FBR96    
+    !
+    !             4.5   friction velocity which contains correction du to rain            
+    !
+    ZUSTAR2(J)= - (ZTAU(J) + ZTAUR(J)) / PRHOA(J)
+    PUSTAR(J) =  SQRT(ZUSTAR2(J))
+    !
+    !             4.6   Total surface fluxes
+    !           
+    PSFTH (J) =  ZHF(J) + ZRF(J)
+    PSFTQ (J) =  ZEF(J) / ZLV(J)
+    ! 
+  ENDIF
+ENDDO 
+
+
+coeffs = [PCD,&
+       PCE,&
+       PCH]
+                     
+!-------------------------------------------------------------------------------
+!
+!       5.  FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS 
+!           -----------
+!       5.1    Richardson number
+!             
+!
+!------------STOP LA --------------------
+!ZDIRCOSZW(:) = 1.
+! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,ZTA,ZQASAT,&
+!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI   )  
+!!
+!!       5.2     Aerodynamical conductance and resistance
+!!             
+!ZAC(:) = PCH(:)*ZVMOD(:)
+!PRESA(:) = 1. / MAX(ZAC(:),XSURF_EPSILON)
+!
+!!       5.3 Z0 and Z0H over sea
+!!
+!PZ0SEA(:) =  ZCHARN(:) * ZUSTAR2(:) / XG + XVZ0CM * PCD(:) / PCDN(:)
+!!
+!!PZ0HSEA(:) = PZ0SEA(:)
+!!
+!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',1,ZHOOK_HANDLE)
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE COARE30_FLUX_CNRM
Index: /LMDZ6/trunk/libf/phylmd/coare_cp.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/coare_cp.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/coare_cp.F90	(revision 4722)
@@ -0,0 +1,252 @@
+real function psit_30(zet)
+
+  real, intent(in) :: zet
+
+  if(zet<0) then
+     x=(1.-(15*zet))**.5 
+     psik=2*log((1+x)/2) 
+     x=(1.-(34.15*zet))**.3333 
+     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 
+     f=zet*zet/(1+zet*zet) 
+     psit_30=(1-f)*psik+f*psic   
+
+  else
+     c=min(50.,.35*zet) 
+     psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
+  endif
+end FUNCTION psit_30
+
+real function psiuo(zet)
+
+  real, intent(in) :: zet
+
+  if (zet<0) then 
+
+     x=(1.-15.*zet)**.25 
+     psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) 
+     x=(1.-10.15*zet)**.3333 
+     psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 
+     f=zet*zet/(1+zet*zet) 
+     psiuo=(1-f)*psik+f*psic                                                
+  else
+     c=min(50.,.35*zet) 
+     psiuo=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
+  endif
+END FUNCTION psiuo
+
+
+real function gen_q1(ts,p,dq)
+
+  real, intent(in) :: ts,p,dq
+  real es,qs
+
+  es = 6.112*exp(17.502*(ts-273.15)/(ts-32.18))*.98*(1.0007+3.46e-8*p)
+  qs = es*62.197/(p-37.8*es)
+
+  q1 = dq + qs  
+
+end function gen_q1
+
+
+
+
+subroutine coare_cp(du,dt,dq,&
+     t,q,&
+     zu,zt,zq,&
+     p,zi,&
+     nits,&
+     coeffs,rugosm,rugosh)
+  !     zo,tau,hsb,hlb,&
+  !     var)
+  !version with shortened iteration  modified Rt and Rq
+  !uses wave information wave period in s and wave ht in m
+  !no wave, standard coare 2.6 charnock:  jwave=0 
+  !Oost et al.  zo=50/2/pi L (u*/c)**4.5 if jwave=1
+  !taylor and yelland  zo=1200 h*(L/h)**4.5 jwave=2
+  !
+  USE MODD_CSTS,       ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
+                            XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, &
+                            XP00
+
+  IMPLICIT NONE
+
+  real, intent(in) :: du,dt,dq,t,q
+  real, intent(in) :: zu,zt,zq,p,zi
+  integer, intent(in) :: nits
+  ! real, dimension (nits), intent(out) :: zo,tau,hsb,hlb
+  ! real, dimension(nits,3), intent(out) :: var
+  real, dimension(3), intent(out) :: coeffs
+  real, intent(out) :: rugosm
+  real, intent(out) :: rugosh
+
+  real, parameter ::  beta=1.2, von=.4, fdg = 1. ,&
+       tdk = 273.16, pi = 3.141593, grav = 9.82,&
+       rgas = 287.1
+
+  integer, dimension(3) :: shape_input
+
+  real bf, cc, rhoa, visa,&
+       u10, ut, uts, ut0, ug,&
+       cd10, ch10, ct, ct10,&
+       charn, ribu,&
+       l, l10, &
+       ribcu,rr,&
+       zet,zetu, zom10, zoh10, zot, zoq, zot10,&
+       cd, ch, le, cpa, cpv,&
+       usr,qsr,tsr,&
+       zom, t_c
+
+!, ZTWAVE, ZHWAVE, ZCWAVE, ZLWAVE
+
+
+
+
+  real old_usr, old_tsr, old_qsr,tmp
+
+  real, external :: psit_30, psiuo, grv
+  integer i,j,k
+!---------------- Rajout pour prendre en compte différent Z0 --------------------------------!
+!  INTEGER :: NGRVWAVES        ! Pour le choix du z0
+!  NGRVWAVES = 2
+!---------------------------------------------------------------------------------------------
+!-------------------- Attention Modif réalisée pas SURRR -------------------------------------
+!---------------------------------------------------------------------------------------------
+
+  Ribcu=-zu/zi/.004/Beta**3 
+  cpa=1004.67 
+
+  t_c = t - 273.15
+
+
+  Le=(2.501-.00237*(t_c-dt))*1e6 
+
+
+  visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) 
+
+
+  cpv=cpa*(1+0.84*Q) 
+  rhoa=P/(Rgas*t*(1+0.61*Q)) 
+
+
+
+
+
+
+  ug=.2 
+
+  ut=sqrt(du*du+ug*ug) 
+  ut=MAX(ut , 0.1 * MIN(10.,zu) )
+
+  u10=ut*log(10/1e-4)/log(zu/1e-4) 
+  usr=.035*u10                              ! turbulent friction velocity (m/s), including gustiness
+  zom10=0.011*usr*usr/grav+0.11*visa/usr     ! roughness length for u (smith 88)
+  Cd10=(von/log(10/zom10))**2 
+!  zoh10=0.40*visa/usr+1.4e-5
+!  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca 
+  Ch10=0.00115 
+  Ct10=Ch10/sqrt(Cd10) 
+  zot10=10/exp(von/Ct10)                    ! roughness length for t
+  Cd=(von/log(zu/zom10))**2 
+  Ct=von/log(zt/zot10) 
+  CC=von*Ct/Cd 
+
+  ut0 = ut
+
+
+  ut = ut0
+  Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 
+
+  if (Ribu .LT. 0) then 
+     zetu=CC*Ribu/(1+Ribu/Ribcu) 
+  else 
+     zetu=CC*Ribu*(1+27/9*Ribu/CC) 
+  endif
+
+  L10=zu/zetu 
+
+  ! if (zetu .GT. 50) then 
+  !    nits=1 
+  ! endif
+
+  usr=ut*von/(log(zu/zom10)-psiuo(zetu))
+  tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) 
+  qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) 
+  !  tkt=.001
+
+  ! charnock constant - lin par morceau - constant
+  if ( ut <= 10. ) then
+     charn=0.011 
+  else
+     if (ut .GT. 18) then
+        charn=0.018
+     else
+        charn=0.011+(ut-10)/(18-10)*(0.018-0.011) 
+     endif
+  endif
+
+!  ZHWAVE = 0.018*ut*ut*(1.+0.015*ut)
+!  ZTWAVE = 0.729*ut
+!  ZCWAVE = XG*ZTWAVE/(2.*pi)
+!  ZLWAVE = ZTWAVE*ZCWAVE
+
+
+  !***************  bulk loop ************
+  do i=1, nits 
+
+     zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q)
+
+!     IF (NGRVWAVES==0) THEN
+!       zom = charn*usr*usr/XG + 0.11*visa/usr !Smith 1988
+!     ELSE IF (NGRVWAVES==1) THEN
+!       zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 &
+!               + 0.11*visa/usr                       !Oost et al. 2002  
+!     ELSE IF (NGRVWAVES==2) THEN
+!       zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 &
+!               + 0.11*visa/usr                       !Taulor and Yelland 2001  
+!     ENDIF 
+
+     zom=charn*usr*usr/grav+0.11*visa/usr 
+
+     rr=zom*usr/visa 
+     L=zu/zet 
+     zoq=min(1.15e-4,5.5e-5/rr**.6)  ! a modifier
+     zot=zoq                         ! a modifier
+!     zot=0.40*visa/usr+1.4e-5 
+
+     old_usr = usr
+     old_tsr = tsr
+     old_qsr = qsr
+
+     usr=ut*von/(log(zu/zom)-psiuo(zu/L)) 
+     tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) 
+     qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) 
+
+     Bf=-grav/t*usr*(tsr+.61*t*qsr) 
+
+     if (Bf .GT. 0) then 
+        ug=Beta*(Bf*zi)**.333 
+     else 
+        ug=.2 
+     endif
+
+     ut=sqrt(du*du+ug*ug) 
+     ut=MAX(ut , 0.1 * MIN(10.,zu) )
+
+
+  enddo !bulk iter loop
+
+  ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress
+  ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr 
+  ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr 
+
+  tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du
+
+  coeffs = [tmp**2,&
+       von*fdg/(log(zt/zot)-psit_30(zt/L)) * tmp,&
+       von*fdg/(log(zq/zoq)-psit_30(zq/L)) * tmp]
+
+  rugosm = zom
+  rugosh = zot
+
+
+end subroutine coare_cp
Index: /LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4721)
+++ /LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4722)
@@ -243,4 +243,8 @@
     LOGICAL, SAVE :: ok_new_lscp_omp
     LOGICAL, SAVE :: ok_icefra_lscp_omp
+    !rajout de choix_bulk et nit_bulk par Olivier Torres
+    INTEGER,SAVE  :: choix_bulk_omp
+    INTEGER,SAVE  :: nit_bulk_omp
+    INTEGER,SAVE  :: kz0_omp
     LOGICAL, SAVE :: ok_bs_omp, ok_rad_bs_omp
 
@@ -936,4 +940,31 @@
     nbapp_rad_omp = 12
     CALL getin('nbapp_rad',nbapp_rad_omp)
+
+    !rajout Olivier Torres
+    !Config  Key  = choix_bulk
+    !Config  Desc = choix de la formulation bulk a prendre dans clcdrag au-dessus de l'ocean
+    !Config  Def  = 0
+    !Config         0 -> originale (lmdz/Louis 79)
+    !Config         1 -> COARE
+    !Config         2 -> CORE-"pure" (cf. Large)
+    !Config         3 -> CORE-"mixte" (avec z_0 et C_T^N donnees par Smith 88)
+    choix_bulk_omp = 0
+    CALL getin('choix_bulk',choix_bulk_omp)
+
+    !Config  Key  = nit_bulk
+    !Config  Desc = choix du nombre d'it de pt fixe dans la bulk
+    !Config  Def  = 5
+    nit_bulk_omp = 1
+    CALL getin('nit_bulk',nit_bulk_omp)
+
+    !Config  Key  = kz0
+    !Config  Desc = choix de la formulation z0 pour la bulk ECUME
+    !Config  Def  = 1
+    !Config         0 -> ARPEGE formulation
+    !Config         1 -> Smith Formulation
+    !Config         2 -> Direct computation using the stability functions
+    kz0_omp = 0
+    CALL getin('kz0',kz0_omp)
+
 
     !Config  Key  = iflag_con
@@ -2496,4 +2527,8 @@
     var_fco2_land_cor = var_fco2_land_cor_omp
     dms_cycle_cpl = dms_cycle_cpl_omp
+    !rajout Olivier Torres
+    kz0=kz0_omp
+    choix_bulk = choix_bulk_omp
+    nit_bulk = nit_bulk_omp
 
     ! Test of coherence between type_ocean and version_ocean
@@ -2847,4 +2882,8 @@
     WRITE(lunout,*) ' buf_sph_pol = ', buf_sph_pol
     WRITE(lunout,*) ' buf_siz_pol= ', buf_siz_pol
+    !rajout Olivier Torres
+    write(lunout,*) 'choix_bulk = ', choix_bulk
+    write(lunout,*) 'nit_bulk = ', nit_bulk
+    write(lunout,*) 'kz0 = ', kz0
 
     !$OMP END MASTER
Index: /LMDZ6/trunk/libf/phylmd/ecumev6_flux.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/ecumev6_flux.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/ecumev6_flux.F90	(revision 4722)
@@ -0,0 +1,777 @@
+!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
+!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!SFX_LIC for details. version 1.
+!     #########
+    SUBROUTINE ECUMEV6_FLUX(PZ0SEA,PTA,PSST,PQA,PQSAT,PVMOD, &
+                            PZREF,PSSS,PUREF,PPS,PPA,OPRECIP,OPWEBB,        &
+                            PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,        &
+                            PRI,PRESA,PRAIN,PZ0HSEA,OPERTFLUX,coeffs   )
+!###############################################################################
+!!
+!!****  *ECUMEV6_FLUX*
+!!
+!!    PURPOSE
+!!    -------
+!       Calculate the surface turbulent fluxes of heat, moisture, and momentum 
+!       over sea surface + corrections due to rainfall & Webb effect.
+!!
+!!**  METHOD
+!!    ------
+!       The estimation of the transfer coefficients relies on the iterative 
+!       computation of the scaling parameters U*/Teta*/q*. The convergence is
+!       supposed to be reached in NITERFL iterations maximum.
+!       Neutral transfer coefficients for momentum/temperature/humidity
+!       are computed as a function of the 10m-height neutral wind speed using
+!       the ECUME_V6 formulation based on the multi-campaign (POMME,FETCH,CATCH,
+!       SEMAPHORE,EQUALANT) ALBATROS dataset.
+!!
+!!    EXTERNAL
+!!    --------
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!
+!!    REFERENCE
+!!    ---------
+!!      Fairall et al (1996), JGR, 3747-3764
+!!      Gosnell et al (1995), JGR, 437-442
+!!      Fairall et al (1996), JGR, 1295-1308
+!!
+!!    AUTHOR
+!!    ------
+!!      C. Lebeaupin  *Météo-France* (adapted from S. Belamari's code)
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original     15/03/2005
+!!      Modified        01/2006  C. Lebeaupin (adapted from  A. Pirani's code)
+!!      Modified        08/2009  B. Decharme: limitation of Ri
+!!      Modified        09/2012  B. Decharme: CD correction
+!!      Modified        09/2012  B. Decharme: limitation of Ri in surface_ri.F90
+!!      Modified        10/2012  P. Le Moigne: extra inputs for FLake use
+!!      Modified        06/2013  B. Decharme: bug in z0 (output) computation 
+!!      Modified        12/2013  S. Belamari: ZRF computation updated:
+!!                                1. ZP00/PPA in ZDWAT, ZLVA in ZDQSDT/ZBULB/ZRF
+!!                                2. ZDWAT/ZDTMP in ZBULB/ZRF (Gosnell et al 95)
+!!                                3. cool skin correction included
+!!      Modified        01/2014  S. Belamari: salinity impact on latent heat of
+!!                                vaporization of seawater included
+!!      Modified        01/2014  S. Belamari: new formulation for pure water
+!!                                specific heat (ZCPWA)
+!!      Modified        01/2014  S. Belamari: 4 choices for PZ0SEA computation
+!!      Modified        12/2015  S. Belamari: ECUME now provides parameterisations
+!!                                for:  U10n*sqrt(CDN)          instead of CDN
+!!                                      U10n*CHN/sqrt(CDN)         "       CHN
+!!                                      U10n*CEN/sqrt(CDN)         "       CEN
+!!      Modified        01/2016  S. Belamari: New ECUME formulation
+!!
+!!      To be done:
+!!      include gustiness computation following Mondon & Redelsperger (1998)
+!!!
+!-------------------------------------------------------------------------------
+!!
+!!    MODIFICATIONS RELATED TO SST CORRECTION COMPUTATION
+!!    ---------------------------------------------------
+!!      Modified        09/2013  S. Belamari: use 0.98 for the ocean emissivity
+!!                                following up to date satellite measurements in
+!!                                the 8-14 μm range (obtained values range from
+!!                                0.98 to 0.99).
+!!!
+!-------------------------------------------------------------------------------
+!
+!       0.   DECLARATIONS
+!            ------------
+!
+USE dimphy
+USE indice_sol_mod
+USE MODD_CSTS,             ONLY : XPI, XDAY, XKARMAN, XG, XP00, XSTEFAN, XRD, XRV,   &
+                                  XCPD, XCPV, XCL, XTT, XLVTT
+
+
+!USE MODD_SURF_PAR,         ONLY : XUNDEF
+!USE MODD_SURF_ATM,         ONLY : XVCHRNK, XVZ0CM
+!USE MODD_REPROD_OPER,      ONLY : CCHARNOCK
+!
+!USE MODE_THERMOS
+!USE MODI_WIND_THRESHOLD
+!USE MODI_SURFACE_RI
+!
+!USE YOMHOOK,   ONLY : LHOOK,   DR_HOOK
+!USE PARKIND1,  ONLY : JPRB
+!
+!USE MODI_ABOR1_SFX
+!
+IMPLICIT NONE
+!
+!       0.1. Declarations of arguments
+!
+REAL, DIMENSION(klon), INTENT(IN)    :: PVMOD      ! module of wind at atm level (m/s)
+REAL, DIMENSION(klon), INTENT(IN)    :: PTA        ! air temperature at atm level (K)
+REAL, DIMENSION(klon), INTENT(IN)    :: PQA        ! air spec. hum. at atm level (kg/kg)
+REAL, DIMENSION(klon), INTENT(IN)    :: PQSAT      ! sea surface spec. hum. (kg/kg)
+REAL, DIMENSION(klon), INTENT(IN)    :: PPA        ! air pressure at atm level (Pa)
+!REAL, DIMENSION(:), INTENT(IN)    :: PRHOA      ! air density at atm level (kg/m3)
+!REAL, DIMENSION(:), INTENT(IN)    :: PEXNA      ! Exner function at atm level
+REAL, DIMENSION(klon), INTENT(IN)    :: PUREF      ! atm level for wind (m)
+REAL, DIMENSION(klon), INTENT(IN)    :: PZREF      ! atm level for temp./hum. (m)
+REAL, DIMENSION(klon), INTENT(IN)    :: PSSS       ! Sea Surface Salinity (g/kg)
+REAL, DIMENSION(klon), INTENT(IN)    :: PPS        ! air pressure at sea surface (Pa)
+!REAL, DIMENSION(:), INTENT(IN)    :: PEXNS      ! Exner function at sea surface
+!REAL, DIMENSION(:), INTENT(IN)    :: PPERTFLUX  ! stochastic flux perturbation pattern
+! for correction
+!REAL,               INTENT(IN)    :: PICHCE    !
+LOGICAL,            INTENT(IN)    :: OPRECIP   !
+LOGICAL,            INTENT(IN)    :: OPWEBB    !
+LOGICAL,            INTENT(IN)    :: OPERTFLUX
+REAL, DIMENSION(klon), INTENT(IN)    :: PRAIN     ! precipitation rate (kg/s/m2)
+!
+!INTEGER,            INTENT(IN)    :: KZ0
+!
+REAL, DIMENSION(klon), INTENT(INOUT) :: PSST       ! Sea Surface Temperature (K)
+REAL, DIMENSION(klon), INTENT(INOUT) :: PZ0SEA     ! roughness length over sea
+REAL, DIMENSION(klon), INTENT(OUT)   :: PZ0HSEA    ! heat roughness length over sea
+
+! surface fluxes : latent heat, sensible heat, friction fluxes
+REAL, DIMENSION(klon), INTENT(OUT)   :: PUSTAR     ! friction velocity (m/s)
+REAL, DIMENSION(klon), INTENT(OUT)   :: PSFTH      ! heat flux (W/m2)
+REAL, DIMENSION(klon), INTENT(OUT)   :: PSFTQ      ! water flux (kg/m2/s)
+
+! diagnostics
+REAL, DIMENSION(klon), INTENT(OUT)   :: PCD        ! transfer coef. for momentum
+REAL, DIMENSION(klon), INTENT(OUT)   :: PCH        ! transfer coef. for temperature
+REAL, DIMENSION(klon), INTENT(OUT)   :: PCE        ! transfer coef. for humidity
+REAL, DIMENSION(klon), INTENT(OUT)   :: PCDN       ! neutral coef. for momentum
+REAL, DIMENSION(klon), INTENT(OUT)   :: PRI        ! Richardson number
+REAL, DIMENSION(klon), INTENT(OUT)   :: PRESA      ! aerodynamical resistance
+real, dimension(3), intent(out)   :: coeffs
+
+!       0.2. Declarations of local variables
+!
+! specif SB
+INTEGER, DIMENSION(SIZE(PTA))     :: JCV        ! convergence index
+INTEGER, DIMENSION(SIZE(PTA))     :: JITER      ! nb of iterations to converge
+!rajout
+REAL, DIMENSION(SIZE(PTA))        :: PEXNA      ! Exner function at atm level
+REAL, DIMENSION(SIZE(PTA))       :: PEXNS      ! Exner function at atm level
+!
+REAL, DIMENSION(SIZE(PTA))        :: ZTAU       ! momentum flux (N/m2)
+REAL, DIMENSION(SIZE(PTA))        :: ZHF        ! sensible heat flux (W/m2)
+REAL, DIMENSION(SIZE(PTA))        :: ZEF        ! latent heat flux (W/m2)
+REAL, DIMENSION(SIZE(PTA))        :: ZTAUR      ! momentum flx due to rain (N/m2)
+REAL, DIMENSION(SIZE(PTA))        :: ZRF        ! sensible flx due to rain (W/m2)
+REAL, DIMENSION(SIZE(PTA))        :: ZEFWEBB    ! Webb corr. on latent flx (W/m2)
+
+REAL, DIMENSION(SIZE(PTA))        :: ZVMOD      ! wind intensity at atm level (m/s)
+REAL, DIMENSION(SIZE(PTA))        :: ZQSATA     ! sat.spec.hum. at atm level (kg/kg)
+REAL, DIMENSION(SIZE(PTA))        :: ZLVA       ! vap.heat of pure water at atm level (J/kg)
+REAL, DIMENSION(SIZE(PTA))        :: ZLVS       ! vap.heat of seawater at sea surface (J/kg)
+REAL, DIMENSION(SIZE(PTA))        :: ZCPA       ! specif.heat moist air (J/kg/K)
+REAL, DIMENSION(SIZE(PTA))        :: ZVISA      ! kinemat.visc. of dry air (m2/s)
+REAL, DIMENSION(SIZE(PTA))        :: ZDU        ! U   vert.grad. (real atm)
+REAL, DIMENSION(SIZE(PTA))        :: ZDT,ZDQ    ! T,Q vert.grad. (real atm)
+REAL, DIMENSION(SIZE(PTA))        :: ZDDU       ! U   vert.grad. (real atm + gust)
+REAL, DIMENSION(SIZE(PTA))        :: ZDDT,ZDDQ  ! T,Q vert.grad. (real atm + WL/CS)
+REAL, DIMENSION(SIZE(PTA))        :: ZUSR       ! velocity scaling param. (m/s)
+                                                ! =friction velocity
+REAL, DIMENSION(SIZE(PTA))        :: ZTSR       ! temperature scaling param. (K)
+REAL, DIMENSION(SIZE(PTA))        :: ZQSR       ! humidity scaling param. (kg/kg)
+REAL, DIMENSION(SIZE(PTA))        :: ZDELTAU10N,ZDELTAT10N,ZDELTAQ10N
+                                                ! U,T,Q vert.grad. (10m, neutral atm)
+REAL, DIMENSION(SIZE(PTA))        :: ZUSR0,ZTSR0,ZQSR0    ! ITERATIVE PROCESS
+REAL, DIMENSION(SIZE(PTA))        :: ZDUSTO,ZDTSTO,ZDQSTO ! ITERATIVE PROCESS
+REAL, DIMENSION(SIZE(PTA))        :: ZPSIU,ZPSIT! PSI funct for U, T/Q (Z0 comp)
+REAL, DIMENSION(SIZE(PTA))        :: ZCHARN     ! Charnock parameter   (Z0 comp)
+
+REAL, DIMENSION(SIZE(PTA))        :: ZUSTAR2    ! square of friction velocity
+REAL, DIMENSION(SIZE(PTA))        :: ZAC        ! aerodynamical conductance
+REAL, DIMENSION(SIZE(PTA))        :: ZDIRCOSZW  ! orography slope cosine
+                                                ! (=1 on water!)
+REAL, DIMENSION(SIZE(PTA))        :: ZPARUN,ZPARTN,ZPARQN ! neutral parameter for U,T,Q
+
+!-- rajout pour la pression saturante
+REAL, DIMENSION(SIZE(PPA))        :: ZFOES                                 ! [OPWEBB]
+REAL, DIMENSION(SIZE(PPA))        :: ZWORK1
+REAL, DIMENSION(SIZE(PPA))        :: ZWORK2  
+REAL, DIMENSION(SIZE(PPS))        :: ZWORK1A
+REAL, DIMENSION(SIZE(PPS))        :: ZWORK2A
+!#####################
+
+REAL, DIMENSION(0:5)              :: ZCOEFU,ZCOEFT,ZCOEFQ
+
+!--------- Modif Olive -----------------
+REAL, DIMENSION(SIZE(PTA))        :: PRHOA
+REAL, PARAMETER                   :: XUNDEF = 1.E+20
+
+
+REAL       :: XVCHRNK = 0.021
+REAL       :: XVZ0CM = 1.0E-5 
+!REAL       :: XRIMAX
+
+CHARACTER  :: CCHARNOCK = 'NEW'
+
+
+!--------------------------------------
+
+
+! local constants
+LOGICAL :: OPCVFLX              ! to force convergence
+INTEGER :: NITERMAX             ! nb of iterations to get free convergence
+INTEGER :: NITERSUP             ! nb of additional iterations if OPCVFLX=.TRUE.
+INTEGER :: NITERFL              ! maximum number of iterations
+REAL    :: ZETV,ZRDSRV          ! thermodynamic constants
+REAL    :: ZSQR3
+REAL    :: ZLMOMIN,ZLMOMAX      ! min/max value of Obukhovs stability param. z/l
+REAL    :: ZBTA,ZGMA            ! parameters of the stability functions
+REAL    :: ZDUSR0,ZDTSR0,ZDQSR0 ! maximum gap for USR/TSR/QSR between 2 steps
+REAL    :: ZP00                 ! [OPRECIP] - water vap. diffusiv.ref.press.(Pa)
+REAL    :: ZUTU,ZUTT,ZUTQ       ! U10n threshold in ECUME parameterisation
+REAL    :: ZCDIRU,ZCDIRT,ZCDIRQ ! coef directeur pour fonction affine U,T,Q
+REAL    :: ZORDOU,ZORDOT,ZORDOQ ! ordonnee a l'origine pour fonction affine U,T,Q
+
+INTEGER :: JJ                                   ! for ITERATIVE PROCESS
+INTEGER :: JLON,JK
+REAL    :: ZLMOU,ZLMOT                          ! Obukhovs param. z/l for U, T/Q
+REAL    :: ZPSI_U,ZPSI_T                        ! PSI funct. for U, T/Q
+REAL    :: Z0TSEA,Z0QSEA                        ! roughness length for T, Q
+REAL    :: ZCHIC,ZCHIK,ZPSIC,ZPSIK,ZLOGUS10,ZLOGTS10
+REAL    :: ZTAC,ZCPWA,ZDQSDT,ZDWAT,ZDTMP,ZBULB  ! [OPRECIP]
+REAL    :: ZWW                                  ! [OPWEBB]
+
+
+INTEGER :: PREF             ! reference pressure for exner function
+REAL, DIMENSION(klon)    :: PQSATA      ! sea surface spec. hum. (kg/kg)
+
+REAL    :: qsat_seawater2,qsat_seawater
+
+INCLUDE "YOMCST.h"
+INCLUDE "clesphys.h"
+
+!REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+!-------------------------------------------------------------------------------
+!----------------------- Modif Olive calcul de PRHOA ---------------------------
+
+!write(*,*) "PZ0SEA ",PZ0SEA
+!write(*,*) "PTA ",PTA
+!write(*,*) "PSST ",PSST
+!write(*,*) "PQA ",PQA
+!write(*,*) "PVMOD ",PVMOD
+!write(*,*) "PZREF ",PZREF
+!write(*,*) "PUREF ",PUREF
+!write(*,*) "PPS ",PPS
+!write(*,*) "PPA ",PPA
+!write(*,*) "OPRECIP ",OPRECIP
+!write(*,*) "PZ0HSEA ",PZ0HSEA
+!write(*,*) "PRAIN ",PRAIN
+
+
+PRHOA(:) = PPS(:) / (287.1 * PTA(:) * (1.+.61*PQA(:)))
+!write(*,*) "klon klon ",klon,PTA
+!write(*,*) "PRHOA ",SIZE(PRHOA),PRHOA
+
+PREF = 100900.                    ! = 1000 hPa
+
+!PEXNA = (PPA/PPS)**RKAPPA
+!PEXNS = (PPS/PPS)**RKAPPA
+
+PEXNA = (PPA/PREF)**(RD/RCPD)
+PEXNS = (PPS/PREF)**(RD/RCPD)
+
+!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',0,ZHOOK_HANDLE)
+!
+ZDUSR0   = 1.E-06
+ZDTSR0   = 1.E-06
+ZDQSR0   = 1.E-09
+!
+NITERMAX = 5
+NITERSUP = 5
+OPCVFLX  = .TRUE.
+!
+NITERFL = NITERMAX
+IF (OPCVFLX) NITERFL = NITERMAX+NITERSUP
+!
+ZCOEFU = (/ 1.00E-03, 3.66E-02, -1.92E-03, 2.32E-04, -7.02E-06,  6.40E-08 /)
+ZCOEFT = (/ 5.36E-03, 2.90E-02, -1.24E-03, 4.50E-04, -2.06E-05,       0.0 /)
+ZCOEFQ = (/ 1.00E-03, 3.59E-02, -2.87E-04,      0.0,       0.0,       0.0 /)
+!
+ZUTU = 40.0
+ZUTT = 14.4
+ZUTQ = 10.0
+!
+ZCDIRU = ZCOEFU(1) + 2.0*ZCOEFU(2)*ZUTU + 3.0*ZCOEFU(3)*ZUTU**2   &
+                   + 4.0*ZCOEFU(4)*ZUTU**3 + 5.0*ZCOEFU(5)*ZUTU**4
+ZCDIRT = ZCOEFT(1) + 2.0*ZCOEFT(2)*ZUTT + 3.0*ZCOEFT(3)*ZUTT**2   &
+                   + 4.0*ZCOEFT(4)*ZUTT**3
+ZCDIRQ = ZCOEFQ(1) + 2.0*ZCOEFQ(2)*ZUTQ
+!
+ZORDOU = ZCOEFU(0) + ZCOEFU(1)*ZUTU + ZCOEFU(2)*ZUTU**2 + ZCOEFU(3)*ZUTU**3   &
+                   + ZCOEFU(4)*ZUTU**4 + ZCOEFU(5)*ZUTU**5
+ZORDOT = ZCOEFT(0) + ZCOEFT(1)*ZUTT + ZCOEFT(2)*ZUTT**2 + ZCOEFT(3)*ZUTT**3   &
+                   + ZCOEFT(4)*ZUTT**4
+ZORDOQ = ZCOEFQ(0) + ZCOEFQ(1)*ZUTQ + ZCOEFQ(2)*ZUTQ**2
+!
+!-------------------------------------------------------------------------------
+!
+!       1.   AUXILIARY CONSTANTS & ARRAY INITIALISATION BY UNDEFINED VALUES.
+!       --------------------------------------------------------------------
+!
+ZDIRCOSZW(:) = 1.0
+!
+ZETV    = XRV/XRD-1.0   !~0.61 (cf Liu et al. 1979)
+ZRDSRV  = XRD/XRV       !~0.622
+ZSQR3   = SQRT(3.0)
+ZLMOMIN = -200.0
+ZLMOMAX = 0.25
+ZBTA    = 16.0
+ZGMA    = 7.0           !initially =4.7, modified to 7.0 following G. Caniaux
+!
+ZP00    = 1013.25E+02
+!
+PCD  = XUNDEF
+PCH  = XUNDEF
+PCE  = XUNDEF
+PCDN = XUNDEF
+ZUSR = XUNDEF
+ZTSR = XUNDEF
+ZQSR = XUNDEF
+ZTAU = XUNDEF
+ZHF  = XUNDEF
+ZEF  = XUNDEF
+!
+PSFTH  = XUNDEF
+PSFTQ  = XUNDEF
+PUSTAR = XUNDEF
+PRESA  = XUNDEF
+PRI    = XUNDEF
+!
+ZTAUR   = 0.0
+ZRF     = 0.0
+ZEFWEBB = 0.0
+!
+!-------------------------------------------------------------------------------
+!
+!       2.   INITIALISATIONS BEFORE ITERATIVE LOOP.
+!       -------------------------------------------
+!
+!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))    !set a minimum value to wind
+ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) )    !set a minimum value to wind
+
+write(*,*) "ZVMOD ",SIZE(ZVMOD)
+
+!
+!       2.0. Radiative fluxes - For warm layer & cool skin
+!
+!       2.0b. Warm Layer correction
+!
+!       2.1. Specific humidity at saturation
+!
+WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF)
+  PQSATA = QSAT_SEAWATER2 (PSST(:),PPS(:),PSSS(:))    !at sea surface
+ELSEWHERE
+  PQSATA (:) = QSAT_SEAWATER (PSST(:),PPS(:))            !at sea surface
+ENDWHERE
+
+!ZQSATA(:) = QSAT(PTA(:),PPA(:))                         !at atm level
+
+!### OLIVIER POUR PRESSION SATURANTE #####
+!-------------------------------------------------------------------------------
+!
+ZFOES  = 1 !PSAT(PT(:))
+ZFOES  = 0.98*ZFOES
+!
+ZWORK1    = ZFOES/PPS
+ZWORK2    = XRD/XRV
+
+ZWORK1A    = ZFOES/PPA
+ZWORK2A    = XRD/XRV
+
+!write(*,*) "ZFOES ",ZFOES
+!write(*,*) "PPS ",PPS
+!write(*,*) "ZWORK1 ",ZWORK1
+!write(*,*) "XRD ",XRD
+!write(*,*) "XRV ",XRV
+!write(*,*) "PPA ",PPA
+!write(*,*) "ZWORK1A ",ZWORK1A
+
+write(*,*) "PQSAT : ",PQSAT
+write(*,*) "PQSATA : ",PQSATA
+
+
+!
+!*       2.    COMPUTE SATURATION HUMIDITY
+!              ---------------------------
+!
+!PQSAT  = ZWORK2*ZWORK1 / (1.+(ZWORK2-1.)*ZWORK1)
+!ZQSATA = ZWORK2A*ZWORK1A / (1.+(ZWORK2A-1.)*ZWORK1A)
+ZQSATA = QSAT_SEAWATER (PTA(:),PPA(:))            !at sea surface
+
+
+!
+!       2.2. Gradients at the air-sea interface
+!
+ZDU(:) = ZVMOD(:)               !one assumes u is measured / sea surface current
+ZDT(:) = PTA(:)/PEXNA(:)-PSST(:)/PEXNS(:)
+ZDQ(:) = PQA(:)-PQSATA(:)
+
+write(*,*) "PQA ",PQA(:)
+write(*,*) "PQSAT",PQSAT(:)
+write(*,*) "ZDQ",ZDQ(:)
+!
+!       2.3. Latent heat of vaporisation
+!
+ZLVA(:) = XLVTT+(XCPV-XCL)*(PTA (:)-XTT)                !of pure water at atm level
+ZLVS(:) = XLVTT+(XCPV-XCL)*(PSST(:)-XTT)                !of pure water at sea surface
+
+
+
+write(*,*) "ZLVA ",ZLVA
+write(*,*) "ZLVS ",ZLVS
+
+
+WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF)
+  ZLVS(:) = ZLVS(:)*(1.0-1.00472E-3*PSSS(:))            !of seawater at sea surface
+ENDWHERE
+!
+!       2.4. Specific heat of moist air (Businger 1982)
+!
+!ZCPA(:) = XCPD*(1.0+(XCPV/XCPD-1.0)*PQA(:))
+ZCPA(:) = XCPD
+!
+!       2.4b Kinematic viscosity of dry air (Andreas 1989, CRREL Rep. 89-11)
+!
+ZVISA(:) = 1.326E-05*(1.0+6.542E-03*(PTA(:)-XTT)+8.301E-06*(PTA(:)-XTT)**2   &
+           -4.84E-09*(PTA(:)-XTT)**3)
+!
+!       2.4c Coefficients for warm layer and/or cool skin correction
+!
+!       2.5. Initial guess
+!
+ZDDU(:) = ZDU(:)
+ZDDT(:) = ZDT(:)
+ZDDQ(:) = ZDQ(:)
+ZDDU(:) = SIGN(MAX(ABS(ZDDU(:)),10.0*ZDUSR0),ZDDU(:))
+ZDDT(:) = SIGN(MAX(ABS(ZDDT(:)),10.0*ZDTSR0),ZDDT(:))
+ZDDQ(:) = SIGN(MAX(ABS(ZDDQ(:)),10.0*ZDQSR0),ZDDQ(:))
+
+write(*,*) "ZDDU ",ZDDU
+write(*,*) "ZDDQ ",ZDDQ
+write(*,*) "ZDDT ",ZDDT
+
+!
+JCV (:) = -1
+ZUSR(:) = 0.04*ZDDU(:)
+ZTSR(:) = 0.04*ZDDT(:)
+ZQSR(:) = 0.04*ZDDQ(:)
+ZDELTAU10N(:) = ZDDU(:)
+ZDELTAT10N(:) = ZDDT(:)
+ZDELTAQ10N(:) = ZDDQ(:)
+JITER(:) = 99
+!
+! In the following, we suppose that Richardson number PRI < XRIMAX
+! If not true, Monin-Obukhov theory can't (and therefore shouldn't) be applied !
+!-------------------------------------------------------------------------------
+!
+!       3.   ITERATIVE LOOP TO COMPUTE U*, T*, Q*.
+!       ------------------------------------------
+!
+DO JJ=1,NITERFL
+  DO JLON=1,SIZE(PTA)
+!
+  IF (JCV(JLON) == -1) THEN
+    ZUSR0(JLON)=ZUSR(JLON)
+    ZTSR0(JLON)=ZTSR(JLON)
+    ZQSR0(JLON)=ZQSR(JLON)
+    IF (JJ == NITERMAX+1 .OR. JJ == NITERMAX+NITERSUP) THEN
+      ZDELTAU10N(JLON) = 0.5*(ZDUSTO(JLON)+ZDELTAU10N(JLON))    !forced convergence
+      ZDELTAT10N(JLON) = 0.5*(ZDTSTO(JLON)+ZDELTAT10N(JLON))
+      ZDELTAQ10N(JLON) = 0.5*(ZDQSTO(JLON)+ZDELTAQ10N(JLON))
+      IF (JJ == NITERMAX+NITERSUP) JCV(JLON)=3
+    ENDIF
+    ZDUSTO(JLON) = ZDELTAU10N(JLON)
+    ZDTSTO(JLON) = ZDELTAT10N(JLON)
+    ZDQSTO(JLON) = ZDELTAQ10N(JLON)
+!
+!       3.1. Neutral parameter for wind speed (ECUME_V6 formulation)
+!
+    IF (ZDELTAU10N(JLON) <= ZUTU) THEN
+      ZPARUN(JLON) = ZCOEFU(0) + ZCOEFU(1)*ZDELTAU10N(JLON)      &
+                               + ZCOEFU(2)*ZDELTAU10N(JLON)**2   &
+                               + ZCOEFU(3)*ZDELTAU10N(JLON)**3   &
+                               + ZCOEFU(4)*ZDELTAU10N(JLON)**4   &
+                               + ZCOEFU(5)*ZDELTAU10N(JLON)**5
+    ELSE
+      ZPARUN(JLON) = ZCDIRU*(ZDELTAU10N(JLON)-ZUTU) + ZORDOU
+    ENDIF
+    PCDN(JLON) = (ZPARUN(JLON)/ZDELTAU10N(JLON))**2
+!
+!       3.2. Neutral parameter for temperature (ECUME_V6 formulation)
+!
+    IF (ZDELTAU10N(JLON) <= ZUTT) THEN
+      ZPARTN(JLON) = ZCOEFT(0) + ZCOEFT(1)*ZDELTAU10N(JLON)      &
+                               + ZCOEFT(2)*ZDELTAU10N(JLON)**2   &
+                               + ZCOEFT(3)*ZDELTAU10N(JLON)**3   &
+                               + ZCOEFT(4)*ZDELTAU10N(JLON)**4
+    ELSE
+      ZPARTN(JLON) = ZCDIRT*(ZDELTAU10N(JLON)-ZUTT) + ZORDOT
+    ENDIF
+!
+!       3.3. Neutral parameter for humidity (ECUME_V6 formulation)
+!
+    IF (ZDELTAU10N(JLON) <= ZUTQ) THEN
+      ZPARQN(JLON) = ZCOEFQ(0) + ZCOEFQ(1)*ZDELTAU10N(JLON)      &
+                               + ZCOEFQ(2)*ZDELTAU10N(JLON)**2
+    ELSE
+      ZPARQN(JLON) = ZCDIRQ*(ZDELTAU10N(JLON)-ZUTQ) + ZORDOQ
+    ENDIF
+!
+!       3.4. Scaling parameters U*, T*, Q*
+!
+    ZUSR(JLON) = ZPARUN(JLON)
+    ZTSR(JLON) = ZPARTN(JLON)*ZDELTAT10N(JLON)/ZDELTAU10N(JLON)
+    ZQSR(JLON) = ZPARQN(JLON)*ZDELTAQ10N(JLON)/ZDELTAU10N(JLON)
+!
+!       3.4b Gustiness factor (Deardorff 1970)
+!
+!       3.4c Cool skin correction
+!
+!       3.5. Obukhovs stability param. z/l following Liu et al. (JAS, 1979)
+!
+! For U
+    ZLMOU = PUREF(JLON)*XG*XKARMAN*(ZTSR(JLON)/PTA(JLON)   &
+            +ZETV*ZQSR(JLON)/(1.0+ZETV*PQA(JLON)))/ZUSR(JLON)**2
+! For T/Q
+    ZLMOT = ZLMOU*(PZREF(JLON)/PUREF(JLON))
+    ZLMOU = MAX(MIN(ZLMOU,ZLMOMAX),ZLMOMIN)
+    ZLMOT = MAX(MIN(ZLMOT,ZLMOMAX),ZLMOMIN)
+!
+!       3.6. Stability function psi (see Liu et al, 1979 ; Dyer and Hicks, 1970)
+!            Modified to include convective form following Fairall (unpublished)
+!
+! For U
+    IF (ZLMOU == 0.0) THEN
+      ZPSI_U = 0.0
+    ELSEIF (ZLMOU > 0.0) THEN
+      ZPSI_U = -ZGMA*ZLMOU
+    ELSE
+      ZCHIK  = (1.0-ZBTA*ZLMOU)**0.25
+      ZPSIK  = 2.0*LOG((1.0+ZCHIK)/2.0)  &
+                +LOG((1.0+ZCHIK**2)/2.0) &
+                -2.0*ATAN(ZCHIK)+0.5*XPI
+      ZCHIC  = (1.0-12.87*ZLMOU)**(1.0/3.0)       !for very unstable conditions
+      ZPSIC  = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0)   &
+                -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) &
+                +XPI/ZSQR3
+      ZPSI_U = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOU**2) !match Kansas & free-conv. forms
+    ENDIF
+    ZPSIU(JLON) = ZPSI_U
+! For T/Q
+    IF (ZLMOT == 0.0) THEN
+      ZPSI_T = 0.0
+    ELSEIF (ZLMOT > 0.0) THEN
+      ZPSI_T = -ZGMA*ZLMOT
+    ELSE
+      ZCHIK  = (1.0-ZBTA*ZLMOT)**0.25
+      ZPSIK  = 2.0*LOG((1.0+ZCHIK**2)/2.0)
+      ZCHIC  = (1.0-12.87*ZLMOT)**(1.0/3.0)       !for very unstable conditions
+      ZPSIC  = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0)   &
+                -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) &
+                +XPI/ZSQR3
+      ZPSI_T = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOT**2) !match Kansas & free-conv. forms
+    ENDIF
+    ZPSIT(JLON) = ZPSI_T
+!
+!       3.7. Update air-sea gradients
+!
+    ZDDU(JLON) = ZDU(JLON)
+    ZDDT(JLON) = ZDT(JLON)
+    ZDDQ(JLON) = ZDQ(JLON)
+    ZDDU(JLON) = SIGN(MAX(ABS(ZDDU(JLON)),10.0*ZDUSR0),ZDDU(JLON))
+    ZDDT(JLON) = SIGN(MAX(ABS(ZDDT(JLON)),10.0*ZDTSR0),ZDDT(JLON))
+    ZDDQ(JLON) = SIGN(MAX(ABS(ZDDQ(JLON)),10.0*ZDQSR0),ZDDQ(JLON))
+    ZLOGUS10   = LOG(PUREF(JLON)/10.0)
+    ZLOGTS10   = LOG(PZREF(JLON)/10.0)
+    ZDELTAU10N(JLON) = ZDDU(JLON)-ZUSR(JLON)*(ZLOGUS10-ZPSI_U)/XKARMAN
+    ZDELTAT10N(JLON) = ZDDT(JLON)-ZTSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN
+    ZDELTAQ10N(JLON) = ZDDQ(JLON)-ZQSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN
+    ZDELTAU10N(JLON) = SIGN(MAX(ABS(ZDELTAU10N(JLON)),10.0*ZDUSR0),   &
+                            ZDELTAU10N(JLON))
+    ZDELTAT10N(JLON) = SIGN(MAX(ABS(ZDELTAT10N(JLON)),10.0*ZDTSR0),   &
+                            ZDELTAT10N(JLON))
+    ZDELTAQ10N(JLON) = SIGN(MAX(ABS(ZDELTAQ10N(JLON)),10.0*ZDQSR0),   &
+                            ZDELTAQ10N(JLON))
+!
+!       3.8. Test convergence for U*, T*, Q*
+!
+    IF (ABS(ZUSR(JLON)-ZUSR0(JLON)) < ZDUSR0 .AND.   &
+        ABS(ZTSR(JLON)-ZTSR0(JLON)) < ZDTSR0 .AND.   &
+        ABS(ZQSR(JLON)-ZQSR0(JLON)) < ZDQSR0) THEN
+      JCV(JLON) = 1                                     !free convergence
+      IF (JJ >= NITERMAX+1) JCV(JLON) = 2               !leaded convergence
+    ENDIF
+    JITER(JLON) = JJ
+  ENDIF
+!
+  ENDDO
+ENDDO
+!
+!-------------------------------------------------------------------------------
+!
+!       4.   COMPUTATION OF TURBULENT FLUXES AND EXCHANGE COEFFICIENTS.
+!       ---------------------------------------------------------------
+!
+DO JLON=1,SIZE(PTA)
+!
+!       4.1. Surface turbulent fluxes
+!            (ATM CONV.: ZTAU<<0 ; ZHF,ZEF<0 if atm looses heat)
+!
+  ZTAU(JLON) = -PRHOA(JLON)*ZUSR(JLON)**2
+  ZHF(JLON)  = -PRHOA(JLON)*ZCPA(JLON)*ZUSR(JLON)*ZTSR(JLON)
+  ZEF(JLON)  = -PRHOA(JLON)*ZLVS(JLON)*ZUSR(JLON)*ZQSR(JLON)
+
+  write(*,*) "ZTAU = ",ZTAU(JLON)
+  write(*,*) "SENS = ",ZHF(JLON)
+  write(*,*) "LAT = ",ZEF(JLON)
+
+!
+!       4.2. Exchange coefficients PCD, PCH, PCE
+!
+  PCD(JLON) = (ZUSR(JLON)/ZDDU(JLON))**2
+  PCH(JLON) = (ZUSR(JLON)*ZTSR(JLON))/(ZDDU(JLON)*ZDDT(JLON))
+  PCE(JLON) = (ZUSR(JLON)*ZQSR(JLON))/(ZDDU(JLON)*ZDDQ(JLON))
+
+  write(*,*) "ZUSR = ",ZUSR(JLON)
+  write(*,*) "ZTSR = ",ZTSR(JLON)
+  write(*,*) "ZQSR = ",ZQSR(JLON)
+
+!
+!       4.3. Stochastic perturbation of turbulent fluxes
+!
+!  IF( OPERTFLUX )THEN
+!    ZTAU(JLON) = ZTAU(JLON)* ( 1. + PPERTFLUX(JLON) / 2. )
+!    ZHF (JLON) = ZHF(JLON)*  ( 1. + PPERTFLUX(JLON) / 2. )
+!    ZEF (JLON) = ZEF(JLON)*  ( 1. + PPERTFLUX(JLON) / 2. )
+!  ENDIF
+!
+ENDDO
+!
+!-------------------------------------------------------------------------------
+!
+!       5.   COMPUTATION OF FLUX CORRECTIONS DUE TO RAINFALL.
+!            (ATM conv: ZRF<0 if atm. looses heat, ZTAUR<<0)
+!       -----------------------------------------------------
+!
+IF (OPRECIP) THEN
+  DO JLON=1,SIZE(PTA)
+!
+!       5.1. Momentum flux due to rainfall (ZTAUR, N/m2)
+!
+! See pp3752 in FBR96.
+    ZTAUR(JLON) = -0.85*PRAIN(JLON)*PVMOD(JLON)
+!
+!       5.2. Sensible heat flux due to rainfall (ZRF, W/m2)
+!
+! See Eq.12 in GoF95 with ZCPWA as specific heat of water at atm level (J/kg/K),
+! ZDQSDT from Clausius-Clapeyron relation, ZDWAT as water vapor diffusivity 
+! (Eq.13-3 of Pruppacher and Klett, 1978), ZDTMP as heat diffusivity, and ZBULB
+! as wet-bulb factor (Eq.11 in GoF95).
+!
+    ZTAC   = PTA(JLON)-XTT
+    ZCPWA  = 4217.51 -3.65566*ZTAC +0.1381*ZTAC**2       &
+              -2.8309E-03*ZTAC**3 +3.42061E-05*ZTAC**4   &
+              -2.18107E-07*ZTAC**5 +5.74535E-10*ZTAC**6
+    ZDQSDT = (ZLVA(JLON)*ZQSATA(JLON))/(XRV*PTA(JLON)**2)
+    ZDWAT  = 2.11E-05*(ZP00/PPA(JLON))*(PTA(JLON)/XTT)**1.94
+    ZDTMP  = (1.0+3.309E-03*ZTAC-1.44E-06*ZTAC**2)   &
+              *0.02411/(PRHOA(JLON)*ZCPA(JLON))
+    ZBULB  = 1.0/(1.0+ZDQSDT*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP))
+    ZRF(JLON) = PRAIN(JLON)*ZCPWA*ZBULB*((PSST(JLON)-PTA(JLON))   &
+                +(PQSATA(JLON)-PQA(JLON))*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP))
+!
+  ENDDO
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+!       6.   COMPUTATION OF WEBB CORRECTION TO LATENT HEAT FLUX (ZEFWEBB, W/m2).
+!       ------------------------------------------------------------------------
+!
+! See Eq.21 and Eq.22 in FBR96.
+IF (OPWEBB) THEN
+  DO JLON=1,SIZE(PTA)
+    ZWW = (1.0+ZETV)*(ZUSR(JLON)*ZQSR(JLON))   &
+           +(1.0+(1.0+ZETV)*PQA(JLON))*(ZUSR(JLON)*ZTSR(JLON))/PTA(JLON)
+    ZEFWEBB(JLON) = -PRHOA(JLON)*ZLVS(JLON)*ZWW*PQA(JLON)
+  ENDDO
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+!       7.   FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS. 
+!       ---------------------------------------------------------------
+!
+!       7.1. Richardson number
+!
+! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,PTA,PQA,   &
+!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI)
+!
+!       7.2. Friction velocity which contains correction due to rain
+!
+ZUSTAR2(:) = -(ZTAU(:)+ZTAUR(:))/PRHOA(:)       !>>0 as ZTAU<<0 & ZTAUR<=0
+!
+IF (OPRECIP) THEN
+  PCD(:) = ZUSTAR2(:)/ZDDU(:)**2
+ENDIF
+!
+PUSTAR(:) = SQRT(ZUSTAR2(:))                    !>>0
+!
+!       7.3. Aerodynamical conductance and resistance
+!
+ZAC  (:) = PCH(:)*ZDDU(:)
+PRESA(:) = 1.0/ZAC(:)
+!
+!       7.4. Total surface fluxes
+!
+PSFTH(:) =  ZHF(:)+ZRF(:)
+PSFTQ(:) = (ZEF(:)+ZEFWEBB(:))/ZLVS(:)
+!
+!       7.5. Charnock number
+!
+IF (CCHARNOCK == 'OLD') THEN
+  ZCHARN(:) = XVCHRNK
+ELSE            !modified for moderate wind speed as in COARE3.0
+  ZCHARN(:) = MIN(0.018,MAX(0.011,0.011+(0.007/8.0)*(ZDDU(:)-10.0)))
+ENDIF
+!
+!       7.6. Roughness lengths Z0 and Z0H over sea
+!
+!IF (KZ0 == 0) THEN      ! ARPEGE formulation
+!  PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + XVZ0CM*PCD(:)/PCDN(:)
+!  PZ0HSEA(:) = PZ0SEA (:)
+!ELSEIF (KZ0 == 1) THEN  ! Smith (1988) formulation
+!  PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + 0.11*ZVISA(:)/PUSTAR(:)
+!  PZ0HSEA(:) = PZ0SEA (:)
+!ELSEIF (KZ0 == 2) THEN  ! Direct computation using the stability functions
+!  DO JLON=1,SIZE(PTA)
+!    PZ0SEA (JLON) = PUREF(JLON)/EXP(XKARMAN*ZDDU(JLON)/PUSTAR(JLON)+ZPSIU(JLON))
+!    Z0TSEA        = PZREF(JLON)/EXP(XKARMAN*ZDDT(JLON)/ZTSR  (JLON)+ZPSIT(JLON))
+!    Z0QSEA        = PZREF(JLON)/EXP(XKARMAN*ZDDQ(JLON)/ZQSR  (JLON)+ZPSIT(JLON))
+!    PZ0HSEA(JLON) = 0.5*(Z0TSEA+Z0QSEA)
+!  ENDDO
+!ENDIF
+
+write(*,*) "JLON ",JLON
+write(*,*) "PTA ",klon,PTA
+write(*,*) "PCD ",SIZE(PCD),PCD
+write(*,*) "PCQ ",SIZE(PCE),PCE
+write(*,*) "PCH ",SIZE(PCH),PCH
+
+coeffs = [PCD,&
+       PCE,&
+       PCH]
+!
+!
+!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',1,ZHOOK_HANDLE)
+!
+!-------------------------------------------------------------------------------
+   END SUBROUTINE ECUMEV6_FLUX
Index: /LMDZ6/trunk/libf/phylmd/ini_csts.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/ini_csts.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/ini_csts.F90	(revision 4722)
@@ -0,0 +1,187 @@
+!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
+!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!SFX_LIC for details. version 1.
+!     #########
+      SUBROUTINE INI_CSTS 
+!     ##################
+!
+!!****  *INI_CSTS * - routine to initialize the module MODD_CST
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this routine is to initialize  the physical constants
+!     stored in  module MODD_CST.
+!      
+!
+!!**  METHOD
+!!    ------
+!!      The physical constants are set to their numerical values 
+!!     
+!!
+!!    EXTERNAL
+!!    --------
+!!      FMLOOK : to retrieve logical unit number associated to a file
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_CST     : contains physical constants
+!!
+!!    REFERENCE
+!!    ---------
+!!      Book2 of the documentation (module MODD_CST, routine INI_CSTS)
+!!      
+!!
+!!    AUTHOR
+!!    ------
+!!      V. Ducrocq       * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    18/05/94 
+!!      J. Stein    02/01/95  add the volumic mass of liquid water
+!!      J.-P. Pinty 13/12/95  add the water vapor pressure over solid ice
+!!      J. Stein    29/06/97  add XTH00
+!!      V. Masson   05/10/98  add XRHOLI
+!!      C. Mari     31/10/00  add NDAYSEC
+!!      V. Masson   01/03/03  add XCONDI
+!!      A. Voldoire 01/12/09  add XTTSI, XICEC, XTTS for ESM
+!!      J. Escobar  28/03/2014 for pb with emissivity/aerosol reset XSURF_TINY=1.0e-80 in real8 case
+!!
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_CSTS
+!
+!
+!USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
+!USE PARKIND1  ,ONLY : JPRB
+!
+!USE MODI_INI_CTURBS
+!
+!USE MODI_INI_OCEAN_CSTS
+!
+!USE MODI_INI_SURF_CSTS
+!
+IMPLICIT NONE
+!  
+!-------------------------------------------------------------------------------
+!
+!*       1.     FUNDAMENTAL CONSTANTS
+!               ---------------------
+!
+
+!REAL(KIND=JPRB) :: ZHOOK_HANDLE
+
+!IF (LHOOK) CALL DR_HOOK('INI_CSTS',0,ZHOOK_HANDLE)
+
+!#ifdef SFX_MNH
+!#ifdef MNH_MPI_DOUBLE_PRECISION
+XSURF_TINY    = 1.0e-80
+!#else
+!XSURF_TINY    = TINY    (XSURF_TINY    )
+!#endif
+!#else
+XSURF_TINY    = 1.0e-80
+!#endif
+XSURF_TINY_12 = SQRT    (XSURF_TINY    )
+XSURF_EPSILON = EPSILON (XSURF_EPSILON ) * 10.0
+
+XPI         = 2.*ASIN(1.)
+XKARMAN     = 0.4
+XBOLTZ      = 1.380658E-23
+XLIGHTSPEED = 299792458.
+XPLANCK     = 6.6260755E-34
+XAVOGADRO   = 6.0221367E+23
+!
+!-------------------------------------------------------------------------------
+!
+!*       2.     ASTRONOMICAL CONSTANTS
+!               ----------------------
+!
+XDAY   = 86400.
+XSIYEA = 365.25*XDAY*2.*XPI/ 6.283076
+XSIDAY = XDAY/(1.+XDAY/XSIYEA)
+XOMEGA = 2.*XPI/XSIDAY
+NDAYSEC = 24*3600 ! Number of seconds in a day
+!
+!-------------------------------------------------------------------------------!
+!
+!
+!*       3.     TERRESTRIAL GEOIDE CONSTANTS
+!               ----------------------------
+!
+XRADIUS = 6371229.
+XG      = 9.80665
+!
+!-------------------------------------------------------------------------------
+!
+!*       4.     REFERENCE PRESSURE
+!               -------------------
+!
+XP00 = 1.E5
+XTH00 = 300.
+!-------------------------------------------------------------------------------
+!
+!*       5.     RADIATION CONSTANTS
+!               -------------------
+!
+!JUAN OVERFLOW XSTEFAN = 2.* XPI**5 * XBOLTZ**4 / (15.* XLIGHTSPEED**2 * XPLANCK**3)
+XSTEFAN = ( 2.* XPI**5 / 15. ) * ( (XBOLTZ / XPLANCK)* XBOLTZ ) * (XBOLTZ/(XLIGHTSPEED*XPLANCK))**2 
+XI0     = 1370.
+!
+!-------------------------------------------------------------------------------
+!
+!*       6.     THERMODYNAMIC CONSTANTS
+!               -----------------------
+!
+XMD    = 28.9644E-3
+XMV    = 18.0153E-3
+XRD    = XAVOGADRO * XBOLTZ / XMD
+XRV    = XAVOGADRO * XBOLTZ / XMV
+XCPD   = 7.* XRD /2.
+XCPV   = 4.* XRV
+XRHOLW = 1000.
+XRHOLI = 917.
+XCONDI = 2.22
+XCL    = 4.218E+3
+XCI    = 2.106E+3
+XTT    = 273.16
+XTTSI  = XTT - 1.8
+XICEC  = 0.5
+XTTS   = XTT*(1-XICEC) + XTTSI*XICEC
+XLVTT  = 2.5008E+6
+XLSTT  = 2.8345E+6
+XLMTT  = XLSTT - XLVTT
+XESTT  = 611.14
+XGAMW  = (XCL - XCPV) / XRV
+XBETAW = (XLVTT/XRV) + (XGAMW * XTT)
+XALPW  = LOG(XESTT) + (XBETAW /XTT) + (XGAMW *LOG(XTT))
+XGAMI  = (XCI - XCPV) / XRV
+XBETAI = (XLSTT/XRV) + (XGAMI * XTT)
+XALPI  = LOG(XESTT) + (XBETAI /XTT) + (XGAMI *LOG(XTT))
+!
+!-------------------------------------------------------------------------------
+!
+!*       7.     TURBULENCE CONSTANTS
+!               --------------------
+!
+! CALL INI_CTURBS
+!-------------------------------------------------------------------------------
+!
+!*       8.     OCEAN CONSTANTS
+!               ---------------
+!
+! CALL INI_OCEAN_CSTS
+!
+!*       9.     SURFACE CONSTANTS
+!               -----------------
+!
+! CALL INI_SURF_CSTS
+!IF (LHOOK) CALL DR_HOOK('INI_CSTS',1,ZHOOK_HANDLE)
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE INI_CSTS 
Index: /LMDZ6/trunk/libf/phylmd/modd_csts.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/modd_csts.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/modd_csts.F90	(revision 4722)
@@ -0,0 +1,91 @@
+!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
+!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!SFX_LIC for details. version 1.
+!     ###############
+      MODULE MODD_CSTS      
+!     ###############
+!
+!!****  *MODD_CSTS* - declaration of Physic constants 
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this declarative module is to declare  the 
+!     Physics constants.    
+!
+!!
+!!**  IMPLICIT ARGUMENTS
+!!    ------------------
+!!      None 
+!!
+!!    REFERENCE
+!!    ---------
+!!          
+!!    AUTHOR
+!!    ------
+!!      V. Ducrocq   *Meteo France*
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    16/05/94  
+!!      J. Stein    02/01/95  add xrholw                    
+!!      J.-P. Pinty 13/12/95  add XALPI,XBETAI,XGAMI
+!!      J. Stein    25/07/97  add XTH00                    
+!!      V. Masson   05/10/98  add XRHOLI
+!!      C. Mari     31/10/00  add NDAYSEC
+!!      J. Escobar     06/13  add XSURF_TIMY XSURF_TIMY_12 XSURF_EPSILON for REAL*4
+!-------------------------------------------------------------------------------
+!
+!*       0.   DECLARATIONS
+!             ------------
+!
+IMPLICIT NONE 
+REAL,SAVE :: XPI                ! Pi
+!
+REAL,SAVE :: XDAY,XSIYEA,XSIDAY ! day duration, sideral year duration,
+                                ! sideral day duration
+!
+REAL,SAVE :: XKARMAN            ! von karman constant
+REAL,SAVE :: XLIGHTSPEED        ! light speed
+REAL,SAVE :: XPLANCK            ! Planck constant
+REAL,SAVE :: XBOLTZ             ! Boltzman constant 
+REAL,SAVE :: XAVOGADRO          ! Avogadro number
+!
+REAL,SAVE :: XRADIUS,XOMEGA     ! Earth radius, earth rotation
+REAL,SAVE :: XG                 ! Gravity constant
+!
+REAL,SAVE :: XP00               ! Reference pressure
+!
+REAL,SAVE :: XSTEFAN,XI0        ! Stefan-Boltzman constant, solar constant
+!
+REAL,SAVE :: XMD,XMV            ! Molar mass of dry air and molar mass of vapor
+REAL,SAVE :: XRD,XRV            ! Gaz constant for dry air, gaz constant for vapor
+REAL,SAVE :: XCPD,XCPV          ! Cpd (dry air), Cpv (vapor)
+REAL,SAVE :: XRHOLW             ! Volumic mass of liquid water
+REAL,SAVE :: XCL,XCI            ! Cl (liquid), Ci (ice)
+REAL,SAVE :: XTT                ! Triple point temperature
+REAL,SAVE :: XTTSI              ! Temperature of ice fusion over salty sea
+REAL,SAVE :: XTTS               ! Equivalent temperature of ice fusion over a mixed of sea and sea-ice
+REAL,SAVE :: XICEC              ! Threshold fraction over which the tile is considered as only covered with ice
+REAL,SAVE :: XLVTT              ! Vaporization heat constant
+REAL,SAVE :: XLSTT              ! Sublimation heat constant
+REAL,SAVE :: XLMTT              ! Melting heat constant
+REAL,SAVE :: XESTT              ! Saturation vapor pressure  at triple point
+                                ! temperature  
+REAL,SAVE :: XALPW,XBETAW,XGAMW ! Constants for saturation vapor 
+                                !  pressure  function 
+REAL,SAVE :: XALPI,XBETAI,XGAMI ! Constants for saturation vapor
+                                !  pressure  function over solid ice
+REAL, SAVE        :: XTH00      ! reference value  for the potential
+                                ! temperature
+REAL,SAVE :: XRHOLI             ! Volumic mass of ice
+REAL,SAVE :: XCONDI             ! thermal conductivity of ice (W m-1 K-1)
+!
+INTEGER, SAVE :: NDAYSEC        ! Number of seconds in a day
+!
+REAL,SAVE     :: XSURF_TINY          ! minimum real on this machine
+REAL,SAVE     :: XSURF_TINY_12       ! sqrt(minimum real on this machine)
+REAL,SAVE     :: XSURF_EPSILON       ! minimum space with 1.0
+!
+END MODULE MODD_CSTS
+
Index: /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4721)
+++ /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4722)
@@ -348,5 +348,5 @@
     REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro  ! rugosity length
     REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosine of solar zenith angle
-    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
+    REAL, DIMENSION(klon),        INTENT(INOUT)     :: rain_f  ! rain fall
     REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
     REAL, DIMENSION(klon),        INTENT(IN)        :: bs_f  ! blowing snow fall
@@ -1071,5 +1071,5 @@
 !albedo SB <<<
     yrain_f = 0.0 ; ysnow_f = 0.0  ; ybs_f=0.0  ; yfder = 0.0     ; ysolsw = 0.0
-    ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yu1 = 0.0    
+    ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yz0h_oupas = 0.0 ; yu1 = 0.0    
     yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0     ; yqbs1 = 0.0
     ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
@@ -1598,7 +1598,7 @@
         ENDDO
         CALL cdrag(knon, nsrf, &
-            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
+            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, &
             yts, yqsurf, yz0m, yz0h, yri0, 0, &
-            ycdragm, ycdragh, zri1, pref )
+            ycdragm, ycdragh, zri1, pref, rain_f, zxtsol, ypplay(:,1))
 
 ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
@@ -1632,7 +1632,7 @@
 
             CALL cdrag(knon, nsrf, &
-            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
+            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),s_pblh_x,&
             yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
-            ycdragm_x, ycdragh_x, zri1_x, pref_x )
+            ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, zxtsol, ypplay(:,1) )
 
 ! --- special Dice. JYG+MPL 25112013
@@ -1659,7 +1659,7 @@
         ENDDO
         CALL cdrag(knon, nsrf, &
-            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
+            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,&
             yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
-            ycdragm_w, ycdragh_w, zri1_w, pref_w )
+            ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, zxtsol, ypplay(:,1) )
 !
 !!!bug !!        zgeo1(:) = wake_s(:)*zgeo1_w(:) + (1.-wake_s(:))*zgeo1_x(:)
@@ -2086,5 +2086,5 @@
                yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
                yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
-               yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
+               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
           ENDIF
           
@@ -3072,5 +3072,5 @@
        IF (iflag_split .eq.0) THEN
         IF (iflag_new_t2mq2m==1) THEN
-         CALL stdlevvarn(klon, knon, nsrf, zxli, &
+           CALL stdlevvarn(klon, knon, nsrf, zxli, &
             uzon, vmer, tair1, qair1, zgeo1, &
             tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
@@ -3081,5 +3081,5 @@
             uzon, vmer, tair1, qair1, zgeo1, &
             tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
-            yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
+            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
         ENDIF
        ELSE  !(iflag_split .eq.0)
@@ -3099,9 +3099,9 @@
             uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
             tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
-            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
+            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol)
         CALL stdlevvar(klon, knon, nsrf, zxli, &
             uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
             tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
-            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
+            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol)
         ENDIF
 !!!
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4721)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4722)
@@ -1259,4 +1259,7 @@
     REAL, dimension(klon,klev) :: t_env,q_env
 
+    REAL, dimension(klon) :: pr_et
+    REAL, dimension(klon) :: w_et, jlr_g_c, jlr_g_s
+
     REAL pi
     INTEGER ieru
@@ -2804,4 +2807,15 @@
        ELSE IF (iflag_gusts==2) THEN
           gustiness(1:klon)=f_gust_bl*ale_bl_stat(1:klon)+f_gust_wk*ale_wake(1:klon)
+       !!!! modif olivier torres
+       ELSE IF (iflag_gusts==3) THEN
+          w_et=wstar(1,3)
+          jlr_g_s=(0.65*w_et)**2
+          pr_et=rain_con*8640
+          jlr_g_c = (((19.8*(pr_et(1:klon)**2))/(1.5+pr_et(1:klon)+pr_et(1:klon)**2))**(0.4))**2
+          gustiness(1:klon)=jlr_g_c+jlr_g_s
+!!       write(*,*) "rain ",pr_et
+!!       write(*,*) "jlr_g_c",jlr_g_c
+!!       write(*,*) "wstar",wstar(1,3)
+!!       write(*,*) "jlr_g_s",jlr_g_s
           ! ELSE IF (iflag_gusts==2) THEN
           !    do i = 1, klon
Index: /LMDZ6/trunk/libf/phylmd/qsat_seawater.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/qsat_seawater.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/qsat_seawater.F90	(revision 4722)
@@ -0,0 +1,119 @@
+!-------------------------------------------------------------------------------
+!
+!     ######################################
+REAL      FUNCTION QSAT_SEAWATER(PT,PP) 
+!     ######################################
+!
+!!****  *QSATW * - function to compute saturation vapor humidity from
+!!                 temperature
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this function is to compute the saturation vapor 
+!     pressure from temperature over saline seawater
+!      
+!
+!!**  METHOD
+!!    ------
+!!       Given temperature T (PT), the saturation vapor pressure es(T)
+!!    (FOES(PT)) is computed by integration of the Clapeyron equation
+!!    from the triple point temperature Tt (XTT) and the saturation vapor 
+!!    pressure of the triple point es(Tt) (XESTT), i.e  
+!!    The reduction due to salinity is compute with the factor 0.98 (reduction of 2%)
+!!     
+!!         es(T)= 0.98*EXP( alphaw - betaw /T - gammaw Log(T) )
+!!  
+!!     with :
+!!       alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt) 
+!!       betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!!       gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!!      Then, the specific humidity at saturation is deduced.
+!!  
+!!
+!!    EXTERNAL
+!!    --------
+!!      NONE
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_CST : comtains physical constants
+!!        XALPW   : Constant for saturation vapor pressure function
+!!        XBETAW  : Constant for saturation vapor pressure function
+!!        XGAMW   : Constant for saturation vapor pressure function  
+!!      
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH 
+!!      Zeng, X., Zhao, M., and Dickinson, R. E., 1998 : Intercomparaison of bulk
+!!      aerodynamic algorithm for the computation of sea surface fluxes using
+!!      TOGA COARE and TAO data. Journal of Climate, vol 11, n°10, pp 2628--2644
+!!
+!!
+!!    AUTHOR
+!!    ------
+!!      C. Lebeaupin    * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    6/04/2005 
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_CSTS
+USE dimphy
+USE indice_sol_mod
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of arguments and results
+!
+!
+
+REAL, DIMENSION(klon), INTENT(IN)                :: PT     ! Temperature
+                                                        ! (Kelvin)
+REAL, DIMENSION(klon), INTENT(IN)                :: PP     ! Pressure
+                                                        ! (Pa)
+REAL, DIMENSION(SIZE(PT))                        :: PQSAT  ! saturation vapor 
+                                                        ! specific humidity
+                                                        ! with respect to
+                                                        ! water (kg/kg)
+
+!
+!*       0.2   Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT))                   :: ZFOES  ! saturation vapor 
+                                                        ! pressure
+                                                        ! (Pascal) 
+!
+REAL, DIMENSION(SIZE(PT))                   :: ZWORK1
+REAL                                        :: ZWORK2
+!REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!-------------------------------------------------------------------------------
+!
+!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',0,ZHOOK_HANDLE)
+!
+!ZFOES  = 1 !PSAT(PT(:))
+!ZFOES  = 0.98*ZFOES
+!ZFOES(:) = ZFOES(:)*1013.25E+02             !convert from atm to Pa
+ZFOES = 0.98*EXP( XALPW - XBETAW/PT - XGAMW*LOG(PT)  )
+! vapor pressure reduction of 2% over saline seawater could have a significant 
+! impact on the computation of surface latent heat flux under strong wind 
+! conditions (Zeng et al, 1998). 
+!
+ZWORK1 = ZFOES/PP
+ZWORK2    = XRD/XRV
+!
+!*       2.    COMPUTE SATURATION HUMIDITY
+!              ---------------------------
+!
+PQSAT = ZWORK2*ZWORK1 / (1.+(ZWORK2-1.)*ZWORK1)
+!
+!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',1,ZHOOK_HANDLE)
+!-------------------------------------------------------------------------------
+!
+END FUNCTION QSAT_SEAWATER
+!
+!-------------------------------------------------------------------------------
Index: /LMDZ6/trunk/libf/phylmd/qsat_seawater2.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/qsat_seawater2.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/qsat_seawater2.F90	(revision 4722)
@@ -0,0 +1,104 @@
+!-------------------------------------------------------------------------------
+!
+!     ######################################
+REAL      FUNCTION QSAT_SEAWATER2(PT,PP,PSSS) 
+!     ######################################
+!
+!!****  *QSATW * - function to compute saturation vapor humidity from
+!!                 temperature
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this function is to compute the saturation vapor 
+!     pressure from temperature over saline seawater
+!      
+!
+!!**  METHOD
+!!    ------
+!!       Given temperature T (PT) and salinity S (PSSS), the saturation vapor 
+!!    pressure es(T,S) (FOES(PT,PSSS)) is computed following Weiss and Price
+!!    (1980).
+!!
+!!      Then, the specific humidity at saturation is deduced.
+!!  
+!!
+!!    EXTERNAL
+!!    --------
+!!      NONE
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_CST : contains physical constants
+!!      
+!!    REFERENCE
+!!    ---------
+!!      Weiss, R.F., and Price, B.A., 1980 : Nitrous oxide solubility in water
+!!      and seawater. Marine Chemistry, n°8, pp 347-359.
+!!
+!!
+!!    AUTHOR
+!!    ------
+!!      S. Belamari     * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    19/03/2014 
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_CSTS, ONLY : XRD, XRV
+USE dimphy
+USE indice_sol_mod
+!
+IMPLICIT NONE
+!
+!*       0.1   Declarations of arguments and results
+!
+!
+REAL, DIMENSION(klon), INTENT(IN)                :: PT     ! Temperature
+                                                        ! (Kelvin)
+REAL, DIMENSION(klon), INTENT(IN)                :: PP     ! Pressure
+                                                        ! (Pascal)
+REAL, DIMENSION(klon), INTENT(IN)                :: PSSS   ! Salinity
+                                                        ! (g/kg)
+REAL, DIMENSION(SIZE(PT))                   :: PQSATA  ! saturation vapor
+                                                        ! specific humidity
+                                                        ! with respect to
+                                                        ! water (kg/kg)
+!
+!*       0.2   Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT))                   :: ZFOES  ! saturation vapor
+                                                        ! pressure
+                                                        ! (Pascal)
+!
+REAL, DIMENSION(SIZE(PT))                   :: ZWORK1
+REAL                                        :: ZWORK2
+!REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!-------------------------------------------------------------------------------
+!
+!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW2_1D',0,ZHOOK_HANDLE)
+!
+!*       1.    COMPUTE SATURATION VAPOR PRESSURE
+!              ---------------------------------
+!
+ZFOES(:) = EXP( 24.4543 -67.4509*(100.0/PT(:)) -4.8489*LOG(PT(:)/100.0)   &
+                -5.44E-04*(PSSS(:)/1.00472) ) !see Sharqawy et al (2010) Eq32 p368
+ZFOES(:) = ZFOES(:)*1013.25E+02             !convert from atm to Pa
+!
+ZWORK1(:) = ZFOES(:)/PP(:)
+ZWORK2    = XRD/XRV
+!
+!*       2.    COMPUTE SATURATION SPECIFIC HUMIDITY
+!              ------------------------------------
+!
+PQSATA(:) = ZWORK2*ZWORK1(:) / (1.0+(ZWORK2-1.0)*ZWORK1(:))
+!
+!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW2_1D',1,ZHOOK_HANDLE)
+!-------------------------------------------------------------------------------
+!
+END FUNCTION QSAT_SEAWATER2
+!
+!-------------------------------------------------------------------------------
Index: /LMDZ6/trunk/libf/phylmd/qsatseaw_1D.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/qsatseaw_1D.F90	(revision 4722)
+++ /LMDZ6/trunk/libf/phylmd/qsatseaw_1D.F90	(revision 4722)
@@ -0,0 +1,113 @@
+!-------------------------------------------------------------------------------
+!
+!     ######################################
+REAL       FUNCTION QSATSEAW_1D(PT,PP) 
+!     ######################################
+!
+!!****  *QSATW * - function to compute saturation vapor humidity from
+!!                 temperature
+!!
+!!    PURPOSE
+!!    -------
+!       The purpose of this function is to compute the saturation vapor 
+!     pressure from temperature over saline seawater
+!      
+!
+!!**  METHOD
+!!    ------
+!!       Given temperature T (PT), the saturation vapor pressure es(T)
+!!    (FOES(PT)) is computed by integration of the Clapeyron equation
+!!    from the triple point temperature Tt (XTT) and the saturation vapor 
+!!    pressure of the triple point es(Tt) (XESTT), i.e  
+!!    The reduction due to salinity is compute with the factor 0.98 (reduction of 2%)
+!!     
+!!         es(T)= 0.98*EXP( alphaw - betaw /T - gammaw Log(T) )
+!!  
+!!     with :
+!!       alphaw (XALPW) = LOG(es(Tt))+ betaw/Tt + gammaw Log(Tt) 
+!!       betaw (XBETAW) = Lv(Tt)/Rv + gammaw Tt
+!!       gammaw (XGAMW) = (Cl -Cpv) /Rv
+!!
+!!      Then, the specific humidity at saturation is deduced.
+!!  
+!!
+!!    EXTERNAL
+!!    --------
+!!      NONE
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!      Module MODD_CST : comtains physical constants
+!!        XALPW   : Constant for saturation vapor pressure function
+!!        XBETAW  : Constant for saturation vapor pressure function
+!!        XGAMW   : Constant for saturation vapor pressure function  
+!!      
+!!    REFERENCE
+!!    ---------
+!!      Book2 of documentation of Meso-NH 
+!!      Zeng, X., Zhao, M., and Dickinson, R. E., 1998 : Intercomparaison of bulk
+!!      aerodynamic algorithm for the computation of sea surface fluxes using
+!!      TOGA COARE and TAO data. Journal of Climate, vol 11, n�10, pp 2628--2644
+!!
+!!
+!!    AUTHOR
+!!    ------
+!!      C. Lebeaupin    * Meteo France *
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    6/04/2005 
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODD_CSTS
+USE dimphy
+USE indice_sol_mod
+
+IMPLICIT NONE
+
+!*       0.1   Declarations of arguments and results
+!
+!
+REAL, DIMENSION(klon), INTENT(IN)                :: PT     ! Temperature
+                                                        ! (Kelvin)
+REAL, DIMENSION(klon), INTENT(IN)                :: PP     ! Pressure
+                                                        ! (Pa)
+REAL, DIMENSION(SIZE(PT))                        :: PQSAT  ! saturation vapor 
+                                                        ! specific humidity
+                                                        ! with respect to
+                                                        ! water (kg/kg)
+!
+!*       0.2   Declarations of local variables
+!
+REAL, DIMENSION(SIZE(PT))                   :: ZFOES  ! saturation vapor 
+                                                        ! pressure
+                                                        ! (Pascal) 
+!
+INTEGER                         :: JJ   ! loop index
+!REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!-------------------------------------------------------------------------------
+!
+!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',0,ZHOOK_HANDLE)
+!DO JJ = 1, SIZE(PT)
+!*       1.    COMPUTE SATURATION VAPOR PRESSURE
+!              ---------------------------------
+!
+ZFOES = 0.98*EXP( XALPW - XBETAW/PT - XGAMW*LOG(PT)  )
+! vapor pressure reduction of 2% over saline seawater could have a significant 
+! impact on the computation of surface latent heat flux under strong wind 
+! conditions (Zeng et al, 1998). 
+!
+!*       2.    COMPUTE SATURATION HUMIDITY
+!              ---------------------------
+!
+PQSAT = XRD/XRV*ZFOES/PP /(1.+(XRD/XRV-1.)*ZFOES/PP)  
+!
+!ENDDO
+!IF (LHOOK) CALL DR_HOOK('MODE_THERMOS:QSATSEAW_1D',1,ZHOOK_HANDLE)
+!-------------------------------------------------------------------------------
+!
+END FUNCTION QSATSEAW_1D
+!------------------------------------
Index: /LMDZ6/trunk/libf/phylmd/screenc_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/screenc_mod.F90	(revision 4721)
+++ /LMDZ6/trunk/libf/phylmd/screenc_mod.F90	(revision 4722)
@@ -18,5 +18,5 @@
                          ts, qsurf, z0m, z0h, psol, &
                          ustar, testar, qstar, okri, ri1, &
-                         pref, delu, delte, delq)
+                         pref, delu, delte, delq, s_pblh, prain, tsol, pat1)
       IMPLICIT NONE
 !-----------------------------------------------------------------------
@@ -60,6 +60,13 @@
       REAL, dimension(klon), intent(in) :: speed, temp, q_zref
       REAL, intent(in) :: zref
-      REAL, dimension(klon), intent(in) :: ts, qsurf, z0m, z0h, psol
-      REAL, dimension(klon), intent(in) :: ustar, testar, qstar, ri1         
+      REAL, dimension(klon), intent(IN) :: ts
+      REAL, dimension(klon), intent(in) :: qsurf, psol
+      REAL, dimension(klon), intent(inout):: z0m, z0h
+      REAL, dimension(klon), intent(in) :: ustar, testar, qstar, ri1   
+
+      REAL, dimension(klon), intent(inout) :: s_pblh    
+      REAL, dimension(klon), intent(in) :: prain    
+      REAL, dimension(klon), intent(in) :: tsol    
+      REAL, DIMENSION(klon), INTENT(IN)    :: pat1 !pression premier lev      
 !
       REAL, dimension(klon), intent(out) :: pref, delu, delte, delq 
@@ -88,7 +95,7 @@
       CALL cdrag (knon, nsrf, &
                     speed, temp, q_zref, gref, &
-                    psol, ts, qsurf, z0m, z0h, &
+                    psol, s_pblh, ts, qsurf, z0m, z0h, &
                     zri_zero,0,                &
-                    cdram, cdrah, zri1, pref)
+                    cdram, cdrah, zri1, pref, prain, tsol, pat1)
       DO i = 1, knon
         IF(ok_prescr_ust) THEN
@@ -114,5 +121,5 @@
                          cdrm, cdrh,  okri, &
                          ri1, iri1, &
-                         pref, delm, delh, zri1)
+                         pref, delm, delh, zri1, s_pblh, prain, tsol, pat1)
       IMPLICIT NONE
 !-----------------------------------------------------------------------
@@ -156,6 +163,11 @@
       REAL, dimension(klon), intent(in) :: speed, temp, q_zref
       REAL, intent(in) :: zref
-      REAL, dimension(klon), intent(in) :: ts, qsurf, z0m, z0h, psol
+      REAL, dimension(klon), intent(in) :: ts, qsurf, psol
+      REAL, dimension(klon), intent(inout) :: z0m, z0h
       REAL, dimension(klon), intent(in) :: cdrm, cdrh, ri1         
+      REAL, dimension(klon), intent(inout) :: s_pblh    
+      REAL, dimension(klon), intent(in) :: prain    
+      REAL, dimension(klon), intent(in) :: tsol    
+      REAL, DIMENSION(klon), INTENT(IN) :: pat1 !pression premier lev      
       INTEGER, INTENT(IN)  :: iri1 ! Richardson de la 1ere couche
 !
@@ -180,7 +192,7 @@
       CALL cdrag(knon, nsrf, &
                     speed, temp, q_zref, gref, &
-                    psol, ts, qsurf, z0m, z0h, &
+                    psol, s_pblh, ts, qsurf, z0m, z0h, &
                     ri1, iri1, &
-                    cdram, cdrah, zri1, pref)
+                    cdram, cdrah, zri1, pref, prain, tsol, pat1)
       DO i = 1, knon
         delm(i) = sqrt(cdrm(i))/sqrt(cdram(i))
Index: /LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90	(revision 4721)
+++ /LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90	(revision 4722)
@@ -19,5 +19,5 @@
                            u1, v1, t1, q1, z1, &
                            ts1, qsurf, z0m, z0h, psol, pat1, &
-                           t_2m, q_2m, t_10m, q_10m, u_10m, ustar)
+                           t_2m, q_2m, t_10m, q_10m, u_10m, ustar, s_pblh, prain, tsol)
       IMPLICIT NONE
 !-------------------------------------------------------------------------
@@ -61,9 +61,13 @@
       LOGICAL, intent(in) :: zxli
       REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1
-      REAL, dimension(klon), intent(in) :: qsurf, z0m, z0h
+      REAL, dimension(klon), intent(in) :: qsurf
+      REAL, dimension(klon), intent(inout) :: z0m, z0h
       REAL, dimension(klon), intent(in) :: psol, pat1
 !
       REAL, dimension(klon), intent(out) :: t_2m, q_2m, ustar
       REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m
+      REAL, DIMENSION(klon), INTENT(INOUT) :: s_pblh
+      REAL, DIMENSION(klon), INTENT(IN) :: prain
+      REAL, DIMENSION(klon), INTENT(IN) :: tsol
 !-------------------------------------------------------------------------
       include "flux_arp.h"
@@ -120,7 +124,7 @@
       CALL cdrag(knon, nsrf, &
  &                   speed, t1, q1, z1, &
- &                   psol, ts1, qsurf, z0m, z0h, &
+ &                   psol, s_pblh, ts1, qsurf, z0m, z0h, &
  &                   zri_zero, 0, &
- &                   cdram, cdrah, zri1, pref)
+ &                   cdram, cdrah, zri1, pref, prain, tsol, pat1)
 
 ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
@@ -178,5 +182,5 @@
  &                   ts1, qsurf, z0m, z0h, psol, &           
  &                   ustar, testar, qstar, okri, ri1, &
- &                   pref, delu, delte, delq) 
+ &                   pref, delu, delte, delq, s_pblh ,prain, tsol, pat1) 
 !
         DO i = 1, knon
@@ -280,5 +284,5 @@
  &                   ts1, qsurf, z0m, z0h, psol, &
  &                   ustar, testar, qstar, okri, ri1, &
- &                   pref, delu, delte, delq)
+ &                   pref, delu, delte, delq, s_pblh ,prain, tsol, pat1)
 !
         DO i = 1, knon
@@ -357,6 +361,7 @@
       INTEGER, intent(in) :: klon, knon, nsrf
       LOGICAL, intent(in) :: zxli
-      REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, z1, ts1
-      REAL, dimension(klon), intent(in) :: qsurf, z0m, z0h
+      REAL, dimension(klon), intent(in) :: u1, v1, t1, q1, ts1, z1
+      REAL, dimension(klon), intent(inout) :: z0m, z0h
+      REAL, dimension(klon), intent(in) :: qsurf
       REAL, dimension(klon), intent(in) :: psol, pat1
 !
@@ -371,4 +376,7 @@
       REAL, dimension(klon) :: cdmn2m, cdhn2m, fm2m, fh2m
       REAL, dimension(klon) :: ri2m_new
+      REAL, DIMENSION(klon) :: s_pblh
+      REAL, DIMENSION(klon) :: prain
+      REAL, DIMENSION(klon) :: tsol
 !-------------------------------------------------------------------------
       include "flux_arp.h"
@@ -444,7 +452,7 @@
       CALL cdrag(knon, nsrf, &
  &                   speed, t1, q1, z1, &
- &                   psol, ts1, qsurf, z0m, z0h, &
+ &                   psol, s_pblh, ts1, qsurf, z0m, z0h, &
  &                   zri_zero, 0, &
- &                   cdram, cdrah, zri1, pref)
+ &                   cdram, cdrah, zri1, pref, prain, tsol, pat1)
 
 !
@@ -468,5 +476,6 @@
  &                   cdram, cdrah,  okri, &
  &                   ri1, 1, &
- &                   pref_new, delm_new, delh_new, ri2m)
+ &                   pref_new, delm_new, delh_new, ri2m, &
+ &                   s_pblh, prain, tsol, pat1      )
 !
        DO i = 1, knon
@@ -535,5 +544,6 @@
  &                   cdram, cdrah,  okri, &
  &                   ri1, 0, &
- &                   pref, delm, delh, ri2m)
+ &                   pref, delm, delh, ri2m, &
+ &                   s_pblh, prain, tsol, pat1      )
 !
         DO i = 1, knon
@@ -614,5 +624,6 @@
  &                   cdram, cdrah,  okri, &
  &                   ri1, 1, &
- &                   pref_new, delm_new, delh_new, ri10m)
+ &                   pref_new, delm_new, delh_new, ri10m, &
+ &                   s_pblh, prain, tsol, pat1      )
 !
        DO i = 1, knon
@@ -671,5 +682,6 @@
  &                   cdram, cdrah,  okri, &
  &                   ri1, 0, &
- &                   pref, delm, delh, ri10m)
+ &                   pref, delm, delh, ri10m, &
+ &                   s_pblh, prain, tsol, pat1      )
 !
         DO i = 1, knon
