Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/cdrag_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/cdrag_mod.F90	(revision 3823)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/cdrag_mod.F90	(revision 3823)
@@ -0,0 +1,822 @@
+!
+MODULE cdrag_mod
+!
+! This module contains some procedures for calculation of the cdrag
+! coefficients for turbulent diffusion at surface
+!
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!
+!r original routine svn3623
+!
+ SUBROUTINE cdrag(knon,  nsrf,   &
+     speed, t1,    q1,    zgeop1, &
+     psol,  tsurf, qsurf, z0m, z0h,  &
+     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)        :: 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 momentum
+  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
+
+! 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)
+
+
+! 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 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/branches/IPSL-CM6A-MR/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/pbl_surface_mod.F90	(revision 3822)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/pbl_surface_mod.F90	(revision 3823)
@@ -23,4 +23,7 @@
   USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
   USE coef_diff_turb_mod,  ONLY : coef_diff_turb
+  USE ioipsl_getin_p_mod,  ONLY : getin_p
+  USE cdrag_mod
+  USE stdlevvar_mod
   USE wx_pbl_mod,          ONLY : wx_pbl_init, wx_pbl_final, &
 !!                                  wx_pbl_fuse_no_dts, wx_pbl_split_no_dts, &
@@ -42,4 +45,6 @@
   INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
   !$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
+  integer, save :: iflag_new_t2mq2m
+  !$OMP THREADPRIVATE(iflag_new_t2mq2m)
 !FC
 !  integer, save :: iflag_frein
@@ -177,5 +182,5 @@
        alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
        alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
-       zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
+       zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout, &
        d_t,       d_q,       d_u,      d_v, d_t_diss, & 
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
@@ -394,4 +399,5 @@
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
+    INTEGER, DIMENSION(klon, 6),  INTENT(OUT)       :: zn2mout ! nb of times t2m is out of the [tsol,temp]
     REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature 
@@ -537,4 +543,6 @@
     REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
     REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
+    INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
+    INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout, n2mout_x, n2mout_w
     REAL, DIMENSION(klon)              :: yustar
     REAL, DIMENSION(klon)              :: ywstar
@@ -838,4 +846,9 @@
 
     IF (first_call) THEN
+
+       iflag_new_t2mq2m=0
+       CALL getin_p('iflag_new_t2mq2m',iflag_new_t2mq2m)
+       print*,'pbl_iflag_new_t2mq2m=',iflag_new_t2mq2m
+
        print*,'PBL SURFACE AVEC GUSTINESS'
        first_call=.FALSE.
@@ -903,4 +916,5 @@
  zxfluxlat(:)=0.
  zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
+ zn2mout(:,:)=0 ;
  d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
  zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
@@ -1821,9 +1835,15 @@
 
           ! Calculate the temperature et relative humidity at 2m and the wind at 10m 
+          IF (iflag_new_t2mq2m==1) THEN
+           CALL stdlevvarn(klon, knon, is_ter, zxli, &
+               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
+               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
+               yt2m, yq2m, yt10m, yq10m, yu10m, yn2mout(:, nsrf, :))
+          ELSE
           CALL stdlevvar(klon, knon, is_ter, zxli, &
                yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
                yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
                yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
-          
+          ENDIF
        ENDIF
 
@@ -2614,9 +2634,26 @@
 !!! jyg le 07/02/2012
        IF (iflag_split .eq.0) THEN
+        IF (iflag_new_t2mq2m==1) THEN
+         CALL stdlevvarn(klon, knon, nsrf, zxli, &
+            uzon, vmer, tair1, qair1, zgeo1, &
+            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
+            yt2m, yq2m, yt10m, yq10m, yu10m, yn2mout(:, nsrf, :))
+        ELSE
         CALL stdlevvar(klon, knon, nsrf, zxli, &
             uzon, vmer, tair1, qair1, zgeo1, &
             tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
             yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
+        ENDIF
        ELSE  !(iflag_split .eq.0)
+        IF (iflag_new_t2mq2m==1) THEN
+         CALL stdlevvarn(klon, knon, nsrf, zxli, &
+            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, yn2mout_x(:, nsrf, :))
+         CALL stdlevvarn(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, yn2mout_w(:, nsrf, :))
+        ELSE
         CALL stdlevvar(klon, knon, nsrf, zxli, &
             uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
@@ -2627,4 +2664,5 @@
             tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
             yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
+        ENDIF
 !!!
        ENDIF  ! (iflag_split .eq.0)
@@ -2640,4 +2678,9 @@
           u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
           v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
+!
+          DO k = 1, 6
+           n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
+          END DO
+!
         ENDDO
        ELSE  !(iflag_split .eq.0)
@@ -2650,4 +2693,9 @@
           u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
           v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
+!
+          DO k = 1, 6
+           n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
+          END DO
+!
         ENDDO
         DO j=1, knon
@@ -2663,4 +2711,9 @@
           u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
           v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
+!
+          DO k = 1, 6
+           n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
+          END DO
+!
         ENDDO
 !!!
@@ -2919,5 +2972,5 @@
 !
     zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
-    zt2m(:) = 0.0    ; zq2m(:) = 0.0 
+    zt2m(:) = 0.0    ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
     zustar(:)=0.0 ; zu10m(:) = 0.0   ; zv10m(:) = 0.0
     s_pblh(:) = 0.0  ; s_plcl(:) = 0.0 
@@ -2972,4 +3025,9 @@
           zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
           zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
+!
+          DO k = 1, 6
+           zn2mout(i,k)  = zn2mout(i,k)  + n2mout(i,nsrf,k)  * pctsrf(i,nsrf)
+          ENDDO
+!
           zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
           wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
@@ -3003,4 +3061,9 @@
           zt2m(i)  = zt2m(i)  + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)
           zq2m(i)  = zq2m(i)  + q2m_x(i,nsrf)  * pctsrf(i,nsrf)
+!
+          DO k = 1, 6
+           zn2mout(i,k)  = zn2mout(i,k)  + n2mout_x(i,nsrf,k)  * pctsrf(i,nsrf)
+          ENDDO
+!
           zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
           wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
@@ -3081,5 +3144,5 @@
     DO nsrf = 1, nbsrf
        DO i = 1, klon
-          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
+          zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
           zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
        ENDDO
Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_local_var_mod.F90	(revision 3822)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_local_var_mod.F90	(revision 3823)
@@ -380,10 +380,4 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: t2m_min_mon, t2m_max_mon
 !$OMP THREADPRIVATE(t2m_min_mon, t2m_max_mon)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zq2m_cor, zt2m_cor
-!$OMP THREADPRIVATE(zq2m_cor, zt2m_cor)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zu10m_cor, zv10m_cor
-!$OMP THREADPRIVATE(zu10m_cor, zv10m_cor)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zrh2m_cor, zqsat2m_cor
-!$OMP THREADPRIVATE(zrh2m_cor, zqsat2m_cor)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: weak_inversion
 !$OMP THREADPRIVATE(weak_inversion)
@@ -475,4 +469,7 @@
       REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: p_tropopause, z_tropopause, t_tropopause
 !$OMP THREADPRIVATE(p_tropopause, z_tropopause, t_tropopause)
+
+      INTEGER,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zn2mout
+!$OMP THREADPRIVATE(zn2mout)
 
 #ifdef CPP_StratAer
@@ -760,6 +757,4 @@
       ALLOCATE(zt2m_min_mon(klon), zt2m_max_mon(klon))
       ALLOCATE(t2m_min_mon(klon), t2m_max_mon(klon))
-      ALLOCATE(zq2m_cor(klon), zt2m_cor(klon), zu10m_cor(klon), zv10m_cor(klon))
-      ALLOCATE(zrh2m_cor(klon), zqsat2m_cor(klon))
       ALLOCATE(sens(klon), flwp(klon), fiwp(klon))
       ALLOCATE(alp_bl_conv(klon), alp_bl_det(klon))
@@ -831,4 +826,6 @@
       ALLOCATE (z_tropopause(klon))
       ALLOCATE (t_tropopause(klon))
+
+      ALLOCATE(zn2mout(klon,6))
 
 #ifdef CPP_StratAer
@@ -1057,6 +1054,4 @@
       DEALLOCATE(zt2m_min_mon, zt2m_max_mon)
       DEALLOCATE(t2m_min_mon, t2m_max_mon)
-      DEALLOCATE(zq2m_cor, zt2m_cor, zu10m_cor, zv10m_cor)
-      DEALLOCATE(zrh2m_cor, zqsat2m_cor)
       DEALLOCATE(sens, flwp, fiwp)
       DEALLOCATE(alp_bl_conv,alp_bl_det)
@@ -1119,4 +1114,6 @@
       DEALLOCATE (z_tropopause)
       DEALLOCATE (t_tropopause)
+
+      DEALLOCATE(zn2mout)
 
 #ifdef CPP_StratAer
Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3822)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3823)
@@ -271,4 +271,18 @@
       ctrl_out((/ 10, 6, 10, 10, 10, 10, 11, 11, 11, 11/), &
     't2m_sic', "Temp 2m "//clnsurf(4), "K", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE :: o_nt2mout = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
+    'nt2mout', 'Nbt2m out of range complete computation', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_nq2mout = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
+    'nq2mout', 'Nbq2m out of range complete computation', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_nu2mout = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
+    'nu2mout', 'Nbu2m out of range complete computation', '-', (/ ('', i=1, 10) /))
+
+  TYPE(ctrl_out), SAVE :: o_nt2moutfg = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
+    'nt2moutfg', 'Nbt2m out of range complete/fgRi1 computation', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_nq2moutfg = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
+    'nq2moutfg', 'Nbq2m out of range complete/fgRi1 computation', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_nu2moutfg = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
+    'nu2moutfg', 'Nbu2m out of range complete/fgRi1 computation', '-', (/ ('', i=1, 10) /))
 
   TYPE(ctrl_out), SAVE :: o_gusts = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11, 11/), &
Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_output_write_mod.F90	(revision 3822)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/phys_output_write_mod.F90	(revision 3823)
@@ -38,4 +38,7 @@
          o_t2m, o_t2m_min, o_t2m_max, &
          o_t2m_min_mon, o_t2m_max_mon, &
