Changeset 8055

Timestamp:

2017-05-20T13:49:38+02:00 (8 years ago)

Author:

gm

Message:

wrk_OMP: 2nd step: add OMP processor distribution in ZDF package

Location:

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM

Files:

• NEMO/OPA_SRC/ZDF/zdfddm.F90 (modified) (9 diffs)
• NEMO/OPA_SRC/ZDF/zdfevd.F90 (modified) (8 diffs)
• NEMO/OPA_SRC/ZDF/zdfgls.F90 (modified) (59 diffs)
• NEMO/OPA_SRC/ZDF/zdfiwm.F90 (modified) (13 diffs)
• NEMO/OPA_SRC/ZDF/zdfphy.F90 (modified) (7 diffs)
• NEMO/OPA_SRC/ZDF/zdfric.F90 (modified) (6 diffs)
• NEMO/OPA_SRC/ZDF/zdfsh2.F90 (modified) (5 diffs)
• NEMO/OPA_SRC/ZDF/zdfswm.F90 (modified) (4 diffs)
• NEMO/OPA_SRC/ZDF/zdftke.F90 (modified) (33 diffs)
• NEMO/OPA_SRC/domain_substitute.h90 (added)
• NEMO/OPA_SRC/step.F90 (modified) (4 diffs)
• SETTE/sette.sh (modified) (5 diffs)

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

@@ -31 +31 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "domain_substitute.h90"   
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO Consortium (2014)
@@ -39 +39 @@
 CONTAINS
 
-   SUBROUTINE zdf_ddm( kt )
+   SUBROUTINE zdf_ddm( ARG_2D, kt, p_avm, p_avt, p_avs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_ddm  ***
@@ -69 +69 @@
       !! References :   Merryfield et al., JPO, 29, 1124-1142, 1999.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
+      INTEGER, INTENT(in   ) ::   ARG_2D   ! inner domain start-end i-indices
+      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step indexocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm   !  Kz on momentum    (w-points)
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avt   !  Kz on temperature (w-points)
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avs   !  Kz on salinity    (w-points)
       !
       INTEGER  ::   ji, jj , jk     ! dummy loop indices
@@ -77 +81 @@
       REAL(wp) ::   zavft, zavfs    !   -      -
       REAL(wp) ::   zavdt, zavds    !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
+      REAL(wp), DIMENSION(WRK_2D) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
       !!----------------------------------------------------------------------
       !
@@ -91 +95 @@
 !!gm                            and many acces in memory
          
-         DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s])
-            DO ji = 1, jpi
+         DO jj = k_Jstr, k_Jend        !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
+            DO ji = k_Jstr, k_Iend
                zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
 !!gm please, use e3w_n below 
@@ -109 +113 @@
          END DO
 
-         DO jj = 1, jpj                                     ! indicators:
-            DO ji = 1, jpi
+         DO jj = k_Jstr, k_Jend          !==  indicators  ==!
+            DO ji = k_Jstr, k_Iend
                ! stability indicator: msks=1 if rn2>0; 0 elsewhere
                IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
@@ -135 +139 @@
          END DO
          ! mask zmsk in order to have avt and avs masked
-         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk)
+         zmsks(WRK_2D) = zmsks(WRK_2D) * wmask(WRK_2D,jk)
 
 
@@ -141 +145 @@
          ! ------------------
          ! Constant eddy coefficient: reset to the background value
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Jstr, k_Iend
                zinr = 1._wp / zrau(ji,jj)
                ! salt fingering
@@ -154 +158 @@
                   &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
                ! add to the eddy viscosity coef. previously computed
-               avs(ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
-               avt(ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
-               avm(ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
+               p_avs(ji,jj,jk) = p_avt(ji,jj,jk) + zavfs + zavds
+               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zavft + zavdt
+               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )
             END DO
          END DO

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90

@@ -31 +31 @@
    PUBLIC   zdf_evd    ! called by step.F90
 
+   !! * Substitutions
+#  include "domain_substitute.h90"   
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -38 +40 @@
 CONTAINS
 
-   SUBROUTINE zdf_evd( kt )
+   SUBROUTINE zdf_evd( ARG_2D, kt, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_evd  ***
@@ -55 +57 @@
       !! ** Action  :   avt, avm   enhanced where static instability occurs
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
+      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time-step indexocean time step
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_evd, zavm_evd
+!!OMP issue to be solved with outputs...
+!      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_evd, zavm_evd
+!!OMP end
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_evd')
       !
+!!OMP  this is to be suppressed
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
@@ -69 +76 @@
          IF(lwp) WRITE(numout,*)
       ENDIF
+!!OMP end
       !
       !
-      zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application
+!!OMP issue with outputs
+!      zavt_evd(:,:,:) = p_avt(:,:,:)           ! set avt prior to evd application
+!!OMP end
       !
       SELECT CASE ( nn_evdm )
@@ -77 +87 @@
       CASE ( 1 )           !==  enhance tracer & momentum Kz  ==!   (if rn2<-1.e-12)
          !
-         zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application
+!!OMP issue with outputs
+!         zavm_evd(:,:,:) = p_avm(:,:,:)           ! set avm prior to evd application
+!!OMP end
          !
-!! change last digits results
-!         WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 ) THEN
-!            avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
-!            avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
+!! change last digits results ????
+!         WHERE( MAX( rn2(WRK_2D,2:jpkm1), rn2b(WRK_2D,2:jpkm1) )  <= -1.e-12 ) THEN
+!            p_avt(WRK_2D,2:jpkm1) = rn_evd * wmask(WRK_2D,2:jpkm1)
+!            p_avm(WRK_2D,2:jpkm1) = rn_evd * wmask(WRK_2D,2:jpkm1)
 !         END WHERE
          !
          DO jk = 1, jpkm1 
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Jstr, k_Iend
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
-                     avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
-                     avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
+                     p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
+                     p_avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
                   ENDIF
                END DO
@@ -96 +108 @@
          END DO 
          !
-         zavm_evd(:,:,:) = avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd
-         CALL iom_put( "avm_evd", zavm_evd )              ! output this change
+!!OMP issue with outputs
+!         zavm_evd(:,:,:) = p_avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd
+!         CALL iom_put( "avm_evd", zavm_evd )              ! output this change
+!!gm end
          !
       CASE DEFAULT         !==  enhance tracer Kz  ==!   (if rn2<-1.e-12) 
 !! change last digits results
-!         WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 )
-!            avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
+!         WHERE( MAX( rn2(WRK_2D,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 )
+!            avt(WRK_2D,2:jpkm1) = rn_evd * wmask(WRK_2D,2:jpkm1)
 !         END WHERE
 
@@ -109 +123 @@
                DO ji = 2, jpim1
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
-                     avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
+                     p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
                END DO
             END DO
@@ -116 +130 @@
       END SELECT 
       !
-      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd
-      CALL iom_put( "avt_evd", zavt_evd )              ! output this change
-      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd )
+!!OMP issue with outputs
+!      zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd
+!      CALL iom_put( "avt_evd", zavt_evd )              ! output this change
+!      IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd )
+!!OMP end
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_evd')

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

@@ -19 +19 @@
    USE domvvl         ! ocean space and time domain : variable volume layer
    USE zdf_oce        ! ocean vertical physics
-   USE zdfsh2         ! vertical physics: shear production term of TKE
    USE zdfbfr         ! bottom friction (only for rn_bfrz0)
    USE sbc_oce        ! surface boundary condition: ocean
@@ -103 +102 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "domain_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -123 +122 @@
 
 
-   SUBROUTINE zdf_gls( kt )
+   SUBROUTINE zdf_gls( ARG_2D, kt, psh2, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_gls  ***
@@ -130 +129 @@
       !!              coefficients using the GLS turbulent closure scheme.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt ! ocean time step
+      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step
+      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   psh2           ! shear production term
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
+      !
       INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
       REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars
@@ -137 +140 @@
       REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
       REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hmxl_b      ! mixing length at time before
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   psi         ! psi at time now
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a    ! element of the first  matrix diagonal
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zstt, zstm  ! stability function on tracer and momentum
+      REAL(wp), DIMENSION(WRK_2D) ::   zdep
+      REAL(wp), DIMENSION(WRK_2D) ::   zkar
+      REAL(wp), DIMENSION(WRK_2D) ::   zflxs       ! Turbulence fluxed induced by internal waves 
+      REAL(wp), DIMENSION(WRK_2D) ::   zhsro       ! Surface roughness (surface waves)
+      REAL(wp), DIMENSION(WRK_3D) ::   eb          ! tke at time before
+      REAL(wp), DIMENSION(WRK_3D) ::   hmxl_b      ! mixing length at time before
+      REAL(wp), DIMENSION(WRK_3D) ::   eps         ! dissipation rate
+      REAL(wp), DIMENSION(WRK_3D) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
+      REAL(wp), DIMENSION(WRK_3D) ::   psi         ! psi at time now
+      REAL(wp), DIMENSION(WRK_3D) ::   z_elem_a    ! element of the first  matrix diagonal
+      REAL(wp), DIMENSION(WRK_3D) ::   z_elem_b    ! element of the second matrix diagonal
+      REAL(wp), DIMENSION(WRK_3D) ::   z_elem_c    ! element of the third  matrix diagonal
+      REAL(wp), DIMENSION(WRK_3D) ::   zstt, zstm  ! stability function on tracer and momentum
       !!--------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('zdf_gls')
       !
-      CALL wrk_alloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
-      CALL wrk_alloc( jpi,jpj,jpk,   eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
-      CALL wrk_alloc( jpi,jpj,jpk,   zstt, zstm )
-
       ! Preliminary computing
 
       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
-      ! restore before values computed GLS alone
-      avt(:,:,:) = avt_k(:,:,:)
-      avm(:,:,:) = avm_k(:,:,:)
 
       ! Compute surface and bottom friction at T-points
-      DO jj = 2, jpjm1          
-         DO ji = fs_2, fs_jpim1   ! vector opt.         
+      DO jj = k_Jstr, k_Jend 
+         DO ji = k_Istr, k_Iend
             !
             ! surface friction
@@ -186 +181 @@
       SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==!
       CASE ( 0 )                          ! Constant roughness          
-         zhsro(:,:) = rn_hsro
+         zhsro(WRK_2D) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
-         zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro)
+         zhsro(WRK_2D) = MAX(rsbc_zs1 * ustars2(WRK_2D), rn_hsro)
       CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
 !!gm         zcof = 2._wp * 0.6_wp / 28._wp
-!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustars2(:,:),rsmall) )  )      ! Wave age (eq. 10)
-         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
-         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
+!!gm         zdep(WRK_2D)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustars2(WRK_2D),rsmall) )  )      ! Wave age (eq. 10)
+         zdep(WRK_2D)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(WRK_2D),rsmall))))             ! Wave age (eq. 10)
+         zhsro(WRK_2D) = MAX(rsbc_zs2 * ustars2(WRK_2D) * zdep(WRK_2D)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
 !!gm  BUG missing a multiplicative coefficient....
-         zhsro(:,:) = hsw(:,:)
+         zhsro(WRK_2D) = hsw(WRK_2D)
       END SELECT
-
-      !                             !==  Compute shear and dissipation rate  ==!
-      CALL zdf_sh2( shear )
-
-
+      !
       DO jk = 2, jpkm1              !==  Compute dissipation rate  ==!
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
             END DO
@@ -212 +203 @@
 
       ! Save tke at before time step
-      eb    (:,:,:) = en    (:,:,:)
-      hmxl_b(:,:,:) = hmxl_n(:,:,:)
+      eb    (WRK_2D,:) = en    (WRK_2D,:)
+      hmxl_b(WRK_2D,:) = hmxl_n(WRK_2D,:)
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1 
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
                   zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
@@ -242 +233 @@
 
       DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               !
-               buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
+               !
+               buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction
                !
                diss = eps(ji,jj,jk)                         ! dissipation
                !
-               dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
-               !
-               zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk)          ! production term
+               dir = 0.5_wp + SIGN( 0.5_wp, psh2(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               !
+               zesh2 = dir*(psh2(ji,jj,jk)+buoy)+(1._wp-dir)*psh2(ji,jj,jk)          ! production term
                zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
                !
@@ -269 +260 @@
                zcof = rfact_tke * tmask(ji,jj,jk)
                !                                               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk  ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
                !                                               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
                !                                               ! diagonal
                z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
@@ -280 +271 @@
       END DO
       !
-      z_elem_b(:,:,jpk) = 1._wp
+      z_elem_b(WRK_2D,jpk) = 1._wp
       !
       ! Set surface condition on zwall_psi (1 at the bottom)
-      zwall_psi(:,:, 1 ) = zwall_psi(:,:,2)
-      zwall_psi(:,:,jpk) = 1._wp
+      zwall_psi(WRK_2D, 1 ) = zwall_psi(WRK_2D,2)
+      zwall_psi(WRK_2D,jpk) = 1._wp
       !
       ! Surface boundary condition on tke
@@ -293 +284 @@
       CASE ( 0 )             ! Dirichlet case
       ! First level
-      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
-      z_elem_a(:,:,1) = en(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      en(WRK_2D,1) = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(WRK_2D,1) = MAX(en(WRK_2D,1), rn_emin) 
+      z_elem_a(WRK_2D,1) = en(WRK_2D,1)
+      z_elem_c(WRK_2D,1) = 0._wp
+      z_elem_b(WRK_2D,1) = 1._wp
       ! 
       ! One level below
-      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) &
-         &               / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
-      en(:,:,2) = MAX(en(:,:,2), rn_emin )
-      z_elem_a(:,:,2) = 0._wp 
-      z_elem_c(:,:,2) = 0._wp
-      z_elem_b(:,:,2) = 1._wp
+      en(WRK_2D,2) = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1 * ((zhsro(WRK_2D)+gdepw_n(WRK_2D,2)) &
+         &               / zhsro(WRK_2D) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
+      en(WRK_2D,2) = MAX(en(WRK_2D,2), rn_emin )
+      z_elem_a(WRK_2D,2) = 0._wp 
+      z_elem_c(WRK_2D,2) = 0._wp
+      z_elem_b(WRK_2D,2) = 1._wp
       !
       !
@@ -311 +302 @@
       !
       ! Dirichlet conditions at k=1
-      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
-      z_elem_a(:,:,1) = en(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      en(WRK_2D,1)       = rc02r * ustars2(WRK_2D) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(WRK_2D,1)       = MAX(en(WRK_2D,1), rn_emin)      
+      z_elem_a(WRK_2D,1) = en(WRK_2D,1)
+      z_elem_c(WRK_2D,1) = 0._wp
+      z_elem_b(WRK_2D,1) = 1._wp
       !
       ! at k=2, set de/dz=Fw
       !cbr
-      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-      z_elem_a(:,:,2) = 0._wp
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
-      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
-          &                       * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf)
-
-      en(:,:,2) = en(:,:,2) + zflxs(:,:)/e3w_n(:,:,2)
+      z_elem_b(WRK_2D,2) = z_elem_b(WRK_2D,2) +  z_elem_a(WRK_2D,2) ! Remove z_elem_a from z_elem_b
+      z_elem_a(WRK_2D,2) = 0._wp
+      zkar(WRK_2D)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(WRK_2D,1)/zhsro(WRK_2D)) ))
+      zflxs(WRK_2D)      = rsbc_tke2 * ustars2(WRK_2D)**1.5_wp * zkar(WRK_2D) &
+          &                       * ((zhsro(WRK_2D)+gdept_n(WRK_2D,1)) / zhsro(WRK_2D) )**(1.5_wp*ra_sf)
+
+      en(WRK_2D,2) = en(WRK_2D,2) + zflxs(WRK_2D)/e3w_n(WRK_2D,2)
       !
       !
