Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/acama_gwd_rando_m.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/acama_gwd_rando_m.F90	(revision 3191)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/acama_gwd_rando_m.F90	(revision 3200)
@@ -20,4 +20,7 @@
     use dimphy, only: klon, klev
     use assert_m, only: assert
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE vertical_layers_mod, ONLY : presnivs
+
     include "YOMCST.h"
     include "clesphys.h"
@@ -111,4 +114,27 @@
     REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
     REAL BVSEC ! Security to avoid negative BVF
+
+    REAL, DIMENSION(klev+1) ::HREF
+    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
+    LOGICAL, SAVE :: firstcall = .TRUE.
+  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
+
+    CHARACTER (LEN=20) :: modname='flott_gwd_rando'
+    CHARACTER (LEN=80) :: abort_message
+
+
+
+  IF (firstcall) THEN
+    ! Cle introduite pour resoudre un probleme de non reproductibilite
+    ! Le but est de pouvoir tester de revenir a la version precedenete
+    ! A eliminer rapidement
+    CALL getin_p('gwd_reproductibilite_mpiomp',gwd_reproductibilite_mpiomp)
+    IF (NW+4*(NA-1)+NA>=KLEV) THEN
+       abort_message = 'NW+3*NA>=KLEV Probleme pour generation des ondes'
+       CALL abort_physic (modname,abort_message,1)
+    ENDIF
+    firstcall=.false.
+!    CALL iophys_ini
+  ENDIF
 
     !-----------------------------------------------------------------
@@ -205,11 +231,22 @@
     ! Launching altitude
 
+    IF (gwd_reproductibilite_mpiomp) THEN
+       ! Reprend la formule qui calcule PH en fonction de PP=play
+       DO LL = 2, KLEV
+          HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1))) / 2.)
+       end DO
+       HREF(KLEV + 1) = 0.
+       HREF(1) = 2. * presnivs(1) - HREF(2)
+    ELSE
+       HREF(1:KLEV)=PH(KLON/2,1:KLEV)
+    ENDIF
+
     LAUNCH=0
     LTROP =0
     DO LL = 1, KLEV
-       IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XLAUNCH) LAUNCH = LL
+       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
     ENDDO
     DO LL = 1, KLEV
-       IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XTROP) LTROP = LL
+       IF (HREF(LL) / HREF(1) > XTROP) LTROP = LL
     ENDDO
 
@@ -293,8 +330,5 @@
 
     JW = 0
-    DO JP = 1, NP
-       DO JK = 1, NK
-          DO JO = 1, NO
-             JW = JW + 1
+    DO JW = 1, NW
              ! Angle
              DO II = 1, KLON
@@ -340,6 +374,4 @@
                 ! RUW0(JW, II) = RUWFRT
              ENDDO
-          end DO
-       end DO
     end DO
 
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/flott_gwd_rando_m.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/flott_gwd_rando_m.F90	(revision 3191)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/flott_gwd_rando_m.F90	(revision 3200)
@@ -18,4 +18,7 @@
       use dimphy, only: klon, klev
       use assert_m, only: assert
+      USE ioipsl_getin_p_mod, ONLY : getin_p
+      USE vertical_layers_mod, ONLY : presnivs
+
       include "YOMCST.h"
       include "clesphys.h"
@@ -103,7 +106,31 @@
     REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
     REAL PSEC ! Security to avoid division by 0 pressure
-    REAL PHM1(KLON, KLEV + 1) ! 1/Press at 1/2 levels
     REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
     REAL BVSEC ! Security to avoid negative BVF
+    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
+
+    REAL, DIMENSION(klev+1) ::HREF
+
+    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.
+    LOGICAL, SAVE :: firstcall = .TRUE.
+  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)
+
+    CHARACTER (LEN=20) :: modname='flott_gwd_rando'
+    CHARACTER (LEN=80) :: abort_message
+
+
+
+  IF (firstcall) THEN
+    ! Cle introduite pour resoudre un probleme de non reproductibilite
+    ! Le but est de pouvoir tester de revenir a la version precedenete
+    ! A eliminer rapidement
+    CALL getin_p('gwd_reproductibilite_mpiomp',gwd_reproductibilite_mpiomp)
+    IF (NW+3*NA>=KLEV) THEN
+       abort_message = 'NW+3*NA>=KLEV Probleme pour generation des ondes'
+       CALL abort_physic (modname,abort_message,1)
+    ENDIF
+    firstcall=.false.
+  ENDIF
+
 
     !-----------------------------------------------------------------
@@ -156,4 +183,5 @@
     ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ
 
+IF (1==0) THEN
     !ONLINE
         call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &
@@ -167,4 +195,5 @@
           "FLOTT_GWD_RANDO klev")
     !END ONLINE
+ENDIF
 
     IF(DELTAT < DTIME)THEN
@@ -183,20 +212,29 @@
     DO LL = 2, KLEV
        PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
-       PHM1(:, LL) = 1. / PH(:, LL) 
-    end DO
-
+    end DO
     PH(:, KLEV + 1) = 0. 
-    PHM1(:, KLEV + 1) = 1. / PSEC
     PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
 
     ! Launching altitude
