Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/ARCH/CMCC/arch-ifort_athena_xios_omp.fcm
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/ARCH/CMCC/arch-ifort_athena_xios_omp.fcm	(revision 8279)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/ARCH/CMCC/arch-ifort_athena_xios_omp.fcm	(revision 8279)
@@ -0,0 +1,62 @@
+# mpi ifort compiler options for ATHENA using XIOS parallel writer server
+#
+#
+# NCDF_HOME   root directory containing lib and include subdirectories for netcdf4
+# HDF5_HOME   root directory containing lib and include subdirectories for HDF5
+# XIOS_HOME   root directory containing lib for XIOS
+# OASIS_HOME  root directory containing lib for OASIS
+#
+# NCDF_INC    netcdf4 include file
+# NCDF_LIB    netcdf4 library
+# XIOS_INC    xios include file    (taken into accound only if key_iomput is activated)
+# XIOS_LIB    xios library         (taken into accound only if key_iomput is activated)
+# OASIS_INC   oasis include file   (taken into accound only if key_oasis3 is activated)
+# OASIS_LIB   oasis library        (taken into accound only if key_oasis3 is activated)
+#
+# FC          Fortran compiler command
+# FCFLAGS     Fortran compiler flags
+# FFLAGS      Fortran 77 compiler flags
+# LD          linker
+# LDFLAGS     linker flags, e.g. -L<lib dir> if you have libraries
+# FPPFLAGS    pre-processing flags
+# AR          assembler
+# ARFLAGS     assembler flags
+# MK          make
+# USER_INC    complete list of include files
+# USER_LIB    complete list of libraries to pass to the linker
+# CC          C compiler used to compile conv for AGRIF
+# CFLAGS      compiler flags used with CC
+#
+# Note that:
+#  - unix variables "$..." are accpeted and will be evaluated before calling fcm.
+#  - fcm variables are starting with a % (and not a $)
+#
+
+# required modules
+# module load  INTEL/intel_xe_2013 NETCDF/netcdf-4.3_parallel NETCDF/parallel-netcdf-1.7.0 HDF5/hdf5-1.8.11_parallel
+
+# NETCDF and PNETCDF should be set automatically when loading modules.
+# The following environment variables must be set by the user.
+#export XIOS=/users/home/models/nemo/xios
+#export HDF5=/users/home/opt/hdf5/hdf5-1.8.11_parallel
+
+%NCDF_INC            -I${NETCDF}/include -I${PNETCDF}/include
+%NCDF_LIB            -L${NETCDF}/lib -lnetcdff -lnetcdf -L${PNETCDF}/lib -lpnetcdf
+%HDF5_INC            -I${HDF5}/include
+%HDF5_LIB            -L${HDF5}/lib -lhdf5_hl -lhdf5
+%XIOS_INC            -I${XIOS}/inc
+%XIOS_LIB            -L${XIOS}/lib -lxios
+%CPP                 cpp
+%FC                  mpiifort
+%FCFLAGS             -openmp -r8 -O3 -xHost -fp-model source -traceback
+%FFLAGS              %FCFLAGS
+%LD                  mpiifort
+%FPPFLAGS            -P -C -traditional
+%LDFLAGS             -openmp -lstdc++ -lz -lgpfs -lcurl
+%AR                  ar 
+%ARFLAGS             -r
+%MK                  gmake
+%USER_INC            %XIOS_INC %NCDF_INC %HDF5_INC
+%USER_LIB            %XIOS_LIB %NCDF_LIB %HDF5_LIB
+%CC                  icc
+%CFLAGS              -O0
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 8279)
@@ -255,4 +255,8 @@
 /
 !-----------------------------------------------------------------------
+&namomp        !   Open Multi-Processing
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namctl        !   Control prints & Benchmark
 !-----------------------------------------------------------------------
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8279)
@@ -12,5 +12,5 @@
 !!              7 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf)
 !!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_iwm)
-!!              9 - miscellaneous    (nammpp, namctl)
+!!              9 - miscellaneous    (nammpp, namomp, namctl)
 !!             10 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb, namsto)
 !!             11 - Obs & Assim      (namobs, nam_asminc)
@@ -984,4 +984,5 @@
 !!======================================================================
 !!   nammpp            Massively Parallel Processing                    ("key_mpp_mpi")
+!!   namomp            Open Multi-Processing
 !!   namctl            Control prints 
 !!   namsto            Stochastic parametrization of EOS
@@ -999,4 +1000,8 @@
    jpnij       =   0       !  jpnij  number of local domains (set automatically if < 1)
 /
