--- trunk/phylmd/Interface_surf/yamada4.f	2018/07/24 15:22:48	286
+++ trunk/phylmd/Interface_surf/yamada4.f	2018/09/06 14:25:07	303
@@ -8,18 +8,19 @@
 
 contains
 
-  SUBROUTINE yamada4(dt, zlev, zlay, u, v, teta, q2, km, kn, ustar)
+  SUBROUTINE yamada4(zlev, zlay, u, v, teta, q2, coefm, coefh, ustar)
 
     ! From LMDZ4/libf/phylmd/yamada4.F, version 1.1 2004/06/22 11:45:36
 
+    ! Library:
+    use nr_util, only: assert, assert_eq
+
+    use comconst, only: dtphys
     USE conf_phys_m, ONLY: iflag_pbl
     USE dimphy, ONLY: klev
-    use nr_util, only: assert, assert_eq
     USE suphec_m, ONLY: rg
 
-    REAL, intent(in):: dt ! pas de temps
-
-    REAL zlev(:, :) ! (knon, klev + 1)
+    REAL, intent(in):: zlev(:, :) ! (knon, klev + 1)
     ! altitude \`a chaque niveau (interface inf\'erieure de la couche de
     ! m\^eme indice)
 
@@ -39,17 +40,18 @@
     ! En entr\'ee : la valeur au d\'ebut du pas de temps ; en sortie : la
     ! valeur \`a la fin du pas de temps.
 
-    REAL, intent(out):: km(:, 2:) ! (knon, 2:klev)
+    REAL, intent(out):: coefm(:, 2:) ! (knon, 2:klev)
     ! diffusivit\'e turbulente de quantit\'e de mouvement (au bas de
     ! chaque couche) (en sortie : la valeur \`a la fin du pas de temps)
 
-    REAL, intent(out):: kn(:, 2:) ! (knon, 2:klev)
+    REAL, intent(out):: coefh(:, 2:) ! (knon, 2:klev)
     ! diffusivit\'e turbulente des scalaires (au bas de chaque couche)
     ! (en sortie : la valeur \`a la fin du pas de temps)
 
     real, intent(in):: ustar(:) ! (knon)
 
     ! Local:
+    
     integer knon
     real kmin, qmin
     real pblhmin(size(ustar)), coriol(size(ustar)) ! (knon)
@@ -58,9 +60,7 @@
     REAL unsdzdec(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
     real delta(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
     real aa(size(zlev, 1), size(zlev, 2)) ! (knon, klev + 1)
-    real aa1
     logical:: first = .true.
-    integer:: ipas = 0
     integer ig, k
     real ri
     real, dimension(size(zlev, 1), size(zlev, 2)):: rif, sm ! (knon, klev + 1)
@@ -82,19 +82,13 @@
 
     call assert(any(iflag_pbl == [6, 8, 9]), "yamada4 iflag_pbl")
     knon = assert_eq([size(zlev, 1), size(zlay, 1), size(u, 1), size(v, 1), &
-         size(teta, 1), size(ustar), size(q2, 1), size(km, 1), size(kn, 1)], &
-         "yamada4 knon")
+         size(teta, 1), size(ustar), size(q2, 1), size(coefm, 1), &
+         size(coefh, 1)], "yamada4 knon")
     call assert(klev == [size(zlev, 2) - 1, size(zlay, 2), size(u, 2), &
-         size(v, 2), size(teta, 2), size(q2, 2) - 1, size(km, 2) + 1, &
-         size(kn, 2) + 1], "yamada4 klev")
-
-    ipas = ipas + 1
+         size(v, 2), size(teta, 2), size(q2, 2) - 1, size(coefm, 2) + 1, &
+         size(coefh, 2) + 1], "yamada4 klev")
 
     ! les increments verticaux
-    DO ig = 1, knon
-       ! alerte: zlev n'est pas declare a klev + 1
-       zlev(ig, klev + 1) = zlay(ig, klev) + (zlay(ig, klev) - zlev(ig, klev))
-    ENDDO
 
     DO k = 1, klev
        DO ig = 1, knon
@@ -212,16 +206,15 @@
     else if (iflag_pbl >= 8) then
        ! Yamada 2.5 a la Didi
 
-       ! Calcul de l, km, au pas precedent
+       ! Calcul de l, coefm, au pas precedent
        do k = 2, klev
           do ig = 1, knon
              delta(ig, k) = q2(ig, k)/(l(ig, k)**2*sm(ig, k))
              if (delta(ig, k).lt.1.e-20) then
                 delta(ig, k) = 1.e-20
              endif
-             km(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k)
-             aa1 = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)
-             aa(ig, k) = aa1*dt/(delta(ig, k)*l(ig, k))
+             coefm(ig, k) = l(ig, k)*sqrt(q2(ig, k))*sm(ig, k)
+             aa(ig, k) = (m2(ig, k)*(1.-rif(ig, k))-delta(ig, k)/b1)*dtphys/(delta(ig, k)*l(ig, k))
              qpre = sqrt(q2(ig, k))
              if (iflag_pbl == 8) then
                 if (aa(ig, k).gt.0.) then
@@ -246,8 +239,8 @@
     do k = 2, klev
        do ig = 1, knon
           zq = sqrt(q2(ig, k))
-          km(ig, k) = l(ig, k)*zq*sm(ig, k)
-          kn(ig, k) = km(ig, k)*alpha(ig, k)
+          coefm(ig, k) = l(ig, k)*zq*sm(ig, k)
+          coefh(ig, k) = coefm(ig, k)*alpha(ig, k)
        enddo
     enddo
 
@@ -270,9 +263,9 @@
           else
              kmin = -1. ! kmin n'est utilise que pour les SL stables.
           endif
-          if (kn(ig, k).lt.kmin.or.km(ig, k).lt.kmin) then
-             kn(ig, k) = kmin
-             km(ig, k) = kmin
+          if (coefh(ig, k).lt.kmin.or.coefm(ig, k).lt.kmin) then
+             coefh(ig, k) = kmin
+             coefm(ig, k) = kmin
              ! la longueur de melange est suposee etre l = kap z
              ! K = l q Sm d'ou q2 = (K/l Sm)**2
              q2(ig, k) = (qmin/sm(ig, k))**2
@@ -286,7 +279,7 @@
 
   !*******************************************************************
 
-  real function frif(ri)
+  pure real function frif(ri)
 
     real, intent(in):: ri
 
@@ -296,7 +289,7 @@
 
   !*******************************************************************
 
-  real function falpha(ri)
+  pure real function falpha(ri)
 
     real, intent(in):: ri
 
@@ -306,7 +299,7 @@
 
   !*******************************************************************
 
-  real function fsm(ri)
+  pure real function fsm(ri)
 
     real, intent(in):: ri
 
@@ -316,7 +309,7 @@
 
   !*******************************************************************
 
-  real function fl(zzz, zl0, zq2, zn2)
+  pure real function fl(zzz, zl0, zq2, zn2)
 
     real, intent(in):: zzz, zl0, zq2, zn2