Index: LMDZ5/branches/testing/libf/phylmd/cdrag.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cdrag.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/cdrag.F90	(revision 2886)
@@ -2,22 +2,26 @@
 !$Id$
 !
- SUBROUTINE cdrag( knon,  nsrf,   &
+ SUBROUTINE cdrag(knon,  nsrf,   &
      speed, t1,    q1,    zgeop1, &
      psol,  tsurf, qsurf, z0m, z0h,  &
-     pcfm,  pcfh,  zri,   pref )
+     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).   
-!
-! Modified histroy: 
-!   27-Jan-2014: richardson number inconsistant between
-!   coefcdrag.F90 and clcdrag.F90, Fuxing WANG wrote this subroutine
-!   by merging coefcdrag and clcdrag.
+!         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:
@@ -28,4 +32,7 @@
 !     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.
@@ -34,103 +41,121 @@
 !     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
-!
-! knon----input-I- nombre de points pour un type de surface
-! nsrf----input-I- indice pour le type de surface; voir indicesol.h
-! speed---input-R- module du vent au 1er niveau du modele
-! t1------input-R- temperature de l'air au 1er niveau du modele
-! q1------input-R- humidite de l'air au 1er niveau du modele
-! zgeop---input-R- geopotentiel au 1er niveau du modele
-! psol----input-R- pression au sol
-! tsurf---input-R- temperature de l'air a la surface 
-! qsurf---input-R- humidite de l'air a la surface
-! z0m, z0h---input-R- rugosite
-!! u1, v1 are removed, speed is used. Fuxing WANG, 04/03/2015,
-!! u1------input-R- vent zonal au 1er niveau du modele
-!! v1------input-R- vent meridien au 1er niveau du modele
-!
-! pcfm---output-R- cdrag pour le moment 
-! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible
-! zri----output-R- Richardson number
-! pref---output-R- pression au niveau zgeop/RG
-!
-! Parameters: 
-! ckap-----Karman constant
-! cb,cc,cd-parameters in Louis et al., 1982
 ! ================================================================= 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
+  
+  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)        :: psol  ! pression au sol
-  REAL, DIMENSION(klon), INTENT(IN)        :: t1    ! temperature at 1st level
-  REAL, DIMENSION(klon), INTENT(IN)        :: q1    ! humidity at 1st level
   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) 
-!  paprs, pplay u1, v1: to be deleted
-!  they were in the old clcdrag. Fuxing WANG, 04/03/2015
-!  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
-!  REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
-!  REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1
+  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)       :: pcfm  ! Drag coefficient for heat flux
-  REAL, DIMENSION(klon), INTENT(OUT)       :: pcfh  ! Drag coefficient for momentum
+  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
 
-! Parametres local
-  INTEGER                                  :: ng_q1    ! Number of grids that q1 < 0.0
-  INTEGER                                  :: ng_qsurf ! Number of grids that qsurf < 0.0
-!  zgeop1, psol: to be deleted, they are inputs now. Fuxing WANG, 04/03/2015
-!  REAL, DIMENSION(klon)                    :: zgeop1! geopotentiel au 1er niveau du modele
-!  REAL, DIMENSION(klon)                    :: psol  ! pression au sol
-!
-! ================================================================= c
-!
+! Variables Locales
+!******************************************************************
+ 
+
   INCLUDE "YOMCST.h"
   INCLUDE "YOETHF.h"
-!  INCLUDE "indicesol.h"
   INCLUDE "clesphys.h"