+!-----------------------------------------------------------------------
+&namomp        !   Open Multi-Processing
+!-----------------------------------------------------------------------
+   jtni        =   1       !  jtni   number of openMP threads following i
 !-----------------------------------------------------------------------
 &namctl        !   Control prints 
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/DOM/dom_omp.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/DOM/dom_omp.F90	(revision 8279)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/DOM/dom_omp.F90	(revision 8279)
@@ -0,0 +1,112 @@
+MODULE dom_omp
+   !!======================================================================
+   !!                      ***  MODULE  domain_omp  ***
+   !! OpenMP horizontal decoposition : OpenMP threads domain decomposition 
+   !!======================================================================
+   !! History :  4.0  !  2017-05  (S. Mocavero - F. Mele)  original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tile_dec_2D  : initialize local threads indices
+   !!----------------------------------------------------------------------
+
+  USE par_oce        ! ocean space and time domain
+  USE in_out_manager ! I/O manager
+
+  IMPLICIT NONE
+
+  INTEGER         :: jtni 
+  INTEGER         :: jtnj
+  INTEGER         :: jtnij
+
+  INTEGER         :: tnarea
+  INTEGER         :: tnareai
+  INTEGER         :: tnareaj
+
+  INTEGER,SAVE    :: tnldi,tnldj
+  INTEGER,SAVE    :: tnlei,tnlej
+!$OMP THREADPRIVATE(tnldi,tnldj,tnlei,tnlej)
+
+CONTAINS
+  SUBROUTINE tile_dec_2D
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE tile_dec_2D  ***
+      !!                    
+      !! ** Purpose :   read namomp namelists and initialize local threads indices
+      !!
+      !! ** input   : - namomp namelist
+      !!----------------------------------------------------------------------
+      INTEGER :: ii,jj
+      INTEGER :: tnlci,tnlcj
+      INTEGER :: omp_get_thread_num
+      INTEGER :: omp_get_num_threads
+      INTEGER :: ios
+      !
+      NAMELIST/namomp/ jtni
+      !!----------------------------------------------------------------------
+      REWIND( numnam_cfg )
+      READ  ( numnam_cfg, namomp, IOSTAT = ios )
+
+      CALL OMP_SET_DYNAMIC(.FALSE.)
+
+!$OMP PARALLEL PRIVATE(tnarea,tnareai,tnareaj,ii,jj,tnlci,tnlcj)
+!$OMP SINGLE
+      jtnij = omp_get_num_threads()
+        
+      jtnj = jtnij / jtni          ! threads number along domain latitude
+      IF((jtni.GT.jtnij).OR.(MOD(jtnij,jtni).NE.0)) THEN
+         WRITE(*,*) ' tile_dec_2D: Error while setting OpenMP decomposition' 
+         WRITE(*,*) ' -> disabled OpenMP decomposition on  direction' 
+         jtni = 1
+         jtnj = jtnij
+      ENDIF
+!$OMP END SINGLE
+
+      tnarea = omp_get_thread_num()
+      tnareai = MOD(tnarea, jtni)
+      tnareaj = tnarea / jtni
+
+      tnlei=0
+      DO ii=0,tnareai
+         tnldi=tnlei+1
+         tnlci=jpi/jtni
+         IF (MOD(jpi,jtni)>ii) tnlci=tnlci+1
+         tnlei=tnldi+tnlci-1
+      ENDDO
+
+      tnlej=0
+      DO jj=0,tnareaj
+         tnldj=tnlej+1
+         tnlcj=jpj/jtnj
+         IF (MOD(jpj,jtnj)>jj) tnlcj=tnlcj+1
+         tnlej=tnldj+tnlcj-1
+      ENDDO
+
+      IF(tnareai == 0) THEN
+#if defined key_vectopt_loop
+         tnldi = 1
+#else
+         tnldi = 2
+#endif
+      END IF
+
+      IF(tnareaj == 0) THEN
+         tnldj = 2
+      END IF
+
+      IF(tnareai == jtni-1) THEN
+#if defined key_vectopt_loop
+         tnlei = jpi
+#else
+         tnlei = jpim1
+#endif
+      ENDIF
+
+      IF(tnareaj == jtnj-1) THEN
+         tnlej = jpjm1
+      ENDIF
+!$OMP END PARALLEL
+
+  END SUBROUTINE tile_dec_2D
+
+END MODULE dom_omp
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 8279)
@@ -95,6 +95,6 @@
 !!gm                            and many acces in memory
          
-         DO jj = k_Jstr, k_Jend        !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
-            DO ji = k_Jstr, k_Iend
+         DO jj = tnldj, tnlej        !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
+            DO ji = tnldi, tnlei
                zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
 !!gm please, use e3w_n below 
@@ -113,6 +113,6 @@
          END DO
 
-         DO jj = k_Jstr, k_Jend          !==  indicators  ==!
-            DO ji = k_Jstr, k_Iend
+         DO jj = tnldj, tnlej          !==  indicators  ==!
+            DO ji = tnldi, tnlei
                ! stability indicator: msks=1 if rn2>0; 0 elsewhere
                IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
@@ -145,6 +145,6 @@
          ! ------------------
          ! Constant eddy coefficient: reset to the background value
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Jstr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                zinr = 1._wp / zrau(ji,jj)
                ! salt fingering
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 8279)
@@ -98,6 +98,6 @@
          !
          DO jk = 1, jpkm1 
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Jstr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
                      p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8279)