@@ -338 +329 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -359 +350 @@
       CASE ( 1 )             ! Neumman boundary condition
          !                      
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -382 +373 @@
       !
       DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
             END DO
@@ -389 +380 @@
       END DO
       DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
             END DO
@@ -396 +387 @@
       END DO
       DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
             END DO
@@ -403 +394 @@
       END DO
       !                                            ! set the minimum value of tke 
-      en(:,:,:) = MAX( en(:,:,:), rn_emin )
+      en(WRK_3D) = MAX( en(WRK_3D), rn_emin )
 
       !!----------------------------------------!!
@@ -415 +406 @@
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
                END DO
@@ -424 +415 @@
       CASE( 1 )               ! k-eps
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   psi(ji,jj,jk)  = eps(ji,jj,jk)
                END DO
@@ -433 +424 @@
       CASE( 2 )               ! k-w
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
                END DO
@@ -442 +433 @@
       CASE( 3 )               ! generic
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
                END DO
@@ -459 +450 @@
 
       DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                !
                ! psi / k
@@ -471 +462 @@
                !
                ! shear prod. - stratif. destruction
-               prod = rpsi1 * zratio * shear(ji,jj,jk)
+               prod = rpsi1 * zratio * psh2(ji,jj,jk)
                !
                ! stratif. destruction
-               buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )
+               buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
                !
                ! shear prod. - stratif. destruction
@@ -487 +478 @@
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
                !                                               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk  ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
                !                                               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
                !                                               ! diagonal
                z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
@@ -498 +489 @@
       END DO
       !
-      z_elem_b(:,:,jpk) = 1._wp
+      z_elem_b(WRK_2D,jpk) = 1._wp
 
       ! Surface boundary condition on psi
@@ -508 +499 @@
       !
       ! Surface value
-      zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      zdep(WRK_2D)       = zhsro(WRK_2D) * rl_sf ! Cosmetic
+      psi (WRK_2D,1)     = rc0**rpp * en(WRK_2D,1)**rmm * zdep(WRK_2D)**rnn * tmask(WRK_2D,1)
+      z_elem_a(WRK_2D,1) = psi(WRK_2D,1)
+      z_elem_c(WRK_2D,1) = 0._wp
+      z_elem_b(WRK_2D,1) = 1._wp
       !
       ! One level below
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) )))
-      zdep(:,:)       = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:)
-      psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,2) = 0._wp
-      z_elem_c(:,:,2) = 0._wp
-      z_elem_b(:,:,2) = 1._wp
+      zkar(WRK_2D)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(WRK_2D,2)/zhsro(WRK_2D) )))
+      zdep(WRK_2D)       = (zhsro(WRK_2D) + gdepw_n(WRK_2D,2)) * zkar(WRK_2D)
+      psi (WRK_2D,2)     = rc0**rpp * en(WRK_2D,2)**rmm * zdep(WRK_2D)**rnn * tmask(WRK_2D,1)
+      z_elem_a(WRK_2D,2) = 0._wp
+      z_elem_c(WRK_2D,2) = 0._wp
+      z_elem_b(WRK_2D,2) = 1._wp
       ! 
       !
@@ -526 +517 @@
       !
       ! Surface value: Dirichlet
-      zdep(:,:)       = zhsro(:,:) * rl_sf
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      zdep(WRK_2D)       = zhsro(WRK_2D) * rl_sf
+      psi (WRK_2D,1)     = rc0**rpp * en(WRK_2D,1)**rmm * zdep(WRK_2D)**rnn * tmask(WRK_2D,1)
+      z_elem_a(WRK_2D,1) = psi(WRK_2D,1)
+      z_elem_c(WRK_2D,1) = 0._wp
+      z_elem_b(WRK_2D,1) = 1._wp
       !
       ! Neumann condition at k=2
-      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-      z_elem_a(:,:,2) = 0._wp
+      z_elem_b(WRK_2D,2) = z_elem_b(WRK_2D,2) +  z_elem_a(WRK_2D,2) ! Remove z_elem_a from z_elem_b
+      z_elem_a(WRK_2D,2) = 0._wp
       !
       ! Set psi vertical flux at the surface:
-      zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
-      zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
-      zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
-      zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
-             & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.)
-      zflxs(:,:) = zdep(:,:) * zflxs(:,:)
-      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
+      zkar(WRK_2D) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(WRK_2D,1)/zhsro(WRK_2D) )) ! Lengh scale slope
+      zdep(WRK_2D) = ((zhsro(WRK_2D) + gdept_n(WRK_2D,1)) / zhsro(WRK_2D))**(rmm*ra_sf)
+      zflxs(WRK_2D) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(WRK_2D))*(1._wp + rsbc_tke1*zdep(WRK_2D))**(2._wp*rmm/3._wp-1_wp)
+      zdep(WRK_2D) =  rsbc_psi1 * (zwall_psi(WRK_2D,1)*p_avm(WRK_2D,1)+zwall_psi(WRK_2D,2)*p_avm(WRK_2D,2)) * &
+             & ustars2(WRK_2D)**rmm * zkar(WRK_2D)**rnn * (zhsro(WRK_2D) + gdept_n(WRK_2D,1))**(rnn-1.)
+      zflxs(WRK_2D) = zdep(WRK_2D) * zflxs(WRK_2D)
+      psi(WRK_2D,2) = psi(WRK_2D,2) + zflxs(WRK_2D) / e3w_n(WRK_2D,2)
       !   
       !
