Index: LMDZ4/trunk/libf/phylmd/climb_wind_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/climb_wind_mod.F90	(revision 793)
+++ LMDZ4/trunk/libf/phylmd/climb_wind_mod.F90	(revision 1067)
@@ -1,4 +1,2 @@
-!
-! $Header$
 !
 MODULE climb_wind_mod
@@ -21,9 +19,13 @@
   REAL, DIMENSION(:,:), ALLOCATABLE  :: Ccoef_V, Dcoef_V
   !$OMP THREADPRIVATE(Ccoef_V, Dcoef_V)
+  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_U, Bcoef_U
+  !$OMP THREADPRIVATE(Acoef_U, Bcoef_U)
+  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_V, Bcoef_V
+  !$OMP THREADPRIVATE(Acoef_V, Bcoef_V)
   LOGICAL                            :: firstcall=.TRUE.
   !$OMP THREADPRIVATE(firstcall)
 
   
-  PUBLIC :: climb_wind_down, calcul_wind_flux, climb_wind_up
+  PUBLIC :: climb_wind_down, climb_wind_up
 
 CONTAINS
@@ -62,4 +64,7 @@
     IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Dcoef_V',1)
 
+    ALLOCATE(Acoef_U(klon), Bcoef_U(klon), Acoef_V(klon), Bcoef_V(klon), STAT=ierr)
+    IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_U and Bcoef_U, ierr=', ierr
+
     firstcall=.FALSE.
 
@@ -68,5 +73,6 @@
 !****************************************************************************************
 !
-  SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old)
+  SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old, &
+       Acoef_U_out, Acoef_V_out, Bcoef_U_out, Bcoef_V_out)
 !
 ! This routine calculates for the wind components u and v,
@@ -88,4 +94,11 @@
     REAL, DIMENSION(klon,klev), INTENT(IN)   :: v_old
 
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_U_out
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_V_out
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_U_out
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_V_out
+
 ! Local variables
 !****************************************************************************************
@@ -102,26 +115,14 @@
 !
 !****************************************************************************************
-
 ! - Define alpha (alf1 and alf2) 
     alf1(:) = 1.0
     alf2(:) = 1.0 - alf1(:)
 
-! - Calculte the wind components for the first layer
-    u1lay(1:knon) = u_old(1:knon,1)*alf1(1:knon) + u_old(1:knon,2)*alf2(1:knon)
-    v1lay(1:knon) = v_old(1:knon,1)*alf1(1:knon) + v_old(1:knon,2)*alf2(1:knon)    
-
-! - Calculate K
+! - Calculate the coefficients K
     Kcoefm(:,:) = 0.0
-    DO i = 1, knon
-       Kcoefm(i,1) = coef_in(i,1) * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))* &
-            pplay(i,1)/(RD*temp(i,1)) 
-       Kcoefm(i,1) = Kcoefm(i,1) * dtime*RG
-    END DO
-
     DO k = 2, klev
        DO i=1,knon
-          Kcoefm(i,k) = coef_in(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &
+          Kcoefm(i,k) = coef_in(i,k)*RG*RG*dtime/(pplay(i,k-1)-pplay(i,k)) &
                *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
-          Kcoefm(i,k) = Kcoefm(i,k) * dtime*RG
        END DO
     END DO
@@ -130,10 +131,20 @@
     CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
          u_old(:,:), alf1(:), alf2(:),  &
-         Ccoef_U(:,:), Dcoef_U(:,:))
+         Ccoef_U(:,:), Dcoef_U(:,:), Acoef_U(:), Bcoef_U(:))
 
 ! - Calculate the coefficients C and D, component "v"
     CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
          v_old(:,:), alf1(:), alf2(:),  &
-         Ccoef_V(:,:), Dcoef_V(:,:))
+         Ccoef_V(:,:), Dcoef_V(:,:), Acoef_V(:), Bcoef_V(:))
+
+!****************************************************************************************
+! 6)
+! Return the first layer in output variables
+!
+!****************************************************************************************
+    Acoef_U_out = Acoef_U
+    Bcoef_U_out = Bcoef_U
+    Acoef_V_out = Acoef_V
+    Bcoef_V_out = Bcoef_V
 
   END SUBROUTINE climb_wind_down