@@ -159,10 +159,10 @@
       ! Preliminary computing
 
-      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
+      ustars2(WRK_2D) = 0._wp   ;   ustarb2(WRK_2D) = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
 
       ! Compute surface and bottom friction at T-points
-      DO jj = k_Jstr, k_Jend 
-         DO ji = k_Istr, k_Iend
+      DO jj = tnldj, tnlej 
+         DO ji = tnldi, tnlei
             !
             ! surface friction
@@ -195,6 +195,6 @@
       !
       DO jk = 2, jpkm1              !==  Compute dissipation rate  ==!
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
             END DO
@@ -208,6 +208,6 @@
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   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) )
@@ -233,6 +233,6 @@
 
       DO jk = 2, jpkm1
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                !
                buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)       ! stratif. destruction
@@ -329,6 +329,6 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -350,6 +350,6 @@
       CASE ( 1 )             ! Neumman boundary condition
          !                      
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -373,6 +373,6 @@
       !
       DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                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
@@ -380,6 +380,6 @@
       END DO
       DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                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
@@ -387,6 +387,6 @@
       END DO
       DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                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
@@ -406,6 +406,6 @@
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
                END DO
@@ -415,6 +415,6 @@
       CASE( 1 )               ! k-eps
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   psi(ji,jj,jk)  = eps(ji,jj,jk)
                END DO
@@ -424,6 +424,6 @@
       CASE( 2 )               ! k-w
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
                END DO
@@ -433,6 +433,6 @@
       CASE( 3 )               ! generic
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
                END DO
@@ -450,6 +450,6 @@
 
       DO jk = 2, jpkm1
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                !
                ! psi / k
@@ -548,6 +548,6 @@
          !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -569,6 +569,6 @@
       CASE ( 1 )             ! Neumman boundary condition
          !                      
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
@@ -600,6 +600,6 @@
       !
       DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                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
@@ -607,6 +607,6 @@
       END DO
       DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                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
@@ -614,6 +614,6 @@
       END DO
       DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                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
@@ -628,6 +628,6 @@
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
          DO jk = 1, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   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
@@ -637,6 +637,6 @@
       CASE( 1 )               ! k-eps
          DO jk = 1, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   eps(ji,jj,jk) = psi(ji,jj,jk)
                END DO
@@ -646,6 +646,6 @@
       CASE( 2 )               ! k-w
          DO jk = 1, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
                END DO
@@ -658,6 +658,6 @@
          zex2  = -1._wp / rnn
          DO jk = 1, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
                END DO
@@ -670,6 +670,6 @@
       ! --------------------------------------------------
       DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                ! limitation
                eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
@@ -690,6 +690,6 @@
       CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   ! zcof =  l²/q²
                   zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
@@ -711,6 +711,6 @@
       CASE ( 2, 3 )               ! Canuto stability functions
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend 
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej 
+               DO ji = tnldi, tnlei
                   ! zcof =  l²/q²
                   zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
@@ -745,6 +745,6 @@
 
 !!gm should be done for ISF (top boundary cond.)
-      DO jj = k_Jstr, k_Jend 
-         DO ji = k_Istr, k_Iend
+      DO jj = tnldj, tnlej 
+         DO ji = tnldi, tnlei
             zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
          END DO
@@ -754,6 +754,6 @@
       ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
       DO jk = 1, jpk
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
                zav           = zsqen * zstt(ji,jj,jk)
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90	(revision 8279)
@@ -150,6 +150,6 @@
       !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth,
       !                                                 using an exponential decay from the seafloor.
-      DO jj = k_Jstr, k_Jend
-         DO ji = k_Istr, k_Iend     ! part independent of the level
+      DO jj = tnldj, tnlej
+         DO ji = tnldi, tnlei     ! 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,6 +158,6 @@
       END DO
       DO jk = 2, jpkm1              ! complete with the level-dependent part
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                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) )      &
@@ -185,6 +185,6 @@
          !
          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
+!         DO jj = tnldj, tnlej
+!            DO ji = tnldi, tnlei
 !               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
 !            END DO
@@ -203,6 +203,6 @@
          !
          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
+!         DO jj = tnldj, tnlej
+!            DO ji = tnldi, tnlei
 !               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
 !            END DO
@@ -224,6 +224,6 @@
 !         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
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                zztemp = zfact(ji,jj) + e3w_n(ji,jj,jk) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  )
                zfact(ji,jj)    = zztemp
@@ -234,6 +234,6 @@
       !
       DO jk = 2, jpkm1
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
                   &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
@@ -245,5 +245,5 @@
       zweight(WRK_3D) = 0._wp
       DO jk = 2, jpkm1
