Index: /NEMO/trunk/tests/DOME/EXPREF/1_namelist_cfg
===================================================================
--- /NEMO/trunk/tests/DOME/EXPREF/1_namelist_cfg	(revision 14253)
+++ /NEMO/trunk/tests/DOME/EXPREF/1_namelist_cfg	(revision 14254)
@@ -18,4 +18,15 @@
 /
 !-----------------------------------------------------------------------
+&namusr_def    !   User defined :   OVERFLOW configuration
+!-----------------------------------------------------------------------
+   !                       !  type of vertical coordinate
+   ln_zco      = .false.      ! z-coordinate
+   ln_zps      = .true.       ! z-partial-step coordinate
+   ln_sco      = .false.      ! s-coordinate
+   rn_dx       =   5000.   !  horizontal resolution   [meters]
+   rn_dz       =     60.   !  vertical   resolution   [meters]
+   rn_f0       =  1.e-4    !  coriolis [s-1]
+/
+!-----------------------------------------------------------------------
 &namrun        !   parameters of the run
 !-----------------------------------------------------------------------
@@ -105,5 +116,5 @@
 !-----------------------------------------------------------------------
    ln_traadv_fct = .true. !  FCT scheme
-      nn_fct_h   =  4            !  =2/4, horizontal 2nd / 4th order
+      nn_fct_h   =  2            !  =2/4, horizontal 2nd / 4th order
       nn_fct_v   =  2            !  =2/4, vertical   2nd / COMPACT 4th order
 /
Index: /NEMO/trunk/tests/DOME/EXPREF/anim_dome.py
===================================================================
--- /NEMO/trunk/tests/DOME/EXPREF/anim_dome.py	(revision 14254)
+++ /NEMO/trunk/tests/DOME/EXPREF/anim_dome.py	(revision 14254)
@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+""" Draw a plot (or and animation)
+    of passive tracer at the bottom
+    in the DOME experiment
+"""
+
+from netCDF4 import Dataset
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.animation import FuncAnimation
+# from mpl_toolkits.axes_grid1.inset_locator import inset_axes
+#
+# Parent grid data:
+ncid = Dataset('DOME_grid_T.nc')
+lon0 = ncid.variables['nav_lon'][:, :]
+lat0 = ncid.variables['nav_lat'][:, :]
+work = ncid.variables['btra'][:, :, :]
+zos = ncid.variables['zos'][:, :, :]
+ncid.close()
+#
+(jpt0, jpj0, jpi0) = np.shape(work)
+mask = np.where((zos == 0.), True, False)
+tra0 = np.ma.array(work, mask=mask, hard_mask=True)
+
+#
+# Child grid data:
+ncid = Dataset('1_DOME_grid_T.nc')
+lon1 = ncid.variables['nav_lon'][:, :]
+lat1 = ncid.variables['nav_lat'][:, :]
+work = ncid.variables['btra'][:, :, :]
+zos = ncid.variables['zos'][:, :, :]
+sp = ncid.variables['Agrif_sponge'][:, :, :]
+ncid.close()
+#
+(jpt1, jpj1, jpi1) = np.shape(work)
+mask = np.where((zos == 0.), True, False)
+tra1 = np.ma.array(work, mask=mask, hard_mask=True)
+
+# Get resolution in km:
+res0 = np.abs(lon0[0, 1] - lon0[1, 0])
+res1 = np.abs(lon1[0, 1] - lon1[0, 0])
+
+# Shift lon, lat in order to have a pixel centred around the right location:
+lon0 = lon0 - res0/2.
+lat0 = lat0 - res0/2.
+lon1 = lon1 - res1/2.
+lat1 = lat1 - res1/2.
+
+# Indexes to skip ghost zone:
+nghost = 3 # + int(res0/res1)*2 - 1
+imin = nghost + 1
+jmin = nghost + 1
+imax = jpi1 - nghost - 1
+jmax = jpj1 - nghost - 1
+#
+fig, ax = plt.subplots()
+fig.set_tight_layout(True)
+mycmap = plt.cm.nipy_spectral
+plt.set_cmap(mycmap)
+i = 0 
+label = ' Bottom tracer concentration day: {00}'.format(i+1)
+pcol = plt.pcolor(lon0, lat0, np.squeeze(tra0[i, :, :]),
+                  vmin=0.01, vmax=1., edgecolor='0.9', cmap=mycmap)
+pcol = plt.pcolor(lon1[jmin:jmax, imin:imax], lat1[jmin:jmax, imin:imax], np.squeeze(tra1[i, jmin:jmax, imin:imax]),
+                  vmin=0.01, vmax=1., edgecolor='0.4', cmap=mycmap)
+# plt.contour(lon1[jmin:jmax, imin:imax], lat1[jmin:jmax, imin:imax], np.squeeze(sp[i, jmin:jmax, imin:imax]))
+
+plt.axis('scaled')
+# plt.xlim((-800., 250.))
+# plt.ylim((-500.,  50.))
+plt.xlim((-1400., 200.))
+plt.ylim((-300.,  50.))
+plt.gcf().set_size_inches(6, 3)
+
+plt.ylabel('Y (km)', fontsize=14)
+plt.xlabel('X (km)', fontsize=14)
+cbar = plt.colorbar(pcol, shrink=0.5)
+cbar.ax.tick_params(labelsize=14)
+
+plt.title(label)
+
+
+def update(t):
+    txt = ' Bottom tracer concentration day: {00}'.format(t+1)
+    print('Process day', t+1)
+    pcol = plt.pcolor(lon0, lat0, np.squeeze(tra0[t, :, :]),
+                      vmin=0.01, vmax=1., edgecolor='0.9', cmap=mycmap)
+    pcol = plt.pcolor(lon1[jmin:jmax, imin:imax], lat1[jmin:jmax, imin:imax], np.squeeze(tra1[t, jmin:jmax, imin:imax]),
+                      vmin=0.01, vmax=1., edgecolor='0.4', cmap=mycmap)
+    plt.title(txt)
+    return pcol
+
+
+anim = FuncAnimation(fig, update, frames=np.arange(1, 39), interval=150)
+anim.save('DOME_anim_btra.gif', dpi=80, writer='imagemagick')
+# plt.savefig('DOME_bottom_tracer.png')
+# plt.show()
Index: /NEMO/trunk/tests/DOME/EXPREF/file_def_nemo-oce.xml
===================================================================
--- /NEMO/trunk/tests/DOME/EXPREF/file_def_nemo-oce.xml	(revision 14253)
+++ /NEMO/trunk/tests/DOME/EXPREF/file_def_nemo-oce.xml	(revision 14254)
@@ -8,6 +8,6 @@
     -->
     