@@ -557 +548 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -578 +569 @@
       CASE ( 1 )             ! Neumman boundary condition
          !                      
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -597 +588 @@
                ! Set psi vertical flux at the bottom:
                zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1)
-               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   &
+               zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
                   &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
                psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1)
@@ -609 +600 @@
       !
       DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
             END DO
@@ -616 +607 @@
       END DO
       DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
             END DO
@@ -623 +614 @@
       END DO
       DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
             END DO
@@ -637 +628 @@
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin)
                END DO
@@ -646 +637 @@
       CASE( 1 )               ! k-eps
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   eps(ji,jj,jk) = psi(ji,jj,jk)
                END DO
@@ -655 +646 @@
       CASE( 2 )               ! k-w
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
                END DO
@@ -667 +658 @@
          zex2  = -1._wp / rnn
          DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
                END DO
@@ -679 +670 @@
       ! --------------------------------------------------
       DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                ! limitation
                eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
@@ -699 +690 @@
       CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   ! zcof =  l²/q²
                   zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
@@ -720 +711 @@
       CASE ( 2, 3 )               ! Canuto stability functions
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend 
+               DO ji = k_Istr, k_Iend
                   ! zcof =  l²/q²
                   zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
@@ -730 +721 @@
                   gh = gh * rf6
                   ! Gm =  M²l²/q² Shear number
-                  shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall )
+                  shr = psh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
                   gm = MAX( shr * zcof , 1.e-10 )
                   gm = gm * rf6
@@ -751 +742 @@
       ! Lines below are useless if GOTM style Dirichlet conditions are used
 
-      zstm(:,:,1) = zstm(:,:,2)
+      zstm(WRK_2D,1) = zstm(WRK_2D,2)
 
 !!gm should be done for ISF (top boundary cond.)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = k_Jstr, k_Jend 
+         DO ji = k_Istr, k_Iend
             zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
          END DO
@@ -763 +754 @@
       ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
       DO jk = 1, jpk
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
                zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
                zav           = zsqen * zstt(ji,jj,jk)
-               avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
+               p_avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
                zav           = zsqen * zstm(ji,jj,jk)
-               avm(ji,jj,jk) = MAX( zav, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
+               p_avm(ji,jj,jk) = MAX( zav, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
             END DO
          END DO
       END DO
-      avt(:,:,1)  = 0._wp
-!!gm I'm sure this is needed to compute the shear term at next time-step
-      CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
+      p_avt(WRK_2D,1)  = 0._wp
       !
       IF(ln_ctl) THEN
@@ -781 +770 @@
          CALL prt_ctl( tab3d_1=avm, clinfo1=' gls  - m: ', ovlap=1, kdim=jpk )
       ENDIF
-      !
-      avt_k(:,:,:) = avt(:,:,:)      !==  store avt, avm values computed by GLS only  ==!
-      avm_k(:,:,:) = avm(:,:,:)
-      !
-      IF( lrst_oce )   CALL gls_rst( kt, 'WRITE' )     ! write the GLS info in the restart file
-      !
-      CALL wrk_dealloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
-      CALL wrk_dealloc( jpi,jpj,jpk,   eb, hmxl_b, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
-      CALL wrk_dealloc( jpi,jpj,jpk,   zstt, zstm )
       !
       IF( nn_timing == 1 )   CALL timing_stop('zdf_gls')
@@ -837 +817 @@
       IF(lwp) THEN                     !* Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
+         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
@@ -855 +835 @@
          WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
          WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0
+         WRITE(numout,*)
       ENDIF
 
@@ -861 +842 @@
 
       !                                !* Check of some namelist values
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )
-      IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )
-      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )
-      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
+      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
+      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
+      IF( nn_z0_met  < 0   .OR. nn_z0_met    > 3 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )
+      IF( nn_z0_met == 3  .AND. .NOT.ln_wave     )   CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )
+      IF( nn_stab_func < 0 .OR. nn_stab_func > 3 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )
+      IF( nn_clos      < 0 .OR. nn_clos      > 3 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
 
       SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
@@ -872 +853 @@
       CASE( 0 )                              ! k-kl  (Mellor-Yamada)
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
+         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
+         IF(lwp) WRITE(numout,*)
          rpp     = 0._wp
          rmm     = 1._wp
@@ -890 +872 @@
       CASE( 1 )                              ! k-eps
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
+         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen'
+         IF(lwp) WRITE(numout,*)
          rpp     =  3._wp
          rmm     =  1.5_wp
@@ -908 +891 @@
       CASE( 2 )                              ! k-omega
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
+         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen'
+         IF(lwp) WRITE(numout,*)
          rpp     = -1._wp
          rmm     =  0.5_wp
@@ -926 +910 @@
       CASE( 3 )                              ! generic
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
+         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen'
+         IF(lwp) WRITE(numout,*)
          rpp     = 2._wp
          rmm     = 1._wp
@@ -949 +934 @@
       CASE ( 0 )                             ! Galperin stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin'
          rc2     =  0._wp
          rc3     =  0._wp
@@ -961 +946 @@
       CASE ( 1 )                             ! Kantha-Clayson stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson'
          rc2     =  0.7_wp
          rc3     =  0.2_wp
@@ -973 +958 @@
       CASE ( 2 )                             ! Canuto A stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A'
          rs0 = 1.5_wp * rl1 * rl5*rl5
          rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
@@ -997 +982 @@
       CASE ( 3 )                             ! Canuto B stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B'
          rs0 = 1.5_wp * rm1 * rm5*rm5
          rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
@@ -1052 +1037 @@
       IF(lwp) THEN                     !* Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'Limit values'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) 'Parameter  m = ',rmm
-         WRITE(numout,*) 'Parameter  n = ',rnn
-         WRITE(numout,*) 'Parameter  p = ',rpp
-         WRITE(numout,*) 'rpsi1   = ',rpsi1
-         WRITE(numout,*) 'rpsi2   = ',rpsi2
-         WRITE(numout,*) 'rpsi3m  = ',rpsi3m
-         WRITE(numout,*) 'rpsi3p  = ',rpsi3p
-         WRITE(numout,*) 'rsc_tke = ',rsc_tke
-         WRITE(numout,*) 'rsc_psi = ',rsc_psi
-         WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0
-         WRITE(numout,*) 'rc0     = ',rc0
+         WRITE(numout,*) '   Limit values :'
+         WRITE(numout,*) '      Parameter  m = ',rmm
+         WRITE(numout,*) '      Parameter  n = ',rnn
+         WRITE(numout,*) '      Parameter  p = ',rpp
+         WRITE(numout,*) '      rpsi1   = ',rpsi1
+         WRITE(numout,*) '      rpsi2   = ',rpsi2
+         WRITE(numout,*) '      rpsi3m  = ',rpsi3m
+         WRITE(numout,*) '      rpsi3p  = ',rpsi3p
+         WRITE(numout,*) '      rsc_tke = ',rsc_tke
+         WRITE(numout,*) '      rsc_psi = ',rsc_psi
+         WRITE(numout,*) '      rsc_psi0 = ',rsc_psi0
+         WRITE(numout,*) '      rc0     = ',rc0
          WRITE(numout,*)
-         WRITE(numout,*) 'Shear free turbulence parameters:'
-         WRITE(numout,*) 'rcm_sf  = ',rcm_sf
-         WRITE(numout,*) 'ra_sf   = ',ra_sf
-         WRITE(numout,*) 'rl_sf   = ',rl_sf
-         WRITE(numout,*)
+         WRITE(numout,*) '   Shear free turbulence parameters:'
+         WRITE(numout,*) '      rcm_sf  = ',rcm_sf
+         WRITE(numout,*) '      ra_sf   = ',ra_sf
+         WRITE(numout,*) '      rl_sf   = ',rl_sf
       ENDIF
 
@@ -1090 +1073 @@
       !
       !                                !* Wall proximity function
-      zwall (:,:,:) = 1._wp * tmask(:,:,:)
+      zwall(:,:,:) = 1._wp * tmask(:,:,:)
 
       !                                !* read or initialize all required files  
@@ -1133 +1116 @@
                CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n )
             ELSE                        
-               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxl_n computed by iterative loop'
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values'
                en    (:,:,:) = rn_emin
                hmxl_n(:,:,:) = 0.05_wp
-               avt_k (:,:,:) = avt(:,:,:)
-               avm_k (:,:,:) = avm(:,:,:)
-               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
+               ! avt_k, avm_k already set to the background value in zdf_phy_init
             ENDIF
          ELSE                                   !* Start from rest
-            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxl_n by background values'
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values'
             en    (:,:,:) = rn_emin
             hmxl_n(:,:,:) = 0.05_wp
-            DO jk = 1, jpk
-               avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk)
-               avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk)
-            END DO
+            ! avt_k, avm_k already set to the background value in zdf_phy_init
          ENDIF
          !

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90

@@ -35 +35 @@
    PUBLIC   zdf_iwm         ! called in step module 
    PUBLIC   zdf_iwm_init    ! called in nemogcm module 
-   PUBLIC   zdf_iwm_alloc   ! called in nemogcm module
 
    !                       !!* Namelist  namzdf_iwm : internal wave-driven mixing *
@@ -56 +55 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "domain_substitute.h90"   
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2016)
@@ -78 +77 @@
 
 
-   SUBROUTINE zdf_iwm( kt )
+   SUBROUTINE zdf_iwm( ARG_2D, kt, p_avm, p_avt, p_avs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_iwm  ***
@@ -127 +126 @@
       !! References :  de Lavergne et al. 2015, JPO; 2016, in prep.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step 
+      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points)
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points)
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   ztpc         ! scalar workspace
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid)
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid)
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter
+      REAL(wp) ::   zztemp       ! scalar workspace
+      REAL(wp), DIMENSION(WRK_2D) ::  zfact     ! Used for vertical structure
+      REAL(wp), DIMENSION(WRK_2D) ::  zhdep     ! Ocean depth
+      REAL(wp), DIMENSION(WRK_3D) ::  zwkb      ! WKB-stretched height above bottom
+      REAL(wp), DIMENSION(WRK_3D) ::  zweight   ! Weight for high mode vertical distribution
+      REAL(wp), DIMENSION(WRK_3D) ::  znu_t     ! Molecular kinematic viscosity (T grid)
+      REAL(wp), DIMENSION(WRK_3D) ::  znu_w     ! Molecular kinematic viscosity (W grid)
+      REAL(wp), DIMENSION(WRK_3D) ::  zReb      ! Turbulence intensity parameter
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('zdf_iwm')
       !