-!
-! Quelques constantes et options:
-!!$PB      REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2
-  REAL, PARAMETER :: CKAP=0.40, CB=5.0, CC=5.0, CD=5.0, CEPDU2 = (0.1)**2
-!
-! Variables locales :
-  INTEGER               :: i
-  REAL                  :: zdu2, ztsolv
-  REAL                  :: ztvd, zscf
-  REAL                  :: zucf, zcr
-  REAL                  :: friv, frih
-  REAL, DIMENSION(klon) :: zcfm1, zcfm2 ! Drag coefficient for momentum
-  REAL, DIMENSION(klon) :: zcfh1, zcfh2 ! Drag coefficient for heat flux
-  LOGICAL, PARAMETER    :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li
-  REAL, DIMENSION(klon) :: zcdn_m, zcdn_h         ! Drag coefficient in neutral conditions
+
+
+  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
-!
-! Fonctions thermodynamiques et fonctions d'instabilite
-  REAL                  :: fsta, fins, x
-  fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
-  fins(x) = SQRT(1.0-18.0*x)
+
+  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
-! Fuxing WANG, 04/03/2015, delete the calculation of zgeop1 
-! (le geopotentiel du premier couche de modele).
-! zgeop1 is an input ivariable in this subroutine.
-!  DO i = 1, knon
-!     zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) &
-!          * (paprs(i,1)-pplay(i,1))
-!  END DO
-! ================================================================= 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
@@ -154,80 +179,233 @@
   ENDIF
 
-! Calculer le frottement au sol (Cdrag)
-  DO i = 1, knon
-!------------------------------------------------------------
-! u1, v1 are replaced by speed. Fuxing WANG, 04/03/2015,
-!     zdu2 = MAX(CEPDU2,u1(i)**2+v1(i)**2)
-!------------------------------------------------------------
+
+
+!=============================================================================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)
-!     psol(i) = paprs(i,1)
      pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* &
-                 (1.+ RETV * max(q1(i),0.0))))  ! negative q1 set to zero
-!------------ the old calculations in clcdrag----------------
-!     ztsolv = tsurf(i) * (1.0+RETV*qsurf(i)) 
-!     ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) &
-!          *(1.+RETV*q1(i)) 
-!------------------------------------------------------------
-! Fuxing WANG, 04/03/2015, in this revised version, 
-! the negative qsurf and q1 are set to zero (as in coefcdrag)
-     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
+                 (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*q1(i))) &
+          *(1.+RETV*max(q1(i),0.0))                      ! negative q1 set to zero
      zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
 
 