-    <file_definition type="one_file" name="@expname@" sync_freq="12h" min_digits="4">
-      <file_group id="12h" output_freq="12h"  output_level="10" enabled=".TRUE." >  <!-- 12h files -->  
+    <file_definition type="one_file" name="@expname@" sync_freq="1d" min_digits="4">
+      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE." >  <!-- 1d files -->  
  
 	<file id="file1" name_suffix="_grid_T" description="ocean T grid variables" >
Index: /NEMO/trunk/tests/DOME/EXPREF/namelist_cfg
===================================================================
--- /NEMO/trunk/tests/DOME/EXPREF/namelist_cfg	(revision 14253)
+++ /NEMO/trunk/tests/DOME/EXPREF/namelist_cfg	(revision 14254)
@@ -3,5 +3,15 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !
-!
+!-----------------------------------------------------------------------
+&namusr_def    !   User defined :   OVERFLOW configuration
+!-----------------------------------------------------------------------
+   !                       !  type of vertical coordinate
+   ln_zco      = .false.      ! z-coordinate
+   ln_zps      = .true.       ! z-partial-step coordinate
+   ln_sco      = .false.      ! s-coordinate
+   rn_dx       =   5000.   !  horizontal resolution   [meters]
+   rn_dz       =     60.   !  vertical   resolution   [meters]
+   rn_f0       =  1.e-4    !  coriolis [s-1]
+/
 !-----------------------------------------------------------------------
 &namrun        !   parameters of the run
@@ -92,7 +102,7 @@
 !-----------------------------------------------------------------------
 &namtra_adv    !   advection scheme for tracer
