Index: LMDZ5/trunk/libf/phylmd/cdrag.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cdrag.F90	(revision 2232)
+++ LMDZ5/trunk/libf/phylmd/cdrag.F90	(revision 2232)
@@ -0,0 +1,225 @@
+!
+!$Id: cdrag.F90 ???? 2015-??-?? ??:??:??Z ? $
+!
+ SUBROUTINE cdrag( knon,  nsrf,   &
+     speed, t1,    q1,    zgeop1, &
+     psol,  tsurf, qsurf, rugos,  &
+     pcfm,  pcfh,  zri,   pref )
+
+  USE dimphy
+  USE indice_sol_mod
+  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.
+!
+! 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.
+!   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.
+!
+! ================================================================= 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
+! rugos---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
+!
+!
+! Parametres d'entree
+!***************************************************************** 
+  INTEGER, INTENT(IN)                      :: knon, nsrf
+  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)        :: rugos ! 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
+
+! 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)       :: 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
+!
+  INCLUDE "YOMCST.h"
+  INCLUDE "YOETHF.h"
+!  INCLUDE "indicesol.h"
+  INCLUDE "clesphys.h"
+  INCLUDE "iniprint.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         ! Drag coefficient in neutral conditions
+!
+! 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)
+
+! ================================================================= 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
+!
+! Fuxing WANG, 04/03/2015
+! To check if there are negative q1, qsurf values.
+  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) 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_gcm(modname,abort_message,1)
+  ENDIF
+  IF (ng_qsurf.GT.0) 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_gcm(modname,abort_message,1)
+  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)
+!------------------------------------------------------------
+     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*q1(i))) &
+          *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero
+     zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd)
+     zcdn(i) = (CKAP/LOG(1.+zgeop1(i)/(RG*rugos(i))))**2
+
+     IF (zri(i) .GT. 0.) THEN      ! situation stable
+        zri(i) = MIN(20.,zri(i))
+        IF (.NOT.zxli) THEN
+           zscf = SQRT(1.+CD*ABS(zri(i)))
+           friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min)
+           zcfm1(i) = zcdn(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(i) * frih
+           IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn(i) * frih
+!!$ PB
+           pcfm(i) = zcfm1(i)
+           pcfh(i) = zcfh1(i)
+        ELSE
+           pcfm(i) = zcdn(i)* fsta(zri(i))
+           pcfh(i) = zcdn(i)* fsta(zri(i))
+        ENDIF
+     ELSE                          ! situation instable
+        IF (.NOT.zxli) THEN
+           zucf = 1./(1.+3.0*CB*CC*zcdn(i)*SQRT(ABS(zri(i)) &
+                *(1.0+zgeop1(i)/(RG*rugos(i)))))
+           zcfm2(i) = zcdn(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min)
+!!$ PB     zcfh2(i) = zcdn(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min)
+           zcfh2(i) = f_cdrag_ter*zcdn(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(i)* fins(zri(i))
+           pcfh(i) = zcdn(i)* fins(zri(i))
+        ENDIF
+! cdrah sur l'ocean cf. Miller et al. (1992)
+        zcr = (0.0016/(zcdn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.)
+        IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn(i)*(1.0+zcr**1.25)**(1./1.25)
+     ENDIF
+  END DO
+
+! ================================================================= 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/trunk/libf/phylmd/clcdrag.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/clcdrag.F90	(revision 2231)
+++ LMDZ5/trunk/libf/phylmd/clcdrag.F90	(revision 2232)
@@ -61,4 +61,9 @@
   REAL, DIMENSION(klon) :: zgeop1       ! geopotentiel au 1er niveau du modele
   LOGICAL, PARAMETER    :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li
+
+  CHARACTER (LEN=80) :: abort_message
+  CHARACTER (LEN=20) :: modname = 'clcdrag'
+
+
 !
 ! Fonctions thermodynamiques et fonctions d'instabilite
@@ -66,4 +71,9 @@
   fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
   fins(x) = SQRT(1.0-18.0*x)
+
+  abort_message='obsolete, remplace par cdrag, use at you own risk'
+  CALL abort_gcm(modname,abort_message,1)
+
+
 
 ! ================================================================= c
Index: LMDZ5/trunk/libf/phylmd/coefcdrag.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/coefcdrag.F90	(revision 2231)
+++ LMDZ5/trunk/libf/phylmd/coefcdrag.F90	(revision 2232)
@@ -64,4 +64,11 @@
       REAL, dimension(klon) :: trm0, trm1
 
+  CHARACTER (LEN=80) :: abort_message
+  CHARACTER (LEN=20) :: modname = 'coefcdra'
+
+
+!
+
+
 !-------------------------------------------------------------------------
       REAL :: fsta, fins, x
@@ -69,4 +76,8 @@
       fins(x) = SQRT(1.0-18.0*x)
 !-------------------------------------------------------------------------
+
+  abort_message='obsolete, remplace par cdrag, use at you own risk'
+  CALL abort_gcm(modname,abort_message,1)
+
 !
       DO i = 1, knon
Index: LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 2231)
+++ LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 2232)
@@ -553,5 +553,5 @@
     REAL, DIMENSION(klon)              :: tair1, qair1, tairsol
     REAL, DIMENSION(klon)              :: psfce, patm
-    REAL, DIMENSION(klon)              :: qairsol, zgeo1
+    REAL, DIMENSION(klon)              :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015
     REAL, DIMENSION(klon)              :: rugo1
     REAL, DIMENSION(klon)              :: yfluxsens
@@ -714,8 +714,8 @@
     REAL, DIMENSION(klon)       :: ytrmb3_w
 !
-    REAL, DIMENSION(klon)              :: uzon_x, vmer_x
+    REAL, DIMENSION(klon)              :: uzon_x, vmer_x, speed_x, zri1_x, pref_x !speed_x, zri1_x, pref_x, added by Fuxing WANG, 04/03/2015
     REAL, DIMENSION(klon)              :: zgeo1_x, tair1_x, qair1_x, tairsol_x
 !
-    REAL, DIMENSION(klon)              :: uzon_w, vmer_w
+    REAL, DIMENSION(klon)              :: uzon_w, vmer_w, speed_w, zri1_w, pref_w !speed_w, zri1_w, pref_w, added by Fuxing WANG, 04/03/2015
     REAL, DIMENSION(klon)              :: zgeo1_w, tair1_w, qair1_w, tairsol_w
 
@@ -1362,8 +1362,20 @@
 !!!
 !!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
-        CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
-            yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
+! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
+!       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
+!           yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
+!           yts, yqsurf, yrugos, &
+!           ycdragm, ycdragh )
+! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
+        DO i = 1, klon
+           zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
+                * (ypaprs(i,1)-ypplay(i,1))
+           speed(i) = SQRT(yu(i,1)**2+yv(i,1)**2)
+        END DO
+        CALL cdrag(knon, nsrf, &
+            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
             yts, yqsurf, yrugos, &
-            ycdragm, ycdragh )
+            ycdragm, ycdragh, zri1, pref )
+
 ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
      IF (ok_prescr_ust) then