-! Coefficients CD neutres pour m et h : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h))
+! 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)))
-     zcdn_m(i) = zzzcd*zzzcd
-     zcdn_h(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i))))
-
-     IF (zri(i) .GT. 0.) THEN      ! situation stable
+     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))
-        IF (.NOT.zxli) THEN
+
+       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)))
-           friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
-           zcfm1(i) = zcdn_m(i) * friv
-           frih = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min )
-!!$ PB     zcfh1(i) = zcdn(i) * frih
-!!$ PB     zcfh1(i) = f_cdrag_stable * zcdn(i) * frih
-           zcfh1(i) = f_cdrag_ter * zcdn_h(i) * frih
-           IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn_h(i) * frih
-!!$ PB
-           pcfm(i) = zcfm1(i)
-           pcfh(i) = zcfh1(i)
-        ELSE
-           pcfm(i) = zcdn_m(i)* fsta(zri(i))
-           pcfh(i) = zcdn_h(i)* fsta(zri(i))
-        ENDIF
-     ELSE                          ! situation instable
-        IF (.NOT.zxli) THEN
-           zucf = 1./(1.+3.0*CB*CC*zcdn_m(i)*SQRT(ABS(zri(i)) &
-                *(1.0+zgeop1(i)/(RG*z0m(i)))))
-           zcfm2(i) = zcdn_m(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
-!!$ PB     zcfh2(i) = zcdn_h(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
-           zcfh2(i) = f_cdrag_ter*zcdn_h(i)*amax1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min)
-           pcfm(i) = zcfm2(i)
-           pcfh(i) = zcfh2(i)
-        ELSE
-           pcfm(i) = zcdn_m(i)* fins(zri(i))
-           pcfh(i) = zcdn_h(i)* fins(zri(i))
-        ENDIF
-        IF(iflag_gusts==0) THEN
-! cdrah sur l'ocean cf. Miller et al. (1992) - only active when gustiness parameterization is not active
-           zcr = (0.0016/(zcdn_m(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
-           IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn_h(i)*(1.0+zcr**1.25)**(1./1.25)
-        ENDIF
-     ENDIF
-  END DO
-
+           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
      
-  ! IM cf JLD : on seuille cdrag_m et cdrag_h
-  IF (nsrf == is_oce) THEN
-     DO i=1,knon
-        pcfm(i)=MIN(pcfm(i),cdmmax)
-        pcfh(i)=MIN(pcfh(i),cdhmax)
-     END DO
-  END IF
+ 
 
 END SUBROUTINE cdrag
+
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Index: LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 2886)
@@ -182,10 +182,10 @@
       CALL histwrite3d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr(:,icl,:),nvert,icl)
    enddo
-     CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
-#endif
-
+#endif
+
+   CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
    do k=1,PARASOL_NREFL
     do ip=1, Npoints
-     if (stlidar%cldlayer(ip,4).gt.0.01) then
+     if (stlidar%cldlayer(ip,4).gt.0.01.and.stlidar%parasolrefl(ip,k).ne.missing_val) then
        parasolcrefl(ip,k)=(stlidar%parasolrefl(ip,k)-0.03*(1.-stlidar%cldlayer(ip,4)))/ &
                             stlidar%cldlayer(ip,4)
@@ -196,4 +196,5 @@
      endif
     enddo
+   
    enddo
    CALL histwrite3d_cosp(o_Ncrefl,Ncref,nvertp)
Index: LMDZ5/branches/testing/libf/phylmd/cpl_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cpl_mod.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/cpl_mod.F90	(revision 2886)
@@ -52,4 +52,6 @@
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp
   !$OMP THREADPRIVATE(cpl_windsp)
+  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_sens_rain, cpl_sens_snow
+  !$OMP THREADPRIVATE(cpl_sens_rain, cpl_sens_snow)
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_taumod
   !$OMP THREADPRIVATE(cpl_taumod)
@@ -90,4 +92,6 @@
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_windsp2D
   !$OMP THREADPRIVATE(cpl_windsp2D)
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: cpl_sens_rain2D, cpl_sens_snow2D
+  !$OMP THREADPRIVATE(cpl_sens_rain2D, cpl_sens_snow2D)
   REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: cpl_atm_co22D
   !$OMP THREADPRIVATE(cpl_atm_co22D)
@@ -169,4 +173,8 @@
     sum_error = sum_error + error
     ALLOCATE(cpl_taumod(klon,2), stat = error)
+    sum_error = sum_error + error
+    ALLOCATE(cpl_sens_rain(klon,2), stat = error)
+    sum_error = sum_error + error
+    ALLOCATE(cpl_sens_snow(klon,2), stat = error)
     sum_error = sum_error + error
     ALLOCATE(cpl_rriv2D(nbp_lon,jj_nb), stat=error)
@@ -531,5 +539,6 @@
   SUBROUTINE cpl_send_ocean_fields(itime, knon, knindex, &
        swdown, lwdown, fluxlat, fluxsens, &
-       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy, windsp)
+       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy, windsp,&
+       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
 !
 ! This subroutine cumulates some fields for each time-step during a coupling 
@@ -552,4 +561,6 @@
     REAL, DIMENSION(klon), INTENT(IN)       :: evap, tsurf, fder, albsol
     REAL, DIMENSION(klon), INTENT(IN)       :: taux, tauy, windsp
+    REAL, DIMENSION(klon), INTENT(IN)       :: sens_prec_liq, sens_prec_sol
+    REAL, DIMENSION(klon), INTENT(IN)       :: lat_prec_liq, lat_prec_sol
 
 ! Local variables
@@ -583,4 +594,6 @@
        cpl_tauy(1:knon,cpl_index) = 0.0
        cpl_windsp(1:knon,cpl_index) = 0.0
+       cpl_sens_rain(1:knon,cpl_index) = 0.0
+       cpl_sens_snow(1:knon,cpl_index) = 0.0
        cpl_taumod(1:knon,cpl_index) = 0.0
        IF (carbon_cycle_cpl) cpl_atm_co2(1:knon,cpl_index) = 0.0
@@ -614,4 +627,8 @@
        cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) + &
             windsp(ig)      / REAL(nexca)
+       cpl_sens_rain(ig,cpl_index) = cpl_sens_rain(ig,cpl_index) + &
+            sens_prec_liq(ig)      / REAL(nexca)
+       cpl_sens_snow(ig,cpl_index) = cpl_sens_snow(ig,cpl_index) + &
+            sens_prec_sol(ig)      / REAL(nexca)
        cpl_taumod(ig,cpl_index) =   cpl_taumod(ig,cpl_index) + &
           SQRT ( taux(ig)*taux(ig)+tauy(ig)*tauy(ig) ) / REAL (nexca)
@@ -654,4 +671,8 @@
           sum_error = sum_error + error
           ALLOCATE(cpl_windsp2D(nbp_lon,jj_nb), stat=error)
+          sum_error = sum_error + error
+          ALLOCATE(cpl_sens_rain2D(nbp_lon,jj_nb,2), stat=error)
+          sum_error = sum_error + error
+          ALLOCATE(cpl_sens_snow2D(nbp_lon,jj_nb,2), stat=error)
           sum_error = sum_error + error
           ALLOCATE(cpl_taumod2D(nbp_lon,jj_nb,2), stat=error)
@@ -706,4 +727,10 @@
             knon, knindex)
 