@@ -141,12 +152,12 @@
 !****************************************************************************************
 !
-  SUBROUTINE calc_coef(knon, Kcoef, dels, X, alfa1, alfa2, Ccoef, Dcoef)
-!
-! Find the coefficients C and D in fonction of alfa, K and dels
+  SUBROUTINE calc_coef(knon, Kcoef, delp, X, alfa1, alfa2, Ccoef, Dcoef, Acoef, Bcoef)
+!
+! Find the coefficients C and D in fonction of alfa, K and delp
 !
 ! Input arguments
 !****************************************************************************************
     INTEGER, INTENT(IN)                      :: knon
-    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, dels
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
     REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
     REAL, DIMENSION(klon), INTENT(IN)        :: alfa1, alfa2
@@ -154,4 +165,5 @@
 ! Output arguments
 !****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
     REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
   
@@ -161,7 +173,9 @@
     REAL                                     :: buf
 
+    INCLUDE "YOMCST.h"
 !****************************************************************************************
 ! 
-! Niveau au sommet, k=klev
+
+! Calculate coefficients C and D at top level, k=klev
 !
     Ccoef(:,:) = 0.0
@@ -169,37 +183,29 @@
 
     DO i = 1, knon
-       buf = dels(i,klev) + Kcoef(i,klev)
-
-       Ccoef(i,klev) = X(i,klev)*dels(i,klev)/buf 
+       buf = delp(i,klev) + Kcoef(i,klev)
+
+       Ccoef(i,klev) = X(i,klev)*delp(i,klev)/buf 
        Dcoef(i,klev) = Kcoef(i,klev)/buf
     END DO
     
 ! 
-! Niveau  (klev-1) <= k <= 2
+! Calculate coefficients C and D at top level (klev-1) <= k <= 2
 !
     DO k=(klev-1),2,-1
        DO i = 1, knon
-          buf = dels(i,k) + Kcoef(i,k) + &
-               Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
+          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
           
-          Ccoef(i,k) = X(i,k)*dels(i,k) + &
-               Kcoef(i,k+1)*Ccoef(i,k+1)
-          Ccoef(i,k) = Ccoef(i,k)/buf         
+          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1))/buf
           Dcoef(i,k) = Kcoef(i,k)/buf
        END DO
     END DO
 
-! 
-! Niveau k=1
+!
+! Calculate coeffiecent A and B at surface
 !
     DO i = 1, knon
-       buf = dels(i,1) + &
-            (alfa1(i) + alfa2(i)*Dcoef(i,2)) * Kcoef(i,1) + &
-            (1.-Dcoef(i,2))*Kcoef(i,2)
-       
-       Ccoef(i,1) = X(i,1)*dels(i,1) + &
-            (Kcoef(i,2)-Kcoef(i,1)*alfa2(i)) * Ccoef(i,2)
-       Ccoef(i,1) = Ccoef(i,1)/buf
-       Dcoef(i,1) = Kcoef(i,1)/buf
+       buf = delp(i,1) + Kcoef(i,2)*(1-Dcoef(i,2))
+       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*Ccoef(i,2))/buf
+       Bcoef(i) = -RG/buf
     END DO
 
@@ -208,71 +214,14 @@
 !****************************************************************************************
 !