-         zweight(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * hbot_iwm(:,:)     &
+         zweight(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * hbot_iwm(WRK_2D)  &
             &               * (  EXP( -zwkb(WRK_2D,jk  ) / hbot_iwm(WRK_2D) )  &
             &                  - EXP( -zwkb(WRK_2D,jk-1) / hbot_iwm(WRK_2D) )  )
@@ -256,6 +256,6 @@
       !
       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
+!      DO jj = tnldj, tnlej
+!         DO ji = tnldi, tnlei
 !            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
 !         END DO
@@ -290,6 +290,6 @@
       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 = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   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) )
@@ -310,6 +310,6 @@
 !!gm used of glosum 3D....
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   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)
@@ -317,5 +317,5 @@
             END DO
          END DO
-         IF( lk_mpp )   CALL mpp_sum( zztemp )
+!         IF( lk_mpp )   CALL mpp_sum( zztemp )
          zztemp = rau0 * zztemp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
          !
@@ -335,6 +335,6 @@
       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 = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   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   )    )   &
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 8279)
@@ -32,4 +32,7 @@
    USE lbclnk         ! lateral boundary conditions
    USE lib_mpp        ! distribued memory computing
+#if defined key_omp
+   USE dom_omp        ! omp threads decomposition of MPI subdomain
+#endif
 
    IMPLICIT NONE
@@ -47,4 +50,5 @@
 
    LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC))
+   INTEGER :: ARG_2D
 
    !! * Substitutions
@@ -216,31 +220,28 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER ::   ji, jj, jk   ! dummy loop indice
+      INTEGER :: 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 (kt==nn_it000) THEN
+#ifdef key_omp
+         CALL tile_dec_2D
+#else
 #if defined key_vectopt_loop
-      k_Istr = 1   ;   k_Iend = jpi
+         tnldi = 1   ;   tnlei = jpi
 #else
-      k_Istr = 2   ;   k_Iend = jpim1
+         tnldi = 2   ;   tnlei = jpim1
 #endif
-      k_Jstr = 2   ;   k_Jend = jpjm1
-      !
-!!OMP <<<<===   end
-      
+         tnldj = 2   ;   tnlej = jpjm1
+#endif
+      END IF
+      !
+      CALL zdf_bfr( kt )                        !* bottom friction (if quadratic)
+      !
+!$OMP PARALLEL PRIVATE(zsh2)
       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
@@ -248,4 +249,5 @@
 
       !
+!$OMP BARRIER
       SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
       CASE( np_RIC )   ;   CALL zdf_ric( ARG_2D, kt, gdept_n, zsh2, avm_k, avt_k )         ! Richardson number dependent Kz
@@ -257,4 +259,6 @@
 !!gm         avm_k(WRK_3D) = rn_avm0 * wmask(WRK_3D)
       END SELECT
+      !
+      DEALLOCATE( zsh2 )
       !  
       !                          !==  ocean Kz  ==!   (avt, avs, avm)
@@ -271,5 +275,5 @@
       !
       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
@@ -282,5 +286,6 @@
       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)
-
+!$OMP END PARALLEL
+      !
 
       !                                         !* Lateral boundary conditions (sign unchanged)
@@ -301,6 +306,4 @@
       ENDIF
       !
-      DEALLOCATE( zsh2 )
-      !
    END SUBROUTINE zdf_phy
 
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 8279)
@@ -164,6 +164,6 @@
       !                       !==  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)
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei        ! coefficient = F(richardson number) (avm-weighted Ri)
                zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , p_avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
                zav   = rn_avmri * zcfRi**nn_ric
@@ -180,6 +180,6 @@
       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
+         DO jj = tnldj, tnlej     !* Ekman depth
+            DO ji = tnldi, tnlei
                zustar = SQRT( taum(ji,jj) * r1_rau0 )
                zed    = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )     ! Ekman depth
@@ -188,6 +188,6 @@
          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
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   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)
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfsh2.F90	(revision 8279)
@@ -62,6 +62,6 @@
       !
       DO jk = 2, jpkm1
-         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
+         DO jj = tnldj-1, tnlej      !* 2 x shear production at uw- and vw-points (energy conserving form)
+            DO ji = tnldi-1, tnlei
                zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
                   &         * (   pun(ji,jj,jk-1) -   pun(ji,jj,jk) ) &
@@ -72,6 +72,6 @@
             END DO
          END DO
-         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)
+         DO jj = tnldj, tnlej        !* shear production at w-point
+            DO ji = tnldi, tnlei           ! 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) )   )
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90	(revision 8279)
@@ -64,6 +64,6 @@
       zcoef = 1._wp * 0.353553_wp
       DO jk = 2, jpkm1
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk)
                !
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 8279)
@@ -220,4 +220,5 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('tke_tke')
+
       !
       zbbrau = rn_ebb / rau0       ! Local constant initialisation
@@ -231,12 +232,12 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF ( ln_isfcav ) THEN
-         DO jj = k_Jstr, k_Jend            ! en(mikt(ji,jj))   = rn_emin
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej            ! en(mikt(ji,jj))   = rn_emin
+            DO ji = tnldi, tnlei
                en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
             END DO
          END DO
       END IF