+       CALL gath2cpl(cpl_sens_rain(:,cpl_index), cpl_sens_rain2D(:,:,cpl_index), &
+            knon, knindex)
+
+       CALL gath2cpl(cpl_sens_snow(:,cpl_index), cpl_sens_snow2D(:,:,cpl_index), &
+            knon, knindex)
+
        CALL gath2cpl(cpl_taumod(:,cpl_index), cpl_taumod2D(:,:,cpl_index), &
             knon, knindex)
@@ -722,5 +749,6 @@
        pctsrf, lafin, rlon, rlat, &
        swdown, lwdown, fluxlat, fluxsens, &
-       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy)
+       precip_rain, precip_snow, evap, tsurf, fder, albsol, taux, tauy,&
+       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
 !
 ! This subroutine cumulates some fields for each time-step during a coupling 
@@ -746,4 +774,6 @@
     REAL, DIMENSION(klon), INTENT(IN)       :: albsol, taux, tauy
     REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
+    REAL, DIMENSION(klon), INTENT(IN)       :: sens_prec_liq, sens_prec_sol
+    REAL, DIMENSION(klon), INTENT(IN)       :: lat_prec_liq, lat_prec_sol
     LOGICAL, INTENT(IN)                     :: lafin
 
@@ -778,4 +808,6 @@
        cpl_taux(1:knon,cpl_index) = 0.0
        cpl_tauy(1:knon,cpl_index) = 0.0
+       cpl_sens_rain(1:knon,cpl_index) = 0.0
+       cpl_sens_snow(1:knon,cpl_index) = 0.0
        cpl_taumod(1:knon,cpl_index) = 0.0
     ENDIF
@@ -806,4 +838,8 @@
        cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) + &
             tauy(ig)        / REAL(nexca)     
+       cpl_sens_rain(ig,cpl_index) = cpl_sens_rain(ig,cpl_index) + &
+            sens_prec_liq(ig)      / REAL(nexca)
+       cpl_sens_snow(ig,cpl_index) = cpl_sens_snow(ig,cpl_index) + &
+            sens_prec_sol(ig)      / REAL(nexca)
        cpl_taumod(ig,cpl_index) = cpl_taumod(ig,cpl_index) + &
             SQRT ( taux(ig)*taux(ig)+tauy(ig)*tauy(ig) ) / REAL(nexca) 
