Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90	(revision 13852)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ZDF/zdftke.F90	(revision 13873)
@@ -239,40 +239,10 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      zd_up(ji,jj,1) = 0._wp
-      zd_lw(ji,jj,1) = 1._wp
-
-      IF ( cpl_phioc .and. ln_phioc )  THEN
-         SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
-
-         CASE ( 0 ) ! Dirichlet BC
-
-            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
-               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
-               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
-               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
-                                                   ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)
-            END_2D
-
-         CASE ( 1 ) ! Neumann BC
-            DO_2D( 0, 0, 0, 0 )
-               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
-               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 )/e3w(ji,jj,2,Kmm)
-               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
-               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
-               zd_lw(ji,jj,2) = 0._wp
-               zdiag(ji,jj,1) = 1._wp
-            END_2D
-         
-         END SELECT
-
-      ELSE  ! TKE Dirichlet boundary condition (without wave coupling)
-
-         DO_2D( 0, 0, 0, 0 )       ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-            zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) ! choose to keep coherence with former estimation of
-                                               ! zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)
-         END_2D
-
-      ENDIF
+      DO_2D( 0, 0, 0, 0 )
+         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+         zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
+         zd_lw(ji,jj,1) = 1._wp  
+         zd_up(ji,jj,1) = 0._wp
+      END_2D
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -425,14 +395,32 @@
       END_3D
       !
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                            ! choose to keep coherence with previous estimation of
-      !                            !  Surface boundary condition on tke
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ![EC] should it be removed?
-      IF ( ln_isfcav ) THEN
-         DO_2D( 1, 1, 1, 1 )       ! en(mikt(ji,jj))   = rn_emin
-            en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1)
-         END_2D
-      END IF
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Surface boundary condition on tke if
+      !                     !  coupleing with waves
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      IF ( cpl_phioc .and. ln_phioc )  THEN
+         SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves 
+
+         CASE ( 0 ) ! Dirichlet BC
+            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
+               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
+            END_2D
+
+         CASE ( 1 ) ! Neumann BC
+            DO_2D( 0, 0, 0, 0 )
+               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
+               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
+               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
+               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
+               zdiag(ji,jj,1) = 1._wp
+               zd_lw(ji,jj,2) = 0._wp
+            END_2D
+
+         END SELECT
+
+      ENDIF
       !
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