-      DO jj = k_Jstr, k_Jend            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-         DO ji = k_Istr, k_Iend
+      DO jj = tnldj, tnlej            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+         DO ji = tnldi, tnlei
             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
          END DO
@@ -254,6 +255,6 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
-!!    DO jj = k_Jstr, k_Jend
-!!       DO ji = k_Jstr, k_Iend
+!!    DO jj = tnldj, tnlej
+!!       DO ji = tnldj, tnlei
 !!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
 !!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
@@ -281,6 +282,6 @@
          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 = 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)
+            DO jj = tnldj, tnlej               ! Last w-level at which zpelc>=0.5*us*us 
+               DO ji = tnldi, tnlei            !      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,6 +290,6 @@
          END DO
          !                               ! finite LC depth
-         DO jj = k_Jstr, k_Jend 
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej 
+            DO ji = tnldi, tnlei
                zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj))
             END DO
@@ -296,6 +297,6 @@
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
          DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
                   !                                           ! vertical velocity due to LC
@@ -320,6 +321,6 @@
       IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   !                             ! local Richardson number
                   zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
@@ -332,6 +333,6 @@
       !         
       DO jk = 2, jpkm1           !* Matrix and right hand side in en
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                zcof   = zfact1 * tmask(ji,jj,jk)
                !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
@@ -356,30 +357,30 @@
       !                          !* 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 = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                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 = k_Jstr, k_Jend                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-         DO ji = k_Istr, k_Iend
+      DO jj = tnldj, tnlej                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+         DO ji = tnldi, tnlei
             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 = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                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 = k_Jstr, k_Jend                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-         DO ji = k_Istr, k_Iend
+      DO jj = tnldj, tnlej                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+         DO ji = tnldi, tnlei
             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 = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                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
@@ -387,6 +388,6 @@
       END DO
       DO jk = 2, jpkm1                             ! set the minimum value of tke
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
             END DO
@@ -402,6 +403,6 @@
       IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   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)
@@ -410,6 +411,6 @@
          END DO
       ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                jk = nmln(ji,jj)
                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
@@ -419,6 +420,6 @@
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
@@ -433,4 +434,5 @@
       ENDIF
       !
+
       IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
       !
@@ -498,6 +500,6 @@
       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
          zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
             END DO
@@ -508,6 +510,6 @@
       !
       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
@@ -527,6 +529,6 @@
       CASE ( 0 )           ! bounded by the distance to surface and bottom
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   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) )
@@ -540,6 +542,6 @@
       CASE ( 1 )           ! bounded by the vertical scale factor
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
@@ -551,6 +553,6 @@
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
          DO jk = 2, jpkm1         ! from the surface to the bottom :
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
                END DO
@@ -558,6 +560,6 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
                   zmxlm(ji,jj,jk) = zemxl
@@ -569,6 +571,6 @@
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
          DO jk = 2, jpkm1         ! from the surface to the bottom : lup
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
                END DO
@@ -576,6 +578,6 @@
          END DO
          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
                END DO
@@ -583,6 +585,6 @@
          END DO
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
                   zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
@@ -600,6 +602,6 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
-         DO jj = k_Jstr, k_Jend
-            DO ji = k_Istr, k_Iend
+         DO jj = tnldj, tnlej
+            DO ji = tnldi, tnlei
                zsqen = SQRT( en(ji,jj,jk) )
                zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
@@ -614,6 +616,6 @@
       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
          DO jk = 2, jpkm1
-            DO jj = k_Jstr, k_Jend
-               DO ji = k_Istr, k_Iend
+            DO jj = tnldj, tnlej
+               DO ji = tnldi, tnlei
                   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
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/domain_substitute.h90
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/domain_substitute.h90	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/domain_substitute.h90	(revision 8279)
@@ -4,9 +4,9 @@
    !! ** purpose :   local Sub-Domain: arguments, dimension & semi-column loops
    !!----------------------------------------------------------------------
-#define ARG_2D   k_Istr, k_Iend, k_Jstr, k_Jend
+#define ARG_2D   tnldi, tnlei, tnldj, tnlej
 
-#define WRK_2D             k_Istr:k_Iend,k_Jstr:k_Jend
-#define WRK_2De(s,e)       k_Istr s:k_Iend e,k_Jstr s:k_Jend e
+#define WRK_2D             tnldi:tnlei,tnldj:tnlej
+#define WRK_2De(s,e)       tnldi s:tnlei e,tnldj s:tnlej e
 