@@ -839,4 +875,8 @@
           sum_error = sum_error + error
           ALLOCATE(cpl_windsp2D(nbp_lon,jj_nb), stat=error)
+          sum_error = sum_error + error
+          ALLOCATE(cpl_sens_rain2D(nbp_lon,jj_nb,2), stat=error)
+          sum_error = sum_error + error
+          ALLOCATE(cpl_sens_snow2D(nbp_lon,jj_nb,2), stat=error)
           sum_error = sum_error + error
           ALLOCATE(cpl_taumod2D(nbp_lon,jj_nb,2), stat=error)
@@ -889,4 +929,10 @@
 
        CALL gath2cpl(cpl_tauy(:,cpl_index), cpl_tauy2D(:,:,cpl_index), &
+            knon, knindex)
+
+       CALL gath2cpl(cpl_sens_rain(:,cpl_index), cpl_sens_rain2D(:,:,cpl_index), &
+            knon, knindex)
+
+       CALL gath2cpl(cpl_sens_snow(:,cpl_index), cpl_sens_snow2D(:,:,cpl_index), &
             knon, knindex)
 
@@ -1078,4 +1124,8 @@
     tab_flds(:,:,ids_nsfice) = cpl_nsol2D(:,:,2)
     tab_flds(:,:,ids_dflxdt) = cpl_fder2D(:,:,2)
+    tab_flds(:,:,ids_qraioc) = cpl_sens_rain2D(:,:,1)
+    tab_flds(:,:,ids_qsnooc) = cpl_sens_snow2D(:,:,1)
+    tab_flds(:,:,ids_qraiic) = cpl_sens_rain2D(:,:,2)
+    tab_flds(:,:,ids_qsnoic) = cpl_sens_snow2D(:,:,2)
     
     IF (version_ocean=='nemo') THEN
@@ -1279,4 +1329,7 @@
     DEALLOCATE(cpl_taux2D, cpl_tauy2D, cpl_windsp2D, cpl_taumod2D, stat=error )
     sum_error = sum_error + error
+    DEALLOCATE(cpl_sens_rain2D, cpl_sens_snow2D, stat=error)
+    sum_error = sum_error + error
+
     
     IF (carbon_cycle_cpl) THEN
Index: LMDZ5/branches/testing/libf/phylmd/ener_conserv.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/ener_conserv.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/ener_conserv.F90	(revision 2886)
@@ -24,4 +24,5 @@
 USE phys_local_var_mod, ONLY : d_t_eva,d_t_lsc,d_q_eva,d_q_lsc
 USE phys_output_var_mod, ONLY : bils_ec,bils_ech,bils_tke,bils_kinetic,bils_enthalp,bils_latent,bils_diss
+USE add_phys_tend_mod, ONLY : fl_cor_ebil 
 
 IMPLICIT none
@@ -59,9 +60,13 @@
    DO k = 1, klev
    DO i = 1, klon
-      ZRCPD = RCPD*(1.0+RVTMP2*(pqn(i,k)+pqln(i,k)+pqsn(i,k)))
-      d_t_ec(i,k)=0.5/ZRCPD &
- &      *(puo(i,k)**2+pvo(i,k)**2-pun(i,k)**2-pvn(i,k)**2)
-      ENDDO
-      ENDDO
+     IF (fl_cor_ebil .GT. 0) then
+       ZRCPD = RCPD*(1.0+RVTMP2*(pqn(i,k)+pqln(i,k)+pqsn(i,k)))
+     ELSE
+       ZRCPD = RCPD*(1.0+RVTMP2*pqn(i,k))
+     ENDIF
+     d_t_ec(i,k)=0.5/ZRCPD &
+ &     *(puo(i,k)**2+pvo(i,k)**2-pun(i,k)**2-pvn(i,k)**2)
+   ENDDO
+   ENDDO
 !-jld ec_conser
 
Index: LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90	(revision 2886)
@@ -368,4 +368,9 @@
                 / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
          end if