-  SUBROUTINE calcul_wind_flux(knon, dtime, &
-       flux_u, flux_v)
-
-    INCLUDE "YOMCST.h"
-
-! Input arguments
-!****************************************************************************************
-    INTEGER, INTENT(IN)                  :: knon
-    REAL, INTENT(IN)                     :: dtime
-
-! Output arguments
-!****************************************************************************************
-    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u
-    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v
-
-! Local variables
-!****************************************************************************************
-    INTEGER                              :: i
-    REAL, DIMENSION(klon)                :: u0, v0  ! u and v at niveau 0
-    REAL, DIMENSION(klon)                :: u1, v1  ! u and v at niveau 1
-    REAL, DIMENSION(klon)                :: u2, v2  ! u and v at niveau 2
-    
-
-!****************************************************************************************
-! Les vents de surface sont supposes nuls
-!
-!****************************************************************************************
-    u0(:) = 0.0
-    v0(:) = 0.0
-
-!****************************************************************************************
-! On calcule les vents du couhes 1 et 2 recurviement
-!
-!****************************************************************************************
-    DO i = 1, knon
-       u1(i) = Ccoef_U(i,1) + Dcoef_U(i,1)*u0(i)
-       v1(i) = Ccoef_V(i,1) + Dcoef_V(i,1)*v0(i)
-       u2(i) = Ccoef_U(i,2) + Dcoef_U(i,2)*u1(i)
-       v2(i) = Ccoef_V(i,2) + Dcoef_V(i,2)*v1(i)
-    END DO
-
-!****************************************************************************************
-! On calcule le flux
-!
-!****************************************************************************************
-    flux_u(:) = 0.0
-    flux_v(:) = 0.0
-
-    DO i=1,knon
-       flux_u(i) = Kcoefm(i,1)/RG/dtime * (u1(i)*alf1(i) + u2(i)*alf2(i) - u0(i)) 
-       flux_v(i) = Kcoefm(i,1)/RG/dtime * (v1(i)*alf1(i) + v2(i)*alf2(i) - v0(i)) 
-    END DO
-
-  END SUBROUTINE calcul_wind_flux
-!
-!****************************************************************************************
-!
-  SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, &
+
+  SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, flx_u1, flx_v1,  &
        flx_u_new, flx_v_new, d_u_new, d_v_new)
 !
-! Diffuse the wind components from the surface and up to the top layer. Coefficents
-! C and D are known from before. The values for U and V at surface are supposed to be 
-! zero (this could be modified).
-!
-! u(k) = Cu(k) + Du(k)*u(k-1)
-! v(k) = Cv(k) + Dv(k)*v(k-1)
-! [1 <= k <= klev]
+! Diffuse the wind components from the surface layer and up to the top layer. 
+! Coefficents A, B, C and D are known from before. Start values for the diffusion are the
+! momentum fluxes at surface.
+!
+! u(k=1) = A + B*flx*dtime
+! u(k)   = C(k) + D(k)*u(k-1)  [2 <= k <= klev]
 !
 !****************************************************************************************
@@ -285,4 +234,5 @@
     REAL, DIMENSION(klon,klev), INTENT(IN)  :: u_old
     REAL, DIMENSION(klon,klev), INTENT(IN)  :: v_old
+    REAL, DIMENSION(klon), INTENT(IN)       :: flx_u1, flx_v1 ! momentum flux
 
 ! Output arguments
@@ -294,5 +244,4 @@
 !****************************************************************************************
     REAL, DIMENSION(klon,klev)              :: u_new, v_new
-    REAL, DIMENSION(klon)                   :: u0, v0
     INTEGER                                 :: k, i
     
@@ -300,12 +249,8 @@
 !****************************************************************************************
 
-! Niveau 0
-    u0(1:knon) = 0.0
-    v0(1:knon) = 0.0
-
 ! Niveau 1
     DO i = 1, knon
-       u_new(i,1) = Ccoef_U(i,1) + Dcoef_U(i,1) * u0(i)
-       v_new(i,1) = Ccoef_V(i,1) + Dcoef_V(i,1) * v0(i)
+       u_new(i,1) = Acoef_U(i) + Bcoef_U(i)*flx_u1(i)*dtime
+       v_new(i,1) = Acoef_V(i) + Bcoef_V(i)*flx_v1(i)*dtime
     END DO
 
@@ -329,11 +274,6 @@
     flx_v_new(:,:) = 0.0
 
-! Niveau 1
-    DO i = 1, knon
-       flx_u_new(i,1) = Kcoefm(i,1)/RG/dtime * &
-            (u_new(i,1)*alf1(i)+u_new(i,2)*alf2(i) - u0(i)) 
-       flx_v_new(i,1) = Kcoefm(i,1)/RG/dtime * &
-            (v_new(i,1)*alf1(i)+v_new(i,2)*alf2(i) - v0(i)) 
-    END DO
+    flx_u_new(1:knon,1)=flx_u1(1:knon)
+    flx_v_new(1:knon,1)=flx_v1(1:knon)
 
 ! Niveau 2->klev