-#define WRK_3D             k_Istr:k_Iend,k_Jstr:k_Jend,1:jpk
-#define WRK_3De(s,e,t,b)   k_Istr s:k_Iend e,k_Jstr s:k_Jend e,1 t:jpk t
+#define WRK_3D             tnldi:tnlei,tnldj:tnlej,1:jpk
+#define WRK_3De(s,e,t,b)   tnldi s:tnlei e,tnldj s:tnlej e,1 t:jpk t
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-ifort_athena_xios_omp
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-ifort_athena_xios_omp	(revision 8279)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-ifort_athena_xios_omp	(revision 8279)
@@ -0,0 +1,112 @@
+#!/bin/bash 
+#!
+#BSUB -q poe_short
+#BSUB -n TOTAL_NPROCS
+#BSUB -a poe
+#BSUB -R "span[ptile=NPROC_NODE]"
+#BSUB -J NEMO_SETTE
+#BSUB -o poe.stdout.%J
+#BSUB -e poe.stderr.%J
+#BSUB -x
+
+ 
+
+###############################################################
+# Test specific settings. Do not hand edit these lines; the fcm_job.sh script will set these
+# (via sed operating on this template job file). 
+#
+  OCEANCORES=NPROCS
+  XIOS_NUMPROCS=NXIOPROCS
+  export SETTE_DIR=DEF_SETTE_DIR
+
+  export MP_TASK_AFFINITY=mcm
+  export OMP_NUM_THREADS=OMP_NTHR
+  export OMP_STACKSIZE="500M"
+
+###############################################################
+#
+# set up mpp computing environment
+#
+# Local settings for machine IBM Power6 (VARGAS at IDRIS France)
+#
+
+# export MPIRUN="mpirun -n $OCEANCORES"
+
+module load  INTEL/intel_xe_2013 NETCDF/netcdf-4.3_parallel NETCDF/parallel-netcdf-1.3.1 HDF5/hdf5-1.8.11_parallel
+export MPIRUN="mpirun.lsf"
+
+export LD_LIBRARY_PATH=${XIOS}/lib:${LD_LIBRARY_PATH}
+
+XIOS_SERVER_PATHNAME="${XIOS}/bin/xios_server.exe"
+
+#
+# load sette functions (only post_test_tidyup needed)
+#
+  . ${SETTE_DIR}/all_functions.sh
+
+
+# Don't remove neither change the following line 
+# BODY
+
+#
+# These variables are needed by post_test_tidyup function in all_functions.sh
+#
+  export EXE_DIR=DEF_EXE_DIR
+  export INPUT_DIR=DEF_INPUT_DIR
+  export CONFIG_DIR=DEF_CONFIG_DIR
+  export NEMO_VALIDATION_DIR=DEF_NEMO_VALIDATION
+  export NEW_CONF=DEF_NEW_CONF
+  export CMP_NAM=DEF_CMP_NAM
+  export TEST_NAME=DEF_TEST_NAME
+#
+# end of set up
+
+
+###############################################################
+#
+# change to the working directory 
+#
+cd ${EXE_DIR}
+
+#
+  echo Running on host `hostname`
+  echo Time is `date`
+  echo Directory is `pwd`
+# 
+#  Run the parallel MPI executable 
+#
+#  echo "Running time ${MPIRUN} " $OCEANCORES " ./opa"
+#
+
+  if [ MPI_FLAG == "yes" ]; then
+     if [ $XIOS_NUMPROCS -eq 0 ]; then
+       time ${MPIRUN} ./opa
+     else
+
+       export MP_PGMMODEL=mpmd
+       export MP_CMDFILE=xioscmdfile
+       rm xioscmdfile &> /dev/null
+       touch xioscmdfile
+       k=0
+       while (( k < $OCEANCORES )) ; do
+          echo ./opa >> xioscmdfile
+          k=`expr $k + 1`
+       done
+       k=0
+       while (( k < $XIOS_NUMPROCS )) ; do
+          echo $XIOS_SERVER_PATHNAME >> xioscmdfile
+          k=`expr $k + 1`
+       done
+
+       time ${MPIRUN}
+     fi
+  else
+     time ./opa
+  fi
+#
+  post_test_tidyup
+
+# END_BODY
+# Don't remove neither change the previous line 
+
+  exit
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/prepare_job.sh
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/prepare_job.sh	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/prepare_job.sh	(revision 8279)
@@ -94,4 +94,5 @@
 JOB_FILE=$5
 NXIO_PROC=$6
+NOMP_THR=$7
 
 # export EXE_DIR. This directory is used to execute model 
@@ -293,4 +294,7 @@
                                 echo NB_PROC_NODE ${NB_PROC_NODE}
                                 ;;
+                        ifort_athena_*)
+                                NB_PROC_NODE=$(( 16 / NOMP_THR ))
+                                ;;
 			*)
 				NB_NODES=${NB_PROC}
@@ -305,4 +309,5 @@
              -e"s/TOTAL_NPROCS/${TOTAL_NPROCS}/" \
              -e"s/NPROCS/${NB_PROC}/" \
+             -e"s/OMP_NTHR/${NOMP_THR}/" \
              -e"s/NXIOPROCS/${NXIO_PROC}/" \
              -e"s:DEF_SETTE_DIR:${SETTE_DIR}:" -e"s:DEF_INPUT_DIR:${INPUT_DIR}:" \