-      CALL wrk_alloc( jpi,jpj,       zfact, zhdep )
-      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
-
       !                          ! ----------------------------- !
       !                          !  Internal wave-driven mixing  !  (compute zav_wave)
@@ -151 +150 @@
       !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth,
       !                                                 using an exponential decay from the seafloor.
-      DO jj = 1, jpj                ! part independent of the level
-         DO ji = 1, jpi
+      DO jj = k_Jstr, k_Jend
+         DO ji = k_Istr, k_Iend     ! part independent of the level
             zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
             zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
@@ -158 +157 @@
          END DO
       END DO
-
       DO jk = 2, jpkm1              ! complete with the level-dependent part
-         emix_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      &
-            &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   &
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
+               emix_iwm(ji,jj,jk) = zfact(ji,jj) * wmask(ji,jj,jk)                              & 
+                  &   * (  EXP( ( gde3w_n(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )      &
+                  &      - EXP( ( gde3w_n(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )  )   &
+                  &   / ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )
+
 !!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point
-            &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
+!!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all
+
+            END DO
+         END DO
       END DO
 
       !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying 
       !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
-      !
+      !                                          ! (NB: N2 is masked, so no use of wmask here)
       SELECT CASE ( nn_zpyc )
       !
       CASE ( 1 )               ! Dissipation scales as N (recommended)
          !
-         zfact(:,:) = 0._wp
-         DO jk = 2, jpkm1              ! part independent of the level
-            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
-         END DO
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
-            END DO
-         END DO
+         zfact(WRK_2D) = 0._wp         ! part independent of the level
+         DO jk = 2, jpkm1
+            zfact(WRK_2D) = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  )
+         END DO
+         !
+         WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) )
+!         DO jj = k_Jstr, k_Jend
+!            DO ji = k_Istr, k_Iend
+!               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+!            END DO
+!         END DO
+         !                             ! complete with the level-dependent part
+         DO jk = 2, jpkm1
+            emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zfact(WRK_2D) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  )
+         END DO
+         !
+      CASE ( 2 )               ! Dissipation scales as N^2
+         !
+         zfact(WRK_2D) = 0._wp
+         DO jk = 2, jpkm1              ! part independent of the level 
+            zfact(WRK_2D) = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * MAX( 0._wp, rn2(WRK_2D,jk) )
+         END DO
+         !
+         WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) )
+!         DO jj = k_Jstr, k_Jend
+!            DO ji = k_Istr, k_Iend
+!               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+!            END DO
+!         END DO
          !
          DO jk = 2, jpkm1              ! complete with the level-dependent part
-            emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
-         END DO
-         !
-      CASE ( 2 )               ! Dissipation scales as N^2
-         !
-         zfact(:,:) = 0._wp
-         DO jk = 2, jpkm1              ! part independent of the level
-            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
-         END DO
-         !
-         DO jj= 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
-            END DO
-         END DO
-         !
-         DO jk = 2, jpkm1              ! complete with the level-dependent part
-            emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+            emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zfact(WRK_2D) * MAX( 0._wp, rn2(WRK_2D,jk) )
          END DO
          !
@@ -210 +218 @@
       !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
       !
-      zwkb(:,:,:) = 0._wp
-      zfact(:,:) = 0._wp
+      zwkb (WRK_3D) = 0._wp
+      zfact(WRK_2D) = 0._wp
       DO jk = 2, jpkm1
-         zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
-         zwkb(:,:,jk) = zfact(:,:)
+!!gm  this is potentially dangerous 
+!         zfact(WRK_2D)    = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  )
+!         zwkb (WRK_2D,jk) = zfact(WRK_2D)
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
+               zztemp = zfact(ji,jj) + e3w_n(ji,jj,jk) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  )
+               zfact(ji,jj)    = zztemp
+               zwkb (ji,jj,jk) = zztemp
+            END DO
+         END DO         
       END DO
       !
       DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
-                  &                                     * tmask(ji,jj,jk) / zfact(ji,jj)
+                  &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
             END DO
          END DO
       END DO
-      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1)
-      !
-      zweight(:,:,:) = 0._wp
+      zwkb(WRK_2D,1) = zhdep(WRK_2D) * wmask(WRK_2D,1)
+      !
+      zweight(WRK_3D) = 0._wp
       DO jk = 2, jpkm1
-         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk)                    &
-            &   * (  EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) )  )
-      END DO
-      !
-      zfact(:,:) = 0._wp
+         zweight(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * hbot_iwm(:,:)     &
+            &               * (  EXP( -zwkb(WRK_2D,jk  ) / hbot_iwm(WRK_2D) )  &
+            &                  - EXP( -zwkb(WRK_2D,jk-1) / hbot_iwm(WRK_2D) )  )
+      END DO
+      !
+      zfact(WRK_2D) = 0._wp
       DO jk = 2, jpkm1              ! part independent of the level
-         zfact(:,:) = zfact(:,:) + zweight(:,:,jk)
-      END DO
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
-         END DO
-      END DO
+         zfact(WRK_2D) = zfact(WRK_2D) + zweight(WRK_2D,jk)
+      END DO
+      !
+      WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = ebot_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) )
+!      DO jj = k_Jstr, k_Jend
+!         DO ji = k_Istr, k_Iend
+!            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+!         END DO
+!      END DO
       !
       DO jk = 2, jpkm1              ! complete with the level-dependent part
-         emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
-            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
-      END DO
-      !
+         emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zweight(WRK_2D,jk) * zfact(WRK_2D) * wmask(WRK_2D,jk)   &
+            &                                / ( gde3w_n(WRK_2D,jk) - gde3w_n(WRK_2D,jk-1) )
+!!gm  use of e3t_n just above?
+      END DO
+      !
+!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s
       ! Calculate molecular kinematic viscosity
-      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  &
-         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0
+      znu_t(WRK_3D) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(WRK_3D,jp_tem)             &
+         &                        + 0.00694_wp * tsn(WRK_3D,jp_tem) * tsn(WRK_3D,jp_tem)   &
+         &                        + 0.02305_wp * tsn(WRK_3D,jp_sal)  ) * tmask(WRK_3D) * r1_rau0
       DO jk = 2, jpkm1
-         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
-      END DO
+         znu_w(WRK_2D,jk) = 0.5_wp * ( znu_t(WRK_2D,jk-1) + znu_t(WRK_2D,jk) ) * wmask(WRK_2D,jk)
+      END DO
+!!gm end
       !
       ! Calculate turbulence intensity parameter Reb
       DO jk = 2, jpkm1
-         zReb(:,:,jk) = emix_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
+         zReb(WRK_2D,jk) = emix_iwm(WRK_2D,jk) / MAX( 1.e-20_wp, znu_w(WRK_2D,jk) * rn2(WRK_2D,jk) )
       END DO
       !
       ! Define internal wave-induced diffusivity
       DO jk = 2, jpkm1
-         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
+         zav_wave(WRK_2D,jk) = znu_w(WRK_2D,jk) * zReb(WRK_2D,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
       END DO
       !
       IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
          DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
-            DO jj = 1, jpj
-               DO ji = 1, jpi
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
                   IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
                      zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
@@ -281 +303 @@
       !
       DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s
-         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk)
+         zav_wave(WRK_2D,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(WRK_2D,jk) ), 1.e-2_wp  ) * wmask(WRK_2D,jk)
       END DO
       !
       IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
-         ztpc = 0._wp
+         zztemp = 0._wp
 !!gm used of glosum 3D....
          DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   &
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  zztemp = zztemp + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   &
                      &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
          END DO
-         IF( lk_mpp )   CALL mpp_sum( ztpc )
-         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
+         IF( lk_mpp )   CALL mpp_sum( zztemp )
+         zztemp = rau0 * zztemp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
          !
          IF(lwp) THEN
@@ -303 +325 @@
             WRITE(numout,*) '~~~~~~~ '
             WRITE(numout,*)
-            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW'
+            WRITE(numout,*) '      Total power consumption by av_wave =  ', zztemp * 1.e-12_wp, 'TW'
          ENDIF
       ENDIF
@@ -313 +335 @@
       IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
          DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
-            DO jj = 1, jpj
-               DO ji = 1, jpi
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
                   zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
                       &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
@@ -323 +345 @@
          CALL iom_put( "av_ratio", zav_ratio )
          DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
-            avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)
-            avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
-            avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+            p_avs(WRK_2D,jk) = p_avs(WRK_2D,jk) + zav_wave(WRK_2D,jk) * zav_ratio(WRK_2D,jk)
+            p_avt(WRK_2D,jk) = p_avt(WRK_2D,jk) + zav_wave(WRK_2D,jk)
+            p_avm(WRK_2D,jk) = p_avm(WRK_2D,jk) + zav_wave(WRK_2D,jk)
          END DO
          !
       ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing
          DO jk = 2, jpkm1
-            avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk)
-            avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
-            avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+            p_avs(WRK_2D,jk) = p_avs(WRK_2D,jk) + zav_wave(WRK_2D,jk)
+            p_avt(WRK_2D,jk) = p_avt(WRK_2D,jk) + zav_wave(WRK_2D,jk)
+            p_avm(WRK_2D,jk) = p_avm(WRK_2D,jk) + zav_wave(WRK_2D,jk)
          END DO
       ENDIF
@@ -341 +363 @@
                                     !  vertical integral of rau0 * Kz * N^2 (pcmap_iwm), energy density (emix_iwm)
       IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