+
+    !Pour revenir a la version non reproductible en changeant le nombre de process
+    IF (gwd_reproductibilite_mpiomp) THEN
+       ! Reprend la formule qui calcule PH en fonction de PP=play
+       DO LL = 2, KLEV
+          HREF(LL) = EXP((LOG(presnivs(LL)) + LOG(presnivs(LL - 1))) / 2.)
+       end DO
+       HREF(KLEV + 1) = 0.
+       HREF(1) = 2. * presnivs(1) - HREF(2)
+    ELSE
+       HREF(1:KLEV)=PH(KLON/2,1:KLEV)
+    ENDIF
 
     LAUNCH=0
     LTROP =0
     DO LL = 1, KLEV
-       IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XLAUNCH) LAUNCH = LL
+       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
     ENDDO
     DO LL = 1, KLEV
-       IF (PH(KLON / 2, LL) / PH(KLON / 2, 1) > XTROP) LTROP = LL
+       IF (HREF(LL) / HREF(1) > XTROP) LTROP = LL
     ENDDO
 
@@ -245,21 +283,20 @@
     ! waves characteristics in an almost stochastic way
 
-    JW = 0
-    DO JP = 1, NP
-       DO JK = 1, NK
-          DO JO = 1, NO
-             JW = JW + 1
+    DO JW = 1, NW
              ! Angle
              DO II = 1, KLON
                 ! Angle (0 or PI so far)
-                ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &
+                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
+                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
+                ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
                      * RPI / 2.
                 ! Horizontal wavenumber amplitude
-                ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
+                ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
                 ! Horizontal phase speed
                 CPHA = 0.
                 DO JJ = 1, NA
+                    RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.)
                     CPHA = CPHA + &
-         CMAX*2.*(MOD(TT(II, JW+3*JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.)
+                    CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.)
                 END DO
                 IF (CPHA.LT.0.)  THEN
@@ -276,7 +313,5 @@
                 RUW0(JW, II) = RUWMAX
              ENDDO
-          end DO
-       end DO
-    end DO
+    ENDDO
 
     ! 4. COMPUTE THE FLUXES
@@ -417,4 +452,5 @@
     ENDDO
 
+
   END SUBROUTINE FLOTT_GWD_RANDO
 
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/pbl_surface_mod.F90	(revision 3191)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/pbl_surface_mod.F90	(revision 3200)
@@ -1266,5 +1266,5 @@
        ENDDO
 !!! jyg le 07/02/2012 et le 10/04/2013
-        DO k = 1, klev
+        DO k = 1, klev+1
           DO j = 1, knon
              i = ni(j)
@@ -1272,5 +1272,10 @@
 !!             ytke(j,k)   = tke(i,k,nsrf)
              ytke(j,k)   = tke_x(i,k,nsrf)
+          ENDDO
+        ENDDO
 !>jyg
+        DO k = 1, klev
+          DO j = 1, knon
+             i = ni(j)
 !FC
              y_treedrg(j,k) =  treedrg(i,k,nsrf)
@@ -2398,5 +2403,5 @@
        IF (iflag_split .eq.0) THEN
         wake_dltke(:,:,nsrf) = 0.
-        DO k = 1, klev
+        DO k = 1, klev+1
            DO j = 1, knon
               i = ni(j)
@@ -2411,5 +2416,5 @@
 
        ELSE  ! (iflag_split .eq.0)
-        DO k = 1, klev
+        DO k = 1, klev+1
           DO j = 1, knon
             i = ni(j)
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3191)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3200)
@@ -4373,5 +4373,5 @@
 
 
-    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pbl_tke)
+    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
 
 
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/tend_to_tke.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/tend_to_tke.F90	(revision 3191)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/tend_to_tke.F90	(revision 3200)
@@ -32,5 +32,5 @@
 !**************************************************************************************
 
- SUBROUTINE tend_to_tke(dt,plev,exner,temp,windu,windv,dt_a,du_a,dv_a,tke)
+ SUBROUTINE tend_to_tke(dt,plev,exner,temp,windu,windv,dt_a,du_a,dv_a,pctsrf,tke)
 
  USE dimphy, ONLY: klon, klev
@@ -53,4 +53,5 @@
   REAL du_a(klon,klev)      ! Zonal wind speed tendency [m/s], grid-cell average or for a one subsurface
   REAL dv_a(klon,klev)      ! Meridional wind speed tendency [m/s], grid-cell average or for a one subsurface
+  REAL pctsrf(klon,nbsrf+1)       ! Turbulent Kinetic energy [m2/s2], grid-cell average or for a subsurface
 
 ! Inputs/Outputs
@@ -119,14 +120,14 @@
 
 
- DO isrf=1,nbsrf
+ DO isrf=1,nsrf
     DO k=1,klev
-       tke(:,k,isrf)= tke(:,k,isrf)+tendu(:,k)+tendv(:,k)+tendt(:,k)
-       tke(:,k,isrf)= max(tke(:,k,isrf),1.e-10)
+       DO i=1,klon
+          IF (pctsrf(i,isrf)>0.) THEN
+            tke(i,k,isrf)= tke(i,k,isrf)+tendu(i,k)+tendv(i,k)+tendt(i,k)
+            tke(i,k,isrf)= max(tke(i,k,isrf),1.e-10)
+          ENDIF
+       ENDDO
     ENDDO
  ENDDO
-
-! dtke_t(:,:)=tendt(:,:)
-! dtke_u(:,:)=tendu(:,:)
-! dtke_v(:,:)=tendv(:,:)
 
 