+        ENDDO
+     ELSE  ! IF(k.LE.klevm1) 
+        DO i = 1, klon 
+           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG 
+           zmqc(i) = 0.
         ENDDO
      ENDIF ! end IF(k.LE.klevm1)
Index: LMDZ5/branches/testing/libf/phylmd/oasis.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/oasis.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/oasis.F90	(revision 2886)
@@ -55,5 +55,9 @@
   INTEGER, PARAMETER :: ids_atmco2 = 24
   INTEGER, PARAMETER :: ids_taumod = 25
-  INTEGER, PARAMETER :: maxsend    = 25  ! Maximum number of fields to send
+  INTEGER, PARAMETER :: ids_qraioc = 26
+  INTEGER, PARAMETER :: ids_qsnooc = 27
+  INTEGER, PARAMETER :: ids_qraiic = 28
+  INTEGER, PARAMETER :: ids_qsnoic = 29
+  INTEGER, PARAMETER :: maxsend    = 29  ! Maximum number of fields to send
   
   ! Id for fields received from ocean
@@ -177,4 +181,8 @@
             infosend(ids_atmco2)%action = .TRUE. ; infosend(ids_atmco2)%name = 'COATMCO2'
         ENDIF
+        infosend(ids_qraioc)%action = .TRUE. ; infosend(ids_qraioc)%name = 'COQRAIOC'
+        infosend(ids_qsnooc)%action = .TRUE. ; infosend(ids_qsnooc)%name = 'COQSNOOC'
+        infosend(ids_qraiic)%action = .TRUE. ; infosend(ids_qraiic)%name = 'COQRAIIC'
+        infosend(ids_qsnoic)%action = .TRUE. ; infosend(ids_qsnoic)%name = 'COQSNOIC'
         
     ELSE IF (version_ocean=='opa8') THEN
Index: LMDZ5/branches/testing/libf/phylmd/ocean_cpl_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/ocean_cpl_mod.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/ocean_cpl_mod.F90	(revision 2886)
@@ -191,7 +191,6 @@
     CALL cpl_send_ocean_fields(itime, knon, knindex, &
          swnet, lwnet, fluxlat, fluxsens, &
-         precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1, windsp&
-         ) 
-!         ,sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+         precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1, windsp,&
+         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
     
 
@@ -353,7 +352,6 @@
        pctsrf, lafin, rlon, rlat, &
        swnet, lwnet, fluxlat, fluxsens, &
-       precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1 &
-       )
-!      ,sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
+       precip_rain, precip_snow, evap, tsurf_new, fder_new, alb1, flux_u1, flux_v1,&
+       sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol)
 
  
Index: LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2886)
@@ -2326,4 +2326,5 @@
 !!!
        IF (iflag_split .eq.0) THEN
+        wake_dltke(:,:,nsrf) = 0.
         DO k = 1, klev
            DO j = 1, knon
Index: LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90	(revision 2886)
@@ -352,4 +352,6 @@
 !!!
 !!!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      LOGICAL, SAVE, ALLOCATABLE :: ptconv(:,:)
+      !$OMP THREADPRIVATE(ptconv)
 !>jyg+nrlmd
   !
@@ -676,4 +678,6 @@
       ALLOCATE(cdragm_x(klon), cdragm_w(klon))
       ALLOCATE(kh(klon), kh_x(klon), kh_w(klon))
+!
+      ALLOCATE(ptconv(klon,klev))
 !
       ALLOCATE(wbeff(klon), convoccur(klon), zmax_th(klon))
@@ -934,4 +938,6 @@
       DEALLOCATE(kh, kh_x, kh_w)
 !
+      DEALLOCATE(ptconv)
+!
       DEALLOCATE(wbeff, convoccur, zmax_th)
       DEALLOCATE(zq2m, zt2m, weak_inversion)
Index: LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2886)
@@ -867,19 +867,26 @@
       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)) /)
 