-         bflx_iwm(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
-         pcmap_iwm(:,:) = 0._wp
          DO jk = 2, jpkm1
-            pcmap_iwm(:,:) = pcmap_iwm(:,:) + e3w_n(:,:,jk) * bflx_iwm(:,:,jk) * wmask(:,:,jk)
-         END DO
-         pcmap_iwm(:,:) = rau0 * pcmap_iwm(:,:)
-         CALL iom_put( "bflx_iwm", bflx_iwm )
+            bflx_iwm(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * zav_wave(WRK_2D,jk)
+         END DO
+         pcmap_iwm(WRK_2D) = 0._wp
+         DO jk = 2, jpkm1
+            pcmap_iwm(WRK_2D) = pcmap_iwm(WRK_2D) + e3w_n(WRK_2D,jk) * bflx_iwm(WRK_2D,jk)
+         END DO
+         pcmap_iwm(WRK_2D) = rau0 * pcmap_iwm(WRK_2D)
+         CALL iom_put( "bflx_iwm" , bflx_iwm  )
          CALL iom_put( "pcmap_iwm", pcmap_iwm )
       ENDIF
@@ -353 +377 @@
       CALL iom_put( "emix_iwm", emix_iwm )
       
-      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep )
-      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
-
       IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
       !

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90

@@ -11 +11 @@
    !!   zdf_phy       : upadate at each time-step the vertical mixing coeff. 
    !!----------------------------------------------------------------------
-   USE par_oce        ! ocean parameters
+   USE oce            ! ocean dynamics and tracers variables
    USE zdf_oce        ! vertical physics: shared variables         
    USE zdfbfr         ! vertical physics: bottom friction
-   USE zdfric         ! vertical physics: Richardson vertical mixing   
+   USE zdfsh2         ! vertical physics: shear production term of TKE
+   USE zdfric         ! vertical physics: RIChardson dependent vertical mixing   
    USE zdftke         ! vertical physics: TKE vertical mixing
    USE zdfgls         ! vertical physics: GLS vertical mixing
@@ -44 +45 @@
    INTEGER, PARAMETER ::   np_TKE = 3   ! Turbulente Kinetic Eenergy closure scheme for Kz
    INTEGER, PARAMETER ::   np_GLS = 4   ! Generic Length Scale closure scheme for Kz
-      
+
+   LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC))
+
+   !! * Substitutions
+#  include "domain_substitute.h90"   
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2017)
@@ -61 +66 @@
       !!                set horizontal shape and vertical profile of background mixing coef.
       !!----------------------------------------------------------------------
+      INTEGER ::   jk            ! dummy loop indices
       INTEGER ::   ioptio, ios   ! local integers
       !!
@@ -90 +96 @@
          WRITE(numout,*) '      vertical closure scheme'
          WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst
-         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfric = ', ln_zdfric
-         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdftke = ', ln_zdftke
-         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfgls = ', ln_zdfgls
+         WRITE(numout,*) '         Richardson number dependent closure     ln_zdfric = ', ln_zdfric
+         WRITE(numout,*) '         Turbulent Kinetic Energy closure (TKE)  ln_zdftke = ', ln_zdftke
+         WRITE(numout,*) '         Generic Length Scale closure (GLS)      ln_zdfgls = ', ln_zdfgls
          WRITE(numout,*) '      convection: '
          WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
@@ -137 +143 @@
       ENDIF
       !
-      !                                   ! set 1st/last level Av to zero one for all
-      avt(:,:,1) = 0._wp   ;   avs(:,:,1) = 0._wp   ;   avm(:,:,1) = 0._wp
+      DO jk = 1, jpk                      ! set turbulent closure Kz to the background value (avt_k, avm_k)
+         avt_k(:,:,jk) = avtb_2d(:,:) * avtb(jk) * wmask (:,:,jk)
+         avm_k(:,:,jk) =                avmb(jk) * wmask (:,:,jk)
+      END DO
+!!gm  to be tested only the 1st & last levels
+!      avt  (:,:, 1 ) = 0._wp   ;   avs(:,:, 1 ) = 0._wp   ;   avm  (:,:, 1 ) = 0._wp
+!      avt  (:,:,jpk) = 0._wp   ;   avs(:,:,jpk) = 0._wp   ;   avm  (:,:,jpk) = 0._wp
+!!gm
+      avt  (:,:,:) = 0._wp   ;   avs(:,:,:) = 0._wp   ;   avm  (:,:,:) = 0._wp
 
       !                          !==  Convection  ==!
@@ -170 +183 @@
       IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' )
       ENDIF
+      !                                ! shear production term flag
+      IF( ln_zdfcst ) THEN   ;   l_zdfsh2 = .FALSE.
+      ELSE                   ;   l_zdfsh2 = .TRUE.
+      ENDIF
 
       !                          !== gravity wave-driven mixing  ==!
@@ -200 +217 @@
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indice
+!!OMP      REAL(wp), DIMENSION(WRK_3D) ::   zsh2   ! shear production
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zsh2   ! shear production
       !! ---------------------------------------------------------------------
+      !      
+!!OMP ===>>>>   Open-MP   to be defined elsewhere   
+      !
+      INTEGER ::   ARG_2D
+      !
+#if defined key_vectopt_loop
+      k_Istr = 1   ;   k_Iend = jpi
+#else
+      k_Istr = 2   ;   k_Iend = jpim1
+#endif
+      k_Jstr = 2   ;   k_Jend = jpjm1
+      !
+!!OMP <<<<===   end
+      
+      ALLOCATE( zsh2(WRK_3D) )
+      
+      
       !
       CALL zdf_bfr( kt )                        !* bottom friction (if quadratic)
-      !                       
+      !
+
+      !                          !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
+      
+      IF( l_zdfsh2 )   &                        !* shear production at w-points (energy conserving form)
+         CALL zdf_sh2( ARG_2D, ub, vb, un, vn, avm_k,   &   ! <<== in  fields
+            &                                   zsh2    )   ! ==>> out field  : shear production
+
+      !
       SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
-      CASE( np_RIC )   ;   CALL zdf_ric( kt )         ! Richardson number dependent Kz
-      CASE( np_TKE )   ;   CALL zdf_tke( kt )         ! TKE closure scheme for Kz
-      CASE( np_GLS )   ;   CALL zdf_gls( kt )         ! GLS closure scheme for Kz
+      CASE( np_RIC )   ;   CALL zdf_ric( ARG_2D, kt, gdept_n, zsh2, avm_k, avt_k )         ! Richardson number dependent Kz
+      CASE( np_TKE )   ;   CALL zdf_tke( ARG_2D, kt         , zsh2, avm_k, avt_k )         ! TKE closure scheme for Kz
+      CASE( np_GLS )   ;   CALL zdf_gls( ARG_2D, kt         , zsh2, avm_k, avt_k )         ! GLS closure scheme for Kz
       CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value)
-         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
-         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
+         ! avt_k and avm_k set one for all at initialisation phase
+!!gm         avt_k(WRK_3D) = rn_avt0 * wmask(WRK_3D)
+!!gm         avm_k(WRK_3D) = rn_avm0 * wmask(WRK_3D)
       END SELECT
+      !  
+      !                          !==  ocean Kz  ==!   (avt, avs, avm)
+      !
+      !                                         !* start from turbulent closure values
+      avt(WRK_2D,2:jpkm1) = avt_k(WRK_2D,2:jpkm1)
+      avm(WRK_2D,2:jpkm1) = avm_k(WRK_2D,2:jpkm1)
       !
       IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths
          DO jk = 2, nkrnf
-            avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk)
+            avt(WRK_2D,jk) = avt(WRK_2D,jk) + 2._wp * rn_avt_rnf * rnfmsk(WRK_2D) * wmask(WRK_2D,jk)
          END DO
       ENDIF
       !
-      IF( ln_zdfevd )   CALL zdf_evd( kt )      !* convection: enhanced vertical eddy diffusivity
+      IF( ln_zdfevd )   CALL zdf_evd( ARG_2D, kt, avm, avt )      !* convection: enhanced vertical eddy diffusivity
       !
       !                                         !* double diffusive mixing
       IF( ln_zdfddm ) THEN                            ! update avt and compute avs
-                        CALL zdf_ddm( kt )
+                        CALL zdf_ddm( ARG_2D, kt, avm, avt, avs )
       ELSE                                            ! same mixing on all tracers
-         avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1)
+         avs(WRK_3D) = avt(WRK_3D)
       ENDIF
       !
       !                                         !* wave-induced mixing 
-      IF( ln_zdfswm )   CALL zdf_swm( kt )            ! surface  wave (Qiao et al. 2004) 
-      IF( ln_zdfiwm )   CALL zdf_iwm( kt )            ! internal wave (de Lavergne et al 2017)
+      IF( ln_zdfswm )   CALL zdf_swm( ARG_2D, kt, avm, avt, avs )   ! surface  wave (Qiao et al. 2004) 
+      IF( ln_zdfiwm )   CALL zdf_iwm( ARG_2D, kt, avm, avt, avs )   ! internal wave (de Lavergne et al 2017)
 
 
       !                                         !* Lateral boundary conditions (sign unchanged)
-      CALL lbc_lnk( avm, 'W', 1. )
-      CALL lbc_lnk( avt, 'W', 1. )
-      CALL lbc_lnk( avs, 'W', 1. )
-      !
-!!gm  ===>>>   TO BE REMOVED 
-      DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
-               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
-!!gm end
-
+      CALL lbc_lnk( avm_k, 'W', 1. )
+      CALL lbc_lnk( avt_k, 'W', 1. )
+      CALL lbc_lnk( avm  , 'W', 1. )
+      CALL lbc_lnk( avt  , 'W', 1. )                  !!gm  a priori only avm_k and avm are required
+      CALL lbc_lnk( avs  , 'W', 1. )                  !!gm  To be tested
+      !
 
       CALL zdf_mxl( kt )                        !* mixed layer depth, and level
 
