Index: /LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90	(revision 4478)
@@ -7,5 +7,5 @@
 subroutine atke_compute_km_kh(ngrid,nlay, &
                         wind_u,wind_v,temp,play,pinterf, &
-                        tke,Km,Kh)
+                        tke,Km_out,Kh_out)
 
 !========================================================================
@@ -16,5 +16,5 @@
 ! collective and collaborative workshop, 
 ! the so-called 'Atelier TKE (ATKE)' with
-! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, 
+! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos, 
 ! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon
 !
@@ -26,6 +26,7 @@
 
 
-USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cn, cinf, rpi
-USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, rg, rd
+USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cn, cinf, rpi, rcpd
+USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd
+USE atke_turbulence_ini_mod, ONLY : viscom, viscoh
 
 implicit none
@@ -47,19 +48,25 @@
 REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
 
-REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: Km   ! Exchange coefficient for momentum at interface between layers
-REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT)    :: Kh   ! Exchange coefficient for heat flux at interface between layers
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers
 
 ! Local variables
-REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange
-REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface
-REAL, DIMENSION(ngrid,nlay+1) :: dz_interf ! distance between two consecutive interfaces
-REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number
-REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number 
-REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum
-REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat
+REAL, DIMENSION(ngrid,nlay+1) :: Km          ! Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay+1) :: Kh          ! Exchange coefficient for heat flux at interface between layers
+REAL, DIMENSION(ngrid,nlay)   :: theta       ! Potential temperature
+REAL, DIMENSION(ngrid,nlay+1) :: l_exchange  ! Length of exchange (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: z_interf    ! Altitude at the interface
+REAL, DIMENSION(ngrid,nlay)   :: z_lay       ! Altitude of layers
+REAL, DIMENSION(ngrid,nlay)   :: dz_interf   ! distance between two consecutive interfaces
+REAL, DIMENSION(ngrid,nlay)   :: dz_lay      ! distance between two layer middles (NB: first and last are half layers)
+REAL, DIMENSION(ngrid,nlay+1) :: Ri          ! Richardson's number (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Prandtl     ! Turbulent Prandtl's number (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
 
 INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
-REAL    :: Ri0        ! parameter for Sm stability function
-
+REAL    :: Ri0,Ri1    ! parameter for Sm stability function and Prandlt
+REAL    :: preff      ! reference pressure for potential temperature calculations
+REAL    :: thetam     ! mean potential temperature at interface
 
 
@@ -68,24 +75,48 @@
 
 DO igrid=1,ngrid
+    dz_interf(igrid,1) = 0.0
     z_interf(igrid,1) = 0.0
-    Ri(igrid,1) = 0.1! TO DO
-    Sm(igrid,1) = 0.0 ! TO DO
+!    Km(igrid,1) = 0.0
+!    Kh(igrid,1) = 0.0
+!    tke(igrid,1) = 0.0
+!    Km(igrid,nlay+1) = 0.0              
+!    Kh(igrid,nlay+1) = 0.0             
+!    tke(igrid,nlay+1) = 0.0
+END DO
+
+! Calculation of potential temperature: (if vapor -> todo virtual potential temperature)
+!=====================================
+
+preff=100000.
+! The result should not depend on the choice of preff
+DO ilay=1,nlay
+     DO igrid = 1, ngrid
+        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
+     END DO
 END DO
 
 
 
-! Calculation of altitude of layers' bottom interfaces:
-!=====================================================
+! Calculation of altitude of layers' middle and bottom interfaces:
+!=================================================================
 
 DO ilay=2,nlay+1
     DO igrid=1,ngrid
-        dz_interf(igrid,ilay) = rg * temp(igrid,ilay) / rg * log(play(igrid,ilay-1)/play(igrid,ilay))
-        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay)
+        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
+        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
     ENDDO
+ENDDO
+
+DO ilay=1,nlay
+   DO igrid=1,ngrid 
+      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
+   ENDDO
 ENDDO
 
 
 ! Computing the mixing length:
-!=============================
+! so far, we have neglected the effect of local stratification
+!==============================================================
+
 
 DO ilay=2,nlay+1
@@ -98,20 +129,22 @@
 !===================================================================
 
-DO ilay=2,nlay+1
+DO ilay=2,nlay
     DO igrid=1,ngrid
-        ! Warning: sure that gradient of temp and wind should be calculated with dz_interf and not with dz_lay?
-        Ri(igrid,ilay) = rg * (temp(igrid,ilay) - temp(igrid,ilay-1)) / dz_interf(igrid,ilay) / & 
-        (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 + & 
-        ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 ) 
-
+        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
+        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
+        Ri(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) / & 
+        MAX(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 
+        ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2,1E-10) 
         ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
-        Ri0=2.*ric*cinf/rpi/cn
+        Ri0=2./rpi*(cinf - cn)*ric/cn
+        ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
+        Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope
 
         IF (Ri(igrid,ilay) < 0.) THEN
-            Sm(igrid,ilay) = 2./rpi * cinf * atan(-Ri(igrid,ilay)/Ri0) + cn
-            Prandtl(igrid,ilay) = 1. + Ri(igrid,ilay) * pr_slope
+            Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
+            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut 
         ELSE
             Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric))