-!-----------------------------------------------------------------------
+!----------------------------------------------------------------------
    ln_traadv_fct = .true. !  FCT scheme
-      nn_fct_h   =  4            !  =2/4, horizontal 2nd / 4th order
+      nn_fct_h   =  2            !  =2/4, horizontal 2nd / 4th order
       nn_fct_v   =  2            !  =2/4, vertical   2nd / COMPACT 4th order
 /
Index: /NEMO/trunk/tests/DOME/MY_SRC/usrdef_hgr.F90
===================================================================
--- /NEMO/trunk/tests/DOME/MY_SRC/usrdef_hgr.F90	(revision 14254)
+++ /NEMO/trunk/tests/DOME/MY_SRC/usrdef_hgr.F90	(revision 14254)
@@ -0,0 +1,127 @@
+MODULE usrdef_hgr
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_hgr  ***
+   !!
+   !!                      ===  DOME configuration  ===
+   !!
+   !! User defined :   mesh and Coriolis parameter of a user configuration
+   !!======================================================================
+   !! History :  NEMO  ! 2017-11  (J. Chanut)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_hgr    : initialize the horizontal mesh for DOME configuration
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE par_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE usrdef_nam, ONLY: rn_dx, rn_dy, rn_f0   ! horizontal resolution in meters
+   !                                             and Coriolis freq. 
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_hgr   ! called by domhgr.F90
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_hgr.F90 13295 2020-07-10 18:24:21Z acc $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_hgr( plamt , plamu , plamv  , plamf  ,   &   ! geographic position (required)
+      &                    pphit , pphiu , pphiv  , pphif  ,   &   !
+      &                    kff   , pff_f , pff_t  ,            &   ! Coriolis parameter  (if domain not on the sphere)
+      &                    pe1t  , pe1u  , pe1v   , pe1f   ,   &   ! scale factors       (required)
+      &                    pe2t  , pe2u  , pe2v   , pe2f   ,   &   !
+      &                    ke1e2u_v      , pe1e2u , pe1e2v     )   ! u- & v-surfaces (if gridsize reduction is used in strait(s))
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE usr_def_hgr  ***
+      !!
+      !! ** Purpose :   user defined mesh and Coriolis parameter
+      !!
+      !! ** Method  :   set all intent(out) argument to a proper value
+      !!                DOME configuration : beta-plance with uniform grid spacing (rn_dx)
+      !!
+      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees) 
+      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane)
+      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
+      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs                     [degrees]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs                      [degrees]
+      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter computed here, =0 otherwise
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point                [1/s]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1t, pe1u, pe1v, pe1f       ! i-scale factors                             [m]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe2t, pe2u, pe2v, pe2f       ! j-scale factors                             [m]
+      INTEGER                 , INTENT(out) ::   ke1e2u_v                     ! =1 u- & v-surfaces computed here, =0 otherwise 
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) ::   zphi0, zlam0 
+      REAL(wp) ::   zti, ztj   ! local scalars
+      !!-------------------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_hgr : DOME configuration bassin'
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '          f-plane with regular grid-spacing'
+      IF(lwp) WRITE(numout,*) '          given by rn_dx and rn_dy' 
+      !
+      !                          
+      ! Position coordinates (in kilometers)
+      !                          ==========
+      zlam0 = -REAL( 0.5 + 1700._wp * 1.e3 / rn_dx)
+      zphi0 = -REAL( 0.5 +  800._wp * 1.e3 / rn_dy)
+#if defined key_agrif
+      IF( .NOT.Agrif_Root() ) THEN 
+         zlam0 = - REAL( 0.5 + 1700._wp * 1.e3 / rn_dx + nbghostcells) &
+               & + REAL((nbghostcells + Agrif_Ix() - 1)*Agrif_irhox())  
+         zphi0 = - REAL( 0.5 +  800._wp * 1.e3 / rn_dy + nbghostcells) &
+               & + REAL((nbghostcells + Agrif_Iy() - 1)*Agrif_irhoy())  
+      ENDIF 
+#endif
+         
+      DO_2D( 1, 1, 1, 1 )
+         zti = REAL( mig0_oldcmp(ji) - 1, wp )   ! start at i=0 in the global grid without halos
+         ztj = REAL( mjg0_oldcmp(jj) - 1, wp )   ! start at j=0 in the global grid without halos
+         
+         plamt(ji,jj) = rn_dx * 1.e-3 * ( zlam0 + zti )
+         plamu(ji,jj) = rn_dx * 1.e-3 * ( zlam0 + zti + 0.5_wp )
+         plamv(ji,jj) = plamt(ji,jj) 
+         plamf(ji,jj) = plamu(ji,jj) 
+         
+         pphit(ji,jj) = rn_dy * 1.e-3 * ( zphi0 + ztj )
+         pphiv(ji,jj) = rn_dy * 1.e-3 * ( zphi0 + ztj + 0.5_wp )
+         pphiu(ji,jj) = pphit(ji,jj) 
+         pphif(ji,jj) = pphiv(ji,jj) 
+      END_2D
+      !     
+      ! Horizontal scale factors (in meters)
+      !                              ======
+      pe1t(:,:) = rn_dx  ;   pe2t(:,:) = rn_dy 
+      pe1u(:,:) = rn_dx  ;   pe2u(:,:) = rn_dy 
+      pe1v(:,:) = rn_dx  ;   pe2v(:,:) = rn_dy 
+      pe1f(:,:) = rn_dx  ;   pe2f(:,:) = rn_dy 
+
+      !                             ! NO reduction of grid size in some straits 
+      ke1e2u_v = 0                  !    ==>> u_ & v_surfaces will be computed in dom_hgr routine
+      pe1e2u(:,:) = 0._wp           !    CAUTION: set to zero to avoid error with some compilers that
+      pe1e2v(:,:) = 0._wp           !             require an initialization of INTENT(out) arguments
+      !
+      !
+      !                       !==  Coriolis parameter  ==!
+      kff = 1                       !  indicate not to compute Coriolis parameter afterward
+      !
+      pff_f(:,:) = rn_f0
+      pff_t(:,:) = rn_f0
+      !
+   END SUBROUTINE usr_def_hgr
+
+   !!======================================================================
+END MODULE usrdef_hgr
Index: /NEMO/trunk/tests/DOME/MY_SRC/usrdef_istate.F90
===================================================================
--- /NEMO/trunk/tests/DOME/MY_SRC/usrdef_istate.F90	(revision 14253)
+++ /NEMO/trunk/tests/DOME/MY_SRC/usrdef_istate.F90	(revision 14254)
@@ -16,4 +16,5 @@
    USE phycst         ! physical constants
    USE eosbn2, ONLY: rn_a0