-!!gm  !==>> to be moved in zdftke & zdfgls    
-                                                      ! write TKE or GLS information in the restart file
-      IF( lrst_oce .AND. ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
-      IF( lrst_oce .AND. ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
-      !      
+
+      IF( lrst_oce ) THEN                       !* write TKE, GLS or RIC fields in the restart file
+         IF( ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
+         IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
+         IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' ) 
+      ENDIF
+      !
+      DEALLOCATE( zsh2 )
+      !
    END SUBROUTINE zdf_phy
 

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

@@ -16 +16 @@
 
    !!----------------------------------------------------------------------
+   !!   zdf_ric_init  : initialization, namelist read, & parameters control
    !!   zdf_ric       : update momentum and tracer Kz from the Richardson number
-   !!   zdf_ric_init  : initialization, namelist read, & parameters control
+   !!   ric_rst       : read/write RIC restart in ocean restart file
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables
    USE zdf_oce        ! vertical physics: variables
-   USE zdfsh2         ! vertical physics: shear production term of TKE
    USE phycst         ! physical constants
    USE sbc_oce,  ONLY :   taum
@@ -28 +28 @@
    !
    USE in_out_manager ! I/O manager
-   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
-   USE lib_mpp        ! MPP library
+   USE iom            ! I/O manager library
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -37 +36 @@
    PRIVATE
 
-   PUBLIC   zdf_ric         ! called by step.F90
-   PUBLIC   zdf_ric_init    ! called by opa.F90
+   PUBLIC   zdf_ric         ! called by zdfphy.F90
+   PUBLIC   ric_rst         ! called by zdfphy.F90
+   PUBLIC   zdf_ric_init    ! called by nemogcm.F90
 
    !                        !!* Namelist namzdf_ric : Richardson number dependent Kz *
@@ -52 +52 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "domain_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -59 +59 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE zdf_ric( kt )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE zdfric  ***
-      !!                    
-      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
-      !!                a function of the local richardson number.
-      !!
-      !! ** Method  :   Local richardson number dependent formulation of the 
-      !!                vertical eddy viscosity and diffusivity coefficients. 
-      !!                The eddy coefficients are given by:
-      !!                    avm = avm0 + avmb
-      !!                    avt = avm0 / (1 + rn_alp*ri)
-      !!                with ri  = N^2 / dz(u)**2
-      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
-      !!                    avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric
-      !!      Where ri is the before local Richardson number,
-      !!            rn_avmri is the maximum value reaches by avm and avt 
-      !!            avmb and avtb are the background (or minimum) values
-      !!            and rn_alp, nn_ric are adjustable parameters.
-      !!      Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
-      !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
-      !!      a numerical threshold is impose on the vertical shear (1.e-20)
-      !!      As second step compute Ekman depth from wind stress forcing
-      !!      and apply namelist provided vertical coeff within this depth.
-      !!      The Ekman depth is:
-      !!              Ustar = SQRT(Taum/rho0)
-      !!              ekd= rn_ekmfc * Ustar / f0
-      !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation
-      !!      of the above equation indicates the value is somewhat arbitrary; therefore
-      !!      we allow the freedom to increase or decrease this value, if the
-      !!      Ekman depth estimate appears too shallow or too deep, respectively.
-      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
-      !!      namelist
-      !!        N.B. the mask are required for implicit scheme, and surface
-      !!      and bottom value already set in zdfphy.F90
-      !!
-      !! ** Action  :   avm, avt  mixing coeff (inner domain values only)
-      !!
-      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
-      !!              PFJ Lermusiaux 2001.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step
-      !!
-      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      LOGICAL  ::   ll_av_weight = .TRUE.      ! local logical
-      REAL(wp) ::   zcfRi, zav, zustar   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zh_ekm   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2     ! 3D     -
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('zdf_ric')
-      !
-      !                       !==  avm and avt = F(Richardson number)  ==!
-      !
-      !                          !* Shear production at uw- and vw-points (energy conserving form)
-      CALL zdf_sh2( zsh2 )
-      !
-      DO jk = 2, jpkm1           !* Vertical eddy viscosity and diffusivity coefficients
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               !                          ! coefficient = F(richardson number) (avm-weighted Ri)
-               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2(ji,jj,jk) + 1.e-20 ) )  )
-               zav   = rn_avmri * zcfRi**nn_ric
-               !
-               !                          ! avm and avt coefficients
-               avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
-               avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-!!gm BUG <<<<====  This param can't work at low latitude 
-!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
-      !
-      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
-         !
-         DO jj = 2, jpjm1        !* Ekman depth
-            DO ji = 2, jpim1
-               zustar = SQRT( taum(ji,jj) * r1_rau0 )
-               zh_ekm(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )     ! Ekman depth
-               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zh_ekm(ji,jj) , rn_mldmax )  )   ! set allowed rang
-            END DO
-         END DO
-         DO jk = 2, jpkm1        !* minimum mixing coeff. within the Ekman layer
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF( gdept_n(ji,jj,jk) < zh_ekm(ji,jj) ) THEN
-                     avm(ji,jj,jk) = MAX(  avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
-                     avt(ji,jj,jk) = MAX(  avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('zdf_ric')
-      !
-   END SUBROUTINE zdf_ric
-
 
    SUBROUTINE zdf_ric_init
@@ -205 +105 @@
       ENDIF
       !
-      DO jk = 1, jpk                 ! Initialization of vertical eddy coef. to the background value
-         avt(:,:,jk) = avtb(jk) * tmask(:,:,jk)
-         avm(:,:,jk) = avmb(jk) * umask(:,:,jk)
+      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE zdf_ric_init
+
+
+   SUBROUTINE zdf_ric( ARG_2D, kt, pdept, p_sh2, p_avm, p_avt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE zdfric  ***
+      !!                    
+      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
+      !!                a function of the local richardson number.
+      !!
+      !! ** Method  :   Local richardson number dependent formulation of the 
+      !!                vertical eddy viscosity and diffusivity coefficients. 
+      !!                The eddy coefficients are given by:
+      !!                    avm = avm0 + avmb
+      !!                    avt = avm0 / (1 + rn_alp*ri)
+      !!                with Ri  = N^2 / dz(u)**2
+      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
+      !!                    avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric
+      !!                where ri is the before local Richardson number,
+      !!                rn_avmri is the maximum value reaches by avm and avt 
+      !!                and rn_alp, nn_ric are adjustable parameters.
+      !!                Typical values : rn_alp=5. and nn_ric=2.
+      !!
+      !!      As second step compute Ekman depth from wind stress forcing
+      !!      and apply namelist provided vertical coeff within this depth.
+      !!      The Ekman depth is:
+      !!              Ustar = SQRT(Taum/rho0)
+      !!              ekd= rn_ekmfc * Ustar / f0
+      !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation
+      !!      of the above equation indicates the value is somewhat arbitrary; therefore
+      !!      we allow the freedom to increase or decrease this value, if the
+      !!      Ekman depth estimate appears too shallow or too deep, respectively.
+      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
+      !!      namelist
+      !!        N.B. the mask are required for implicit scheme, and surface
+      !!      and bottom value already set in zdfphy.F90
+      !!
+      !! ** Action  :   avm, avt  mixing coeff (inner domain values only)
+      !!
+      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
+      !!              PFJ Lermusiaux 2001.
+      !!----------------------------------------------------------------------
+      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time-step
+      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdept          ! depth of t-point  [m]
+      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   ! momentum and tracer Kz (w-points)
+      !!
+      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
+      REAL(wp) ::   zcfRi, zav, zustar, zed     ! local scalars
+      REAL(wp), DIMENSION(WRK_2D) ::   zh_ekm   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('zdf_ric')
+      !
+      !                       !==  avm and avt = F(Richardson number)  ==!
+      DO jk = 2, jpkm1
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Jstr, k_Iend        ! coefficient = F(richardson number) (avm-weighted Ri)
+               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
+               zav   = rn_avmri * zcfRi**nn_ric
+               !                          ! avm and avt coefficients
+               p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
+               p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
+            END DO
+         END DO
       END DO
       !
-   END SUBROUTINE zdf_ric_init
+!!gm BUG <<<<====  This param can't work at low latitude 
+!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
+      !
+      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
+         !
+         DO jj = k_Jstr, k_Jend     !* Ekman depth
+            DO ji = k_Jstr, k_Iend
+               zustar = SQRT( taum(ji,jj) * r1_rau0 )
+               zed    = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )     ! Ekman depth
+               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zed , rn_mldmax )  )      ! set allowed range
+            END DO
+         END DO
+         DO jk = 2, jpkm1           !* minimum mixing coeff. within the Ekman layer
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Jstr, k_Iend
+                  IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN
+                     p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
+                     p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zdf_ric')
+      !
+   END SUBROUTINE zdf_ric
+
+
+   SUBROUTINE ric_rst( kt, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ric_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write TKE file (en) in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!                if the restart does not contain TKE, en is either 
+      !!                set to rn_emin or recomputed 
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !
+      INTEGER ::   jit, jk    ! dummy loop indices
+      INTEGER ::   id1, id2   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+         !                                   ! ---------------
+         !           !* Read the restart file
+         IF( ln_rstart ) THEN
+            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it
+               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k )
+               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
+            ENDIF
+         ENDIF
+         !           !* otherwise Kz already set to the background value in zdf_phy_init
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- ric-rst ----'
+         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
+         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
+         !
+      ENDIF
+      !
+   END SUBROUTINE ric_rst
 
    !!======================================================================

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90

@@ -5 +5 @@
    !!=====================================================================
    !! History :   -   !  2014-10  (A. Barthelemy, G. Madec)  original code
-   !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm
+   !!   NEMO     4.0  !  2017-04  (G. Madec)  remove u-,v-pts avm + ready for coarse-grained Open-MP
    !!----------------------------------------------------------------------
 
@@ -11 +11 @@
    !!   zdf_sh2       : compute mixing the shear production term of TKE
    !!----------------------------------------------------------------------
-   USE oce            ! ocean: shared variables
    USE dom_oce        ! domain: ocean
-   USE zdf_oce        ! vertical physics: variables
    !
    USE in_out_manager ! I/O manager
@@ -23 +21 @@
 
    PUBLIC   zdf_sh2        ! called by zdftke, zdfglf, and zdfric
-   
+
+   !! * Substitutions
+#  include "domain_substitute.h90"   
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2017)
@@ -31 +31 @@
 CONTAINS
 
-   SUBROUTINE zdf_sh2(  psh2  ) 
+   SUBROUTINE zdf_sh2( ARG_2D, pub, pvb, pun, pvn, p_avm, p_sh2  ) 
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_sh2  ***
@@ -46 +46 @@
       !!                NB: wet-point only horizontal averaging of shear
       !!
-      !! ** Action  : - psh2 shear prod. term at w-point (interior ocean domain only)
-      !!
+      !! ** Action  : - p_sh2 shear prod. term at w-point (inner domain only)
+      !!                                                   *****
       !! References :   Bruchard, OM 2002
       !! ---------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   psh2   ! shear production of TKE (w-points)
+      INTEGER                    , INTENT(in   ) ::   ARG_2D               ! inner domain start-end i-indices
+      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pub, pvb, pun, pvn   ! before, now horizontal velocities
+      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm                ! vertical eddy viscosity (w-points)
+      REAL(wp), DIMENSION(WRK_3D), INTENT(  out) ::   p_sh2                ! shear production of TKE (w-points)
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop arguments
-      REAL(wp), DIMENSION(jpi,jpj) ::   zsh2u, zsh2v   ! 2D workspace
+      INTEGER  ::   ji, jj, jk                           ! dummy loop indices
+      REAL(wp), DIMENSION(WRK_2De(-1,)) ::   zsh2u, zsh2v   ! 2D workspace
       !!--------------------------------------------------------------------
+      !
       IF( nn_timing == 1 )  CALL timing_start('zdf_sh2')
       !
       DO jk = 2, jpkm1