@@ -321,4 +326,8 @@
               XC40_METO*)
                     cat run_sette_test.job | sed -e"s/SELECT/${SELECT}/" > run_sette_test1.job
+                    mv run_sette_test1.job run_sette_test.job
+                    ;;
+              ifort_athena_*)
+                    cat run_sette_test.job | sed -e"s/NPROC_NODE/${NB_PROC_NODE}/" > run_sette_test1.job
                     mv run_sette_test1.job run_sette_test.job
                     ;;
Index: /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/sette.sh
===================================================================
--- /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/sette.sh	(revision 8278)
+++ /branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/SETTE/sette.sh	(revision 8279)
@@ -36,4 +36,5 @@
 #                      "yes" to run in MPMD (detached) mode with stand-alone IO servers
 #                      "no"  to run in SPMD (attached) mode without separate IO servers 
+# USING_OMP         : flag to control the use of OpenMP parallelization
 # NUM_XIOSERVERS    : number of stand-alone IO servers to employ
 #                     set to zero if USING_MPMD="no"
@@ -88,11 +89,13 @@
 #
 # Compiler among those in NEMOGCM/ARCH
-COMPILER=X64_ADA
-
-export BATCH_COMMAND_PAR="llsubmit"
+#
+COMPILER=ifort_athena_xios_omp
+
+export BATCH_COMMAND_PAR="bsub < "
 export BATCH_COMMAND_SEQ=$BATCH_COMMAND_PAR
 export INTERACT_FLAG="no"
 export MPIRUN_FLAG="yes"
 export USING_XIOS="yes"
+export USING_OMP="yes"
 #
 export DEL_KEYS="key_iomput"
@@ -119,4 +122,14 @@
    echo "Incompatible choices. MPMD mode requires the XIOS server"
    exit
+fi
+#
+# Settings which control the hybrid parallel execution
+#
+OMP_NTHR=1
+ADD_OMP_KEY=""
+if [ ${USING_OMP} == "yes" ]
+ then
+   OMP_NTHR=4
+   ADD_OMP_KEY="key_omp"
 fi
 
@@ -152,5 +165,10 @@
     export TEST_NAME="LONG"
     cd ${CONFIG_DIR0}
-    . ./makenemo -m ${CMP_NAM} -n GYREPIS_LONG -r GYRE_PISCES -j 8 del_key ${DEL_KEYS}
+    if [ ${USING_OMP} == "yes" ]
+    then
+       . ./makenemo -m ${CMP_NAM} -n GYREPIS_LONG_omp -r GYRE_PISCES -j 8 add_key ${ADD_OMP_KEY} del_key ${DEL_KEYS}
+    else
+       . ./makenemo -m ${CMP_NAM} -n GYREPIS_LONG_mpi -r GYRE_PISCES -j 8 del_key ${DEL_KEYS}
+    fi
     cd ${SETTE_DIR}
     . ./param.cfg
@@ -161,4 +179,6 @@
     if [ -f ${JOB_FILE} ] ; then \rm ${JOB_FILE} ; fi
     cd ${EXE_DIR}  
+    set_namelist namelist_cfg nn_GYRE 12
+    set_namelist namelist_cfg ln_bench .true.
     set_namelist namelist_cfg cn_exp \"GYREPIS_LONG\"
     set_namelist namelist_cfg nn_it000 1
@@ -169,14 +189,15 @@
     set_namelist namelist_cfg ln_read_cfg .false.
     set_namelist namelist_cfg nn_fwb 0
-    set_namelist namelist_cfg jpni 4
-    set_namelist namelist_cfg jpnj 8
-    set_namelist namelist_cfg jpnij 32
-    if [ ${USING_MPMD} == "yes" ] ; then
-       set_xio_using_server iodef.xml true
-    else
-       set_xio_using_server iodef.xml false
-    fi
-    cd ${SETTE_DIR}
-    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS}
+    set_namelist namelist_cfg jpni 8
+    set_namelist namelist_cfg jpnj 4
+    set_namelist namelist_cfg jpnij 32
+    set_namelist namelist_cfg jtni 2 		# to be set when using OpenMP
+    if [ ${USING_MPMD} == "yes" ] ; then
+       set_xio_using_server iodef.xml true
+    else
+       set_xio_using_server iodef.xml false
+    fi
+    cd ${SETTE_DIR}
+    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS} ${OMP_NTHR}
 
     cd ${SETTE_DIR}
@@ -184,4 +205,6 @@
     . ./prepare_exe_dir.sh
     cd ${EXE_DIR}
+    set_namelist namelist_cfg nn_GYRE 12
+    set_namelist namelist_cfg ln_bench .true.
     set_namelist namelist_cfg cn_exp \"GYREPIS_SHORT\"
     set_namelist namelist_cfg nn_it000 61
@@ -194,7 +217,8 @@
     set_namelist namelist_cfg ln_read_cfg .false.
     set_namelist namelist_cfg nn_fwb 0