-            Prandtl(igrid,ilay) = 2./rpi * atan(-Ri(igrid,ilay)) - 1. + pr_asym  ! exp(Ri(igrid,ilay))
+            Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope
         END IF
         
@@ -126,9 +159,9 @@
 
 ! stationary solution neglecting the vertical transport of TKE by turbulence
-    DO ilay=2,nlay+1
+    DO ilay=2,nlay
         DO igrid=1,ngrid
             tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
-            (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 + & 
-            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_interf(igrid,ilay))**2 ) * &
+            (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 
+            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) * &
             (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
         ENDDO
@@ -145,12 +178,16 @@
 ! Computing eddy diffusivity coefficients:
 !========================================
-DO ilay=2,nlay+1
+DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
     DO igrid=1,ngrid
-        Km(igrid,ilay) = l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
-        Kh(igrid,ilay) = l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
+        ! we add the molecular viscosity to Km,h 
+        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
+        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
     END DO
 END DO
 
-
+! for output:
+!===========
+Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
+Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
 
 end subroutine atke_compute_km_kh
Index: /LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90	(revision 4478)
@@ -9,20 +9,24 @@
   real :: kappa = 0.4 ! Von Karman constant
   !$OMP THREADPRIVATE(kappa)
-  real :: l0, ric, ri0, cn, cinf, cepsilon, pr_slope, pr_asym
-  !$OMP THREADPRIVATE(l0, ric, cn, cinf, cepsilon, pr_slope, pr_asym)
+  real :: l0, ric, ri0, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut
+  !$OMP THREADPRIVATE(l0, ric, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut)
   integer :: lunout,prt_level
   !$OMP THREADPRIVATE(lunout,prt_level)
-  real :: rg, rd, rpi
-  !$OMP THREADPRIVATE(rg, rd, rpi)
+  real :: rg, rd, rpi, rcpd
+  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd)
+
+  real :: viscom, viscoh
+  !$OMP THREADPRIVATE(viscom,viscoh)
+
 
 
 CONTAINS
 
-SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in)
+SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in)
 
    USE ioipsl_getin_p_mod, ONLY : getin_p
 
   integer, intent(in) :: lunout_in,prt_level_in
-  real, intent(in) :: rg_in, rd_in, rpi_in
+  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in
 
 
@@ -32,4 +36,8 @@
   rg=rg_in
   rpi=rpi_in
+  rcpd=rcpd_in
+
+  viscom=1.46E-5
+  viscoh=2.06E-5
 
   ! flag that controls options in atke_compute_km_kh
@@ -65,4 +73,9 @@
   CALL getin_p('atke_pr_asym',pr_asym)
 
+  ! value of turbulent prandtl number in neutral conditions (Ri=0)
+  ! to guarantee tke>0 in stationary case, pr_neut has to be >= 1
+  pr_neut=1.0
+  CALL getin_p('atke_pr_neut',pr_neut)
+
    
  RETURN
Index: /LMDZ6/trunk/libf/phylmd/cdrag_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmd/cdrag_mod.F90	(revision 4478)
@@ -16,4 +16,5 @@
      speed, t1,    q1,    zgeop1, &
      psol,  tsurf, qsurf, z0m, z0h,  &
+     ri_in, iri_in, &
      cdm,  cdh,  zri,   pref)
 
@@ -22,4 +23,5 @@
   USE print_control_mod, ONLY: lunout, prt_level
   USE ioipsl_getin_p_mod, ONLY : getin_p
+  USE atke_turbulence_ini_mod, ONLY : ric, cn, cinf, cepsilon, pr_slope, pr_asym, pr_neut
 
   IMPLICIT NONE