-  TYPE(ctrl_out), SAVE, DIMENSION(7) :: o_qSTDlevs     = (/                             &
-      ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q850', "Specific humidity 850hPa", &
-      "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
-      ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q700', "Specific humidity 700hPa", &
-      "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
-      ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q500', "Specific humidity 500hPa", &
-      "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
-      ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q200', "Specific humidity 200hPa", &
-      "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
-      ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q100', "Specific humidity 100hPa", &
-      "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
-      ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q50', "Specific humidity 50hPa",  &
-      "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
-      ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q10', "Specific humidity 10hPa", &
-      "kg/kg", (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)) /)
+  TYPE(ctrl_out), SAVE, DIMENSION(7) :: o_qSTDlevs     = (/ &
+       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q850', &
+       "Specific humidity 850hPa", "kg/kg", &
+       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
+       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q700', &
+       "Specific humidity 700hPa", "kg/kg", &
+       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
+       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q500', &
+       "Specific humidity 500hPa", "kg/kg", &
+       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
+       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q200', &
+       "Specific humidity 200hPa", "kg/kg", &
+       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
+       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q100', &
+       "Specific humidity 100hPa", "kg/kg", &
+       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
+       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q50', &
+       "Specific humidity 50hPa", "kg/kg", &
+       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)), &
+       ctrl_out((/ 1, 7, 7, 10, 10, 10, 11, 11, 11, 11/),'q10', &
+       "Specific humidity 10hPa", "kg/kg", &
+       (/ 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)', 'inst(X)' /)) /)
 
   TYPE(ctrl_out), SAVE, DIMENSION(7) :: o_zSTDlevs   = (/                           &
Index: LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/physiq_mod.F90	(revision 2886)
@@ -156,4 +156,5 @@
 !!!       d_s_the, d_dens_the, &            ! due to thermals
        !                                  
+       ptconv, &
        wbeff, convoccur, zmax_th, &
        sens, flwp, fiwp,  &
@@ -599,4 +600,6 @@
     LOGICAL,SAVE :: ok_adjwk=.FALSE.
     !$OMP THREADPRIVATE(ok_adjwk) 
+    INTEGER,SAVE :: iflag_adjwk=0            !jyg
+    !$OMP THREADPRIVATE(iflag_adjwk)         !jyg
     REAL,SAVE :: oliqmax=999.,oicemax=999.
     !$OMP THREADPRIVATE(oliqmax,oicemax) 
@@ -876,5 +879,5 @@
     save iflag_cld_th
     !$OMP THREADPRIVATE(iflag_cld_th)
-    logical ptconv(klon,klev)
+!IM logical ptconv(klon,klev)  !passe dans phys_local_var_mod
     !IM cf. AM 081204 BEG
     logical ptconvth(klon,klev)
@@ -1201,5 +1204,10 @@
        CALL suphel ! initialiser constantes et parametres phys.
        CALL getin_p('random_notrig_max',random_notrig_max)
-       CALL getin_p('ok_adjwk',ok_adjwk)
+       CALL getin_p('ok_adjwk',ok_adjwk) 
+       IF (ok_adjwk) iflag_adjwk=2  ! for compatibility with older versions
+       ! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
+                      ! 1 => convective adjustment but state variables are unchanged
+                      ! 2 => convective adjustment and state variables are changed
+       CALL getin_p('iflag_adjwk',iflag_adjwk)
        CALL getin_p('oliqmax',oliqmax)
        CALL getin_p('oicemax',oicemax)
@@ -2423,5 +2431,5 @@
        ! after the call to the convective scheme.
        IF (iflag_wake>=1) then
-          IF (ok_adjwk) THEN
+          IF (iflag_adjwk >= 1) THEN
              limbas(:) = 1
              CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
@@ -2441,7 +2449,9 @@
                 ENDDO
              ENDDO
-             CALL add_wake_tend &
+             IF (iflag_adjwk == 2) THEN
+               CALL add_wake_tend &
                  (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy) 
-          ENDIF  ! (ok_adjwk)
+             ENDIF  ! (iflag_adjwk == 2)
+          ENDIF  ! (iflag_adjwk >= 1)
        ENDIF ! (iflag_wake>=1)
        !>jyg
@@ -2525,13 +2535,15 @@
           !    Add the tendency due to the dry adjustment of the wake profile
           IF (iflag_wake>=1) THEN
-             DO k=1,klev
-                DO i=1,klon
-                   ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
-                   fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
-                   d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
-                   d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
-                ENDDO
-             ENDDO
-          ENDIF
+            IF (iflag_adjwk == 2) THEN
+              DO k=1,klev
+                 DO i=1,klon
+                    ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
+                    fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
+                    d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
+                    d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
+                 ENDDO
+              ENDDO
+            ENDIF  ! (iflag_adjwk = 2)
+          ENDIF   ! (iflag_wake>=1)
           !>jyg
           !
Index: LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_init_140gp.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_init_140gp.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_init_140gp.F90	(revision 2886)
@@ -81,10 +81,16 @@
 !  Calculate lookup tables for functions needed in routine TAUMOL (TAUGB2)
 
+! FH 2017/05/03 
+! Ce facteur de correction CORR2 est vraiment bizare parce qu'on
+! impose 1. aux bornes,  en I=1 et I=200 mais la fonction
+! CORE=( 1 - sqrt(i/im) ) / ( 1 - i/im ) = 1/ ( 1 + sqrt(i/im))
+! vaut 1 en i=1 et 1/2 en i=im ...
+
 CORR1(0) = 1.0_JPRB
-CORR1(400) = 1.0_JPRB
+CORR1(200) = 1.0_JPRB
 CORR2(0) = 1.0_JPRB
-CORR2(400) = 1.0_JPRB
-DO I = 1,399
-  Z_FP = 0.0025_JPRB*REAL(I)
+CORR2(200) = 1.0_JPRB
+DO I = 1,199
+  Z_FP = 0.005_JPRB*REAL(I)
   Z_RTFP = SQRT(Z_FP)
   CORR1(I) = Z_RTFP/Z_FP
Index: LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_taumol2.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_taumol2.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/rrtm/rrtm_taumol2.F90	(revision 2886)
@@ -111,4 +111,12 @@
   IFP=MAX(0,IFP)
 
+! FH 2017/05/02
+! Modification parce qu'on avait un plantage sur un cas 1D.
+! C'est evidemment une correction suspecte
+  IF (IFP>200) THEN
+      PRINT*,'WARNING IFP=',IFP,' 2.E2_JPRB*Z_FP+0.5_JPRB avec Z_FP=',Z_FP
+      IFP=200
+  ENDIF
+
   Z_FC00(I_LAY) = P_FAC00(I_LAY) * CORR2(IFP)
   Z_FC10(I_LAY) = P_FAC10(I_LAY) * CORR2(IFP)
@@ -147,4 +155,8 @@
 !---MI 981104        
   IF (IFP <= 0) IFP=0
+  IF (IFP>200) THEN
+      PRINT*,'WARNING IFP=',IFP,' 2.E2_JPRB*Z_FP+0.5_JPRB avec Z_FP=',Z_FP
+      IFP=200
+  ENDIF
 
   Z_FC00(I_LAY) = P_FAC00(I_LAY) * CORR2(IFP)
Index: LMDZ5/branches/testing/libf/phylmd/rrtm/yoerrtbg2.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/rrtm/yoerrtbg2.F90	(revision 2870)
+++ LMDZ5/branches/testing/libf/phylmd/rrtm/yoerrtbg2.F90	(revision 2886)
@@ -13,6 +13,6 @@
 !    -------------------------------------------------------------------
 
-REAL(KIND=JPRB) :: CORR1(0:400)
-REAL(KIND=JPRB) :: CORR2(0:400)
+REAL(KIND=JPRB) :: CORR1(0:200)
+REAL(KIND=JPRB) :: CORR2(0:200)
 
 !     -----------------------------------------------------------------