-         DO jj = 1, jpjm1        !* 2 x shear production at uw- and vw-points (energy conserving form)
-            DO ji = 1, jpim1
-               zsh2u(ji,jj) = ( avm(ji+1,jj,jk) + avm(ji,jj,jk) ) &
-                  &         * (  un(ji,jj,jk-1) -  un(ji,jj,jk) ) &
-                  &         * (  ub(ji,jj,jk-1) -  ub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk) 
-               zsh2v(ji,jj) = ( avm(ji,jj+1,jk) + avm(ji,jj,jk) ) &
-                  &         * (  vn(ji,jj,jk-1) -  vn(ji,jj,jk) ) &
-                  &         * (  vb(ji,jj,jk-1) -  vb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk)
+         DO jj = k_Jstr-1, k_Jend      !* 2 x shear production at uw- and vw-points (energy conserving form)
+            DO ji = k_Istr-1, k_Iend
+               zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
+                  &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) &
+                  &         * (   pub(ji,jj,jk-1) -   pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk)
+               zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) &
+                  &         * (   pvn(ji,jj,jk-1) -   pvn(ji,jj,jk) ) &
+                  &         * (   pvb(ji,jj,jk-1) -   pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk)
             END DO
          END DO
-         DO jj = 2, jpjm1        !* shear production at w-point
-            DO ji = 2, jpim1           ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
-               psh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
-                  &                      + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
+         DO jj = k_Jstr, k_Jend        !* shear production at w-point
+            DO ji = k_Jstr, k_Iend           ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
+               p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
+                  &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
             END DO
          END DO

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90

@@ -27 +27 @@
    PUBLIC zdf_swm_init    ! routine called in zdf_phy_init
 
+   !! * Substitutions
+#  include "domain_substitute.h90"   
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
@@ -34 +36 @@
 CONTAINS
 
-   SUBROUTINE zdf_swm( kt )
+   SUBROUTINE zdf_swm( ARG_2D, kt, p_avm, p_avt, p_avs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE zdf_swm ***
@@ -51 +53 @@
       !! reference : Qiao et al. GRL, 2004
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points)
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points)
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
@@ -59 +64 @@
       zcoef = 1._wp * 0.353553_wp
       DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk)
                !
-               avt(ji,jj,jk) = avt(ji,jj,jk) + zqb
-               avs(ji,jj,jk) = avs(ji,jj,jk) + zqb
-               avm(ji,jj,jk) = avm(ji,jj,jk) + zqb
+               p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zqb
+               p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zqb
+               p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zqb
             END DO
          END DO

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -43 +43 @@
    USE sbc_oce        ! surface boundary condition: ocean
    USE zdf_oce        ! vertical physics: ocean variables
-   USE zdfsh2         ! vertical physics: shear production term of TKE
    USE zdfmxl         ! vertical physics: mixed layer
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl         ! Print control
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! work arrays
-   USE timing         ! Timing
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 #if defined key_agrif
    USE agrif_opa_interp
    USE agrif_opa_update
 #endif
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -87 +87 @@
    REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
 
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
-   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau    ! depth of tke penetration (nn_htau)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl   ! now mixing lenght of dissipation
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr   ! now mixing lenght of dissipation
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "domain_substitute.h90"   
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO Consortium (2015)
@@ -112 +112 @@
 
 
-   SUBROUTINE zdf_tke( kt )
+   SUBROUTINE zdf_tke( ARG_2D, kt, p_sh2, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_tke  ***
@@ -157 +157 @@
       !!              Bruchard OM 2002
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step
+      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
       !!----------------------------------------------------------------------
       !
@@ -165 +168 @@
 #endif
       !
-      avt(:,:,:) = avt_k(:,:,:)     ! restore before Kz computed with TKE alone
-      avm(:,:,:) = avm_k(:,:,:) 
-      !
-      CALL tke_tke                  ! now tke (en)
-      !
-      CALL tke_avn                  ! now avt, avm, dissl
-      !
-      avt_k(:,:,:) = avt(:,:,:)     ! store Kz computed with TKE alone
-      avm_k(:,:,:) = avm(:,:,:) 
+      CALL tke_tke( ARG_2D, gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt )        ! now tke (en)
+      !
+
+!!gm not sure we need this lbc ....      ==>>>  certainly NOT !!!
+      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      CALL tke_avn( ARG_2D, gdepw_n, e3t_n, e3w_n, p_avm, p_avt )        ! now avt, avm, dissl
       !
 #if defined key_agrif
@@ -179 +180 @@
       IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
 #endif      
-      ! 
-      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )
       !
   END SUBROUTINE zdf_tke
 
 
-   SUBROUTINE tke_tke
+   SUBROUTINE tke_tke( ARG_2D, pdepw, p_e3t, p_e3w, p_sh2    &
+      &                                    , p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tke_tke  ***
@@ -200 +200 @@
       !! ** Action  : - en : now turbulent kinetic energy)
       !! ---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
+      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   pdepw          ! vertical eddy viscosity (w-points)
+      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_e3t, p_e3w   ! vertical eddy viscosity (w-points)
+      REAL(wp), DIMENSION(WRK_3D), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   !  vertical eddy viscosity (w-points)
+      !
+      INTEGER ::   ji, jj, jk              ! dummy loop arguments
 !!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! local integers
 !!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          !   -      -
 !!bfr      REAL(wp) ::   zebot                           ! local scalars
-      REAL(wp) ::   zrhoa  = 1.22            ! Air density kg/m3
-      REAL(wp) ::   zcdrag = 1.5e-3          ! drag coefficient
-      REAL(wp) ::   zbbrau, zri              ! local scalars
-      REAL(wp) ::   zfact1, zfact2, zfact3   !   -         -
-      REAL(wp) ::   ztx2  , zty2  , zcof     !   -         -
-      REAL(wp) ::   ztau  , zdif             !   -         -
-      REAL(wp) ::   zus   , zwlc  , zind     !   -         -
-      REAL(wp) ::   zzd_up, zzd_lw           !   -         -
-      INTEGER , DIMENSION(jpi,jpj) ::   imlc
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2
+      REAL(wp)::   zrhoa  = 1.22            ! Air density kg/m3
+      REAL(wp)::   zcdrag = 1.5e-3          ! drag coefficient
+      REAL(wp)::   zbbrau, zri              ! local scalars
+      REAL(wp)::   zfact1, zfact2, zfact3   !   -         -
+      REAL(wp)::   ztx2  , zty2  , zcof     !   -         -
+      REAL(wp)::   ztau  , zdif             !   -         -
+      REAL(wp)::   zus   , zwlc  , zind     !   -         -
+      REAL(wp)::   zzd_up, zzd_lw           !   -         -
+      INTEGER , DIMENSION(WRK_2D) ::   imlc
+      REAL(wp), DIMENSION(WRK_2D) ::   zhlc
+      REAL(wp), DIMENSION(WRK_3D) ::   zpelc, zdiag, zd_up, zd_lw
       !!--------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('tke_tke')
-      !
-      CALL wrk_alloc( jpi,jpj,       zhlc ) 
-      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw ) 
       !
       zbbrau = rn_ebb / rau0       ! Local constant initialisation
@@ -233 +235 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF ( ln_isfcav ) THEN
-         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend            ! en(mikt(ji,jj))   = rn_emin
+            DO ji = k_Istr, k_Iend
                en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
             END DO
          END DO
       END IF
-      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = k_Jstr, k_Jend            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+         DO ji = k_Istr, k_Iend
             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
          END DO
@@ -256 +258 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
-!!    DO jj = 2, jpjm1
-!!       DO ji = fs_2, fs_jpim1   ! vector opt.
+!!    DO jj = k_Jstr, k_Jend
+!!       DO ji = k_Jstr, k_Iend
 !!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
 !!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
@@ -269 +271 @@
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
+      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002)
          !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
          !
          !                        !* total energy produce by LC : cumulative sum over jk
-         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)
+         zpelc(WRK_2D,1) =  MAX( rn2b(WRK_2D,1), 0._wp ) * pdepw(WRK_2D,1) * p_e3w(WRK_2D,1)
          DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk)
+            zpelc(WRK_2D,jk) = zpelc(WRK_2D,jk-1)   &
+               &             + MAX( rn2b(WRK_2D,jk), 0._wp )   &
+               &             * pdepw(WRK_2D,jk) * p_e3w(WRK_2D,jk)
          END DO
          !                        !* finite Langmuir Circulation depth
          zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
-         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
+         imlc(WRK_2D) = mbkt(WRK_2D) + 1       ! Initialization to the number of w ocean point (=2 over land)
          DO jk = jpkm1, 2, -1
-            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
-               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
+            DO jj = k_Jstr, k_Jend               ! Last w-level at which zpelc>=0.5*us*us 
+               DO ji = k_Istr, k_Iend            !      with us=0.016*wind(starting from jpk-1)
                   zus  = zcof * taum(ji,jj)
                   IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
@@ -289 +293 @@
          END DO
          !                               ! finite LC depth
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))
+         DO jj = k_Jstr, k_Jend 
+            DO ji = k_Istr, k_Iend
+               zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj))
             END DO
          END DO
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
          DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
                   zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
                   !                                           ! vertical velocity due to LC
-                  zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) )
-                  zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )
+                  zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) )
+                  zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )
                   !                                           ! TKE Langmuir circulation source term
                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
@@ -318 +322 @@
       !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       !
-      !                           !* Shear production at uw- and vw-points (energy conserving form)
-      CALL zdf_sh2( zsh2 )
-      !
       IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
                   !                             ! local Richardson number
-                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zsh2(ji,jj,jk) + rn_bshear )
+                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
                   !                             ! inverse of Prandtl number
                   apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
@@ -335 +336 @@
       !         
       DO jk = 2, jpkm1           !* Matrix and right hand side in en
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                zcof   = zfact1 * tmask(ji,jj,jk)
                !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
                !                                   ! eddy coefficient (ensure numerical stability)
-               zzd_up = zcof * MAX(   avm(ji,jj,jk+1) +   avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
-                  &          /    ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  )  )
-               zzd_lw = zcof * MAX(   avm(ji,jj,jk  ) +   avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
-                  &          /    ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  )
+               zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
+                  &          /    (  p_e3t(ji,jj,jk  ) * p_e3w(ji,jj,jk  )  )
+               zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
+                  &          /    (  p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk  )  )
                !
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
@@ -350 +351 @@
                !
                !                                   ! right hand side in en
-               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zsh2(ji,jj,jk)                           &   ! shear
-                  &                                 - avt(ji,jj,jk) * rn2(ji,jj,jk)            &   ! stratification
+               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear
+                  &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification
                   &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
                   &                                ) * wmask(ji,jj,jk)