+   USE dom_oce
    !
    USE in_out_manager ! I/O manager
Index: /NEMO/trunk/tests/DOME/MY_SRC/usrdef_nam.F90
===================================================================
--- /NEMO/trunk/tests/DOME/MY_SRC/usrdef_nam.F90	(revision 14254)
+++ /NEMO/trunk/tests/DOME/MY_SRC/usrdef_nam.F90	(revision 14254)
@@ -0,0 +1,127 @@
+MODULE usrdef_nam
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_nam  ***
+   !!
+   !!                      ===  DOME configuration  ===
+   !!
+   !! User defined : set the domain characteristics of a user configuration
+   !!======================================================================
+   !! History :  NEMO ! 2020-12  (J. Chanut)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_nam   : read user defined namelist and set global domain size
+   !!   usr_def_hgr   : initialize the horizontal mesh 
+   !!----------------------------------------------------------------------
+   USE dom_oce
+   USE par_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_nam   ! called by nemogcm.F90
+
+   !                              !!* namusr_def namelist *!!
+   REAL(wp), PUBLIC ::   rn_dx      ! resolution in meters defining the horizontal domain size
+   REAL(wp), PUBLIC ::   rn_dy      ! resolution in meters defining the horizontal domain size
+   REAL(wp), PUBLIC ::   rn_dz      ! vertical resolution 
+   REAL(wp), PUBLIC ::   rn_f0      ! Coriolis frequency 
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_nam.F90 14086 2020-12-04 11:37:14Z cetlod $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_nam( cd_cfg, kk_cfg, kpi, kpj, kpk, kperio )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_nam  ***
+      !!                    
+      !! ** Purpose :   read user defined namelist and define the domain size
+      !!
+      !! ** Method  :   read in namusr_def containing all the user specific namelist parameter
+      !!
+      !!                Here DOME configuration
+      !!
+      !! ** input   : - namusr_def namelist found in namelist_cfg
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*)              , INTENT(out) ::   cd_cfg          ! configuration name
+      INTEGER                       , INTENT(out) ::   kk_cfg          ! configuration resolution
+      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes 
+      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c. 
+      !
+      INTEGER ::   ios          ! Local integer
+      REAL(wp)::   zlx, zly, zh ! Local scalars
+      !!
+      NAMELIST/namusr_def/  ln_zco, ln_zps, ln_sco, rn_dx, rn_dz, rn_f0
+      !!----------------------------------------------------------------------
+      !
+      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist' )
+      !
+      rn_dy = rn_dx
+#if defined key_agrif 
+      ! Domain parameters are taken from parent:
+      IF( .NOT. Agrif_Root() ) THEN
+         rn_dx = Agrif_Parent(rn_dx)/Agrif_Rhox()
+         rn_dy = Agrif_Parent(rn_dy)/Agrif_Rhoy()
+         rn_dz = Agrif_Parent(rn_dz)
+         rn_f0 = Agrif_Parent(rn_f0)
+      ENDIF
+#endif
+      !
+      IF(lwm)   WRITE( numond, namusr_def )
+      !
+      cd_cfg = 'DOME'               ! name & resolution (not used)
+      kk_cfg = nINT( rn_dx )
+      !
+      IF( Agrif_Root() ) THEN       ! Global Domain size:  DOME  global domain is  2000 km x 850 Km x 3600 m
+         kpi = NINT( 2000.e3  / rn_dx ) + 2  
+         kpj = NINT(  850.e3  / rn_dy ) + 2 + 1 
+      ELSE                          ! Global Domain size: add nbghostcells + 1 "land" point on each side
+         kpi  = nbcellsx + 2 * ( nbghostcells + 1 )
+         kpj  = nbcellsy + 2 * ( nbghostcells + 1 )
+!!$         kpi  = nbcellsx + nbghostcells_x   + nbghostcells_x   + 2
+!!$         kpj  = nbcellsy + nbghostcells_y_s + nbghostcells_y_n + 2
+      ENDIF
+      kpk = NINT( 3600._wp / rn_dz ) + 1
+      !
+      zlx = (kpi-2)*rn_dx*1.e-3
+      zly = (kpj-2-1)*rn_dy*1.e-3
+      zh  = (kpk-1)*rn_dz
+      !                             ! Set the lateral boundary condition of the global domain
+      kperio = 0                    ! DOME configuration : closed basin
+      !                             ! control print
+      IF(lwp) THEN
+         WRITE(numout,*) '   '
+         WRITE(numout,*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+         WRITE(numout,*) '   Namelist namusr_def : DOME test case'
+         WRITE(numout,*) '      z-coordinate flag                 ln_zco = ', ln_zco
+         WRITE(numout,*) '      z-partial-step coordinate flag    ln_zps = ', ln_zps 
+         WRITE(numout,*) '      s-coordinate flag                 ln_sco = ', ln_sco  
+         WRITE(numout,*) '      horizontal resolution             rn_dx  = ', rn_dx, ' m'
+         WRITE(numout,*) '      vertical resolution               rn_dz  = ', rn_dz, ' m'
+         WRITE(numout,*) '      resulting global domain size :    Ni0glo = ', kpi
+         WRITE(numout,*) '                                        Nj0glo = ', kpj
+         WRITE(numout,*) '                                        jpkglo = ', kpk
+         WRITE(numout,*) '      DOME domain: '
+         WRITE(numout,*) '         LX [km]: ', zlx
+         WRITE(numout,*) '         LY [km]: ', zly
+         WRITE(numout,*) '          H [m] : ', zh
+         WRITE(numout,*) '      Coriolis frequency                rn_f0 = ', rn_f0, ' s-1'
+         WRITE(numout,*) '   '
+         WRITE(numout,*) '   Lateral boundary condition of the global domain'
+         WRITE(numout,*) '      DOME : closed basin            jperio = ', kperio
+      ENDIF
+      !
+   END SUBROUTINE usr_def_nam
+
+   !!======================================================================
+END MODULE usrdef_nam
Index: /NEMO/trunk/tests/DOME/MY_SRC/usrdef_zgr.F90
===================================================================
--- /NEMO/trunk/tests/DOME/MY_SRC/usrdef_zgr.F90	(revision 14254)
+++ /NEMO/trunk/tests/DOME/MY_SRC/usrdef_zgr.F90	(revision 14254)
@@ -0,0 +1,262 @@
+MODULE usrdef_zgr
+   !!======================================================================
+   !!                     ***  MODULE usrdef_zgr  ***
+   !!
+   !!                       ===  DOME case  ===
+   !!
+   !! user defined :  vertical coordinate system of a user configuration
+   !!======================================================================
+   !! History :  4.0  ! 2020-12  (J. Chanut)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_zgr   : user defined vertical coordinate system (required)
+   !!       zgr_z1d   : reference 1D z-coordinate 
+   !!---------------------------------------------------------------------
+   USE oce            ! ocean variables
+   USE dom_oce ,  ONLY: mi0, mi1          ! ocean space and time domain
+   USE dom_oce ,  ONLY: glamt, gphit      ! ocean space and time domain
+   USE usrdef_nam     ! User defined : namelist variables
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! distributed memory computing library
+   USE timing         ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_zgr   ! called by domzgr.F90
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: usrdef_zgr.F90 14053 2020-12-03 13:48:38Z techene $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS             
+
+   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
+      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
+      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
+      &                    pe3t  , pe3u  , pe3v , pe3f ,               &   ! vertical scale factors
+      &                    pe3w  , pe3uw , pe3vw,                      &   !     -      -      -
+      &                    k_top  , k_bot    )                             ! top & bottom ocean level
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE usr_def_zgr  ***
+      !!
+      !! ** Purpose :   User defined the vertical coordinates
+      !!
+      !!----------------------------------------------------------------------
+      LOGICAL                   , INTENT(in   ) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags ( read in namusr_def )
+      LOGICAL                   , INTENT(  out) ::   ld_isfcav                   ! under iceshelf cavity flag
+      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth        [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors 
+      INTEGER , DIMENSION(:,:)  , INTENT(  out) ::   k_top, k_bot                ! first & last ocean level
+      !
+      INTEGER  ::   ji, jj, jk        ! dummy indices
+      INTEGER  ::   ik                ! local integers
+      REAL(wp) ::   zfact, z1_jpkm1   ! local scalar
+      REAL(wp) ::   ze3min            ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj) ::   zht, zhu, zhv, zhf, z2d   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_zgr : DOME configuration (z(ps)- or s-coordinate closed box ocean without cavities)'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      !
+      !
+      ! type of vertical coordinate
+      ! ---------------------------
+      ! already set in usrdef_nam.F90 by reading the namusr_def namelist except for ISF
+      ld_isfcav = .FALSE.
+      !
+      !
+      ! Build the vertical coordinate system
+      ! ------------------------------------
+      !
+      !                       !==  UNmasked meter bathymetry  ==!
+      !
+      zht(:,:) = MAX( 600._wp,  600._wp - gphit(:,:)*1.e3*0.01 )
+      !
+      ! at u/v/f-point: averaging zht
+      zhu(:,:) = 600_wp ; zhv(:,:) = 600_wp ; zhf(:,:) = 600_wp
+      DO ji = 1, jpim1
+         zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) )
+      END DO
+      DO jj = 1, jpjm1
+         zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) )
+      END DO
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            zhf(ji,jj) = 0.25_wp * (  zht(ji,jj  ) + zht(ji+1,jj  ) &
+                       &            + zht(ji,jj+1) + zht(ji+1,jj+1) ) 
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'usrdef_zgr', zhu, 'U', 1.0_wp, zhv, 'V', 1.0_wp, zhf, 'F', 1.0_wp)      
+      !     
+      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
+      !
+      !
+      !                       !==  top masked level bathymetry  ==!  (all coordinates)
+      !
+      ! no ocean cavities : top ocean level is ONE, except over land
+      ! the ocean basin surrounded by land (1 grid-point) set through lbc_lnk call as jperio=0 
+      z2d(:,:) = 1._wp                    ! surface ocean is the 1st level
+      WHERE (gphit(:,:)>0._wp) z2d(:,:) = 0._wp
+      ! Dig inlet:
+      WHERE ((gphit(:,:)>0._wp).AND.(glamt(:,:)>-50._wp).AND.(glamt(:,:)<50._wp)) z2d(:,:) = 1._wp
+      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )        ! closed basin since jperio = 0 (see userdef_nam.F90)
+      k_top(:,:) = NINT( z2d(:,:) )
+      !
+      !                              
+      !
+      IF ( ld_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate)
+         !
+         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask)
+         !
+         !                                !* terrain-following coordinate with e3.(k)=cst)
+         z1_jpkm1 = 1._wp / REAL( jpkm1 , wp)
+         DO jk = 1, jpk
+            pdept(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp )
+            pdepw(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk-1 , wp )          )
+            pe3t (:,:,jk) = zht(:,:) * z1_jpkm1
+            pe3u (:,:,jk) = zhu(:,:) * z1_jpkm1
+            pe3v (:,:,jk) = zhv(:,:) * z1_jpkm1
+            pe3f (:,:,jk) = zhf(:,:) * z1_jpkm1
+            pe3w (:,:,jk) = zht(:,:) * z1_jpkm1
+            pe3uw(:,:,jk) = zhu(:,:) * z1_jpkm1
+            pe3vw(:,:,jk) = zhv(:,:) * z1_jpkm1
+         END DO      
+      ENDIF
+      !
+      !
+      IF ( ld_zco ) THEN      !==  z-coordinate  ==!   (step-like topography)
+         !
+         !                                !* bottom ocean compute from the depth of grid-points
+         k_bot(:,:) = jpkm1 * k_top(:,:)     ! here use k_top as a land mask
+         DO jk = 1, jpkm1
+            WHERE( pdept_1d(jk) < zht(:,:) .AND. zht(:,:) <= pdept_1d(jk+1) )   k_bot(:,:) = jk * k_top(:,:)
+         END DO
+         !                                !* horizontally uniform coordinate (reference z-co everywhere)
+         DO jk = 1, jpk
+            pdept(:,:,jk) = pdept_1d(jk)
+            pdepw(:,:,jk) = pdepw_1d(jk)
+            pe3t (:,:,jk) = pe3t_1d (jk)
+            pe3u (:,:,jk) = pe3t_1d (jk)
+            pe3v (:,:,jk) = pe3t_1d (jk)
+            pe3f (:,:,jk) = pe3t_1d (jk)
+            pe3w (:,:,jk) = pe3w_1d (jk)
+            pe3uw(:,:,jk) = pe3w_1d (jk)
+            pe3vw(:,:,jk) = pe3w_1d (jk)
+         END DO
+      ENDIF
+      !
+      !
+      IF ( ld_zps ) THEN      !==  zps-coordinate  ==!   (partial bottom-steps)
+         !
+         ze3min = 0.1_wp * rn_dz
+         IF(lwp) WRITE(numout,*) '   minimum thickness of the partial cells = 10 % of e3 = ', ze3min
+         !
+         !
+         !                                !* bottom ocean compute from the depth of grid-points
+         k_bot(:,:) = jpkm1
+         DO jk = jpkm1, 1, -1
+            WHERE( zht(:,:) < pdepw_1d(jk) + ze3min )   k_bot(:,:) = jk-1
+         END DO
+         !
+         !                                !* vertical coordinate system
+         DO jk = 1, jpk                      ! initialization to the reference z-coordinate
+            pdept(:,:,jk) = pdept_1d(jk)
+            pdepw(:,:,jk) = pdepw_1d(jk)
+            pe3t (:,:,jk) = pe3t_1d (jk)
+            pe3u (:,:,jk) = pe3t_1d (jk)
+            pe3v (:,:,jk) = pe3t_1d (jk)
+            pe3f (:,:,jk) = pe3t_1d (jk)
+            pe3w (:,:,jk) = pe3w_1d (jk)
+            pe3uw(:,:,jk) = pe3w_1d (jk)
+            pe3vw(:,:,jk) = pe3w_1d (jk)
+         END DO
+         DO_2D( 1, 1, 1, 1 )
+            ik = k_bot(ji,jj)
+            pdepw(ji,jj,ik+1) = MIN( zht(ji,jj) , pdepw_1d(ik+1) )
+            pe3t (ji,jj,ik  ) = pdepw(ji,jj,ik+1) - pdepw(ji,jj,ik)
+            pe3t (ji,jj,ik+1) = pe3t (ji,jj,ik  ) 
+            !
+            pdept(ji,jj,ik  ) = pdepw(ji,jj,ik  ) + pe3t (ji,jj,ik  ) * 0.5_wp
+            pdept(ji,jj,ik+1) = pdepw(ji,jj,ik+1) + pe3t (ji,jj,ik+1) * 0.5_wp
+            pe3w (ji,jj,ik+1) = pdept(ji,jj,ik+1) - pdept(ji,jj,ik)              ! = pe3t (ji,jj,ik  )
+            pe3w (ji,jj,ik  ) = pdept(ji,jj,ik  ) - pdept(ji,jj,ik-1)            ! st caution ik > 1
+         END_2D         
+         !                                   ! bottom scale factors and depth at  U-, V-, UW and VW-points
+         !                                   ! usually Computed as the minimum of neighbooring scale factors
+         pe3u (:,:,:) = pe3t(:,:,:)          ! HERE DOME configuration : 
+         pe3v (:,:,:) = pe3t(:,:,:)          !    e3 increases with i-index and identical with j-index
+         pe3f (:,:,:) = pe3t(:,:,:)          !    so e3 minimum of (i,i+1) points is (i) point
+         pe3uw(:,:,:) = pe3w(:,:,:)          !    in j-direction e3v=e3t and e3f=e3v
+         pe3vw(:,:,:) = pe3w(:,:,:)          !    ==>>  no need of lbc_lnk calls
+         !      
+      ENDIF
+      !
+   END SUBROUTINE usr_def_zgr
+
+
+   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zgr_z1d  ***
+      !!
+      !! ** Purpose :   set the depth of model levels and the resulting 
+      !!      vertical scale factors.
+      !!
+      !! ** Method  :   1D z-coordinate system (use in all type of coordinate)
+      !!       The depth of model levels is set from dep(k), an analytical function:
+      !!                   w-level: depw_1d  = dep(k)
+      !!                   t-level: dept_1d  = dep(k+0.5)
+      !!       The scale factors are the discrete derivative of the depth:
+      !!                   e3w_1d(jk) = dk[ dept_1d ] 
+      !!                   e3t_1d(jk) = dk[ depw_1d ]
+      !!
+      !!            ===    Here constant vertical resolution   ===
+      !!
+      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
+      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:), INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
+      REAL(wp), DIMENSION(:), INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
+      !
+      INTEGER  ::   jk       ! dummy loop indices
+      REAL(wp) ::   zt, zw   ! local scalar
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN                         ! Parameter print
+         WRITE(numout,*)
+         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
+         WRITE(numout,*) '    ~~~~~~~'
+      ENDIF
+      !
+      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
+      ! ----------------------
+      DO jk = 1, jpk
+         zw = REAL( jk , wp )
+         zt = REAL( jk , wp ) + 0.5_wp
+         pdepw_1d(jk) =    rn_dz *   REAL( jk-1 , wp )
+         pdept_1d(jk) =    rn_dz * ( REAL( jk-1 , wp ) + 0.5_wp )
+         pe3w_1d (jk) =    rn_dz
+         pe3t_1d (jk) =    rn_dz
+      END DO
+      !
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
+         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
+         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
+      ENDIF
+      !
+   END SUBROUTINE zgr_z1d
+   
+   !!======================================================================
+END MODULE usrdef_zgr