+         o_nt2mout, o_nt2moutfg, &
+         o_nq2mout, o_nq2moutfg, &
+         o_nu2mout, o_nu2moutfg, &
          o_q2m, o_ustar, o_u10m, o_v10m, &
          o_wind10m, o_wind10max, o_wind100m, o_gusts, o_sicf, &
@@ -249,6 +252,5 @@
 
     USE phys_local_var_mod, ONLY: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, &
-         zt2m_cor,zq2m_cor,zu10m_cor,zv10m_cor, zrh2m_cor, zqsat2m_cor, &
-         t2m_min_mon, t2m_max_mon, evap, &
+         zn2mout, t2m_min_mon, t2m_max_mon, evap, &
          l_mixmin,l_mix, &
          zu10m, zv10m, zq2m, zustar, zxqsurf, &
@@ -676,7 +678,7 @@
        CALL histwrite_phy(o_slp, slp)
        CALL histwrite_phy(o_tsol, zxtsol)
-       CALL histwrite_phy(o_t2m, zt2m_cor)
-       CALL histwrite_phy(o_t2m_min, zt2m_cor)
-       CALL histwrite_phy(o_t2m_max, zt2m_cor)
+       CALL histwrite_phy(o_t2m, zt2m)
+       CALL histwrite_phy(o_t2m_min, zt2m)
+       CALL histwrite_phy(o_t2m_max, zt2m)
        CALL histwrite_phy(o_t2m_max_mon, t2m_max_mon)
        CALL histwrite_phy(o_t2m_min_mon, t2m_min_mon)
@@ -684,5 +686,47 @@
        IF (vars_defined) THEN
           DO i=1, klon
-             zx_tmp_fi2d(i)=SQRT(zu10m_cor(i)*zu10m_cor(i)+zv10m_cor(i)*zv10m_cor(i))
+             zx_tmp_fi2d(i)=real(zn2mout(i,1))
+          ENDDO
+       ENDIF
+       CALL histwrite_phy(o_nt2mout, zx_tmp_fi2d)
+
+       IF (vars_defined) THEN
+          DO i=1, klon
+             zx_tmp_fi2d(i)=real(zn2mout(i,2))
+          ENDDO
+       ENDIF
+       CALL histwrite_phy(o_nt2moutfg, zx_tmp_fi2d)
+
+       IF (vars_defined) THEN
+          DO i=1, klon
+             zx_tmp_fi2d(i)=real(zn2mout(i,3))
+          ENDDO
+       ENDIF
+       CALL histwrite_phy(o_nq2mout, zx_tmp_fi2d)
+
+       IF (vars_defined) THEN
+          DO i=1, klon
+             zx_tmp_fi2d(i)=real(zn2mout(i,4))
+          ENDDO
+       ENDIF
+       CALL histwrite_phy(o_nq2moutfg, zx_tmp_fi2d)
+
+       IF (vars_defined) THEN
+          DO i=1, klon
+             zx_tmp_fi2d(i)=real(zn2mout(i,5))
+          ENDDO
+       ENDIF
+       CALL histwrite_phy(o_nu2mout, zx_tmp_fi2d)
+
+       IF (vars_defined) THEN
+          DO i=1, klon
+             zx_tmp_fi2d(i)=real(zn2mout(i,6))
+          ENDDO
+       ENDIF
+       CALL histwrite_phy(o_nu2moutfg, zx_tmp_fi2d)
+
+       IF (vars_defined) THEN
+          DO i=1, klon
+             zx_tmp_fi2d(i)=SQRT(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
           ENDDO
        ENDIF
@@ -691,5 +735,5 @@
        IF (vars_defined) THEN
           DO i=1, klon
-             zx_tmp_fi2d(i)=SQRT(zu10m_cor(i)*zu10m_cor(i)+zv10m_cor(i)*zv10m_cor(i))
+             zx_tmp_fi2d(i)=SQRT(zu10m(i)*zu10m(i)+zv10m(i)*zv10m(i))
           ENDDO
        ENDIF
@@ -771,8 +815,9 @@
        ENDIF
        CALL histwrite_phy(o_sicf, zx_tmp_fi2d)
-       CALL histwrite_phy(o_q2m, zq2m_cor)
-       CALL histwrite_phy(o_ustar, zustar)
-       CALL histwrite_phy(o_u10m, zu10m_cor)
-       CALL histwrite_phy(o_v10m, zv10m_cor)
+       CALL histwrite_phy(o_q2m, zq2m)
+       IF (vars_defined) zx_tmp_fi2d = zustar
+       CALL histwrite_phy(o_ustar, zx_tmp_fi2d)
+       CALL histwrite_phy(o_u10m, zu10m)
+       CALL histwrite_phy(o_v10m, zv10m)
 
        IF (vars_defined) THEN
@@ -1412,5 +1457,5 @@
 !       CALL histwrite_phy(o_rh2m_max, zx_tmp_fi2d)
 
-       CALL histwrite_phy(o_qsat2m, zqsat2m_cor)
+       CALL histwrite_phy(o_qsat2m, qsat2m)
        CALL histwrite_phy(o_tpot, tpot)
        CALL histwrite_phy(o_tpote, tpote)
Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/physiq_mod.F90	(revision 3822)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/physiq_mod.F90	(revision 3823)
@@ -119,7 +119,5 @@
        cdragm, cdragh,                   &
        zustar, zu10m, zv10m, rh2m, qsat2m, &
-       zq2m, zt2m, weak_inversion, &
-       zq2m_cor,zt2m_cor,zu10m_cor,zv10m_cor, & ! pour corriger d'un bug
-       zrh2m_cor,zqsat2m_cor, &
+       zq2m, zt2m, zn2mout, weak_inversion, &
        zt2m_min_mon, zt2m_max_mon,   &         ! pour calcul_divers.h
        t2m_min_mon, t2m_max_mon,  &            ! pour calcul_divers.h
@@ -603,9 +601,4 @@
     !$OMP THREADPRIVATE(iflag_alp_wk_cond)
 
-    INTEGER,  SAVE               :: iflag_bug_t2m_ipslcm61=1 !
-    !$OMP THREADPRIVATE(iflag_bug_t2m_ipslcm61)
-    INTEGER,  SAVE               :: iflag_bug_t2m_stab_ipslcm61=-1 !
-    !$OMP THREADPRIVATE(iflag_bug_t2m_stab_ipslcm61)
-
     REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
     REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
@@ -1314,9 +1307,4 @@
        tau_gl=86400.*tau_gl
        WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
-
-       iflag_bug_t2m_ipslcm61 = 1
-       CALL getin_p('iflag_bug_t2m_ipslcm61', iflag_bug_t2m_ipslcm61)
-       iflag_bug_t2m_stab_ipslcm61 = -1
-       CALL getin_p('iflag_bug_t2m_stab_ipslcm61', iflag_bug_t2m_stab_ipslcm61)
 
        CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond) 