@@ -359 +360 @@
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
       DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
             END DO
          END DO
       END DO
-      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = k_Jstr, k_Jend                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+         DO ji = k_Istr, k_Iend
             zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
          END DO
       END DO
       DO jk = 3, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
             END DO
          END DO
       END DO
-      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = k_Jstr, k_Jend                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+         DO ji = k_Istr, k_Iend
             en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
          END DO
       END DO
       DO jk = jpk-2, 2, -1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1    ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
             END DO
@@ -390 +391 @@
       END DO
       DO jk = 2, jpkm1                             ! set the minimum value of tke
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
             END DO
@@ -405 +406 @@
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
                      &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
@@ -413 +414 @@
          END DO
       ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                jk = nmln(ji,jj)
-               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
+               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
                   &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
@@ -422 +423 @@
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
@@ -429 +430 @@
                   zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
-                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
+                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
                      &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
@@ -435 +436 @@
          END DO
       ENDIF
-!!gm not sure we need this lbc ....
-      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
-      !
-      CALL wrk_dealloc( jpi,jpj,       zhlc ) 
-      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
@@ -446 +442 @@
 
 
-   SUBROUTINE tke_avn
+   SUBROUTINE tke_avn( ARG_2D, pdepw, p_e3t, p_e3w, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tke_avn  ***
@@ -480 +476 @@
       !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
       !!----------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! vertical eddy viscosity (w-points)
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w     ! vertical eddy viscosity (w-points)
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   !  vertical eddy viscosity (w-points)
+      !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
       REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
       REAL(wp) ::   zemxl, zemlm, zemlp       !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
+      REAL(wp), DIMENSION(WRK_3D) ::   zmxlm, zmxld
       !!--------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('tke_avn')
 
-      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
-
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                     !  Mixing length
@@ -498 +497 @@
       !
       ! initialisation of interior minimum value (avoid a 2d loop with mikt)
-      zmxlm(:,:,:)  = rmxl_min    
-      zmxld(:,:,:)  = rmxl_min
+      zmxlm(WRK_3D)  = rmxl_min    
+      zmxld(WRK_3D)  = rmxl_min
       !
       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
             END DO
          END DO
       ELSE 
-         zmxlm(:,:,1) = rn_mxl0
+         zmxlm(WRK_2D,1) = rn_mxl0
       ENDIF
       !
       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
@@ -523 +522 @@
       !                     !* Physical limits for the mixing length
       !
-      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value 
-      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
+      zmxld(WRK_2D, 1 ) = zmxlm(WRK_2D,1)   ! surface set to the minimum value 
+      zmxld(WRK_2D,jpk) = rmxl_min          ! last level  set to the minimum value
       !
       SELECT CASE ( nn_mxl )
       !
  !!gm Not sure of that coding for ISF....
-      ! where wmask = 0 set zmxlm == e3w_n
+      ! where wmask = 0 set zmxlm == p_e3w
       CASE ( 0 )           ! bounded by the distance to surface and bottom
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
-                  &            gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) )
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
+                  &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) )
                   ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
-                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
-                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
+                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
+                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
                END DO
             END DO
@@ -545 +544 @@
       CASE ( 1 )           ! bounded by the vertical scale factor
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
                   zmxld(ji,jj,jk) = zemxl
@@ -556 +555 @@
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
          DO jk = 2, jpkm1         ! from the surface to the bottom :
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
                END DO
             END DO
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
                   zmxld(ji,jj,jk) = zemxl
@@ -574 +573 @@
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
          DO jk = 2, jpkm1         ! from the surface to the bottom : lup
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
                END DO
             END DO
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
                END DO
             END DO
          END DO
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
                   zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
                   zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
@@ -605 +604 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+         DO jj = k_Jstr, k_Jend
+            DO ji = k_Istr, k_Iend
                zsqen = SQRT( en(ji,jj,jk) )
                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
-               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
-               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
+               p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
+               p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
                dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
             END DO
@@ -619 +618 @@
       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
          DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
+            DO jj = k_Jstr, k_Jend
+               DO ji = k_Istr, k_Iend
+                  p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
               END DO
             END DO
          END DO
       ENDIF
-!!gm I'm sure this is needed to compute the shear term at next time-step
-      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
-      CALL lbc_lnk( avt, 'W', 1. )      ! Lateral boundary conditions on avt  (sign unchanged)
 
       IF(ln_ctl) THEN
@@ -634 +630 @@
          CALL prt_ctl( tab3d_1=avm, clinfo1=' tke  - m: ', ovlap=1, kdim=jpk )
       ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
@@ -737 +731 @@
       !                               !* set vertical eddy coef. to the background value
       DO jk = 1, jpk
-         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
-         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
+         avt(:,:,jk) = avtb(jk) * wmask(:,:,jk)
+         avm(:,:,jk) = avmb(jk) * wmask(:,:,jk)
       END DO
       dissl(:,:,:) = 1.e-12_wp
@@ -748 +742 @@
 
    SUBROUTINE tke_rst( kt, cdrw )
-     !!---------------------------------------------------------------------
-     !!                   ***  ROUTINE tke_rst  ***
-     !!                     
-     !! ** Purpose :   Read or write TKE file (en) in restart file
-     !!
-     !! ** Method  :   use of IOM library
-     !!                if the restart does not contain TKE, en is either 
-     !!                set to rn_emin or recomputed 
-     !!----------------------------------------------------------------------
-     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
-     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
-     !
-     INTEGER ::   jit, jk              ! dummy loop indices
-     INTEGER ::   id1, id2, id3, id4   ! local integers
-     !!----------------------------------------------------------------------
-     !
-     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
-        !                                   ! ---------------
-        IF( ln_rstart ) THEN                   !* Read the restart file
-           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
-           id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
-           id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
-           id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
-           !
-           IF( id1 > 0 ) THEN                       ! 'en' exists
-              CALL iom_get( numror, jpdom_autoglo, 'en', en )
-              IF( MIN( id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
-                 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k )
-                 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
-                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
-              ELSE                                                 ! one at least array is missing
-                 CALL tke_avn                                          ! compute avt, avm, and dissl (approximation)
-              ENDIF
-           ELSE                                     ! No TKE array found: initialisation
-              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
-              en (:,:,:) = rn_emin * wmask(:,:,:)
-              CALL tke_avn                               ! recompute avt, avm, and dissl (approximation)
-              avt_k(:,:,:) = avt(:,:,:)
-              avm_k(:,:,:) = avm(:,:,:)
-              !
-              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
-              avt_k(:,:,:) = avt(:,:,:)
-              avm_k(:,:,:) = avm(:,:,:)
-           ENDIF
-        ELSE                                   !* Start from rest
-           en(:,:,:) = rn_emin * tmask(:,:,:)
-           DO jk = 1, jpk                           ! set the Kz to the background value
-              avt_k(:,:,jk) = avtb(jk) * wmask (:,:,jk)
-              avm_k(:,:,jk) = avmb(jk) * wmask (:,:,jk)
-           END DO
-        ENDIF
-        !
-     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-        !                                   ! -------------------
-        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
-        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
-        CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
-        CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
-        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
-        !
-     ENDIF
-     !
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE tke_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write TKE file (en) in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!                if the restart does not contain TKE, en is either 
+      !!                set to rn_emin or recomputed 
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !
+      INTEGER ::   jit, jk              ! dummy loop indices
+      INTEGER ::   id1, id2, id3, id4   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
+            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
+            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
+               CALL iom_get( numror, jpdom_autoglo, 'en', en )
+               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k )
+               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
+               CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
+            ELSE                                          ! start TKE from rest
+               IF(lwp) WRITE(numout,*) '   ==>>   previous run without TKE scheme, set en to background values'
+               en(:,:,:) = rn_emin * wmask(:,:,:)
+               ! avt_k, avm_k already set to the background value in zdf_phy_init
+            ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set en to the background value'
+            en(:,:,:) = rn_emin * wmask(:,:,:)
+            ! avt_k, avm_k already set to the background value in zdf_phy_init
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- tke-rst ----'
+         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
+         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
+         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
+         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
+         !
+      ENDIF
+      !
    END SUBROUTINE tke_rst
 

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -36 +36 @@
    !!----------------------------------------------------------------------
    USE step_oce         ! time stepping definition modules
+!!gm to be removed when removing avmu, avmv
+   USE zdf_oce        ! ocean vertical physics variables
+!!gm
    !
    USE iom              ! xIOs server
@@ -102 +105 @@
       IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
 
+CALL FLUSH( numout)
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
@@ -109 +113 @@
       IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
                          CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice)
+CALL FLUSH( numout)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -124 +129 @@
                          CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
                          CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+CALL FLUSH( numout)
 
       !  VERTICAL PHYSICS
                          CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
+CALL FLUSH( numout)
+
+!!gm  ===>>>   TO BE REMOVED   (require changes in zdf_oce, dynzdf(_imp,_exp), dynldf_iso, diawri)
+      DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
+!!gm end
 
       !  LATERAL  PHYSICS

branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/sette.sh

@@ -1098 +1098 @@
     export TEST_NAME="SHORT"
     cd ${CONFIG_DIR0}
-    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_1_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top"
+    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_1_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top"
     cd ${SETTE_DIR}
     . ./param.cfg
@@ -1145 +1145 @@
     export TEST_NAME="SHORT_NOZOOM"
     cd ${CONFIG_DIR0}
-    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top"
+    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top"
     cd ${SETTE_DIR}
     . ./param.cfg
@@ -1182 +1182 @@
     export TEST_NAME="SHORT_NOAGRIF"
     cd ${CONFIG_DIR0}
-    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2_NAG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 del_key "key_zdftmx key_top"
+    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_2_2_NAG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 del_key "key_top"
     cd ${SETTE_DIR}
     . ./param.cfg
@@ -1220 +1220 @@
     export TEST_NAME="LONG"
     cd ${CONFIG_DIR0}
-    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_LONG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top"
+    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_LONG -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top"
     cd ${SETTE_DIR}
     . ./param.cfg
@@ -1319 +1319 @@
     export TEST_NAME="REPRO_4_4"
     cd ${CONFIG_DIR0}
-    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_16 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_zdftmx key_top"
+    . ./makenemo -m ${CMP_NAM} -n ORCA2AGUL_16 -r ORCA2_LIM3_PISCES -d "OPA_SRC LIM_SRC_3 NST_SRC" -j 8 add_key "key_agrif" del_key "key_top"
     cd ${SETTE_DIR}
     . ./param.cfg