@@ -1380,8 +1392,21 @@
         IF (prt_level >=10) print *,'clcdrag -> ycdragh ', ycdragh
        ELSE  !(iflag_split .eq.0)
-        CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
-            yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
+
+! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
+!       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
+!           yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
+!           yts_x, yqsurf, yrugos, &
+!           ycdragm_x, ycdragh_x )
+! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
+        DO i = 1, knon
+           zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
+                * (ypaprs(i,1)-ypplay(i,1))
+           speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2)
+        END DO
+        CALL cdrag(knon, nsrf, &
+            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
             yts_x, yqsurf, yrugos, &
-            ycdragm_x, ycdragh_x )
+            ycdragm_x, ycdragh_x, zri1_x, pref_x )
+
 ! --- special Dice. JYG+MPL 25112013
         IF (ok_prescr_ust) then
Index: LMDZ5/trunk/libf/phylmd/screenc.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/screenc.F90	(revision 2231)
+++ LMDZ5/trunk/libf/phylmd/screenc.F90	(revision 2232)
@@ -66,11 +66,15 @@
 ! Richardson at reference level 
 !
-      CALL coefcdrag (klon, knon, nsrf, zxli, &
+!      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, rugos, &
-                    okri, ri1, &
-                    cdram, cdrah, cdran, zri1, &
-                    pref)
-!
+                    cdram, cdrah, zri1, pref)
       DO i = 1, knon
         delu(i) = ustar(i)/sqrt(cdram(i))
Index: LMDZ5/trunk/libf/phylmd/stdlevvar.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/stdlevvar.F90	(revision 2231)
+++ LMDZ5/trunk/libf/phylmd/stdlevvar.F90	(revision 2232)
@@ -98,8 +98,14 @@
 !
       okri=.FALSE.
-      CALL coefcdrag(klon, knon, nsrf, zxli, &
- &                   speed, t1, q1, z1, psol, &
- &                   ts1, qsurf, rugos, okri, ri1,  &         
- &                   cdram, cdrah, cdran, zri1, pref)            
+!      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, rugos, &
+ &                   cdram, cdrah, zri1, pref)
+
 ! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
      IF (ok_prescr_ust) then
