--- trunk/Sources/phylmd/coefkz.f	2016/03/11 18:47:26	178
+++ trunk/Sources/phylmd/coefkz.f	2017/11/06 17:20:45	229
@@ -4,40 +4,37 @@
 
 contains
 
-  SUBROUTINE coefkz(nsrf, knon, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, &
-       t, q, qsurf, coefm, coefh)
+  SUBROUTINE coefkz(nsrf, paprs, pplay, ksta, ksta_ter, ts, rugos, u, v, t, &
+       q, qsurf, coefm, coefh)
 
     ! Authors: F. Hourdin, M. Forichon, Z. X. Li (LMD/CNRS)
-    ! date: 1993/09/22
+    ! Date: September 22nd, 1993
     ! Objet : calculer le coefficient de frottement du sol ("Cdrag") et les
     ! coefficients d'échange turbulent dans l'atmosphère.
 
-    USE indicesol, ONLY: is_oce
+    use clcdrag_m, only: clcdrag
+    USE conf_phys_m, ONLY: iflag_pbl
     USE dimphy, ONLY: klev, klon
+    USE fcttre, ONLY: foede, foeew
+    USE indicesol, ONLY: is_oce
     USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlstt, rlvtt, rtt
     USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
-    USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
-    USE conf_phys_m, ONLY: iflag_pbl
-    use clcdrag_m, only: clcdrag
-
-    ! Arguments:
 
     integer, intent(in):: nsrf ! indicateur de la nature du sol
-    INTEGER, intent(in):: knon ! nombre de points a traiter
 
-    REAL, intent(in):: paprs(klon, klev+1)
+    REAL, intent(in):: paprs(:, :) ! (klon, klev+1)
     ! pression a chaque intercouche (en Pa)
 
-    real, intent(in):: pplay(klon, klev)
+    real, intent(in):: pplay(:, :) ! (klon, klev)
     ! pression au milieu de chaque couche (en Pa)
 
     REAL, intent(in):: ksta, ksta_ter
-    REAL, intent(in):: ts(klon) ! temperature du sol (en Kelvin)
-    REAL, intent(in):: rugos(klon) ! longeur de rugosite (en m)
-    REAL, intent(in):: u(klon, klev), v(klon, klev) ! wind
-    REAL, intent(in):: t(klon, klev) ! temperature (K)
-    real, intent(in):: q(klon, klev) ! vapeur d'eau (kg/kg)
-    real, intent(in):: qsurf(klon) 
+    REAL, intent(in):: ts(:) ! (knon) temperature du sol (en Kelvin)
+    REAL, intent(in):: rugos(:) ! (klon) longeur de rugosite (en m)
+    REAL, intent(in):: u(:, :), v(:, :) ! (klon, klev) wind
+    REAL, intent(in):: t(:, :) ! (klon, klev) temperature (K)
+    real, intent(in):: q(:, :) ! (klon, klev) vapeur d'eau (kg/kg)
+    real, intent(in):: qsurf(:) ! (knon) 
     REAL, intent(out):: coefm(:, :) ! (knon, klev) coefficient, vitesse
 
     real, intent(out):: coefh(:, :) ! (knon, klev) 
@@ -45,7 +42,9 @@
 
     ! Local:
 
-    INTEGER itop(knon) ! numero de couche du sommet de la couche limite
+    INTEGER knon ! nombre de points a traiter
+    INTEGER itop(size(coefm, 1)) ! (knon) numero de couche du sommet
+                                 ! de la couche limite
 
     ! Quelques constantes et options:
 
@@ -56,7 +55,7 @@
     REAL, PARAMETER:: cd = 5.
     REAL, PARAMETER:: clam = 160.
     REAL, PARAMETER:: ratqs = 0.05 ! largeur de distribution de vapeur d'eau
-    
+
     LOGICAL, PARAMETER:: richum = .TRUE.
     ! utilise le nombre de Richardson humide
 
@@ -86,11 +85,12 @@
     REAL zt, zq, zcvm5, zcor, zqs, zfr, zdqs
     logical zdelta
     REAL z2geomf, zalh2, alm2, zscfh, scfm
-    REAL, PARAMETER:: t_coup = 273.15
     REAL gamt(2:klev) ! contre-gradient pour la chaleur sensible: Kelvin/metre
 
     !--------------------------------------------------------------------
 
+    knon = size(coefm, 1)
+
     ! Prescrire la valeur de contre-gradient
     if (iflag_pbl.eq.1) then
        DO k = 3, klev
@@ -132,8 +132,8 @@
        z1(i) = zgeop(i, 1)
     ENDDO
 
-    CALL clcdrag(klon, knon, nsrf, .false., u1, v1, t1, q1, z1, ts, qsurf, &
-         rugos, coefm(:, 1), coefh(:, 1)) 
+    CALL clcdrag(nsrf, u1, v1, t1, q1, z1, ts, qsurf, rugos, coefm(:, 1), &
+         coefh(:, 1)) 
 
     ! Calculer les coefficients turbulents dans l'atmosphere
 
@@ -150,24 +150,14 @@
 
           ! calculer Qs et dQs/dT:
 
-          IF (thermcep) THEN
-             zdelta = RTT >=zt
-             zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
-                  / (1. + RVTMP2*zq)
-             zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
-             zqs = MIN(0.5, zqs)
-             zcor = 1./(1.-RETV*zqs)
-             zqs = zqs*zcor
-             zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
-          ELSE
-             IF (zt  <  t_coup) THEN
-                zqs = qsats(zt) / pplay(i, k)
-                zdqs = dqsats(zt, zqs)
-             ELSE
-                zqs = qsatl(zt) / pplay(i, k)
-                zdqs = dqsatl(zt, zqs)
-             ENDIF
-          ENDIF
+          zdelta = RTT >=zt
+          zcvm5 = merge(R5IES * RLSTT, R5LES * RLVTT, zdelta) / RCPD &
+               / (1. + RVTMP2*zq)
+          zqs = R2ES * FOEEW(zt, zdelta) / pplay(i, k)
+          zqs = MIN(0.5, zqs)
+          zcor = 1./(1.-RETV*zqs)
+          zqs = zqs*zcor
+          zdqs = FOEDE(zt, zdelta, zcvm5, zqs, zcor)
 
           ! calculer la fraction nuageuse (processus humide):