@@ -84,4 +86,6 @@
   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)        :: 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)
@@ -123,4 +127,6 @@
   REAL, DIMENSION(klon) :: cdmn, cdhn               ! Drag coefficient in neutral conditions
   REAL zzzcd
+  REAL, DIMENSION(klon) :: sm, prandtl              ! Stability function from atke turbulence scheme
+  REAL                  ri0, ri1                    ! to have dimensionless term in sm and prandtl
 
   LOGICAL, SAVE :: firstcall = .TRUE.
@@ -162,5 +168,5 @@
 
   ALPHA=5.0
-  
+
 
 ! ================================================================= c
@@ -224,4 +230,7 @@
           *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
      zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
+     IF (iri_in.EQ.1) THEN
+       zri(i) = ri_in(i)
+     ENDIF
 
 
@@ -280,5 +289,12 @@
            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)
+         
+           ri0=2./rpi*(cinf - cn)*ric/cn
+           ri1=-2./rpi * (pr_asym - pr_neut)/pr_slope
+           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
@@ -376,7 +392,9 @@
           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
@@ -414,409 +432,3 @@
 END SUBROUTINE cdrag
 
-!
- SUBROUTINE cdragn_ri1(knon,  nsrf,   &
-     speed, t1,    q1,    zgeop1, &
-     psol,  tsurf, qsurf, z0m, z0h,  &
-     ri1, iri1, &
-     cdm,  cdh,  zri,   pref)
-
-  USE dimphy
-  USE indice_sol_mod
-  USE print_control_mod, ONLY: lunout, prt_level
-  USE ioipsl_getin_p_mod, ONLY : getin_p
-
-  IMPLICIT NONE
-! ================================================================= c
-!
-! Objet : calcul des cdrags pour le moment (pcfm) et 
-!         les flux de chaleur sensible et latente (pcfh) d'apr??s
-!         Louis 1982, Louis 1979, King et al 2001
-!         ou Zilitinkevich et al 2002  pour les cas stables, Louis 1979
-!         et 1982 pour les cas instables 
-!
-! Modified history: 
-!  writting on the 20/05/2016
-!  modified on the 13/12/2016 to be adapted to LMDZ6
-!
-! References:
-!   Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
-!     atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202. 
-!   Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
-!     operational PBL parametrization at ECMWF'. Workshop on boundary layer
-!     parametrization, November 1981, ECMWF, Reading, England. 
-!     Page: 19. Equations in Table 1.
-!   Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS 
-!    IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
-!    Boundary-Layer Meteorology 72: 331-344
-!   Anton Beljaars. May 1992. The parametrization of the planetary boundary layer.  
-!     European Centre for Medium-Range Weather Forecasts.
-!     Equations: 110-113. Page 40.
-!   Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
-!     model to the parameterization of evaporation from the tropical oceans. J.
-!     Climate, 5:418-434.
-!   King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
-!   to surface and boundary-layer flux parametrizations
-!   QJRMS, 127, pp 779-794
-!
-! ================================================================= c
-! ================================================================= c
-! On choisit le couple de fonctions de correction avec deux flags:
-! Un pour les cas instables, un autre pour les cas stables
-!
-! iflag_corr_insta:
-!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
-!                   2: Louis 1982 
-!                   3: Laurent Li
-!
-! iflag_corr_sta:
-!                   1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
-!                   2: Louis 1982 
-!                   3: Laurent Li
-!                   4: King  2001 (SHARP)
-!                   5: MO 1st order theory (allow collapse of turbulence)
-!            
-!
-!*****************************************************************
-! Parametres d'entree
-!***************************************************************** 
-  
-  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)        :: ri1   ! Richardson 1ere couche 
-  INTEGER, INTENT(IN)                      :: iri1   !
-  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
-
-
-
-! Parametres de sortie
-!******************************************************************
-  REAL, DIMENSION(klon), INTENT(OUT)       :: cdm  ! Drag coefficient for heat flux
-  REAL, DIMENSION(klon), INTENT(OUT)       :: cdh  ! Drag coefficient for momentum
-  REAL, DIMENSION(klon), INTENT(OUT)       :: zri   ! Richardson number
-  REAL, DIMENSION(klon), INTENT(OUT)       :: pref  ! Pression au niveau zgeop/RG
-! Variables Locales
-!******************************************************************
- 
-
-  INCLUDE "YOMCST.h"
-  INCLUDE "YOETHF.h"
-  INCLUDE "clesphys.h"
-
-
-  REAL, PARAMETER       :: CKAP=0.40, CKAPT=0.42
-  REAL                   CEPDU2
-  REAL                   ALPHA
-  REAL                   CB,CC,CD,C2,C3
-  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
-  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
-  REAL zzzcd
-
-  LOGICAL, SAVE :: firstcall = .TRUE.
-  !$OMP THREADPRIVATE(firstcall)
-  INTEGER, SAVE :: iflag_corr_sta
-  !$OMP THREADPRIVATE(iflag_corr_sta)
-  INTEGER, SAVE :: iflag_corr_insta
-  !$OMP THREADPRIVATE(iflag_corr_insta)
-
-!===================================================================c
-! Valeurs numeriques des constantes
-!===================================================================c
-
-
-! Minimum du carre du vent
-
- CEPDU2 = (0.1)**2
-! Louis 1982
-
-  CB=5.0
-  CC=5.0
-  CD=5.0
-
-
-! King 2001
-
-  C2=0.25
-  C3=0.0625
-
-  
-! Louis 1979
-
-  BPRIME=4.7
-  B=9.4
-  
-
-!MO
-
-  ALPHA=5.0
-  
-
-! ================================================================= c
-! Tests avant de commencer
-! Fuxing WANG, 04/03/2015
-! To check if there are negative q1, qsurf values.
-!====================================================================c
-  ng_q1 = 0     ! Initialization
-  ng_qsurf = 0  ! Initialization
-  DO i = 1, knon
-     IF (q1(i).LT.0.0)     ng_q1 = ng_q1 + 1
-     IF (qsurf(i).LT.0.0)  ng_qsurf = ng_qsurf + 1
-  ENDDO
-  IF (ng_q1.GT.0 .and. prt_level > 5) THEN
-      WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
-      WRITE(lunout,*)" The total number of the grids is: ", ng_q1
-      WRITE(lunout,*)" The negative q1 is set to zero "
-!      abort_message="voir ci-dessus"
-!      CALL abort_physic(modname,abort_message,1)
-  ENDIF
-  IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN
-      WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
-      WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf
-      WRITE(lunout,*)" The negative qsurf is set to zero "
-!      abort_message="voir ci-dessus"
-!      CALL abort_physic(modname,abort_message,1)
-  ENDIF
-
-
-
-!=============================================================================c
-! Calcul du cdrag
-!=============================================================================c
-
-! On choisit les fonctions de stabilite utilisees au premier appel
-!**************************************************************************
-  IF (firstcall) THEN
-   iflag_corr_sta=2
-   iflag_corr_insta=2
- 
-   CALL getin_p('iflag_corr_sta',iflag_corr_sta)
-   CALL getin_p('iflag_corr_insta',iflag_corr_insta)
-
-   firstcall = .FALSE.
- ENDIF
-
-!xxxxxxxxxxxxxxxxxxxxxxx
-  DO i = 1, knon        ! Boucle sur l'horizontal
-!xxxxxxxxxxxxxxxxxxxxxxx
-
-
-! calculs preliminaires:
-!***********************
-     
-
-     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 
-     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))) &
-          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
-     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
-
-     IF (iri1.EQ.1) THEN
-       zri(i) = ri1(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 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                          
-
-!'''''''''''''''
-! 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 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
-  END DO   !  Fin de la boucle sur l'horizontal
-!xxxxxxxxxxx
-! ================================================================= c
-     
-END SUBROUTINE cdragn_ri1
-
 END MODULE cdrag_mod
Index: /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 4478)
@@ -876,4 +876,5 @@
     REAL, DIMENSION(klon)              :: yrmu0
     ! Martin
+    REAL, DIMENSIOn(klon)              :: yri0
 
     REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
@@ -1060,4 +1061,5 @@
 
     ytke=0.
+    yri0(:)=0.
 !FC 
     y_treedrg=0.
@@ -1569,5 +1571,5 @@
         CALL cdrag(knon, nsrf, &
             speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
-            yts, yqsurf, yz0m, yz0h, &
+            yts, yqsurf, yz0m, yz0h, yri0, 0, &
             ycdragm, ycdragh, zri1, pref )
 
@@ -1603,5 +1605,5 @@
             CALL cdrag(knon, nsrf, &
             speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
-            yts_x, yqsurf_x, yz0m, yz0h, &
+            yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
             ycdragm_x, ycdragh_x, zri1_x, pref_x )
 
@@ -1630,5 +1632,5 @@
         CALL cdrag(knon, nsrf, &
             speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
-            yts_w, yqsurf_w, yz0m, yz0h, &
+            yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
             ycdragm_w, ycdragh_w, zri1_w, pref_w )
 !
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4478)
@@ -1748,5 +1748,5 @@
        CALL wake_ini(rg,rd,rv,prt_level)
        CALL yamada_ini(klon,lunout,prt_level)
-       CALL atke_ini(prt_level, lunout, RG, RD, RPI)
+       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD)
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
Index: /LMDZ6/trunk/libf/phylmd/screenc_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/screenc_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmd/screenc_mod.F90	(revision 4478)
@@ -70,5 +70,5 @@
 ! Variables locales  
       INTEGER :: i 
-      REAL, dimension(klon) :: cdram, cdrah, cdran, zri1, gref,ycdragm
+      REAL, dimension(klon) :: cdram, cdrah, cdran, zri1, gref,ycdragm,zri_zero
 !
 !------------------------------------------------------------------------- 
@@ -89,4 +89,5 @@
                     speed, temp, q_zref, gref, &
                     psol, ts, qsurf, z0m, z0h, &
+                    zri_zero,0,                &
                     cdram, cdrah, zri1, pref)
       DO i = 1, knon
@@ -177,5 +178,5 @@
 ! Richardson at reference level 
 !
-      CALL cdragn_ri1 (knon, nsrf, &
+      CALL cdrag(knon, nsrf, &
                     speed, temp, q_zref, gref, &
                     psol, ts, qsurf, z0m, z0h, &
Index: /LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmd/stdlevvar_mod.F90	(revision 4478)
@@ -103,5 +103,5 @@
       REAL, dimension(klon) :: te_zref_con, q_zref_con
       REAL, dimension(klon) :: u_zref_c, te_zref_c, temp_c, q_zref_c
-      REAL, dimension(klon) :: ok_pred, ok_corr
+      REAL, dimension(klon) :: ok_pred, ok_corr, zri_zero
 !     REAL, dimension(klon) :: conv_te, conv_q
 !------------------------------------------------------------------------- 
@@ -121,4 +121,5 @@
  &                   speed, t1, q1, z1, &
  &                   psol, ts1, qsurf, z0m, z0h, &
+ &                   zri_zero, 0, &
  &                   cdram, cdrah, zri1, pref)
 
@@ -413,5 +414,5 @@
       REAL, dimension(klon) :: cdrm10m, cdrh10m, ri10m
       REAL, dimension(klon) :: cdmn10m, cdhn10m, fm10m, fh10m
-      REAL, dimension(klon) :: cdn2m, cdn1
+      REAL, dimension(klon) :: cdn2m, cdn1, zri_zero
       REAL :: CEPDUE,zdu2 
       INTEGER :: nzref, nz1
@@ -444,4 +445,5 @@
  &                   speed, t1, q1, z1, &
  &                   psol, ts1, qsurf, z0m, z0h, &
+ &                   zri_zero, 0, &
  &                   cdram, cdrah, zri1, pref)
 
Index: /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 4477)
+++ /LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90	(revision 4478)
@@ -1030,5 +1030,5 @@
     REAL, DIMENSION(klon)              :: yrmu0
     ! Martin
-
+    REAL, DIMENSION(klon)              :: yri0
     REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
          ydser, ydt_ds, ytkt, ytks, ytaur, ysss
@@ -1279,4 +1279,5 @@
 
     ytke=0.
+    yri0(:)=0.
 !FC 
     y_treedrg=0.
@@ -1849,5 +1850,5 @@
         CALL cdrag(knon, nsrf, &
             speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
-            yts, yqsurf, yz0m, yz0h, &
+            yts, yqsurf, yz0m, yz0h, yri0, 0, &
             ycdragm, ycdragh, zri1, pref )
 
@@ -1883,5 +1884,5 @@
             CALL cdrag(knon, nsrf, &
             speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
-            yts_x, yqsurf_x, yz0m, yz0h, &
+            yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
             ycdragm_x, ycdragh_x, zri1_x, pref_x )
 
@@ -1910,5 +1911,5 @@
         CALL cdrag(knon, nsrf, &
             speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
-            yts_w, yqsurf_w, yz0m, yz0h, &
+            yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
             ycdragm_w, ycdragh_w, zri1_w, pref_w )
 !