-    set_namelist namelist_cfg jpni 4
-    set_namelist namelist_cfg jpnj 8
-    set_namelist namelist_cfg jpnij 32
+    set_namelist namelist_cfg jpni 8
+    set_namelist namelist_cfg jpnj 4
+    set_namelist namelist_cfg jpnij 32
+    set_namelist namelist_cfg jtni 2 		# to be set when using OpenMP
     set_namelist namelist_top_cfg ln_rsttr .true.
     set_namelist namelist_top_cfg nn_rsttr 2
@@ -218,5 +242,5 @@
     fi
     cd ${SETTE_DIR}
-    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS}
+    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS} ${OMP_NTHR}
     cd ${SETTE_DIR}
     . ./fcm_job.sh $NPROC ${JOB_FILE} ${INTERACT_FLAG} ${MPIRUN_FLAG}
@@ -227,13 +251,20 @@
     export TEST_NAME="REPRO_4_8"
     cd ${CONFIG_DIR0}
-    . ./makenemo -m ${CMP_NAM} -n GYREPIS_32 -r GYRE_PISCES -j 8 del_key ${DEL_KEYS}
-    cd ${SETTE_DIR}
-    . ./param.cfg
-    . ./all_functions.sh
-    . ./prepare_exe_dir.sh
-    JOB_FILE=${EXE_DIR}/run_job.sh
-    NPROC=32
-    if [ -f ${JOB_FILE} ] ; then \rm ${JOB_FILE} ; fi
-    cd ${EXE_DIR}
+    if [ ${USING_OMP} == "yes" ]
+    then
+       . ./makenemo -m ${CMP_NAM} -n GYREPIS_32_omp -r GYRE_PISCES -j 8 add_key ${ADD_OMP_KEY} del_key ${DEL_KEYS}
+    else
+       . ./makenemo -m ${CMP_NAM} -n GYREPIS_32_mpi -r GYRE_PISCES -j 8 del_key ${DEL_KEYS}
+    fi
+    cd ${SETTE_DIR}
+    . ./param.cfg
+    . ./all_functions.sh
+    . ./prepare_exe_dir.sh
+    JOB_FILE=${EXE_DIR}/run_job.sh
+    NPROC=32
+    if [ -f ${JOB_FILE} ] ; then \rm ${JOB_FILE} ; fi
+    cd ${EXE_DIR}
+    set_namelist namelist_cfg nn_GYRE 12
+    set_namelist namelist_cfg ln_bench .true.
     set_namelist namelist_cfg cn_exp \"GYREPIS_48\"
     set_namelist namelist_cfg nn_it000 1
@@ -247,11 +278,12 @@
     set_namelist namelist_cfg jpnj 8
     set_namelist namelist_cfg jpnij 32
-    if [ ${USING_MPMD} == "yes" ] ; then
-       set_xio_using_server iodef.xml true
-    else
-       set_xio_using_server iodef.xml false
-    fi
-    cd ${SETTE_DIR}
-    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS}
+    set_namelist namelist_cfg jtni 2 		# to be set when using OpenMP
+    if [ ${USING_MPMD} == "yes" ] ; then
+       set_xio_using_server iodef.xml true
+    else
+       set_xio_using_server iodef.xml false
+    fi
+    cd ${SETTE_DIR}
+    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS} ${OMP_NTHR}
     cd ${SETTE_DIR}
     . ./fcm_job.sh $NPROC ${JOB_FILE} ${INTERACT_FLAG} ${MPIRUN_FLAG}
@@ -264,4 +296,6 @@
     if [ -f ${JOB_FILE} ] ; then \rm ${JOB_FILE} ; fi
     cd ${EXE_DIR}
+    set_namelist namelist_cfg nn_GYRE 12
+    set_namelist namelist_cfg ln_bench .true.
     set_namelist namelist_cfg cn_exp \"GYREPIS_84\"
     set_namelist namelist_cfg nn_it000 1
@@ -275,11 +309,12 @@
     set_namelist namelist_cfg jpnj 4
     set_namelist namelist_cfg jpnij 32
-    if [ ${USING_MPMD} == "yes" ] ; then
-       set_xio_using_server iodef.xml true
-    else
-       set_xio_using_server iodef.xml false
-    fi
-    cd ${SETTE_DIR}
-    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS}
+    set_namelist namelist_cfg jtni 2 		# to be set when using OpenMP
+    if [ ${USING_MPMD} == "yes" ] ; then
+       set_xio_using_server iodef.xml true
+    else
+       set_xio_using_server iodef.xml false
+    fi
+    cd ${SETTE_DIR}
+    . ./prepare_job.sh input_GYRE.cfg $NPROC ${TEST_NAME} ${MPIRUN_FLAG} ${JOB_FILE} ${NUM_XIOSERVERS} ${OMP_NTHR}
     cd ${SETTE_DIR}
     . ./fcm_job.sh $NPROC ${JOB_FILE} ${INTERACT_FLAG} ${MPIRUN_FLAG}