@@ -2529,5 +2517,5 @@
                                 !albedo SB <<<
             albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
-            zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
+            zxtsol,    zxfluxlat, zt2m,    qsat2m, zn2mout, &
             d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf, d_t_diss, &
                                 !nrlmd<
@@ -2576,21 +2564,4 @@
 !>jyg
        ENDIF
-
-!add limitation for t,q at and wind at 10m
-        if ( iflag_bug_t2m_ipslcm61 == 0 ) THEN
-          CALL borne_var_surf( klon,klev,nbsrf,                 &
-            iflag_bug_t2m_stab_ipslcm61,                        &
-            t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),    &
-            ftsol,zxqsurf,pctsrf,paprs,                         &
-            t2m, q2m, u10m, v10m,                               &
-            zt2m_cor, zq2m_cor, zu10m_cor, zv10m_cor,           &
-            zrh2m_cor, zqsat2m_cor)
-        ELSE
-          zt2m_cor(:)=zt2m(:)
-          zq2m_cor(:)=zq2m(:)
-          zu10m_cor(:)=zu10m(:)
-          zv10m_cor(:)=zv10m(:)
-          zqsat2m_cor=999.999
-        ENDIF
 
        !---------------------------------------------------------------------
Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/screenc_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/screenc_mod.F90	(revision 3823)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/screenc_mod.F90	(revision 3823)
@@ -0,0 +1,193 @@
+!
+MODULE screenc_mod
+!
+! This module contains some procedures for calculation of the correction
+! of temperature, specific humidity and wind at a reference level
+!
+  USE cdrag_mod
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!
+!r original routine svn3623
+!
+      SUBROUTINE screenc(klon, knon, nsrf, zxli, &
+                         speed, temp, q_zref, zref, &
+                         ts, qsurf, z0m, z0h, psol, &
+                         ustar, testar, qstar, okri, ri1, &
+                         pref, delu, delte, delq)
+      IMPLICIT NONE
+!-----------------------------------------------------------------------
+! 
+! Objet : calcul "correcteur" des anomalies du vent, de la temperature 
+!         potentielle et de l'humidite relative au niveau de reference zref et 
+!         par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) 
+!         a partir des equations de Louis.
+!
+! Reference : Hess, Colman et McAvaney (1995)
+!
+! I. Musat, 01.07.2002
+!-----------------------------------------------------------------------
+!
+! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
+! knon----input-I- nombre de points pour un type de surface
+! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
+! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
+! speed---input-R- module du vent au 1er niveau du modele
+! temp----input-R- temperature de l'air au 1er niveau du modele
+! q_zref--input-R- humidite relative au 1er niveau du modele
+! zref----input-R- altitude de reference
+! ts------input-R- temperature de l'air a la surface
+! qsurf---input-R- humidite relative a la surface
+! z0m, z0h---input-R- rugosite
+! psol----input-R- pression au sol
+! ustar---input-R- facteur d'echelle pour le vent
+! testar--input-R- facteur d'echelle pour la temperature potentielle
+! qstar---input-R- facteur d'echelle pour l'humidite relative
+! okri----input-L- TRUE si on veut tester le nb. Richardson entre la sfce 
+!                  et zref par rapport au Ri entre la sfce et la 1ere couche
+! ri1-----input-R- nb. Richardson entre la surface et la 1ere couche 
+!
+! pref----input-R- pression au niveau de reference
+! delu----input-R- anomalie du vent par rapport au 1er niveau
+! delte---input-R- anomalie de la temperature potentielle par rapport a la surface
+! delq----input-R- anomalie de l'humidite relative par rapport a la surface
+!
+      INTEGER, intent(in) :: klon, knon, nsrf
+      LOGICAL, intent(in) :: zxli, okri 
+      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(out) :: pref, delu, delte, delq 
+!-----------------------------------------------------------------------
+      include "YOMCST.h"
+      include "flux_arp.h"
+!
+! Variables locales  
+      INTEGER :: i 
+      REAL, dimension(klon) :: cdram, cdrah, cdran, zri1, gref,ycdragm
+!
+!------------------------------------------------------------------------- 
+      DO i=1, knon
+        gref(i) = zref*RG
+      ENDDO 
+!
+! Richardson at reference level 
+!
+!      CALL coefcdrag (klon, knon, nsrf, zxli, &
+!                    speed, temp, q_zref, gref, &
+!                    psol, ts, qsurf, rugos, &
+!                    okri, ri1, &
+!                    cdram, cdrah, cdran, zri1, &
+!                    pref)
+! Fuxing WANG, 04/03/2015, replace the coefcdrag by the merged version: cdrag
+      CALL cdrag (knon, nsrf, &
+                    speed, temp, q_zref, gref, &
+                    psol, ts, qsurf, z0m, z0h, &
+                    cdram, cdrah, zri1, pref)
+      DO i = 1, knon
+        IF(ok_prescr_ust) THEN
+! La aussi il faut forcer avec ust (FC + MPL 20160210)
+        ycdragm(i) = ust*ust/(1.+speed(i))/speed(i)
+        cdram=ycdragm
+        delu(i) = ust/sqrt(cdram(i))
+        ELSE
+        delu(i) = ustar(i)/sqrt(cdram(i))
+        ENDIF
+        delte(i)= (testar(i)* sqrt(cdram(i)))/ &
+                   cdrah(i)
+        delq(i)= (qstar(i)* sqrt(cdram(i)))/ &
+                  cdrah(i)
+      ENDDO 
+!
+      RETURN 
+      END SUBROUTINE screenc
+!
+      SUBROUTINE screencn(klon, knon, nsrf, zxli, &
+                         speed, temp, q_zref, zref, &
+                         ts, qsurf, z0m, z0h, psol, &
+                         cdrm, cdrh,  okri, &
+                         ri1, iri1, &
+                         pref, delm, delh, zri1)
+      IMPLICIT NONE
+!-----------------------------------------------------------------------
+! 
+! Objet : calcul "correcteur" des anomalies du vent, de la temperature 
+!         potentielle et de l'humidite relative au niveau de reference zref et 
+!         par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) 
+!         a partir des equations de Louis.
+!
+! Reference : Hess, Colman et McAvaney (1995)
+!
+! I. Musat, 01.07.2002
+!-----------------------------------------------------------------------
+!
+! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
+! knon----input-I- nombre de points pour un type de surface
+! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
+! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
+! speed---input-R- module du vent au 1er niveau du modele
+! temp----input-R- temperature de l'air au 1er niveau du modele
+! q_zref--input-R- humidite relative au 1er niveau du modele
+! zref----input-R- altitude de reference
+! ts------input-R- temperature de l'air a la surface
+! qsurf---input-R- humidite relative a la surface
+! z0m, z0h---input-R- rugosite
+! psol----input-R- pression au sol
+! ustar---input-R- facteur d'echelle pour le vent
+! testar--input-R- facteur d'echelle pour la temperature potentielle
+! qstar---input-R- facteur d'echelle pour l'humidite relative
+! okri----input-L- TRUE si on veut tester le nb. Richardson entre la sfce 
+!                  et zref par rapport au Ri entre la sfce et la 1ere couche
+! ri1-----input-R- nb. Richardson entre la surface et la 1ere couche 
+!
+! pref----input-R- pression au niveau de reference
+! delu----input-R- anomalie du vent par rapport au 1er niveau
+! delte---input-R- anomalie de la temperature potentielle par rapport a la surface
+! delq----input-R- anomalie de l'humidite relative par rapport a la surface
+!
+      INTEGER, intent(in) :: klon, knon, nsrf
+      LOGICAL, intent(in) :: zxli, okri 
+      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) :: cdrm, cdrh, ri1         
+      INTEGER, INTENT(IN)  :: iri1 ! Richardson de la 1ere couche
+!
+      REAL, dimension(klon), intent(out) :: pref, delm, delh, zri1
+!171220     REAL, dimension(klon) :: cdram, cdrah, zri1 
+      REAL, dimension(klon) :: cdram, cdrah
+!-----------------------------------------------------------------------
+      include "YOMCST.h"
+      include "flux_arp.h"
+!
+! Variables locales  
+      INTEGER :: i 
+      REAL, dimension(klon) :: cdran, gref,ycdragm
+!
+!------------------------------------------------------------------------- 
+      DO i=1, knon
+        gref(i) = zref*RG
+      ENDDO 
+!
+! Richardson at reference level 
+!
+      CALL cdragn_ri1 (knon, nsrf, &
+                    speed, temp, q_zref, gref, &
+                    psol, ts, qsurf, z0m, z0h, &
+                    ri1, iri1, &
+                    cdram, cdrah, zri1, pref)
+      DO i = 1, knon
+        delm(i) = sqrt(cdrm(i))/sqrt(cdram(i))
+!verifier que temp est un temperat potentielle.
+        delh(i)= (cdrh(i)* sqrt(cdram(i)))/ &
+                   (cdrah(i)*sqrt(cdrm(i)))
+      ENDDO 
+!
+      RETURN 
+      END SUBROUTINE screencn
+END MODULE screenc_mod
Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/screenp_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/screenp_mod.F90	(revision 3823)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/screenp_mod.F90	(revision 3823)
@@ -0,0 +1,185 @@
+!
+MODULE screenp_mod
+!
+! This module contains some procedures for calculation of the first
+! guess of temperature, specific humidity and wind at a reference level
+! coefficients for turbulent diffusion at surface
+!
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!
+!r original routine
+!
+      SUBROUTINE screenp(klon, knon, nsrf, &
+     &                   speed, tair, qair, &
+     &                   ts, qsurf, rugos, lmon, &
+     &                   ustar, testar, qstar, zref, &
+     &                   delu, delte, delq) 
+      IMPLICIT none
+!-------------------------------------------------------------------------
+!
+! Objet : calcul "predicteur" des anomalies du vent, de la temperature 
+!         potentielle et de l'humidite relative au niveau de reference zref et 
+!         par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) 
+!         a partir des relations de Dyer-Businger.
+!
+! Reference : Hess, Colman et McAvaney (1995)
+!
+! I. Musat, 01.07.2002
+!-------------------------------------------------------------------------
+!
+! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
+! knon----input-I- nombre de points pour un type de surface
+! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
+! speed---input-R- module du vent au 1er niveau du modele
+! tair----input-R- temperature de l'air au 1er niveau du modele
+! qair----input-R- humidite relative au 1er niveau du modele
+! ts------input-R- temperature de l'air a la surface
+! qsurf---input-R- humidite relative a la surface
+! rugos---input-R- rugosite
+! lmon----input-R- longueur de Monin-Obukov
+! ustar---input-R- facteur d'echelle pour le vent
+! testar--input-R- facteur d'echelle pour la temperature potentielle
+! qstar---input-R- facteur d'echelle pour l'humidite relative
+! zref----input-R- altitude de reference
+!
+! delu----input-R- anomalie du vent par rapport au 1er niveau
+! delte---input-R- anomalie de la temperature potentielle par rapport a la surface
+! delq----input-R- anomalie de l'humidite relative par rapport a la surface
+!
+      INTEGER, intent(in) :: klon, knon, nsrf
+      REAL, dimension(klon), intent(in) :: speed, tair, qair
+      REAL, dimension(klon), intent(in) :: ts, qsurf, rugos
+      DOUBLE PRECISION, dimension(klon), intent(in) :: lmon
+      REAL, dimension(klon), intent(in) :: ustar, testar, qstar
+      REAL, intent(in) :: zref
+!
+      REAL, dimension(klon), intent(out) :: delu, delte, delq
+!
+!-------------------------------------------------------------------------
+! Variables locales et constantes :
+      REAL, PARAMETER :: RKAR=0.40
+      INTEGER :: i
+      REAL :: xtmp, xtmp0
+!-------------------------------------------------------------------------
+      DO i = 1, knon
+!
+        IF (lmon(i).GE.0.) THEN
+!
+! STABLE CASE
+!
+          IF (speed(i).GT.1.5.AND.lmon(i).LE.1.0                        &
+     &                      .AND. rugos(i).LE.1.0) THEN
+            delu(i) = (ustar(i)/RKAR)* &
+                      (log(zref/(rugos(i))+1.) + &
+                      min(5.d0, 5.0 *(zref - rugos(i))/lmon(i)))
+            delte(i) = (testar(i)/RKAR)* &
+                       (log(zref/(rugos(i))+1.) + &
+                       min(5.d0, 5.0 * (zref - rugos(i))/lmon(i)))
+            delq(i) = (qstar(i)/RKAR)* &
+                      (log(zref/(rugos(i))+1.) + &
+                      min(5.d0, 5.0 * (zref - rugos(i))/lmon(i)))
+          ELSE
+            delu(i)  = 0.1 * speed(i)
+            delte(i) = 0.1 * (tair(i) - ts(i) )
+            delq(i)  = 0.1 * (max(qair(i),0.0) - max(qsurf(i),0.0))
+          ENDIF
+        ELSE  
+!
+! UNSTABLE CASE
+!
+          IF (speed(i).GT.5.0.AND.abs(lmon(i)).LE.50.0) THEN
+            xtmp = (1. - 16. * (zref/lmon(i)))**(1./4.)
+            xtmp0 = (1. - 16. * (rugos(i)/lmon(i)))**(1./4.)
+            delu(i) = (ustar(i)/RKAR)* &
+                      (log(zref/(rugos(i))+1.) & 
+                      - 2.*log(0.5*(1. + xtmp)) &
+                      + 2.*log(0.5*(1. + xtmp0)) &
+                      - log(0.5*(1. + xtmp*xtmp)) &
+                      + log(0.5*(1. + xtmp0*xtmp0)) &
+                      + 2.*atan(xtmp) - 2.*atan(xtmp0))
+            delte(i) = (testar(i)/RKAR)* &
+                       (log(zref/(rugos(i))+1.) &
+                       - 2.0 * log(0.5*(1. + xtmp*xtmp)) & 
+                       + 2.0 * log(0.5*(1. + xtmp0*xtmp0)))
+            delq(i)  = (qstar(i)/RKAR)* &
+                       (log(zref/(rugos(i))+1.) &
+                       - 2.0 * log(0.5*(1. + xtmp*xtmp)) & 
+                       + 2.0 * log(0.5*(1. + xtmp0*xtmp0)))
+          ELSE
+            delu(i)  = 0.5 * speed(i)
+            delte(i) = 0.5 * (tair(i) - ts(i) )
+            delq(i)  = 0.5 * (max(qair(i),0.0) - max(qsurf(i),0.0))
+          ENDIF
+        ENDIF
+!
+      ENDDO
+      RETURN
+      END SUBROUTINE screenp
+!
+      SUBROUTINE screenpn(klon, knon, nsrf, &
+     &                   speed, tair, qair, &
+     &                   ts, qsurf, rugos, zri1, &
+     &                   zref, &
+     &                   delu, delte, delq) 
+      IMPLICIT none
+!-------------------------------------------------------------------------
+!
+! Objet : calcul "predicteur" des anomalies du vent, de la temperature 
+!         potentielle et de l'humidite relative au niveau de reference zref et 
+!         par rapport au 1er niveau (pour u) ou a la surface (pour theta et q) 
+!         a partir des relations de Dyer-Businger.
+!
+! Reference : Hess, Colman et McAvaney (1995)
+!
+! I. Musat, 01.07.2002
+!-------------------------------------------------------------------------
+!
+! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
+! knon----input-I- nombre de points pour un type de surface
+! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
+! speed---input-R- module du vent au 1er niveau du modele
+! tair----input-R- temperature de l'air au 1er niveau du modele
+! qair----input-R- humidite relative au 1er niveau du modele
+! ts------input-R- temperature de l'air a la surface
+! qsurf---input-R- humidite relative a la surface
+! rugos---input-R- rugosite
+! zref----input-R- altitude de reference
+!
+! delu----input-R- anomalie du vent par rapport au 1er niveau
+! delte---input-R- anomalie de la temperature potentielle par rapport a la surface
+! delq----input-R- anomalie de l'humidite relative par rapport a la surface
+!
+      INTEGER, intent(in) :: klon, knon, nsrf
+      REAL, dimension(klon), intent(in) :: speed, tair, qair
+      REAL, dimension(klon), intent(in) :: ts, qsurf, rugos
+      REAL, dimension(klon), intent(in) :: zri1
+      REAL, intent(in) :: zref
+!
+      REAL, dimension(klon), intent(out) :: delu, delte, delq
+!
+!-------------------------------------------------------------------------
+! Variables locales et constantes :
+      REAL, PARAMETER :: RKAR=0.40
+      INTEGER :: i
+      REAL :: xtmp, xtmp0
+!-------------------------------------------------------------------------
+      DO i = 1, knon
+!
+       IF (zri1(i).GE.0.) THEN
+          delu(i)  = 0.1 * speed(i)
+          delte(i) = 0.1 * (tair(i) - ts(i) )
+          delq(i)  = 0.1 * (max(qair(i),0.0) - max(qsurf(i),0.0))
+       ELSE
+          delu(i)  = 0.5 * speed(i)
+          delte(i) = 0.5 * (tair(i) - ts(i) )
+          delq(i)  = 0.5 * (max(qair(i),0.0) - max(qsurf(i),0.0))
+       ENDIF
+!
+      ENDDO
+      RETURN
+      END SUBROUTINE screenpn
+END MODULE screenp_mod 
Index: LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/stdlevvar_mod.F90
===================================================================
--- LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/stdlevvar_mod.F90	(revision 3823)
+++ LMDZ6/branches/IPSL-CM6A-MR/libf/phylmd/stdlevvar_mod.F90	(revision 3823)
@@ -0,0 +1,727 @@
+!
+MODULE stdlevvar_mod
+!
+! This module contains main procedures for calculation
+! of temperature, specific humidity and wind at a reference level
+!
+  USE cdrag_mod
+  USE screenp_mod
+  USE screenc_mod
+  IMPLICIT NONE
+
+CONTAINS
+!
+!****************************************************************************************
+!
+!r original routine svn3623
+!
+      SUBROUTINE stdlevvar(klon, knon, nsrf, zxli, &
+                           u1, v1, t1, q1, z1, &
+                           ts1, qsurf, z0m, z0h, psol, pat1, &
+                           t_2m, q_2m, t_10m, q_10m, u_10m, ustar)
+      IMPLICIT NONE
+!-------------------------------------------------------------------------
+!
+! Objet : calcul de la temperature et l'humidite relative a 2m et du 
+!         module du vent a 10m a partir des relations de Dyer-Businger et
+!         des equations de Louis.
+!
+! Reference : Hess, Colman et McAvaney (1995)        
+!
+! I. Musat, 01.07.2002
+!
+!AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
+!
+!-------------------------------------------------------------------------
+!
+! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
+! knon----input-I- nombre de points pour un type de surface
+! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
+! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
+! u1------input-R- vent zonal au 1er niveau du modele
+! v1------input-R- vent meridien au 1er niveau du modele
+! t1------input-R- temperature de l'air au 1er niveau du modele
+! q1------input-R- humidite relative au 1er niveau du modele
+! z1------input-R- geopotentiel au 1er niveau du modele
+! ts1-----input-R- temperature de l'air a la surface
+! qsurf---input-R- humidite relative a la surface
+! z0m, z0h---input-R- rugosite
+! psol----input-R- pression au sol
+! pat1----input-R- pression au 1er niveau du modele
+!
+! t_2m---output-R- temperature de l'air a 2m
+! q_2m---output-R- humidite relative a 2m
+! u_10m--output-R- vitesse du vent a 10m
+!AM
+! t_10m--output-R- temperature de l'air a 10m
+! q_10m--output-R- humidite specifique a 10m
+! ustar--output-R- u*
+!
+      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) :: psol, pat1
+!
+      REAL, dimension(klon), intent(out) :: t_2m, q_2m, ustar
+      REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m
+!-------------------------------------------------------------------------
+      include "flux_arp.h"
+      include "YOMCST.h"
+!IM PLUS
+      include "YOETHF.h"
+!
+! Quelques constantes et options:
+!
+! RKAR : constante de von Karman
+      REAL, PARAMETER :: RKAR=0.40
+! niter : nombre iterations calcul "corrector"
+!     INTEGER, parameter :: niter=6, ncon=niter-1
+      INTEGER, parameter :: niter=2, ncon=niter-1
+!
+! Variables locales
+      INTEGER :: i, n
+      REAL :: zref
+      REAL, dimension(klon) :: speed
+! tpot : temperature potentielle
+      REAL, dimension(klon) :: tpot
+      REAL, dimension(klon) :: zri1, cdran
+      REAL, dimension(klon) :: cdram, cdrah
+! ri1 : nb. de Richardson entre la surface --> la 1ere couche
+      REAL, dimension(klon) :: ri1 
+      REAL, dimension(klon) :: testar, qstar
+      REAL, dimension(klon) :: zdte, zdq   
+! lmon : longueur de Monin-Obukhov selon Hess, Colman and McAvaney 
+      DOUBLE PRECISION, dimension(klon) :: lmon
+      DOUBLE PRECISION, parameter :: eps=1.0D-20
+      REAL, dimension(klon) :: delu, delte, delq
+      REAL, dimension(klon) :: u_zref, te_zref, q_zref  
+      REAL, dimension(klon) :: temp, pref
+      LOGICAL :: okri
+      REAL, dimension(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
+!convertgence
+      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) :: conv_te, conv_q
+!------------------------------------------------------------------------- 
+      DO i=1, knon
+       speed(i)=SQRT(u1(i)**2+v1(i)**2)
+       ri1(i) = 0.0
+      ENDDO
+!
+      okri=.FALSE.
+!      CALL coefcdrag(klon, knon, nsrf, zxli, &
+! &                   speed, t1, q1, z1, psol, &
+! &                   ts1, qsurf, rugos, okri, ri1,  &         
+! &                   cdram, cdrah, cdran, zri1, pref)            
+! Fuxing WANG, 04/03/2015, replace the coefcdrag by the merged version: cdrag
+
+      CALL cdrag(knon, nsrf, &
+ &                   speed, t1, q1, z1, &
+ &                   psol, ts1, qsurf, z0m, z0h, &
+ &                   cdram, cdrah, zri1, pref)
+
+! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
+     IF (ok_prescr_ust) then
+      DO i = 1, knon
+       print *,'cdram avant=',cdram(i)
+       cdram(i) = ust*ust/speed(i)/speed(i)
+       print *,'cdram ust speed apres=',cdram(i),ust,speed
+      ENDDO
+     ENDIF
+!
+!---------Star variables----------------------------------------------------
+!
+      DO i = 1, knon
+        ri1(i) = zri1(i)
+        tpot(i) = t1(i)* (psol(i)/pat1(i))**RKAPPA
+        ustar(i) = sqrt(cdram(i) * speed(i) * speed(i))
+        zdte(i) = tpot(i) - ts1(i)
+        zdq(i) = max(q1(i),0.0) - max(qsurf(i),0.0)
+!
+!
+!IM BUG BUG BUG       zdte(i) = max(zdte(i),1.e-10)
+        zdte(i) = sign(max(abs(zdte(i)),1.e-10),zdte(i))
+!
+        testar(i) = (cdrah(i) * zdte(i) * speed(i))/ustar(i)
+        qstar(i) = (cdrah(i) * zdq(i) * speed(i))/ustar(i)
+        lmon(i) = (ustar(i) * ustar(i) * tpot(i))/ &
+ &                (RKAR * RG * testar(i))
+      ENDDO
+!
+!----------First aproximation of variables at zref --------------------------
+      zref = 2.0
+      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
+ &                 ts1, qsurf, z0m, lmon, &
+ &                 ustar, testar, qstar, zref, &
+ &                 delu, delte, delq)
+!
+      DO i = 1, knon
+        u_zref(i) = delu(i)
+        q_zref(i) = max(qsurf(i),0.0) + delq(i)
+        te_zref(i) = ts1(i) + delte(i)
+        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
+        q_zref_p(i) = q_zref(i)
+!       te_zref_p(i) = te_zref(i)
+        temp_p(i) = temp(i)
+      ENDDO
+!
+! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995 
+!
+      DO n = 1, niter
+!
+        okri=.TRUE.
+        CALL screenc(klon, knon, nsrf, zxli, &
+ &                   u_zref, temp, q_zref, zref, &
+ &                   ts1, qsurf, z0m, z0h, psol, &           
+ &                   ustar, testar, qstar, okri, ri1, &
+ &                   pref, delu, delte, delq) 
+!
+        DO i = 1, knon
+          u_zref(i) = delu(i)
+          q_zref(i) = delq(i) + max(qsurf(i),0.0)
+          te_zref(i) = delte(i) + ts1(i) 
+!
+! return to normal temperature
+!
+          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
+!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
+!                 (1 + RVTMP2 * max(q_zref(i),0.0))
+!
+!IM +++
+!         IF(temp(i).GT.350.) THEN
+!           WRITE(*,*) 'temp(i) GT 350 K !!',i,nsrf,temp(i)
+!         ENDIF
+!IM ---
+!
+        IF(n.EQ.ncon) THEN
+          te_zref_con(i) = te_zref(i)
+          q_zref_con(i) = q_zref(i)
+        ENDIF 
+!
+        ENDDO 
+!
+      ENDDO 
+!
+! verifier le critere de convergence : 0.25% pour te_zref et 5% pour qe_zref
+!
+!       DO i = 1, knon
+!         conv_te(i) = (te_zref(i) - te_zref_con(i))/te_zref_con(i)
+!         conv_q(i) = (q_zref(i) - q_zref_con(i))/q_zref_con(i)
+!IM +++
+!         IF(abs(conv_te(i)).GE.0.0025.AND.abs(conv_q(i)).GE.0.05) THEN
+!           PRINT*,'DIV','i=',i,te_zref_con(i),te_zref(i),conv_te(i), &
+!           q_zref_con(i),q_zref(i),conv_q(i)
+!         ENDIF
+!IM ---
+!       ENDDO
+!
+      DO i = 1, knon
+        q_zref_c(i) = q_zref(i)
+        temp_c(i) = temp(i)
+!
+!       IF(zri1(i).LT.0.) THEN
+!         IF(nsrf.EQ.1) THEN
+!           ok_pred(i)=1.
+!           ok_corr(i)=0.
+!         ELSE
+!           ok_pred(i)=0.
+!           ok_corr(i)=1.
+!         ENDIF
+!       ELSE
+!         ok_pred(i)=0.
+!         ok_corr(i)=1.
+!       ENDIF
+!
+        ok_pred(i)=0.
+        ok_corr(i)=1.
+!
+        t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
+        q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
+!IM +++
+!       IF(n.EQ.niter) THEN
+!       IF(t_2m(i).LT.t1(i).AND.t_2m(i).LT.ts1(i)) THEN
+!         PRINT*,' BAD t2m LT ',i,nsrf,t_2m(i),t1(i),ts1(i) 
+!       ELSEIF(t_2m(i).GT.t1(i).AND.t_2m(i).GT.ts1(i)) THEN
+!         PRINT*,' BAD t2m GT ',i,nsrf,t_2m(i),t1(i),ts1(i) 
+!       ENDIF
+!       ENDIF
+!IM ---
+      ENDDO
+!
+!
+!----------First aproximation of variables at zref --------------------------
+!
+      zref = 10.0
+      CALL screenp(klon, knon, nsrf, speed, tpot, q1, &
+ &                 ts1, qsurf, z0m, lmon, &
+ &                 ustar, testar, qstar, zref, &
+ &                 delu, delte, delq)
+!
+      DO i = 1, knon
+        u_zref(i) = delu(i)
+        q_zref(i) = max(qsurf(i),0.0) + delq(i)
+        te_zref(i) = ts1(i) + delte(i)
+        temp(i) = te_zref(i) * (psol(i)/pat1(i))**(-RKAPPA)
+!       temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
+!                 (1 + RVTMP2 * max(q_zref(i),0.0))
+        u_zref_p(i) = u_zref(i)
+      ENDDO
+!
+! Iteration of the variables at the reference level zref : corrector ; see Hess & McAvaney, 1995 
+!
+      DO n = 1, niter
+!
+        okri=.TRUE.
+        CALL screenc(klon, knon, nsrf, zxli, &
+ &                   u_zref, temp, q_zref, zref, &
+ &                   ts1, qsurf, z0m, z0h, psol, &
+ &                   ustar, testar, qstar, okri, ri1, &
+ &                   pref, delu, delte, delq)
+!
+        DO i = 1, knon
+          u_zref(i) = delu(i)
+          q_zref(i) = delq(i) + max(qsurf(i),0.0)
+          te_zref(i) = delte(i) + ts1(i)
+          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
+!         temp(i) = te_zref(i) - (zref* RG)/RCPD/ &
+!                   (1 + RVTMP2 * max(q_zref(i),0.0))
+        ENDDO 
+!
+      ENDDO
+!
+      DO i = 1, knon
+        u_zref_c(i) = u_zref(i)
+!
+        u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
+!
+!AM
+        q_zref_c(i) = q_zref(i)
+        temp_c(i) = temp(i)
+        t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
+        q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
+!MA
+      ENDDO
+! 
+      RETURN
+      END subroutine stdlevvar
+!
+      SUBROUTINE stdlevvarn(klon, knon, nsrf, zxli, &
+                           u1, v1, t1, q1, z1, &
+                           ts1, qsurf, z0m, z0h, psol, pat1, &
+                           t_2m, q_2m, t_10m, q_10m, u_10m, n2mout)
+!
+      USE ioipsl_getin_p_mod, ONLY : getin_p
+      IMPLICIT NONE
+!-------------------------------------------------------------------------
+!
+! Objet : calcul de la temperature et l'humidite relative a 2m et du 
+!         module du vent a 10m a partir des relations de Dyer-Businger et
+!         des equations de Louis.
+!
+! Reference : Hess, Colman et McAvaney (1995)        
+!
+! I. Musat, 01.07.2002
+!
+!AM On rajoute en sortie t et q a 10m pr le calcule d'hbtm2 dans clmain
+!
+!-------------------------------------------------------------------------
+!
+! klon----input-I- dimension de la grille physique (= nb_pts_latitude X nb_pts_longitude)
+! knon----input-I- nombre de points pour un type de surface
+! nsrf----input-I- indice pour le type de surface; voir indice_sol_mod.F90
+! zxli----input-L- TRUE si calcul des cdrags selon Laurent Li
+! u1------input-R- vent zonal au 1er niveau du modele
+! v1------input-R- vent meridien au 1er niveau du modele
+! t1------input-R- temperature de l'air au 1er niveau du modele
+! q1------input-R- humidite relative au 1er niveau du modele
+! z1------input-R- geopotentiel au 1er niveau du modele
+! ts1-----input-R- temperature de l'air a la surface
+! qsurf---input-R- humidite relative a la surface
+! z0m, z0h---input-R- rugosite
+! psol----input-R- pression au sol
+! pat1----input-R- pression au 1er niveau du modele
+!
+! t_2m---output-R- temperature de l'air a 2m
+! q_2m---output-R- humidite relative a 2m
+! u_2m--output-R- vitesse du vent a 2m
+! u_10m--output-R- vitesse du vent a 10m
+!AM
+! t_10m--output-R- temperature de l'air a 10m
+! q_10m--output-R- humidite specifique a 10m
+!
+      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) :: psol, pat1
+!
+      REAL, dimension(klon), intent(out) :: t_2m, q_2m
+      REAL, dimension(klon), intent(out) :: u_10m, t_10m, q_10m
+      INTEGER, dimension(klon, 6), intent(out) :: n2mout
+!
+      REAL, dimension(klon) :: u_2m
+      REAL, dimension(klon) :: cdrm2m, cdrh2m, ri2m
+      REAL, dimension(klon) :: cdram, cdrah, zri1
+      REAL, dimension(klon) :: cdmn1, cdhn1, fm1, fh1
+      REAL, dimension(klon) :: cdmn2m, cdhn2m, fm2m, fh2m
+      REAL, dimension(klon) :: ri2m_new
+!-------------------------------------------------------------------------
+      include "flux_arp.h"
+      include "YOMCST.h"
+!IM PLUS
+      include "YOETHF.h"
+!
+! Quelques constantes et options:
+!
+! RKAR : constante de von Karman
+      REAL, PARAMETER :: RKAR=0.40
+! niter : nombre iterations calcul "corrector"
+!     INTEGER, parameter :: niter=6, ncon=niter-1
+!IM 071020     INTEGER, parameter :: niter=2, ncon=niter-1
+      INTEGER, parameter :: niter=2, ncon=niter
+!     INTEGER, parameter :: niter=6, ncon=niter
+!
+! Variables locales
+      INTEGER :: i, n
+      REAL :: zref
+      REAL, dimension(klon) :: speed
+! tpot : temperature potentielle
+      REAL, dimension(klon) :: tpot
+      REAL, dimension(klon) :: cdran
+! ri1 : nb. de Richardson entre la surface --> la 1ere couche
+      REAL, dimension(klon) :: ri1 
+      DOUBLE PRECISION, parameter :: eps=1.0D-20
+      REAL, dimension(klon) :: delu, delte, delq
+      REAL, dimension(klon) :: delh, delm
+      REAL, dimension(klon) :: delh_new, delm_new
+      REAL, dimension(klon) :: u_zref, te_zref, q_zref  
+      REAL, dimension(klon) :: u_zref_pnew, te_zref_pnew, q_zref_pnew
+      REAL, dimension(klon) :: temp, pref
+      REAL, dimension(klon) :: temp_new, pref_new
+      LOGICAL :: okri
+      REAL, dimension(klon) :: u_zref_p, te_zref_p, temp_p, q_zref_p
+      REAL, dimension(klon) :: u_zref_p_new, te_zref_p_new, temp_p_new, q_zref_p_new
+!convergence
+      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) :: cdrm10m, cdrh10m, ri10m
+      REAL, dimension(klon) :: cdmn10m, cdhn10m, fm10m, fh10m
+      REAL, dimension(klon) :: cdn2m, cdn1
+      REAL :: CEPDUE,zdu2 
+      INTEGER :: nzref, nz1
+      LOGICAL, dimension(klon) :: ok_t2m_toosmall, ok_t2m_toobig
+      LOGICAL, dimension(klon) :: ok_q2m_toosmall, ok_q2m_toobig
+      LOGICAL, dimension(klon) :: ok_u2m_toobig
+      LOGICAL, dimension(klon) :: ok_t10m_toosmall, ok_t10m_toobig
+      LOGICAL, dimension(klon) :: ok_q10m_toosmall, ok_q10m_toobig
+      LOGICAL, dimension(klon) :: ok_u10m_toobig
+      INTEGER, dimension(klon, 6) :: n10mout
+
+!------------------------------------------------------------------------- 
+      CEPDUE=0.1
+!
+! n2mout : compteur des pas de temps ou t2m,q2m ou u2m sont en dehors des intervalles
+!          [tsurf, temp], [qsurf, q1] ou [0, speed]
+!
+      n2mout(:,:)=0
+      
+      DO i=1, knon
+       speed(i)=MAX(SQRT(u1(i)**2+v1(i)**2),CEPDUE)
+       ri1(i) = 0.0
+      ENDDO
+!
+      okri=.FALSE.
+      CALL cdrag(knon, nsrf, &
+ &                   speed, t1, q1, z1, &
+ &                   psol, ts1, qsurf, z0m, z0h, &
+ &                   cdram, cdrah, zri1, pref)
+
+!
+      DO i = 1, knon
+        ri1(i) = zri1(i)
+        tpot(i) = t1(i)* (psol(i)/pat1(i))**RKAPPA
+        zdu2 = MAX(CEPDUE*CEPDUE, speed(i)**2)
+!
+      ENDDO
+!
+!----------First aproximation of variables at zref --------------------------
+      zref = 2.0
+!
+! calcul first-guess en utilisant dans les calculs à 2m
+! le Richardson de la premiere couche atmospherique 
+!
+       CALL screencn(klon, knon, nsrf, zxli, &
+ &                   speed, tpot, q1, zref, &
+ &                   ts1, qsurf, z0m, z0h, psol, &           
+ &                   cdram, cdrah,  okri, &
+ &                   ri1, 1, &
+ &                   pref_new, delm_new, delh_new, ri2m)
+!
+       DO i = 1, knon
+         u_zref(i) = delm_new(i)*speed(i)
+         u_zref_p(i) = u_zref(i)
+         q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
+         &           max(qsurf(i),0.0)*(1-delh_new(i))
+         q_zref_p(i) = q_zref(i)
+         te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
+         te_zref_p(i) = te_zref(i)
+!
+! return to normal temperature
+         temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
+         temp_p(i) = temp(i)
+!
+! compteurs ici
+!
+         ok_t2m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
+         & te_zref(i).LT.ts1(i)
+         ok_t2m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
+         & te_zref(i).GT.ts1(i)
+         ok_q2m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
+         & q_zref(i).LT.qsurf(i)
+         ok_q2m_toobig(i)=q_zref(i).GT.q1(i).AND. &
+         & q_zref(i).GT.qsurf(i)
+         ok_u2m_toobig(i)=u_zref(i).GT.speed(i)
+!
+         IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
+             n2mout(i,1)=n2mout(i,1)+1
+         ENDIF
+         IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
+             n2mout(i,3)=n2mout(i,3)+1
+         ENDIF
+         IF(ok_u2m_toobig(i)) THEN
+             n2mout(i,5)=n2mout(i,5)+1
+         ENDIF
+!
+         IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
+          & ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
+          & ok_u2m_toobig(i)) THEN 
+             delm_new(i)=min(max(delm_new(i),0.),1.)
+             delh_new(i)=min(max(delh_new(i),0.),1.)
+             u_zref(i) = delm_new(i)*speed(i)
+             u_zref_p(i) = u_zref(i)
+             q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
+         &               max(qsurf(i),0.0)*(1-delh_new(i))
+             q_zref_p(i) = q_zref(i)
+             te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
+             te_zref_p(i) = te_zref(i)
+!
+! return to normal temperature
+             temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
+             temp_p(i) = temp(i)
+         ENDIF 
+!
+       ENDDO
+!
+! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995 
+!
+      DO n = 1, niter
+!
+        okri=.TRUE.
+        CALL screencn(klon, knon, nsrf, zxli, &
+ &                   u_zref, temp, q_zref, zref, &
+ &                   ts1, qsurf, z0m, z0h, psol, &
+ &                   cdram, cdrah,  okri, &
+ &                   ri1, 0, &
+ &                   pref, delm, delh, ri2m)
+!
+        DO i = 1, knon
+          u_zref(i) = delm(i)*speed(i)
+          q_zref(i) = delh(i)*max(q1(i),0.0) + &
+          &           max(qsurf(i),0.0)*(1-delh(i))
+          te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
+!
+! return to normal temperature
+          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
+!
+! compteurs ici
+!
+          ok_t2m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
+          & te_zref(i).LT.ts1(i)
+          ok_t2m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
+          & te_zref(i).GT.ts1(i)
+          ok_q2m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
+          & q_zref(i).LT.qsurf(i)
+          ok_q2m_toobig(i)=q_zref(i).GT.q1(i).AND. &
+          & q_zref(i).GT.qsurf(i)
+          ok_u2m_toobig(i)=u_zref(i).GT.speed(i)
+!
+          IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i)) THEN
+              n2mout(i,2)=n2mout(i,2)+1
+          ENDIF
+          IF(ok_q2m_toosmall(i).OR.ok_q2m_toobig(i)) THEN
+              n2mout(i,4)=n2mout(i,4)+1
+          ENDIF
+          IF(ok_u2m_toobig(i)) THEN
+              n2mout(i,6)=n2mout(i,6)+1
+          ENDIF
+!
+          IF(ok_t2m_toosmall(i).OR.ok_t2m_toobig(i).OR. &
+           & ok_q2m_toosmall(i).OR.ok_q2m_toobig(i).OR. &
+           & ok_u2m_toobig(i)) THEN 
+              delm(i)=min(max(delm(i),0.),1.)
+              delh(i)=min(max(delh(i),0.),1.)
+              u_zref(i) = delm(i)*speed(i)
+              q_zref(i) = delh(i)*max(q1(i),0.0) + &
+          &           max(qsurf(i),0.0)*(1-delh(i))
+              te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
+              temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
+          ENDIF
+!
+!
+          IF(n.EQ.ncon) THEN
+            te_zref_con(i) = te_zref(i)
+            q_zref_con(i) = q_zref(i)
+          ENDIF 
+!
+        ENDDO 
+!
+      ENDDO 
+!
+      DO i = 1, knon
+        q_zref_c(i) = q_zref(i)
+        temp_c(i) = temp(i)
+!
+        ok_pred(i)=0.
+        ok_corr(i)=1.
+!
+        t_2m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
+        q_2m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
+!
+        u_zref_c(i) = u_zref(i)
+        u_2m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
+      ENDDO
+!
+!
+!----------First aproximation of variables at zref --------------------------
+!
+      zref = 10.0
+!
+       CALL screencn(klon, knon, nsrf, zxli, &
+ &                   speed, tpot, q1, zref, &
+ &                   ts1, qsurf, z0m, z0h, psol, &           
+ &                   cdram, cdrah,  okri, &
+ &                   ri1, 1, &
+ &                   pref_new, delm_new, delh_new, ri10m)
+!
+       DO i = 1, knon
+         u_zref(i) = delm_new(i)*speed(i)
+         q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
+         &           max(qsurf(i),0.0)*(1-delh_new(i))
+         te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
+         temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
+         u_zref_p(i) = u_zref(i)
+!
+! compteurs ici
+!
+         ok_t10m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
+         & te_zref(i).LT.ts1(i)
+         ok_t10m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
+         & te_zref(i).GT.ts1(i)
+         ok_q10m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
+         & q_zref(i).LT.qsurf(i)
+         ok_q10m_toobig(i)=q_zref(i).GT.q1(i).AND. &
+         & q_zref(i).GT.qsurf(i)
+         ok_u10m_toobig(i)=u_zref(i).GT.speed(i)
+!
+         IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
+          & ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
+          & ok_u10m_toobig(i)) THEN 
+             delm_new(i)=min(max(delm_new(i),0.),1.)
+             delh_new(i)=min(max(delh_new(i),0.),1.)
+             u_zref(i) = delm_new(i)*speed(i)
+             u_zref_p(i) = u_zref(i)
+             q_zref(i) = delh_new(i)*max(q1(i),0.0) + &
+         &               max(qsurf(i),0.0)*(1-delh_new(i))
+             te_zref(i) = delh_new(i)*tpot(i) + ts1(i)*(1-delh_new(i))
+             temp(i) = te_zref(i) * (psol(i)/pref_new(i))**(-RKAPPA)
+         ENDIF 
+!
+       ENDDO
+!
+! Iteration of the variables at the reference level zref : corrector calculation ; see Hess & McAvaney, 1995 
+!
+      DO n = 1, niter
+!
+        okri=.TRUE.
+        CALL screencn(klon, knon, nsrf, zxli, &
+ &                   u_zref, temp, q_zref, zref, &
+ &                   ts1, qsurf, z0m, z0h, psol, &
+ &                   cdram, cdrah,  okri, &
+ &                   ri1, 0, &
+ &                   pref, delm, delh, ri10m)
+!
+        DO i = 1, knon
+          u_zref(i) = delm(i)*speed(i)
+          q_zref(i) = delh(i)*max(q1(i),0.0) + &
+          &           max(qsurf(i),0.0)*(1-delh(i))
+          te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
+!
+! return to normal temperature
+          temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
+!
+! compteurs ici
+!
+          ok_t10m_toosmall(i)=te_zref(i).LT.tpot(i).AND. &
+          & te_zref(i).LT.ts1(i)
+          ok_t10m_toobig(i)=te_zref(i).GT.tpot(i).AND. &
+          & te_zref(i).GT.ts1(i)
+          ok_q10m_toosmall(i)=q_zref(i).LT.q1(i).AND. &
+          & q_zref(i).LT.qsurf(i)
+          ok_q10m_toobig(i)=q_zref(i).GT.q1(i).AND. &
+          & q_zref(i).GT.qsurf(i)
+          ok_u10m_toobig(i)=u_zref(i).GT.speed(i)
+!
+          IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i)) THEN
+              n10mout(i,2)=n10mout(i,2)+1
+          ENDIF
+          IF(ok_q10m_toosmall(i).OR.ok_q10m_toobig(i)) THEN
+              n10mout(i,4)=n10mout(i,4)+1
+          ENDIF
+          IF(ok_u10m_toobig(i)) THEN
+              n10mout(i,6)=n10mout(i,6)+1
+          ENDIF
+!
+          IF(ok_t10m_toosmall(i).OR.ok_t10m_toobig(i).OR. &
+           & ok_q10m_toosmall(i).OR.ok_q10m_toobig(i).OR. &
+           & ok_u10m_toobig(i)) THEN 
+              delm(i)=min(max(delm(i),0.),1.)
+              delh(i)=min(max(delh(i),0.),1.)
+              u_zref(i) = delm(i)*speed(i)
+              q_zref(i) = delh(i)*max(q1(i),0.0) + &
+          &           max(qsurf(i),0.0)*(1-delh(i))
+              te_zref(i) = delh(i)*tpot(i) + ts1(i)*(1-delh(i))
+              temp(i) = te_zref(i) * (psol(i)/pref(i))**(-RKAPPA)
+          ENDIF
+!
+!
+          IF(n.EQ.ncon) THEN
+            te_zref_con(i) = te_zref(i)
+            q_zref_con(i) = q_zref(i)
+          ENDIF 
+!
+        ENDDO 
+!
+      ENDDO 
+!
+      DO i = 1, knon
+        q_zref_c(i) = q_zref(i)
+        temp_c(i) = temp(i)
+!
+        ok_pred(i)=0.
+        ok_corr(i)=1.
+!
+        t_10m(i) = temp_p(i) * ok_pred(i) + temp_c(i) * ok_corr(i)
+        q_10m(i) = q_zref_p(i) * ok_pred(i) + q_zref_c(i) * ok_corr(i)
+!
+        u_zref_c(i) = u_zref(i)
+        u_10m(i) = u_zref_p(i) * ok_pred(i) + u_zref_c(i) * ok_corr(i)
+      ENDDO
+!
+      RETURN
+      END subroutine stdlevvarn
+
+END MODULE stdlevvar_mod
