From 6d4e21f89a3ceee02ff925aff3c24fb3104b3904 Mon Sep 17 00:00:00 2001
From: Diego Bruciaferri <diego.bruciaferri@metoffice.gov.uk>
Date: Tue, 18 Apr 2023 15:51:12 +0000
Subject: [PATCH] Resolve "Introducing Multi-Envelope quasi-Eulerian vertical
coordinates in DOMAINcfg"
---
arch/UKMO/arch-XC40_METO_IFORT_debug.fcm | 67 +
arch/UKMO/arch-vdi_gfortran.fcm | 73 +
tools/DOMAINcfg/namelist_cfg | 84 +-
tools/DOMAINcfg/namelist_ref | 36 +
tools/DOMAINcfg/pysrc/.gitignore | 5 +
.../MEs_4env_2800_r12_r12-r075-r040_v3.inp | 167 +++
.../pysrc/envelopes/generate_envelopes.py | 316 ++++
tools/DOMAINcfg/pysrc/envelopes/lib.py | 244 +++
.../pysrc/local_area/generate_loc_msk.py | 156 ++
tools/DOMAINcfg/pysrc/local_area/lib.py | 463 ++++++
.../pysrc/local_area/loc_area_nordic_ovf.inp | 23 +
tools/DOMAINcfg/pysrc/yml/pyogcm.yml | 33 +
.../DOMAINcfg/src/agrif_recompute_scales.F90 | 2 +-
tools/DOMAINcfg/src/dom_oce.F90 | 19 +-
tools/DOMAINcfg/src/domain.F90 | 6 +-
tools/DOMAINcfg/src/domisf.F90 | 4 +-
tools/DOMAINcfg/src/domloc.F90 | 311 ++++
tools/DOMAINcfg/src/dommes.F90 | 1330 +++++++++++++++++
tools/DOMAINcfg/src/domwri.F90 | 8 +-
tools/DOMAINcfg/src/domzgr.F90 | 197 +--
tools/DOMAINcfg/src/lib_mpp.F90 | 26 +-
21 files changed, 3390 insertions(+), 180 deletions(-)
create mode 100644 arch/UKMO/arch-XC40_METO_IFORT_debug.fcm
create mode 100644 arch/UKMO/arch-vdi_gfortran.fcm
create mode 100644 tools/DOMAINcfg/pysrc/.gitignore
create mode 100644 tools/DOMAINcfg/pysrc/envelopes/MEs_4env_2800_r12_r12-r075-r040_v3.inp
create mode 100755 tools/DOMAINcfg/pysrc/envelopes/generate_envelopes.py
create mode 100644 tools/DOMAINcfg/pysrc/envelopes/lib.py
create mode 100755 tools/DOMAINcfg/pysrc/local_area/generate_loc_msk.py
create mode 100644 tools/DOMAINcfg/pysrc/local_area/lib.py
create mode 100644 tools/DOMAINcfg/pysrc/local_area/loc_area_nordic_ovf.inp
create mode 100644 tools/DOMAINcfg/pysrc/yml/pyogcm.yml
create mode 100644 tools/DOMAINcfg/src/domloc.F90
create mode 100644 tools/DOMAINcfg/src/dommes.F90
diff --git a/arch/UKMO/arch-XC40_METO_IFORT_debug.fcm b/arch/UKMO/arch-XC40_METO_IFORT_debug.fcm
new file mode 100644
index 00000000..1eb02532
--- /dev/null
+++ b/arch/UKMO/arch-XC40_METO_IFORT_debug.fcm
@@ -0,0 +1,67 @@
+# compiler options for Archer CRAY XC-40 (using crayftn compiler)
+#
+# 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_xios is activated)
+# XIOS_LIB xios library (taken into accound only if key_xios 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 $)
+#
+# This arch file depends on loading XIOS-PrgEnv/2.0/24708
+
+%NCDF_HOME $NETCDF_DIR
+%HDF5_HOME $HDF5_DIR
+%XIOS_HOME $xios_path
+%OASIS_HOME $prism_path
+
+%NCDF_INC -I%NCDF_HOME/include -I%HDF5_HOME/include
+%NCDF_LIB -L%HDF5_HOME/lib -L%NCDF_HOME/lib -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 -lz
+
+%XIOS_INC -I%XIOS_HOME/inc
+%XIOS_LIB -L%XIOS_HOME/lib -lxios -lstdc++
+
+%OASIS_INC -I%OASIS_HOME/build/lib/mct -I%OASIS_HOME/build/lib/psmile.MPI1
+%OASIS_LIBDIR -L%OASIS_HOME/lib
+%OASIS_LIB -lpsmile.MPI1 -lmct -lmpeu -lscrip
+
+%CPP cpp
+%FC ftn
+
+%FCFLAGS -r8 -i4 -init=arrays,snan,huge -traceback -debug all -debug inline-debug-info -g -O0 -fp-model strict -check all,noarg_temp_created -fpe-all0 -ftz -ftrapuv
+%FFLAGS %FCFLAGS
+%LD ftn
+%FPPFLAGS -P -E -traditional-cpp
+%LDFLAGS -hbyteswapio
+%AR ar
+%ARFLAGS -r
+%MK gmake
+
+%USER_INC %NCDF_INC %XIOS_INC
+%USER_LIB %NCDF_LIB %XIOS_LIB
+
+%CC cc
+%CFLAGS -O0
+
diff --git a/arch/UKMO/arch-vdi_gfortran.fcm b/arch/UKMO/arch-vdi_gfortran.fcm
new file mode 100644
index 00000000..306129ff
--- /dev/null
+++ b/arch/UKMO/arch-vdi_gfortran.fcm
@@ -0,0 +1,73 @@
+# generic gfortran compiler options for linux
+#
+# 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_xios is activated)
+# XIOS_LIB xios library (taken into accound only if key_xios 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 $)
+#
+# =======================================================================================
+# The following envirnment was used to compile NEMO with this arch file on a UKMO VDI:
+#
+# module purge
+# module load libraries/gcc/6.4.0
+# module load gcc/6.4.0
+# module load mpi/mpich/3.2.1/gnu/6.4.0
+# module load hdf5/1.8.20/gnu/6.4.0
+# module load netcdf/4.6.1/gnu/6.4.0
+# module load xios/r1063/gnu/6.4.0
+# module list
+# export xios_path=/project/ukmo/rhel7/fortran/opt/gfortran/packages/gnu/6.4.0/xios/r1063
+# =======================================================================================
+
+%NCDF_HOME $NETCDF_DIR
+%HDF5_HOME $HDF5_DIR
+%XIOS_HOME $xios_path
+
+%NCDF_INC -I%NCDF_HOME/include -I%HDF5_HOME/include
+%NCDF_LIB -L%NCDF_HOME/lib -lnetcdff -lnetcdf
+%XIOS_INC -I%XIOS_HOME/inc
+%XIOS_LIB -L%XIOS_HOME/lib -lxios -lstdc++
+
+%OASIS_INC -I%OASIS_HOME/build/lib/mct -I%OASIS_HOME/build/lib/psmile.MPI1
+%OASIS_LIB -L%OASIS_HOME/lib -lpsmile.MPI1 -lmct -lmpeu -lscrip
+
+%CPP cpp -Dkey_nosignedzero
+%FC mpif90 -c -cpp
+%FCFLAGS -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer -ffree-line-length-none
+%FFLAGS %FCFLAGS
+%LD mpif90
+%LDFLAGS
+%FPPFLAGS -P -traditional
+%AR ar
+%ARFLAGS rs
+%MK make
+%USER_INC %XIOS_INC %NCDF_INC
+%USER_LIB %XIOS_LIB %NCDF_LIB
+
+%CC cc
+%CFLAGS -O0
diff --git a/tools/DOMAINcfg/namelist_cfg b/tools/DOMAINcfg/namelist_cfg
index cf41b1a1..24e8ac36 100644
--- a/tools/DOMAINcfg/namelist_cfg
+++ b/tools/DOMAINcfg/namelist_cfg
@@ -14,37 +14,44 @@
!-----------------------------------------------------------------------
&namdom ! space and time domain (bathymetry, mesh, timestep)
!-----------------------------------------------------------------------
- ln_read_cfg = .true.
+ ln_read_cfg = .false.
nn_bathy = 1 ! = 0 compute analyticaly
! = 1 read the bathymetry file
! = 2 compute from external bathymetry
! = 3 compute from parent (if "key_agrif")
nn_interp = 1 ! type of interpolation (nn_bathy =2)
- cn_domcfg = 'ORCA_R2_zps_domcfg.nc'
- cn_topo = 'bathymetry_ORCA12_V3.3.nc' ! external topo file (nn_bathy =2)
- cn_bath = 'Bathymetry' ! topo name in file (nn_bathy =2)
- cn_lon = 'nav_lon' ! lon name in file (nn_bathy =2)
- cn_lat = 'nav_lat' ! lat name in file (nn_bathy =2)
+ cn_domcfg = ''
+ cn_fcoord = 'coordinates.nc' ! external coordinates file (jphgr_msh = 0)
+ cn_topo = 'bathy_meter.nc' ! external topo file (nn_bathy =1)
+ cn_bath = 'Bathymetry' ! topo name in file (nn_bathy =1)
+ cn_lon = 'nav_lon' ! lon name in file (nn_bathy =1)
+ cn_lat = 'nav_lat' ! lat name in file (nn_bathy =1)
rn_scale = 1
rn_bathy = 0. ! value of the bathymetry. if (=0) bottom flat at jpkm1
- jphgr_msh = 0 ! type of horizontal mesh
+ nn_msh = 1 ! create (=1) a mesh file or not (=0)
+ jphgr_msh = 0 ! type of horizontal mesh
ppglam0 = 999999.0 ! longitude of first raw and column T-point (jphgr_msh = 1)
ppgphi0 = 999999.0 ! latitude of first raw and column T-point (jphgr_msh = 1)
ppe1_deg = 999999.0 ! zonal grid-spacing (degrees)
ppe2_deg = 999999.0 ! meridional grid-spacing (degrees)
ppe1_m = 999999.0 ! zonal grid-spacing (degrees)
ppe2_m = 999999.0 ! meridional grid-spacing (degrees)
- ppsur = -4762.96143546300 ! ORCA r4, r2 and r05 coefficients
- ppa0 = 255.58049070440 ! (default coefficients)
- ppa1 = 245.58132232490 !
- ppkth = 21.43336197938 !
- ppacr = 3.0 !
+ ppsur = -3958.951371276829 !
+ ppa0 = 103.9530096000000 !
+ ppa1 = 2.415951269000000 !
+ ppkth = 15.35101370000000 !
+ ppacr = 7.0 !
ppdzmin = 999999. ! Minimum vertical spacing
pphmax = 999999. ! Maximum depth
- ldbletanh = .FALSE. ! Use/do not use double tanf function for vertical coordinates
- ppa2 = 999999. ! Double tanh function parameters
- ppkth2 = 999999. !
- ppacr2 = 999999. !
+ ldbletanh = .TRUE. ! Use/do not use double tanf function for vertical coordinates
+ ppa2 = 100.7609285000000 ! Double tanh function parameters
+ ppkth2 = 48.02989372000000 !
+ ppacr2 = 13.0 !
+ rn_atfp = 0.1
+ rn_e3zps_min= 25.0
+ rn_e3zps_rat= 0.2
+ rn_hmin = -8.0
+ rn_rdt = 1350.0
/
!-----------------------------------------------------------------------
&namcfg ! parameters of the configuration
@@ -57,15 +64,15 @@
! ! if ln_e3_dep = T
ln_dept_mid = .false. ! =T : set T points in the middle of cells
! !
- cp_cfg = "orca" ! name of the configuration
- jp_cfg = 2 ! resolution of the configuration
- jpidta = 180 ! 1st lateral dimension ( >= jpi )
- jpjdta = 148 ! 2nd " " ( >= jpj )
- jpkdta = 31 ! number of levels ( >= jpk )
- Ni0glo = 180 ! 1st dimension of global domain --> i =jpidta
- Nj0glo = 148 ! 2nd - - --> j =jpjdta
- jpkglo = 31
- jperio = 4 ! lateral cond. type (between 0 and 6)
+ cp_cfg = "gulf18" ! name of the configuration
+ jp_cfg = 20 ! resolution of the configuration
+ jpidta = 584 ! 1st lateral dimension ( >= jpi )
+ jpjdta = 432 ! 2nd " " ( >= jpj )
+ jpkdta = 52 ! number of levels ( >= jpk )
+ Ni0glo = 584 ! 1st dimension of global domain --> i =jpidta
+ Nj0glo = 432 ! 2nd - - --> j =jpjdta
+ jpkglo = 52
+ jperio = 0 ! lateral cond. type (between 0 and 6)
ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present
! in netcdf input files, as the start j-row for reading
ln_domclo = .false. ! computation of closed sea masks (see namclo)
@@ -75,8 +82,9 @@
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
ln_zco = .false. ! z-coordinate - full steps
- ln_zps = .true. ! z-coordinate - partial steps
+ ln_zps = .false. ! z-coordinate - partial steps
ln_sco = .false. ! s- or hybrid z-s-coordinate
+ ln_mes = .true. ! Multi-Envelope s-coordinate
ln_isfcav = .false. ! ice shelf cavity (T: see namzgr_isf)
/
!-----------------------------------------------------------------------
@@ -88,8 +96,32 @@
!-----------------------------------------------------------------------
/
!-----------------------------------------------------------------------
+&namzgr_mes ! MEs-coordinate
+!-----------------------------------------------------------------------
+ nn_strt = 2 , 1 , 1 , 1 , 1 ! Stretch. funct.: Madec 1996 (0) or
+ ! Song & Haidvogel 1994 (1) or
+ ! Siddorn & Furner 2012 (2)
+ nn_slev = 35 , 17 , 0 , 0 , 0 ! number of s-lev between env(n-1)
+ ! and env(n)
+ rn_e_hc = 50.0 , 0.0 , 0.0 , 0.0 , 0.0 ! critical depth for transition to
+ ! stretch. coord.
+ rn_e_th = 1.0 , 1.0 , 0.0 , 0.0 , 0.0 ! surf. control param.:
+ ! SH94 or MD96: 0<=th<=20
+ ! SF12: thickness surf. cell
+ rn_e_bb = -0.2 , 0.8 , 0.0 , 0.0 , 0.0 ! bot. control param.:
+ ! SH94 or MD96: 0<=bb<=1
+ ! SF12: offset for calculating Zb
+ rn_e_al = 4.4 , 0.0 , 0.0 , 0.0 , 0.0 ! alpha stretching param with SF12
+ rn_e_ba = 0.024, 0.0 , 0.0 , 0.0 , 0.0 ! SF12 bathymetry scaling factor for
+ ! calculating Zb
+ rn_bot_min = 10.0 ! minimum depth of the ocean bottom (>0) (m)
+ rn_bot_max = 7000.0 ! maximum depth of the ocean bottom (= ocean depth) (>0) (m)
+/
+!-----------------------------------------------------------------------
&namclo ! (closed sea : need ln_domclo = .true. in namcfg)
!-----------------------------------------------------------------------
+ rn_shlat = 0.
+ ln_vorlat = .false.
/
!-----------------------------------------------------------------------
&namlbc ! lateral momentum boundary condition (default: NO selection)
diff --git a/tools/DOMAINcfg/namelist_ref b/tools/DOMAINcfg/namelist_ref
index 24dae48f..2d951124 100644
--- a/tools/DOMAINcfg/namelist_ref
+++ b/tools/DOMAINcfg/namelist_ref
@@ -104,7 +104,16 @@
ln_zco = .false. ! z-coordinate - full steps
ln_zps = .false. ! z-coordinate - partial steps
ln_sco = .false. ! s- or hybrid z-s-coordinate
+ ln_mes = .false. ! Multi-Envelope s-coordinate
ln_isfcav = .false. ! ice shelf cavity (T: see namzgr_isf)
+ ln_loczgr = .false. ! to localise or not the chosen vert. coord. system
+ ! If TRUE, the user must provide:
+ ! 1) a bathy_meter.nc file including
+ ! *) s2z_msk: 2D mask for s(=2), s2z(=1) and z(=0) areas
+ ! *) s2z_wgt: 2D field with 1 in s-areas, distance-based
+ ! weights (0<wgt<1) in s2z-areas and 0 elsewhere.
+ ! 2) a domain_cfg_global.nc file including all the geometrical
+ ! information describing the vertical grid used globally.
/
!-----------------------------------------------------------------------
&namzgr_isf ! isf cavity geometry definition (default: OFF)
@@ -146,6 +155,33 @@
!!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
rn_thetb = 1.0 ! bottom control parameter (0<=thetb<= 1)
/
+!-----------------------------------------------------------------------
+&namzgr_mes ! MEs-coordinate (default F)
+!-----------------------------------------------------------------------
+! ! env. 1 ! env. 2 ! env. 3 ! env. 4 ! env. 5 !
+! ! ! ! ! ! !
+ nn_strt = 1 , 1 , 1 , 0 , 0 ! Stretch. funct.: Madec 1996 (0) or
+ ! Song & Haidvogel 1994 (1) or
+ ! Siddorn & Furner 2012 (2)
+ nn_slev = 11 , 10 , 10 , 0 , 0 ! number of s-lev between env(n-1)
+ ! and env(n)
+ rn_e_hc = 20.0 , 10.0 , 0.0 , 0.0 , 0.0 ! critical depth for transition to
+ ! stretch. coord.
+ rn_e_th = 3.9 , 2.98, 2.8 , 0.0 , 0.0 ! surf. control param.:
+ ! SH94 or MD96: 0<=th<=20
+ ! SF12: thickness surf. cell
+ rn_e_bb = 0.77 , 0.15, 0.8 , 0.0 , 0.0 ! bot. control param.:
+ ! SH94 or MD96: 0<=bb<=1
+ ! SF12: offset for calculating Zb
+ rn_e_al = 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ! alpha stretching param with SF12
+ rn_e_ba = 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ! SF12 bathymetry scaling factor for
+ ! calculating Zb
+
+
+ rn_bot_min = 5.0 ! minimum depth of the ocean bottom (>0) (m)
+ rn_bot_max = 2201.5 ! maximum depth of the ocean bottom (= ocean depth) (>0) (m)
+/
+
!-----------------------------------------------------------------------
&namclo ! (closed sea : need ln_domclo = .true. in namcfg) (default: OFF)
!-----------------------------------------------------------------------
diff --git a/tools/DOMAINcfg/pysrc/.gitignore b/tools/DOMAINcfg/pysrc/.gitignore
new file mode 100644
index 00000000..7526cfde
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/.gitignore
@@ -0,0 +1,5 @@
+*.png
+*.out
+*__pycache__
+*Thumbs.db
+*.nc
diff --git a/tools/DOMAINcfg/pysrc/envelopes/MEs_4env_2800_r12_r12-r075-r040_v3.inp b/tools/DOMAINcfg/pysrc/envelopes/MEs_4env_2800_r12_r12-r075-r040_v3.inp
new file mode 100644
index 00000000..ccfb5758
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/envelopes/MEs_4env_2800_r12_r12-r075-r040_v3.inp
@@ -0,0 +1,167 @@
+# ==========================================================================
+# Configuration inupt file to |
+# |
+# 1. define envelopes for generating a multi-envelope (ME) s-coordinate |
+# system as detailed in |
+# Bruciaferri, D., Shapiro, G.I., Wobus, F. 2018, Ocean Dynamics, |
+# https://doi.org/10.1007/s10236-018-1189-x |
+# |
+# 2. localise the ME s-coordinate within another model domain as detailed |
+# in Bruciaferri et al. in preparation (optional). |
+# |
+# ------------------------------------------------------------------------ |
+# |
+# The ME method relies on n envelopes (i.e., arbitrarily defined depth |
+# surfaces) dividing the ocean model vertical domain into n subzones |
+# D_(i), with 1<=i<=n, each one bounded by envelope He_(i-1) at the top |
+# and envelope He_(i) at the bottom, with He_0 = eta (the free surface). |
+# Then: |
+# |
+# (a) FOR ODD i, the transormation from computational space (MEs-space) |
+# to physical space (depth z-space) is |
+# |
+# z = He_0 + hc*s - C(s)*(He_1 - hc - He_0) |
+# |
+# where the depth z and envelopes are downward positive defined, |
+# -1 <= s <= 0, with s(He_0)=0 and s(He_1)=-1 and C(s) is a |
+# stretching function. |
+# |
+# (b) FOR EVEN i: the transormation from MEs-space to z-space is given by |
+# |
+# z = P3(C(s)) |
+# |
+# where P3 is a 3rd order polynomial whose coefficients are computed |
+# locally requiring monotonicity of the transformation and continuity |
+# of its Jacobian and C(s) is stretching function. |
+# |
+# Three options for stretching are given: |
+# *) nn_strt(i) = 0 : Madec et al 1996 cosh/tanh function |
+# *) nn_strt(i) = 1 : Song and Haidvogel 1994 sinh/tanh function |
+# *) nn_strt(i) = 2 : Siddorn and Furner gamma function |
+# |
+#---------------------------------------------------------------------------
+# |
+# SKETCH of the GEOMETRY OF A MEs-COORDINATE SYSTEM |
+# |
+# 3 envelopes are used in this example, such that |
+# 0 < He1 < He2 < He3 : |
+# |
+# === 1st envelope He1 (the shallowest) |
+# ¬¬¬ 2nd envelope He2 |
+# ___ 3rd envelope he3 (the deepest) |
+# --- W-levels |
+# |
+# D1: W-levels marked as D1 belong to the upper |
+# sub-zone D1: |
+# *) The number of discrete levels is controlled |
+# by the nn_slev(1) namelist parameter. |
+# *) The transormation from computational space |
+# to physical space is (a) |
+# |
+# Depth first W-lev: 0 m (surface) |
+# Depth last W-lev: depth of 1st envelope |
+# |
+# D2: W-levels marked as D2 belong to the second |
+# sub-zone D2: |
+# *) The number of discrete levels is controlled |
+# by the nn_slev(2) namelist parameter. |
+# *) The transormation from computational space |
+# to physical space is (b) |
+# |
+# Depth last W-lev: depth of 2nd envelope |
+# |
+# D3: W-levels marked as D3 belong to the third |
+# sub-zone D3: |
+# *) The number of discrete levels is controlled |
+# by the nn_slev(3) namelist parameter. |
+# *) The transormation from computational space |
+# to physical space is (a) |
+# |
+# Depth last W-lev: depth of 3rd envelope |
+# |
+# |~~~~~~~~~~~~~~~~~~~~D1~~~~~~~~~~~~~~~~~~~ SURFACE |
+# | |
+# |--------------------D1------------------- nn_slev(1)=3 |
+# | |
+# |====================D1=================== ENVELOPE 1 |
+# | |
+# |--------------------D2------------------- |
+# | nn_slev(2)=3 |
+# |--------------------D2------------------- |
+# | |
+# |¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬D2¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ ENVELOPE 2 |
+# | |
+# |--------------------D3------------------- nn_slev(3)=2 |
+# | |
+# |____________________D3___________________ ENVELOPE 3 |
+# z \|/ |
+# |
+# ==========================================================================
+
+# A) IINPUT FILES
+
+# Bathymetry of the domain
+bathyFile = '../local_area/bathymetry.loc_area-nord_ovf.dep2800_sig1_stn9_itr1.nc'
+
+# Horizontal grid of the domain
+hgridFile = 'coordinates_eORCA025.nc'
+
+# domain_cfg.nc or mesh_mask.nc file of the external model
+# we want to localise the new ME grid in (optional)
+zgridFile = "domain_cfg_zps.nc"
+
+# B) ENVELOPES geometrical parameters -------------------------------------
+
+# *) e_min_ofs[i] is the offset to be added to the previous envelope env[i-1]
+# in order to calcute the minimum depth of the new envelope.
+#
+# *) To create a flat envelope env[i] with constant depth e_max_dep[i]:
+#
+# e_min_ofs[i] = "flat", e_max_dep[i] > 0
+#
+# *) To create an envelope env[0] which will generate classical s-levels:
+#
+# e_min_ofs[i] > 0, e_max_dep[i] = "max"
+#
+# *) To create a totally flat envelope env[0] which will generate classical
+# z-levels:
+#
+# e_min_ofs[i] = "flat", e_max_dep[i] = "max"
+#
+
+e_min_ofs = [ "flat", 30., 110., "flat"]
+e_max_dep = [ 10. , 500., 2800., 5800.]
+
+# C) ENVELOPES smoothing parameters ------------------------------------------
+
+# For relaxing to geopotential levels
+# near the equator (default=False)
+e_tap_equ = False
+
+# 1. LOCAL SMOOTHING
+
+# List of lists of velocity files to be used for HPGE aware local smoothing.
+# Use an empty list if you don't want to apply any local smoothing to a
+# particular envelope.
+e_loc_vel = [ [],
+ [],
+ ['r12_r12_maximum_hpge.nc',
+ 'r12_r12-r075_maximum_hpge.nc',
+ 'r12_r12-r075-r040_maximum_hpge.nc',
+ 'r12_r12-r075-r040_v2_maximum_hpge.nc'],
+ [] ]
+e_loc_var = [ [], [], ['max_hpge_3','max_hpge_3','max_hpge_3','max_hpge_3'], [] ]
+# List of max spurious currents that will be used as a threshold
+e_loc_vmx = [ [0.], [0.], [0.05,0.05,0.05,0.05], [0.] ]
+
+# List of max slope parameters rx0 for local smoothing.
+e_loc_rmx = [ [0.], [0.], [0.075,0.04,0.04,0.04], [0.] ]
+
+# List of halo for applying the local smoothing.
+e_loc_hal = [ [0], [0], [0,2,2,2], [0] ]
+
+# 2. GLOBAL SMOOTHING
+
+# List of Maximum slope parameter rx0 of each envelope
+e_glo_rmx = [ 0.0, 0.12, 0.12, 0.0 ]
+
diff --git a/tools/DOMAINcfg/pysrc/envelopes/generate_envelopes.py b/tools/DOMAINcfg/pysrc/envelopes/generate_envelopes.py
new file mode 100755
index 00000000..3b4114f2
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/envelopes/generate_envelopes.py
@@ -0,0 +1,316 @@
+#!/usr/bin/env python
+
+# |------------------------------------------------------------|
+# | This module creates general envelope surfaces to be |
+# | used to generate a Localised Multi-Envelope s-coordinate |
+# | vertical grid. |
+# | |
+# | |
+# | -o#&&*''''?d:>b\_ |
+# | _o/"`'' '',, dMF9MMMMMHo_ |
+# | .o&#' `"MbHMMMMMMMMMMMHo. |
+# | .o"" ' vodM*$&&HMMMMMMMMMM?. |
+# | ,' $M&ood,~'`(&##MMMMMMH\ |
+# | / ,MMMMMMM#b?#bobMMMMHMMML |
+# | & ?MMMMMMMMMMMMMMMMM7MMM$R*Hk |
+# | ?$. :MMMMMMMMMMMMMMMMMMM/HMMM|`*L |
+# | | |MMMMMMMMMMMMMMMMMMMMbMH' T, |
+# | $H#: `*MMMMMMMMMMMMMMMMMMMMb#}' `? |
+# | ]MMH# ""*""""*#MMMMMMMMMMMMM' - |
+# | MMMMMb_ |MMMMMMMMMMMP' : |
+# | HMMMMMMMHo `MMMMMMMMMT . |
+# | ?MMMMMMMMP 9MMMMMMMM} - |
+# | -?MMMMMMM |MMMMMMMMM?,d- ' |
+# | :|MMMMMM- `MMMMMMMT .M|. : |
+# | .9MMM[ &MMMMM*' `' . |
+# | :9MMk `MMM#" - |
+# | &M} ` .- |
+# | `&. . |
+# | `~, . ./ |
+# | . _ .- |
+# | '`--._,dd###pp=""' |
+# | |
+# | |
+# | Author: Diego Bruciaferri |
+# | Date and place: 24-06-2021, Met Office, UK |
+# |------------------------------------------------------------|
+
+
+import sys
+import os
+from os.path import isfile, basename, splitext
+import numpy as np
+from matplotlib import pyplot as plt
+import xarray as xr
+from scipy.ndimage import gaussian_filter
+
+from lib import *
+
+# ==============================================================================
+# 1. Checking for input files
+# ==============================================================================
+
+msg_info('GRID ENVELOPES GENERATOR', main=True)
+
+# Checking input file for envelopes
+if len(sys.argv) != 2:
+ raise TypeError(
+ 'You need to specify the absolute path of the envelopes input file'
+ )
+
+# Reading envelopes infos
+env_file = sys.argv[1]
+msg_info('Reading envelopes parameters ....')
+envInfo = read_envInfo(env_file)
+
+# Reading bathymetry and horizontal grid
+msg_info('Reading bathymetry data ... ')
+ds_bathy = xr.open_dataset(envInfo.bathyFile).squeeze()
+
+msg_info('Reading horizontal grid data ... ')
+ds_grid = xr.open_dataset(envInfo.hgridFile)
+glamt = ds_grid["glamt"].squeeze()
+gphit = ds_grid["gphit"].squeeze()
+
+if "nav_lon" in ds_bathy:
+ ds_bathy = ds_bathy.set_coords(["nav_lon","nav_lat"])
+else:
+ ds_bathy["nav_lon"] = glamt
+ ds_bathy["nav_lat"] = gphit
+ ds_bathy = ds_bathy.set_coords(["nav_lon","nav_lat"])
+
+bathy = ds_bathy["Bathymetry"].fillna(0.).squeeze()
+ds_env = ds_bathy.copy()
+
+# Defining local variables -----------------------------------------------------
+
+num_env = len(envInfo.e_min_ofs)
+ni = glamt.shape[1]
+nj = glamt.shape[0]
+
+e_min_ofs = envInfo.e_min_ofs
+e_max_dep = envInfo.e_max_dep
+e_loc_vel = envInfo.e_loc_vel
+e_loc_var = envInfo.e_loc_var
+e_loc_vmx = envInfo.e_loc_vmx
+e_loc_rmx = envInfo.e_loc_rmx
+e_loc_hal = envInfo.e_loc_hal
+e_glo_rmx = envInfo.e_glo_rmx
+e_tap_equ = envInfo.e_tap_equ
+
+tol = 1.e-7
+
+# Computing LSM
+lsm = xr.where(bathy > 0, 1, 0)
+
+#--------------------------------------------------------------------------------------
+# Cretaing envelopes
+for env in range(num_env):
+
+ msg = 'ENVELOPE ' + str(env) + ': min_ofs = '+ str(e_min_ofs[env]) + \
+ ', max_dep = ' + str(e_max_dep[env]) + ', glo_rmx = ' + str(e_glo_rmx[env])
+ msg_info(msg, main=True)
+
+ env_bathy = bathy.copy()
+
+ # =============================================================
+ # 1. GEOMETRY of the envelope
+ #
+ # e_min_ofs[env] > 0 : offset to be added to surf to comupute
+ # minimum depth of envelope env
+ # e_min_ofs[env] = 'flat': flat envelope with constant depth of
+ # e_max_dep[env] m
+ # e_max_dep[env] = 'max' : e_max_dep[env] = np.amax(bathy)
+
+ msg = '1. Computing the geometry of the envelope'
+ msg_info(msg)
+
+ if env == 0:
+ # For the first envelope, the upper envelope
+ # is the ocean surface at rest
+ surf = xr.full_like(bathy,0.,dtype=np.double)
+ else:
+ # For deeper envelopes, the upper
+ # envelope is the previous envelope
+ surf = ds_env["hbatt_"+str(env)]
+
+ if isinstance(e_max_dep[env], str) and e_max_dep[env] == "max":
+ e_max_dep[env] = np.nanmax(bathy)
+
+ # CASE A: flat envelope
+ if isinstance(e_min_ofs[env], str) and e_min_ofs[env] == "flat":
+ msg = 'Generating a flat envelope'
+ msg_info(msg)
+ hbatt = xr.full_like(bathy,e_max_dep[env],dtype=np.double)
+ # CASE B: general envelope
+ else:
+ msg = 'Generating a general envelope'
+ msg_info(msg)
+ hbatt = calc_zenv(env_bathy, surf, e_min_ofs[env], e_max_dep[env])
+
+ # MEs-coordinates are tapered in vicinity of the Equator
+ if e_tap_equ:
+ if (np.nanmax(gphit) * np.nanmin(gphit)) < 0:
+ ztaper = np.exp( -(gphit/8.)**2. )
+ hbatt = np.nanmax(hbatt) * ztaper + hbatt * (1. - ztaper)
+
+ hbatt.plot()
+ plt.show()
+ # =============================================================
+ # 2. SMOOTHING the envelope
+
+ msg = '2. Smoothing the envelope'
+ msg_info(msg)
+
+ # Computing the MB06 Slope Parameter for the raw envelope
+ rmax_raw = np.amax(calc_rmax(hbatt)*lsm).values
+ msg = ' Slope parameter of raw envelope: rmax = ' + str(rmax_raw)
+ msg_info(msg)
+
+ if len(e_loc_vel[env]) == 0:
+ hbatt_smt = hbatt.copy()
+ else:
+ # LOCAL SMOOTHING
+ env_tmp = hbatt.copy()
+ # 1. Generate 2D map for target rmax
+ msk_pge = env_tmp.copy()*0.
+ trg_pge = env_tmp.copy()*0.
+ for m in range(len(e_loc_vel[env])):
+
+ hal = e_loc_hal[env][m]
+
+ filename = e_loc_vel[env][m]
+ varname = e_loc_var[env][m]
+ ds_vel = xr.open_dataset(filename)
+ hpge = ds_vel[varname].data
+
+ nj = hpge.shape[0]
+ ni = hpge.shape[1]
+ for j in range(nj):
+ for i in range(ni):
+ max_hpge = hpge[j,i]
+ if max_hpge >= e_loc_vmx[env][m]:
+ trg_pge[j-hal:j+hal+1,i-hal:i+hal+1] = msk_pge[j,i] + 1
+
+ msk_pge = trg_pge.copy()
+
+ msk_pge.plot.pcolormesh(add_colorbar=True, add_labels=True, \
+ cbar_kwargs=dict(pad=0.15, shrink=1, \
+ label='TOTAL HPGE smoothing mask'))
+ plt.show()
+
+ msg = ' Total number of points with HPGE > ' + str(e_loc_vmx[env][m]) + ' m/s: ' + str(np.nansum(msk_pge.where(msk_pge>0.)))
+ msg_info(msg,)
+
+ # 2. Local smoothing with Martinho & Batteen 2006
+ for m in range(len(e_loc_rmx[env])):
+
+ r0x = e_loc_rmx[env][m]
+ hal = e_loc_hal[env][m]
+ msg = (" Local smoothing of areas of the raw" +\
+ " envelope where HPGE > " + str(e_loc_vmx[env][m]) + " m/s\n"+\
+ " Envelope : " + str(env) + "\n"+\
+ " Local rmax: " + str(r0x) + "\n"+\
+ " Local halo: " + str(hal) + " cells")
+ msg_info(msg)
+
+ msk_smth = msk_pge.where(msk_pge==(m+1))
+
+ msg = ' Total number of points where we target rmax is ' + str(r0x) + ': ' + str(np.nansum(msk_smth))
+ msg_info(msg,)
+
+ msk_smth.plot.pcolormesh(add_colorbar=True, add_labels=True, \
+ cbar_kwargs=dict(pad=0.15, shrink=1, \
+ label='RMAX-HPGE smoothing mask'))
+ plt.show()
+
+ # smoothing with Martinho & Batteen 2006
+ hbatt_smt = smooth_MB06(env_tmp, r0x)
+
+ # applying smoothing only where HPGE are large
+ WRK = hbatt_smt.data
+ TMP = env_tmp.data
+ WRK[np.isnan(msk_smth)] = TMP[np.isnan(msk_smth)]
+ hbatt_smt.data = WRK
+ env_tmp = hbatt_smt.copy()
+
+ ds_env["msk_pge"+str(env+1)] = msk_pge
+
+ # GLOBAL SMOOTHING
+ if e_glo_rmx[env] > 0:
+
+ msg = (' Global smoothing of the raw envelope')
+ msg_info(msg)
+
+ #da_wrk = hbatt_smt.copy()
+ #env_wrk = np.copy(hbatt_smt.data)
+
+ if env == 0:
+ # for the first envelope, we follow NEMO approach
+ # set first land point adjacent to a wet cell to
+ # min_dep as this needs to be included in smoothing
+ cst_lsm = lsm.rolling({dim: 3 for dim in lsm.dims}, min_periods=2).sum()
+ cst_lsm = cst_lsm.shift({dim: -1 for dim in lsm.dims})
+ cst_lsm = (cst_lsm > 0) & (lsm == 0)
+ da_wrk = hbatt_smt.where(cst_lsm == 0, e_min_ofs[env])
+ else:
+ da_wrk = hbatt_smt.copy()
+
+ # smoothing with Martinho & Batteen 2006
+ hbatt_smt = smooth_MB06(da_wrk, e_glo_rmx[env])
+
+
+ # Computing then MB06 Slope Parameter for the smoothed envelope
+ rmax0_smt = calc_rmax(hbatt_smt)*lsm
+ rmax0_smt.plot.pcolormesh(add_colorbar=True, add_labels=True, \
+ cbar_kwargs=dict(pad=0.15, shrink=1, \
+ label='Slope parameter'))
+ plt.show()
+ rmax_smt = np.amax(rmax0_smt).data
+ msg = ' Slope parameter of smoothed envelope: rmax = ' + str(rmax_smt)
+ msg_info(msg)
+
+ # Saving envelope DataArray
+ ds_env["hbatt_"+str(env+1)] = hbatt_smt
+ ds_env["rmax0_"+str(env+1)] = rmax0_smt
+
+# -------------------------------------------------------------------------------------
+# Setting a localised MEs-coord. system if required
+if "s2z_wgt" in ds_env.variables:
+
+ msg = 'SETTING A LOCALISED MEs-coordinates system'
+ msg_info(msg, main=True)
+
+ # Read weights
+ weights = ds_env["s2z_wgt"]
+
+ # Read distribution of levels of global z-coord. grid
+ dsz = xr.open_dataset(envInfo.zgridFile)
+
+ # Computing vertical levels depth
+ # We use e3{t,z}_1d to avoid z-partial steps
+ e3T = dsz.e3t_1d.broadcast_like(dsz.e3t_0).squeeze()
+ e3W = dsz.e3w_1d.broadcast_like(dsz.e3w_0).squeeze()
+ gdepw_0, gdept_0 = e3_to_dep(e3W, e3T)
+
+ # Creating transitioning deeper envelope
+ env = ds_env["hbatt_"+str(num_env)]
+ lev = gdepw_0[{"z":-1}]
+ wrk = xr.full_like(env,None,dtype=np.double)
+ wrk.data = weights * env.data + (1. - weights) * lev.data
+ ds_env["hbatt_"+str(num_env)].data = wrk.data
+
+# -------------------------------------------------------------------------------------
+# Writing the bathy_meter.nc file
+
+msg = 'WRITING the bathy_meter.nc FILE'
+msg_info(msg)
+
+out_name = splitext(basename(envInfo.bathyFile))[0] + "." + splitext(basename(env_file))[0]
+
+print(out_name)
+
+out_file = out_name + ".nc"
+
+ds_env.to_netcdf(out_file)
+
diff --git a/tools/DOMAINcfg/pysrc/envelopes/lib.py b/tools/DOMAINcfg/pysrc/envelopes/lib.py
new file mode 100644
index 00000000..ab9b0b5c
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/envelopes/lib.py
@@ -0,0 +1,244 @@
+#!/usr/bin/env python
+
+from typing import Tuple
+import numpy as np
+import netCDF4 as nc4
+import xarray as xr
+from xarray import Dataset, DataArray
+
+import matplotlib.pyplot as plt
+
+#=======================================================================================
+def read_envInfo(filepath):
+ '''
+ This function reads the parameters which will be used
+ to generate the model grid envelopes.
+
+ If the file is not in the right format and does not contain
+ all the needed infos an error message is returned.
+
+ filepath: string
+ '''
+
+ import importlib.machinery as imp
+
+ loader = imp.SourceFileLoader(filepath,filepath)
+ envInfo = loader.load_module()
+
+ attr = ['bathyFile', 'hgridFile', 'zgridFile', 'e_min_ofs',
+ 'e_max_dep', 'e_loc_vel', 'e_loc_var', 'e_loc_vmx',
+ 'e_loc_rmx', 'e_loc_hal', 'e_glo_rmx', 'e_tap_equ']
+
+ errmsg = False
+ for a in attr:
+ if not hasattr(envInfo,a):
+ errmsg = True
+ attrerr = a
+ else:
+ if getattr(envInfo,a) == "":
+ errmsg = True
+ attrerr = a
+
+ if errmsg:
+ raise AttributeError(
+ 'Attribute ' + attrerr + ' is missing'
+ )
+
+ return envInfo
+
+#=======================================================================================
+def calc_zenv(bathy, surf, offset, max_dep):
+ """
+ Constructs an envelop bathymetry
+ for the sigma levels definition.
+
+ INPUT:
+ *) bathy: the bathymetry field which will be
+ used to compute the envelope.
+ *) surf: 2D field used as upper surface to model
+ the envelope
+ *) offset: offset to be added to surf to calc
+ the minimum depth of the envelope.
+
+ *) max_depth : maximum depth of the envelope.
+ """
+
+ zenv = bathy.copy()
+ env_up = surf.copy()
+
+ # Set maximum depth of the envelope
+ zenv = np.minimum(zenv, max_dep)
+
+ # Set minimum depth of the envelope
+ env_up += offset
+ zenv = np.maximum(zenv, env_up)
+
+ return zenv
+
+#=======================================================================================
+def calc_rmax(depth: DataArray) -> DataArray:
+ """
+ Calculate rmax: measure of steepness
+ This function returns the slope paramater field
+
+ r = abs(Hb - Ha) / (Ha + Hb)
+
+ where Ha and Hb are the depths of adjacent grid cells (Mellor et al 1998).
+
+ Reference:
+ *) Mellor, Oey & Ezer, J Atm. Oce. Tech. 15(5):1122-1131, 1998.
+
+ Parameters
+ ----------
+ depth: DataArray
+ Bottom depth (units: m).
+
+ Returns
+ -------
+ DataArray
+ 2D slope parameter (units: None)
+
+ Notes
+ -----
+ This function uses a "conservative approach" and rmax is overestimated.
+ rmax at T points is the maximum rmax estimated at any adjacent U/V point.
+ """
+ # Mask land
+ depth = depth.where(depth > 0)
+
+ # Loop over x and y
+ both_rmax = []
+ for dim in depth.dims:
+
+ # Compute rmax
+ rolled = depth.rolling({dim: 2}).construct("window_dim")
+ diff = rolled.diff("window_dim").squeeze("window_dim")
+ rmax = np.abs(diff) / rolled.sum("window_dim")
+
+ # Construct dimension with velocity points adjacent to any T point
+ # We need to shift as we rolled twice
+ rmax = rmax.rolling({dim: 2}).construct("vel_points")
+ rmax = rmax.shift({dim: -1})
+
+ both_rmax.append(rmax)
+
+ # Find maximum rmax at adjacent U/V points
+ rmax = xr.concat(both_rmax, "vel_points")
+ rmax = rmax.max("vel_points", skipna=True)
+
+ # Mask halo points
+ for dim in rmax.dims:
+ rmax[{dim: [0, -1]}] = 0
+
+ return rmax.fillna(0)
+
+#=======================================================================================
+def smooth_MB06(
+ depth: DataArray,
+ rmax: float,
+ tol: float = 1.0e-8,
+ max_iter: int = 10_000,
+) -> DataArray:
+ """
+ Direct iterative method of Martinho and Batteen (2006) consistent
+ with NEMO implementation.
+
+ The algorithm ensures that
+
+ H_ij - H_n
+ ---------- < rmax
+ H_ij + H_n
+
+ where H_ij is the depth at some point (i,j) and H_n is the
+ neighbouring depth in the east, west, south or north direction.
+
+ Reference:
+ *) Martinho & Batteen, Oce. Mod. 13(2):166-175, 2006.
+
+ Parameters
+ ----------
+ depth: DataArray
+ Bottom depth.
+ rmax: float
+ Maximum slope parameter allowed
+ tol: float, default = 1.0e-8
+ Tolerance for the iterative method
+ max_iter: int, default = 10000
+ Maximum number of iterations
+
+ Returns
+ -------
+ DataArray
+ Smooth version of the bottom topography with
+ a maximum slope parameter < rmax.
+ """
+
+ # Set scaling factor used for smoothing
+ zrfact = (1.0 - rmax) / (1.0 + rmax)
+
+ # Initialize envelope bathymetry
+ zenv = depth
+
+ for _ in range(max_iter):
+
+ # Initialize lists of DataArrays to concatenate
+ all_ztmp = []
+ all_zr = []
+ for dim in zenv.dims:
+
+ # Shifted arrays
+ zenv_m1 = zenv.shift({dim: -1})
+ zenv_p1 = zenv.shift({dim: +1})
+
+ # Compute zr
+ zr = (zenv_m1 - zenv) / (zenv_m1 + zenv)
+ zr = zr.where((zenv > 0) & (zenv_m1 > 0), 0)
+ for dim_name in zenv.dims:
+ zr[{dim_name: -1}] = 0
+ all_zr += [zr]
+
+ # Compute ztmp
+ zr_p1 = zr.shift({dim: +1})
+ all_ztmp += [zenv.where(zr <= rmax, zenv_m1 * zrfact)]
+ all_ztmp += [zenv.where(zr_p1 >= -rmax, zenv_p1 * zrfact)]
+
+ # Update envelope bathymetry
+ zenv = xr.concat([zenv] + all_ztmp, "dummy_dim").max("dummy_dim")
+
+ # Check target rmax
+ zr = xr.concat(all_zr, "dummy_dim")
+ #print(np.nanmax(np.abs(zr)))
+ if ((np.abs(zr) - rmax) <= tol).all():
+ return zenv
+
+ raise ValueError(
+ "Iterative method did NOT converge."
+ " You might want to increase the number of iterations and/or the tolerance."
+ " The maximum slope parameter rmax is " + str(np.nanmax(np.abs(zr)))
+ )
+
+#=======================================================================================
+def e3_to_dep(e3W: DataArray, e3T: DataArray) -> Tuple[DataArray, ...]:
+
+ gdepT = xr.full_like(e3T, None, dtype=np.double).rename('gdepT')
+ gdepW = xr.full_like(e3W, None, dtype=np.double).rename('gdepW')
+
+ gdepW[{"z":0}] = 0.0
+ gdepT[{"z":0}] = 0.5 * e3W[{"z":0}]
+ for k in range(1, e3W.sizes["z"]):
+ gdepW[{"z":k}] = gdepW[{"z":k-1}] + e3T[{"z":k-1}]
+ gdepT[{"z":k}] = gdepT[{"z":k-1}] + e3W[{"z":k}]
+
+ return tuple([gdepW, gdepT])
+
+#=======================================================================================
+def msg_info(message,main=False):
+
+ if main:
+ print('')
+ print('='*76)
+ print(' '*11 + message)
+ print('='*76)
+ else:
+ print(message)
+ print('')
diff --git a/tools/DOMAINcfg/pysrc/local_area/generate_loc_msk.py b/tools/DOMAINcfg/pysrc/local_area/generate_loc_msk.py
new file mode 100755
index 00000000..396bd6c8
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/local_area/generate_loc_msk.py
@@ -0,0 +1,156 @@
+#!/usr/bin/env python
+
+# |------------------------------------------------------------|
+# | This module creates general envelope surfaces to be |
+# | used to generate a Localised Multi-Envelope s-coordinate |
+# | vertical grid. |
+# | |
+# | |
+# | -o#&&*''''?d:>b\_ |
+# | _o/"`'' '',, dMF9MMMMMHo_ |
+# | .o&#' `"MbHMMMMMMMMMMMHo. |
+# | .o"" ' vodM*$&&HMMMMMMMMMM?. |
+# | ,' $M&ood,~'`(&##MMMMMMH\ |
+# | / ,MMMMMMM#b?#bobMMMMHMMML |
+# | & ?MMMMMMMMMMMMMMMMM7MMM$R*Hk |
+# | ?$. :MMMMMMMMMMMMMMMMMMM/HMMM|`*L |
+# | | |MMMMMMMMMMMMMMMMMMMMbMH' T, |
+# | $H#: `*MMMMMMMMMMMMMMMMMMMMb#}' `? |
+# | ]MMH# ""*""""*#MMMMMMMMMMMMM' - |
+# | MMMMMb_ |MMMMMMMMMMMP' : |
+# | HMMMMMMMHo `MMMMMMMMMT . |
+# | ?MMMMMMMMP 9MMMMMMMM} - |
+# | -?MMMMMMM |MMMMMMMMM?,d- ' |
+# | :|MMMMMM- `MMMMMMMT .M|. : |
+# | .9MMM[ &MMMMM*' `' . |
+# | :9MMk `MMM#" - |
+# | &M} ` .- |
+# | `&. . |
+# | `~, . ./ |
+# | . _ .- |
+# | '`--._,dd###pp=""' |
+# | |
+# | |
+# | Author: Diego Bruciaferri |
+# | Date and place: 24-05-2022, Met Office, UK |
+# |------------------------------------------------------------|
+
+
+import sys
+import os
+from os.path import isfile, basename, splitext
+import numpy as np
+from matplotlib import pyplot as plt
+import xarray as xr
+from scipy.ndimage import gaussian_filter
+
+from lib import *
+
+# ==============================================================================
+# 1. Checking for input files
+# ==============================================================================
+
+msg_info('LOCALISATION AREA GENERATOR', main=True)
+
+# Checking input file for envelopes
+if len(sys.argv) != 2:
+ raise TypeError(
+ 'You need to specify the absolute path of the localisation input file'
+ )
+
+# Reading local area infos
+loc_file = sys.argv[1]
+msg_info('Reading localisation parameters ....')
+locInfo = read_locInfo(loc_file)
+
+bathyFile = locInfo.bathyFile
+hgridFile = locInfo.hgridFile
+loc_isobt = locInfo.loc_isobt
+loc_area = locInfo.loc_area
+s2z_factr = locInfo.s2z_factr
+s2z_sigma = locInfo.s2z_sigma
+s2z_stenc = locInfo.s2z_stenc
+s2z_itera = locInfo.s2z_itera
+
+tol = 1.e-7
+
+# Reading bathymetry and horizontal grid
+msg_info('Reading bathymetry data ... ')
+ds_bathy = xr.open_dataset(bathyFile).squeeze()
+ds_bathy = ds_bathy.fillna(0.)
+msg_info('Reading horizontal grid data ... ')
+ds_grid = xr.open_dataset(hgridFile)
+glamt = ds_grid["glamt"].squeeze()
+gphit = ds_grid["gphit"].squeeze()
+
+if "nav_lon" in ds_bathy:
+ ds_bathy = ds_bathy.set_coords(["nav_lon","nav_lat"])
+else:
+ ds_bathy["nav_lon"] = glamt
+ ds_bathy["nav_lat"] = gphit
+ ds_bathy = ds_bathy.set_coords(["nav_lon","nav_lat"])
+bathy = ds_bathy["Bathymetry"] #.fillna(0.).squeeze()
+
+# Defining local variables -----------------------------------------------------
+
+ni = glamt.shape[1]
+nj = glamt.shape[0]
+ds_loc = ds_bathy.copy()
+
+# Computing LSM
+lsm = xr.where(bathy > 0, 1, 0)
+
+msg = 'Creating areas for localising a vertical coordinate system ... '
+msg_info(msg)
+
+# Create mask identifying zone where
+# we want local coordinates
+msk_loc = generate_loc_area_novf(bathy, loc_isobt)
+ds_loc["loc_area"] = msk_loc
+msk_loc.plot()
+plt.show()
+
+# Creating the transition zone: gaussian filtering
+# for identifying limits of transition zone
+msk_wrk = xr.where(msk_loc==1, 1., 1.+s2z_factr)
+wrk = msk_wrk.copy()
+s2z_trunc = (((s2z_stenc - 1.)/2.)-0.5)/s2z_sigma
+for nit in range(s2z_itera):
+ zwrk_gauss = gaussian_filter(wrk, sigma=s2z_sigma, truncate=s2z_trunc)
+ zwrk_gauss = msk_wrk.where(msk_loc==1, zwrk_gauss)
+ wrk = zwrk_gauss.copy()
+
+# Mask identifying zones
+# global = 0
+# loc to glob = 1
+# local = 2
+msk_zones = xr.full_like(msk_loc, None, dtype=np.double)
+msk = np.ones(shape=msk_zones.values.shape)
+
+diff = np.absolute(zwrk_gauss - msk_wrk)
+values, counts = np.unique(diff, return_counts=True)
+ind = np.argmax(counts)
+print(values[ind]) # prints the most frequent element
+diff = diff.where(diff != values[ind])
+msk[np.isnan(diff)] = 0
+msk[msk_loc==1] = 2
+msk_zones.data = msk
+ds_loc["s2z_msk"] = msk_zones
+
+# 4. Weights to generate transitioning envelope
+msg = 'Computing weights ...'
+msg_info(msg)
+weights = weighting_dist(msk_zones, a=1.7)
+da_wgt = xr.full_like(bathy,None,dtype=np.double)
+da_wgt.data = weights
+ds_loc["s2z_wgt"] = da_wgt
+
+# -------------------------------------------------------------------------------------
+# Writing the bathy_meter.nc file
+
+msg = 'WRITING the bathy_meter.nc FILE'
+msg_info(msg)
+
+out_name = '.dep' + str(int(loc_isobt)) + '_sig' + str(s2z_sigma) + '_stn' + str(s2z_stenc) + '_itr' + str(s2z_itera)
+out_file = "bathymetry.loc_area-" + loc_area + out_name + ".nc"
+ds_loc.to_netcdf(out_file)
diff --git a/tools/DOMAINcfg/pysrc/local_area/lib.py b/tools/DOMAINcfg/pysrc/local_area/lib.py
new file mode 100644
index 00000000..19492fa0
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/local_area/lib.py
@@ -0,0 +1,463 @@
+#!/usr/bin/env python
+
+from typing import Tuple
+import numpy as np
+import netCDF4 as nc4
+import xarray as xr
+from xarray import Dataset, DataArray
+
+#=======================================================================================
+def read_locInfo(filepath):
+ '''
+ This function reads the parameters which will be used
+ to generate the localisation areas.
+
+ If the file is not in the right format and does not contain
+ all the needed infos an error message is returned.
+
+ filepath: string
+ '''
+
+ import importlib.machinery as imp
+
+ loader = imp.SourceFileLoader(filepath,filepath)
+ locInfo = loader.load_module()
+
+ attr = ['bathyFile', 'hgridFile', 'loc_isobt', 'loc_area',
+ 's2z_factr', 's2z_sigma', 's2z_stenc', 's2z_itera' ]
+
+ errmsg = False
+ for a in attr:
+ if not hasattr(locInfo,a):
+ errmsg = True
+ attrerr = a
+ else:
+ if getattr(locInfo,a) == "":
+ errmsg = True
+ attrerr = a
+
+ if errmsg:
+ raise AttributeError(
+ 'Attribute ' + attrerr + ' is missing'
+ )
+
+ return locInfo
+
+#=======================================================================================
+def e3_to_dep(e3W: DataArray, e3T: DataArray) -> Tuple[DataArray, ...]:
+
+ gdepT = xr.full_like(e3T, None, dtype=np.double).rename('gdepT')
+ gdepW = xr.full_like(e3W, None, dtype=np.double).rename('gdepW')
+
+ gdepW[{"z":0}] = 0.0
+ gdepT[{"z":0}] = 0.5 * e3W[{"z":0}]
+ for k in range(1, e3W.sizes["z"]):
+ gdepW[{"z":k}] = gdepW[{"z":k-1}] + e3T[{"z":k-1}]
+ gdepT[{"z":k}] = gdepT[{"z":k-1}] + e3W[{"z":k}]
+
+ return tuple([gdepW, gdepT])
+
+#=======================================================================================
+def msg_info(message,main=False):
+
+ if main:
+ print('')
+ print('='*76)
+ print(' '*11 + message)
+ print('='*76)
+ else:
+ print(message)
+ print('')
+
+#=======================================================================================
+def find_bdy_JI(mask):
+
+ nJ = mask.shape[0]
+ nI = mask.shape[1]
+ JI = []
+ for j in range(1,nJ):
+ for i in range(1,nI):
+ if (mask[j,i] == 0.0) and (np.sum(mask[j-1:j+2,i-1:i+2]) != 0.0):
+ JI.append([j,i])
+ return JI
+
+#=======================================================================================
+def weighting_dist(msk_zones,a):
+ '''
+ msk_zones == 2: s-levels
+ msk_zones == 1: s- to z-levels
+ msk_zones == 0: z-levels
+ '''
+
+ nJ = msk_zones.shape[0]
+ nI = msk_zones.shape[1]
+
+ # Area where we WANT ONLY s-levels
+ msk_inner = np.ones(shape=(nJ,nI))
+ msk_inner[msk_zones==2] = 0
+ # Area where we DO NOT want z-levels
+ msk_outer = np.ones(shape=(nJ,nI))
+ msk_outer[msk_zones==0] = 0
+
+ bdy_inner = find_bdy_JI(msk_inner)
+ bdy_outer = find_bdy_JI(msk_outer)
+
+ w = np.zeros(shape=(nJ,nI))
+ w[msk_outer == 0.0] = 0.0 # where we want only z-levels
+ w[msk_inner == 0.0] = 1.0 # where we want only s-levels
+
+ for j in range(1,nJ):
+ for i in range(1,nI):
+ if msk_outer[j,i]*msk_inner[j,i] != 0.0:
+ d_in = 999999.
+ for kbdy in bdy_inner:
+ jj = kbdy[0]
+ ii = kbdy[1]
+ dist = np.sqrt((j-jj)**2 + (i-ii)**2)
+ if dist < d_in:
+ d_in = dist
+
+ d_out = 999999.
+ for kbdy in bdy_outer:
+ jj = kbdy[0]
+ ii = kbdy[1]
+ dist = np.sqrt((j-jj)**2 + (i-ii)**2)
+ if dist < d_out:
+ d_out = dist
+
+ w[j,i] = 0.5*(np.tanh(a*(2*d_out/(d_in + d_out)-1))/np.tanh(a)+1)
+
+ return w
+
+#=======================================================================================
+def hvrsn_dst(lon1, lat1, lon2, lat2):
+ '''
+ This function calculates the great-circle distance in meters between
+ point1 (lon1,lat1) and point2 (lon2,lat2) using the Haversine formula
+ on a spherical earth of radius 6372.8 km.
+
+ The great-circle distance is the shortest distance over the earth's surface.
+ ( see http://www.movable-type.co.uk/scripts/latlong.html)
+
+ If lon2 and lat2 are 2D matrixes, then dist will be a 2D matrix of distances
+ between all the points in the 2D field and point(lon1,lat1).
+
+ If lon1, lat1, lon2 and lat2 are vectors of size N dist wil be a vector of
+ size N of distances between each pair of points (lon1(i),lat1(i)) and
+ (lon2(i),lat2(i)), with 0 => i > N .
+ '''
+ deg2rad = np.pi / 180.
+ ER = 6372.8 * 1000. # Earth Radius in meters
+
+ dlon = np.multiply(deg2rad, (lon2 - lon1))
+ dlat = np.multiply(deg2rad, (lat2 - lat1))
+
+ lat1 = np.multiply(deg2rad, lat1)
+ lat2 = np.multiply(deg2rad, lat2)
+
+ # Computing the square of half the chord length between the points:
+ a = np.power(np.sin(np.divide(dlat, 2.)),2) + \
+ np.multiply(np.multiply(np.cos(lat1),np.cos(lat2)),np.power(np.sin(np.divide(dlon, 2.)),2))
+
+ # Computing the angular distance in radians between the points
+ angle = np.multiply(2., np.arctan2(np.sqrt(a), np.sqrt(1. -a)))
+
+ # Computing the distance
+ dist = np.multiply(ER, angle)
+
+ return dist
+
+#===================================================================================================
+def get_ij_from_lon_lat(LON, LAT, lon, lat):
+ '''
+ This function finds the closest model
+ grid point i/j to a given lat/lon.
+
+ Syntax:
+ i, j = get_ij_from_lon_lat(LON, LAT, lon, lat)
+
+ LON, LAT: target longitude and latitude
+ lon, lat: 2D arrays of model grid's longitude and latidtude
+ '''
+
+ dist = hvrsn_dst(LON, LAT, lon, lat)
+
+ min_dist = np.amin(dist)
+
+ find_min = np.where(dist == min_dist)
+ sort_j = np.argsort(find_min[1])
+
+ j_indx = find_min[0][sort_j]
+ i_indx = find_min[1][sort_j]
+
+ return j_indx[0], i_indx[0]
+
+# =====================================================================================================
+def bresenham_line(x0, x1, y0, y1):
+ '''
+ point0 = (y0, x0), point1 = (y1, x1)
+
+ It determines the points of an n-dimensional raster that should be
+ selected in order to form a close approximation to a straight line
+ between two points. Taken from the generalised algotihm on
+
+ http://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
+ '''
+
+ steep = abs(y1 - y0) > abs(x1 - x0)
+
+ if steep:
+ # swap(x0, y0)
+ t = y0
+ y0 = x0
+ x0 = t
+ # swap(x1, y1)
+ t = y1
+ y1 = x1
+ x1 = t
+
+ if x0 > x1:
+ # swap(x0, x1)
+ t = x1
+ x1 = x0
+ x0 = t
+ # swap(y0, y1)
+ t = y1
+ y1 = y0
+ y0 = t
+
+ deltax = np.fix(x1 - x0)
+ deltay = np.fix(abs(y1 - y0))
+ error = 0.0
+
+ deltaerr = deltay / deltax
+ y = y0
+
+ if y0 < y1:
+ ystep = 1
+ else:
+ ystep = -1
+
+ c=0
+ pi = np.zeros(shape=[x1-x0+1])
+ pj = np.zeros(shape=[x1-x0+1])
+ for x in np.arange(x0,x1+1) :
+ if steep:
+ pi[c]=y
+ pj[c]=x
+ else:
+ pi[c]=x
+ pj[c]=y
+ error = error + deltaerr
+ if error >= 0.5:
+ y = y + ystep
+ error = error - 1.0
+ c += 1
+
+ return pj, pi
+
+# =====================================================================================================
+
+def get_poly_line_ij(points_i, points_j):
+ '''
+ get_poly_line_ij draw rasterised line between vector-points
+
+ Description:
+ get_poly_line_ij takes a list of points (specified by
+ pairs of indexes i,j) and draws connecting lines between them
+ using the Bresenham line-drawing algorithm.
+
+ Syntax:
+ line_i, line_j = get_poly_line_ij(points_i, points_i)
+
+ Input:
+ points_i, points_j: vectors of equal length of pairs of i, j
+ coordinates that define the line or polyline. The
+ points will be connected in the order they're given
+ in these vectors.
+ Output:
+ line_i, line_j: vectors of the same length as the points-vectors
+ giving the i,j coordinates of the points on the
+ rasterised lines.
+ '''
+ line_i=[]
+ line_j=[]
+
+ line_n=0
+
+ if len(points_i) == 1:
+ line_i = points_i
+ line_j = points_j
+ else:
+ for fi in np.arange(len(points_i)-1):
+ # start point of line
+ i1 = points_i[fi]
+ j1 = points_j[fi]
+ # end point of line
+ i2 = points_i[fi+1]
+ j2 = points_j[fi+1]
+ # 'draw' line from i1,j1 to i2,j2
+ pj, pi = bresenham_line(i1,i2,j1,j2)
+ if pi[0] != i1 or pj[0] != j1:
+ # beginning of line doesn't match end point,
+ # so we flip both vectors
+ pi = np.flipud(pi)
+ pj = np.flipud(pj)
+
+ plen = len(pi)
+
+ for PI in np.arange(plen):
+ line_n = PI
+ if len(line_i) == 0 or line_i[line_n-1] != pi[PI] or line_j[line_n-1] != pj[PI]:
+ line_i.append(int(pi[PI]))
+ line_j.append(int(pj[PI]))
+
+
+ return line_j, line_i
+
+#=======================================================================================
+def floodfill(field,j,i,checkValue,newValue):
+ '''
+ This is a modified version of the original algorithm:
+
+ 1) checkValue is the value we do not want to change,
+ i.e. is the value identifying the boundaries of the
+ region we want to flood.
+ 2) newValue is the new value we want for points whose initial value
+ is not checkValue and is not newValue.
+ N.B. if a point with initial value = to newValue is met, then the
+ flooding stops.
+
+ Example:
+
+ a = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [0, 0, 3, 2, 1, 5, 6, 9, 0],
+ [0, 0, 8, 9, 0, 0, 0, 4, 0],
+ [0, 0, 8, 9, 7, 2, 3, 0, 0],
+ [0, 0, 4, 4, 0, 0, 0, 0, 0],
+ [0, 0, 0, 0, 0, 0, 0, 0, 0]])
+
+ j_start = 3
+ i_start = 4
+ b = com.floodfill(a,j_start,i_start,0,2)
+
+ b = array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [0, 0, 2, 2, 1, 5, 6, 9, 0],
+ [0, 0, 2, 2, 0, 0, 0, 4, 0],
+ [0, 0, 2, 2, 2, 2, 3, 0, 0],
+ [0, 0, 2, 2, 0, 0, 0, 0, 0],
+ [0, 0, 0, 0, 0, 0, 0, 0, 0]])
+
+ '''
+ Field = np.copy(field)
+
+ theStack = [ (j, i) ]
+
+ while len(theStack) > 0:
+ try:
+ j, i = theStack.pop()
+ if Field[j,i] == checkValue:
+ continue
+ if Field[j,i] == newValue:
+ continue
+ Field[j,i] = newValue
+ theStack.append( (j, i + 1) ) # right
+ theStack.append( (j, i - 1) ) # left
+ theStack.append( (j + 1, i) ) # down
+ theStack.append( (j - 1, i) ) # up
+ except IndexError:
+ continue # bounds reached
+
+ return Field
+
+#=======================================================================================
+def get_poly_area_ij(points_i, points_j, a_ji):
+ '''
+ Syntax:
+ area_j, area_i = get_poly_area_ij(points_i, points_j, a_ji)
+
+ Input:
+ points_i, points_j: j,i indexes of the the points
+ defining the polygon
+ a_ji: shape (i.e., (nj,ni)) of the 2d matrix from which points_i
+ and points_j are selected
+ '''
+ [jpj, jpi] = a_ji
+ pnt_i = np.array(points_i)
+ pnt_j = np.array(points_j)
+
+ if (pnt_i[0] == pnt_i[-1]) and (pnt_i[0] == pnt_i[-1]):
+ # polygon is already closed
+ i_indx = np.copy(pnt_i)
+ j_indx = np.copy(pnt_j)
+ else:
+ # close polygon
+ i_indx = np.append( pnt_i, pnt_i[0])
+ j_indx = np.append( pnt_j, pnt_j[0])
+
+ [bound_j, bound_i] = get_poly_line_ij(i_indx, j_indx)
+
+ mask = np.zeros(shape=(jpj,jpi))
+ for n in range(len(bound_i)):
+ i = bound_i[n]
+ j = bound_j[n]
+ mask[j,i] = 1
+ corners_j = [ 0, 0, jpj, jpj ]
+ corners_i = [ 0, jpi, 0, jpi ]
+
+ for n in range(4):
+ j = corners_j[n]
+ i = corners_i[n]
+ if mask[j,i] == 0:
+ fill_i = i
+ fill_j = j
+ break
+
+ mask_filled = floodfill(mask,fill_j,fill_i,1,1)
+ for n in range(len(bound_i)):
+ j = bound_j[n]
+ i = bound_i[n]
+ mask_filled[j,i] = 0
+
+ # the points that are still 0 are within the required area
+ [area_j, area_i] = np.where(mask_filled == 0);
+
+ return area_j, area_i
+
+#=======================================================================================
+def generate_loc_area_novf(bathy, max_dep):
+
+ msk_val = -99
+ s_msk = bathy.where(bathy < max_dep, msk_val)
+
+ # WEST BOUNDARY
+ cond1 = s_msk.nav_lon > -43.
+ cond2 = np.logical_or(s_msk.nav_lon > -38., s_msk.nav_lat > 60.5)
+ s_msk = s_msk.where(np.logical_and(cond1,cond2), msk_val)
+
+ # EAST BOUNDARY
+ cond = s_msk.nav_lon < 0.
+ s_msk = s_msk.where(cond, msk_val)
+ cond1 = np.logical_or(s_msk.nav_lon < -12.5, s_msk.nav_lat > 56.3)
+ s_msk = s_msk.where(cond1, msk_val)
+ cond3 = np.logical_or(s_msk.nav_lon < -5., s_msk.nav_lat > 64)
+ cond4 = bathy > 107.5
+ s_msk = s_msk.where(np.logical_or(cond3,cond4), msk_val)
+
+ # NORTH BOUNDARY
+ cond1 = s_msk.nav_lat < 71.6
+ cond2 = np.logical_or(s_msk.nav_lat < 70., s_msk.nav_lon < -6.)
+ s_msk = s_msk.where(np.logical_and(cond1,cond2), msk_val)
+
+ msk = s_msk.data
+
+ # iceland
+ i_start = 1065
+ j_start = 1010
+
+ msk_fill = floodfill(msk, j_start, i_start, msk_val, -1.)
+ msk_fill *= -1
+ msk_fill[msk_fill != 1] = 0
+ s_msk.data = msk_fill
+
+ return s_msk
diff --git a/tools/DOMAINcfg/pysrc/local_area/loc_area_nordic_ovf.inp b/tools/DOMAINcfg/pysrc/local_area/loc_area_nordic_ovf.inp
new file mode 100644
index 00000000..db5ec53a
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/local_area/loc_area_nordic_ovf.inp
@@ -0,0 +1,23 @@
+# ==========================================================================
+# INPUT FILE TO SETUP A LOCALISED VERTICAL COORDINATE SYSTEM |
+# ==========================================================================
+
+# A) INPUT FILES
+
+# *) Bathymetry of the domain
+bathyFile = 'bathymetry_eORCA025.nc'
+
+# *) Horizontal grid of the domain
+hgridFile = 'coordinates_eORCA025.nc'
+
+# B) LOCALISATION AREA: parameters to identify
+# limits of localisation area
+loc_isobt = 2800.
+loc_area = 'nord_ovf' # stem to be used in the output files name
+
+# c) TRANSITION AREA: parameters for Gaussian filtering
+# to identify limits of transition area
+s2z_factr = 1.e-10
+s2z_sigma = 1
+s2z_stenc = 9
+s2z_itera = 1
diff --git a/tools/DOMAINcfg/pysrc/yml/pyogcm.yml b/tools/DOMAINcfg/pysrc/yml/pyogcm.yml
new file mode 100644
index 00000000..aecc2ab2
--- /dev/null
+++ b/tools/DOMAINcfg/pysrc/yml/pyogcm.yml
@@ -0,0 +1,33 @@
+name: pyogcm
+channels:
+ - conda-forge
+dependencies:
+ - bottleneck
+ - numba
+ - cf_xarray
+ - dask
+ - gsw
+ - netCDF4
+ - pandas
+ - pooch
+ - scipy
+ - xarray
+ - numpy
+ - matplotlib
+ - cmocean
+ - ipython
+ - cartopy
+ - curl
+ - geopy
+ - pip
+ - proj
+ - pyproj
+ - scikit-image
+ - xoak
+ - scikit-learn
+ - pys2index
+ - xgcm
+ - pip:
+ - xnemogcm
+ - git+https://github.com/willirath/xorca.git@master
+ - notebook
diff --git a/tools/DOMAINcfg/src/agrif_recompute_scales.F90 b/tools/DOMAINcfg/src/agrif_recompute_scales.F90
index 2a56d97b..b74fa625 100644
--- a/tools/DOMAINcfg/src/agrif_recompute_scales.F90
+++ b/tools/DOMAINcfg/src/agrif_recompute_scales.F90
@@ -22,7 +22,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: zk ! workspace
- IF ( ln_sco ) RETURN
+ IF ( ln_sco .OR. ln_mes ) RETURN
! Scale factors and depth at U-, V-, UW and VW-points
DO jk = 1, jpk ! initialisation to z-scale factors
diff --git a/tools/DOMAINcfg/src/dom_oce.F90 b/tools/DOMAINcfg/src/dom_oce.F90
index 57ef5d3a..c42e607b 100644
--- a/tools/DOMAINcfg/src/dom_oce.F90
+++ b/tools/DOMAINcfg/src/dom_oce.F90
@@ -213,7 +213,10 @@ MODULE dom_oce
LOGICAL, PUBLIC :: ln_zco !: z-coordinate - full step
LOGICAL, PUBLIC :: ln_zps !: z-coordinate - partial step
LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate
- LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF
+ LOGICAL, PUBLIC :: ln_mes !: Multi-Envelope s-coordinate
+ LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF
+ LOGICAL, PUBLIC :: ln_loczgr !: To localise (.TRUE.) or not (.FALSE.)
+ !: the chosen vertical coordinate system
!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0, e3u_0 , e3v_0 , e3f_0 !: t-,u-,v-,f-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0, e3uw_0, e3vw_0 !: w-,uw-,vw-vert. scale factor [m]
@@ -256,6 +259,13 @@ MODULE dom_oce
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: msk_csgrpglo, msk_csgrprnf, msk_csgrpemp !: closed sea masks
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: tpol, fpol !: north fold mask (jperio= 3 or 4)
+ !
+ ! LOCALISED VERTICAL COORDINATE SYSTEM
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: l2g_msk ! mask to identify areas using
+ ! -) the LOCAL vert. coord. system (=2),
+ ! -) vert. coord. system TRANSITIONING from local to global (=1),
+ ! -) the GLOBAL vert. coord. system (=0)
+ ! If ln_loczgr=.FALSE. then l2g_msk(:,:)=2
!!----------------------------------------------------------------------
!! calendar variables
@@ -366,9 +376,10 @@ CONTAINS
!
ALLOCATE( gdept_1d(jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) ,gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(6) )
!
- ALLOCATE( bathy(jpi,jpj),mbathy(jpi,jpj), tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , &
- & ssmask (jpi,jpj) , ssfmask(jpi,jpj), ssumask(jpi,jpj) , ssvmask(jpi,jpj) , &
- & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj), mbkf(jpi,jpj) , STAT=ierr(7) )
+ ALLOCATE( bathy (jpi,jpj) , mbathy (jpi,jpj) , tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , &
+ & ssmask (jpi,jpj) , ssfmask(jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , &
+ & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , mbkf (jpi,jpj) , &
+ & l2g_msk(jpi,jpj) , STAT=ierr(7) )
!
ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) , &
& risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(8) )
diff --git a/tools/DOMAINcfg/src/domain.F90 b/tools/DOMAINcfg/src/domain.F90
index dfeca554..89cbace6 100644
--- a/tools/DOMAINcfg/src/domain.F90
+++ b/tools/DOMAINcfg/src/domain.F90
@@ -416,9 +416,9 @@ CONTAINS
CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
!
! ! type of vertical coordinate
- IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
- IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
- IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
+ IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
+ IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
+ IF( ln_sco .OR. ln_mes ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
diff --git a/tools/DOMAINcfg/src/domisf.F90 b/tools/DOMAINcfg/src/domisf.F90
index e65f62ef..dbbacaa9 100644
--- a/tools/DOMAINcfg/src/domisf.F90
+++ b/tools/DOMAINcfg/src/domisf.F90
@@ -83,8 +83,8 @@ CONTAINS
!
! 0.1 compatibility option
ierr = 0
- IF ( ln_zco .OR. ln_sco ) ierr = ierr + 1
- IF ( ierr > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
+ IF ( ln_zco .OR. ln_sco .OR. ln_mes ) ierr = ierr + 1
+ IF ( ierr > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco), sigma (ln_sco) or multi-envelope (ln_mes)' )
END SUBROUTINE zgr_isf_nam
diff --git a/tools/DOMAINcfg/src/domloc.F90 b/tools/DOMAINcfg/src/domloc.F90
new file mode 100644
index 00000000..0e88c058
--- /dev/null
+++ b/tools/DOMAINcfg/src/domloc.F90
@@ -0,0 +1,311 @@
+MODULE domloc
+ !!==============================================================================
+ !! *** MODULE domloc ***
+ !! Ocean domain : To localise a vertical coordinate system
+ !!==============================================================================
+ !! History :
+ !! NEMO 4.2 ! 2023-02 (D. Bruciaferri) Original code
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! loc_zgr : localise a vertical coordinate system
+ !!----------------------------------------------------------------------
+ !!
+ !
+ USE dom_oce ! ocean domain
+ USE depth_e3 ! depth <=> e3
+ !
+ USE in_out_manager ! I/O manager
+ USE iom ! I/O library
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
+ USE lib_mpp ! distributed memory computing library
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC loc_zgr ! called by dommes.F90
+ PUBLIC zgr_read ! called by domzgr.F90
+ !
+ ! INPUT GLOBAL VERTICAL COORDINATE SYSTEM
+ !
+ LOGICAL :: ln_zco_in !: z-coordinate - full step
+ LOGICAL :: ln_zps_in !: z-coordinate - partial step
+ LOGICAL :: ln_sco_in !: s-coordinate or hybrid z-s coordinate
+ LOGICAL :: ln_isfcav_in !: presence of ISF
+ INTEGER , ALLOCATABLE, DIMENSION(:,:) :: k_top_in, k_bot_in
+ REAL(wp), ALLOCATABLE, DIMENSION(:) :: e3t_1d_in , e3w_1d_in !: ref. scale factors (m)
+ REAL(wp), ALLOCATABLE, DIMENSION(:) :: gdept_1d_in, gdepw_1d_in !: ref. depth (m)
+ !
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: mbathy_in
+ !
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: e3t_in, e3u_in , e3v_in , e3f_in !: scale factors [m]
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: e3w_in, e3uw_in, e3vw_in !: - -
+ !
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: gdept_in, gdepw_in !: depths [m]
+ !
+ !! * Substitutions
+# include "vectopt_loop_substitute.h90"
+
+CONTAINS
+
+ SUBROUTINE loc_zgr
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE zgr_loc ***
+ !!
+ !! ** Purpose : Wrap the steps to generate a localised vertical
+ !! coordinate systems
+ !!----------------------------------------------------------------------
+ INTEGER :: ji, jj, jk ! dummy loop index
+ INTEGER, DIMENSION(4) :: ierr
+ !-----------------------------------------------------------------------
+ !
+ ! Allocating input grid arrays
+ ALLOCATE(gdept_1d_in(jpk), gdepw_1d_in(jpk), e3t_1d_in (jpk), e3w_1d_in(jpk), STAT=ierr(1))
+ ALLOCATE(k_top_in(jpi,jpj), k_bot_in(jpi,jpj), mbathy_in(jpi,jpj), STAT=ierr(2))
+ ALLOCATE(gdept_in(jpi,jpj,jpk), gdepw_in(jpi,jpj,jpk), e3t_in(jpi,jpj,jpk), e3w_in(jpi,jpj,jpk), STAT=ierr(3))
+ ALLOCATE(e3u_in(jpi,jpj,jpk), e3v_in(jpi,jpj,jpk), e3f_in(jpi,jpj,jpk), e3uw_in(jpi,jpj,jpk), e3vw_in(jpi,jpj,jpk), STAT=ierr(4))
+ IF( MAXVAL(ierr) /= 0 ) CALL ctl_stop( 'STOP', 'loc_zgr: unable to allocate standard ocean arrays' )
+ !
+ ! Reading the input vertical grid that will be used globally
+ ! -----------------------------------------------------------------
+ !
+ CALL zgr_read( ln_zco_in , ln_zps_in , ln_sco_in, ln_isfcav_in, &
+ & gdept_1d_in, gdepw_1d_in, e3t_1d_in, e3w_1d_in , & ! 1D gridpoints depth
+ & gdept_in , gdepw_in , & ! gridpoints depth
+ & e3t_in , e3u_in , e3v_in , e3f_in , & ! vertical scale factors
+ & e3w_in , e3uw_in , e3vw_in , & ! vertical scale factors
+ & k_top_in , k_bot_in )
+ !
+ ! Creating the local coordinate system within the global input vertical grid
+ ! ---------------------------------------------------------------------------
+ !
+ IF ( ln_mes ) CALL loc_zgr_mes
+ ! TO ADD CALL for ln_sco IN THE CASE OF JDHA szt
+ !
+ DEALLOCATE( gdept_1d_in, gdepw_1d_in, e3t_1d_in , e3w_1d_in )
+ DEALLOCATE( k_top_in, k_bot_in, mbathy_in )
+ DEALLOCATE( gdept_in, gdepw_in, e3t_in, e3w_in )
+ DEALLOCATE( e3u_in , e3v_in , e3f_in, e3uw_in, e3vw_in )
+
+ END SUBROUTINE loc_zgr
+
+ SUBROUTINE zgr_read( 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 zgr_read ***
+ !!
+ !! ** Purpose : Read the vertical information in the domain configuration file
+ !!
+ !!----------------------------------------------------------------------
+ LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags
+ 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 vertical scale factors [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 ! - - -
+ INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level
+ !
+ INTEGER :: jk ! dummy loop index
+ INTEGER :: inum ! local logical unit
+ REAL(WP) :: z_zco, z_zps, z_sco, z_cav
+ REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
+ !!----------------------------------------------------------------------
+ !
+ IF(lwp) THEN
+ WRITE(numout,*)
+ WRITE(numout,*) ' zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
+ WRITE(numout,*) ' ~~~~~~~~'
+ ENDIF
+ !
+ CALL iom_open( cn_domcfg, inum )
+ !
+ ! !* type of vertical coordinate
+ CALL iom_get( inum, 'ln_zco' , z_zco )
+ CALL iom_get( inum, 'ln_zps' , z_zps )
+ CALL iom_get( inum, 'ln_sco' , z_sco )
+ IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF
+ IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF
+ IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF
+ !
+ ! !* ocean cavities under iceshelves
+ CALL iom_get( inum, 'ln_isfcav', z_cav )
+ IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF
+ !
+ ! !* vertical scale factors
+ CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate
+ CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d )
+ !
+ CALL iom_get( inum, jpdom_global, 'e3t_0' , pe3t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy ) ! 3D coordinate
+ CALL iom_get( inum, jpdom_global, 'e3u_0' , pe3u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
+ CALL iom_get( inum, jpdom_global, 'e3v_0' , pe3v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
+ CALL iom_get( inum, jpdom_global, 'e3f_0' , pe3f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
+ CALL iom_get( inum, jpdom_global, 'e3w_0' , pe3w , cd_type = 'W', psgn = 1._wp, kfill = jpfillcopy )
+ CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
+ CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
+ !
+ ! !* depths
+ ! !- old depth definition (obsolescent feature)
+ IF( iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0 .AND. &
+ & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0 .AND. &
+ & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0 .AND. &
+ & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0 ) THEN
+ CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', &
+ & ' depths at t- and w-points read in the domain configuration file')
+ CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )
+ CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
+ CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept, kfill = jpfillcopy )
+ CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw, kfill = jpfillcopy )
+ !
+ ELSE !- depths computed from e3. scale factors
+ CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) ! 1D reference depth
+ CALL e3_to_depth( pe3t , pe3w , pdept , pdepw ) ! 3D depths
+ IF(lwp) THEN
+ 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
+ ENDIF
+ !
+ ! !* ocean top and bottom level
+ CALL iom_get( inum, jpdom_global, 'top_level' , z2d ) ! 1st wet T-points (ISF)
+ k_top(:,:) = NINT( z2d(:,:) )
+ CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d ) ! last wet T-points
+ k_bot(:,:) = NINT( z2d(:,:) )
+ !
+ ! reference depth for negative bathy (wetting and drying only)
+ ! IF( ll_wd ) CALL iom_get( inum, 'rn_wd_ref_depth' , ssh_ref )
+ !
+ CALL iom_close( inum )
+ !
+ END SUBROUTINE zgr_read
+
+ SUBROUTINE loc_zgr_mes
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE loc_zgr_mes ***
+ !!
+ !! ** Purpose : Create a local MEs grid within a gloabal grid
+ !! using different vertical coordinates.
+ !!
+ !!----------------------------------------------------------------------
+ INTEGER :: ji, jj, jk ! dummy loop index
+ INTEGER :: inum ! local logical unit
+ REAL :: zwrk
+ REAL(wp), DIMENSION(jpi,jpj) :: l2g_wgt ! weigths for computing model
+ ! levels in transition area
+ !!----------------------------------------------------------------------
+
+ IF(lwp) THEN
+ WRITE(numout,*)
+ WRITE(numout,*) ' loc_zgr_mes: localising ME s-coordinates'
+ WRITE(numout,*) ' ~~~~~~~~'
+ ENDIF
+ !
+ CALL iom_open( 'bathy_meter.nc', inum )
+ CALL iom_get( inum, jpdom_data, 's2z_msk', l2g_msk)
+ CALL iom_get( inum, jpdom_data, 's2z_wgt', l2g_wgt)
+
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ SELECT CASE (INT(l2g_msk(ji,jj)))
+ CASE (0) ! global zps area
+ mbathy (ji,jj ) = mbathy_in(ji,jj)
+ gdept_0(ji,jj,:) = gdept_in(ji,jj,:)
+ gdepw_0(ji,jj,:) = gdepw_in(ji,jj,:)
+ e3t_0 (ji,jj,:) = e3t_in (ji,jj,:)
+ e3w_0 (ji,jj,:) = e3w_in (ji,jj,:)
+ e3u_0 (ji,jj,:) = e3u_in (ji,jj,:)
+ e3v_0 (ji,jj,:) = e3v_in (ji,jj,:)
+ e3f_0 (ji,jj,:) = e3f_in (ji,jj,:)
+ e3uw_0 (ji,jj,:) = e3uw_in (ji,jj,:)
+ e3vw_0 (ji,jj,:) = e3vw_in (ji,jj,:)
+ CASE (1) ! MEs to zps transition area
+ gdept_0(ji,jj,:) = l2g_wgt(ji,jj) * gdept_0(ji,jj,:) + &
+ & ( 1._wp - l2g_wgt(ji,jj) ) * gdept_in(ji,jj,:)
+ gdepw_0(ji,jj,:) = l2g_wgt(ji,jj) * gdepw_0(ji,jj,:) + &
+ & ( 1._wp - l2g_wgt(ji,jj) ) * gdepw_in(ji,jj,:)
+ CASE (2) ! MEs area
+ CYCLE
+ END SELECT
+ END DO
+ END DO
+ !
+ ! e3t, e3w for transition zone
+ ! as finite differences
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ IF ( l2g_msk(ji,jj) == 1._wp ) THEN
+ DO jk = 1,jpkm1
+ e3t_0(ji,jj,jk) = gdepw_0(ji,jj,jk+1) - gdepw_0(ji,jj,jk)
+ e3w_0(ji,jj,jk+1) = gdept_0(ji,jj,jk+1) - gdept_0(ji,jj,jk)
+ ENDDO
+ ! Surface
+ jk = 1
+ e3w_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,1) - gdepw_0(ji,jj,1))
+ !
+ ! Bottom
+ jk = jpk
+ e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk))
+ END IF
+ END DO
+ END DO
+ !
+ ! MBATHY transition zone
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ IF ( l2g_msk(ji,jj) == 1._wp ) THEN
+ DO jk = 1, jpkm1
+ IF( bathy(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
+ END DO
+ END IF
+ END DO
+ END DO
+ !
+ WHERE (bathy(:,:)<=0) mbathy(:,:) = 0
+ !
+ ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0
+ ! for transition zone
+ !
+ DO jj = 1, jpjm1
+ DO ji = 1, jpim1
+ IF ( l2g_msk(ji,jj) == 1._wp ) THEN
+ DO jk = 1, jpk
+
+ zwrk = MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)))
+ e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji ,jj))*e3t_0(ji ,jj,jk) + &
+ MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk)) / zwrk
+
+ zwrk = MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)))
+ e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj ))*e3t_0(ji,jj ,jk) + &
+ MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk)) / zwrk
+
+ zwrk = MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)))
+ e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji ,jj))*e3w_0(ji ,jj,jk) + &
+ MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk)) / zwrk
+
+ zwrk = MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)))
+ e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj ))*e3w_0(ji,jj ,jk) + &
+ MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk)) / zwrk
+
+ zwrk = MAX(1, MIN(1,mbathy(ji ,jj))+MIN(1,mbathy(ji ,jj+1)) &
+ & + MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1)))
+ e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji ,jj ))*e3t_0(ji ,jj ,jk) + &
+ & MIN(1,mbathy(ji+1,jj ))*e3t_0(ji+1,jj ,jk) + &
+ & MIN(1,mbathy(ji+1,jj+1))*e3t_0(ji+1,jj+1,jk) + &
+ & MIN(1,mbathy(ji ,jj+1))*e3t_0(ji ,jj+1,jk)) / &
+ & zwrk
+
+ END DO
+ END IF
+ END DO
+ END DO
+
+ END SUBROUTINE loc_zgr_mes
+
+END MODULE domloc
+
diff --git a/tools/DOMAINcfg/src/dommes.F90 b/tools/DOMAINcfg/src/dommes.F90
new file mode 100644
index 00000000..9d42ea3d
--- /dev/null
+++ b/tools/DOMAINcfg/src/dommes.F90
@@ -0,0 +1,1330 @@
+MODULE dommes
+ !!==============================================================================
+ !! *** MODULE dommes ***
+ !! Ocean domain: Multi-Envelope s-coordinate system (MEs)
+ !!==============================================================================
+ !! History :
+ !! NEMO 4.2 ! 2023-02 (D. Bruciaferri) Original code
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! zgr_mes : build a global or local ME s-coordinate system
+ !! mes_ini : initialise ME vertical grid
+ !! mes_bld : define discrete ME s-levels
+ !!----------------------------------------------------------------------
+ !!
+ USE dom_oce ! ocean domain
+ !
+ USE in_out_manager ! I/O manager
+ USE iom ! I/O library
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
+ USE lib_mpp ! distributed memory computing library
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC zgr_mes ! called by domzgr.F90
+ !
+ CHARACTER(lc) :: ctlmes ! control message error (lc=256)
+ ! * Namelist namzgr_mes *
+ INTEGER :: nn_env ! Number of envelopes used
+ INTEGER, PARAMETER :: max_nn_env = 5 ! Maximum allowed number of envelopes
+ INTEGER :: nn_slev(max_nn_env) ! Array specifing number of levels of each enveloped vertical sub-zone
+ INTEGER :: nn_strt(max_nn_env) ! Array specifing the stretching function
+ ! for each envelope: Madec 1996 (0),
+ ! Song and Haidvogel 1994 (1), Siddorn & Furner 2012 (2)
+ REAL(wp) :: rn_bot_min ! min depth of ocean bottom (>0) (m)
+ REAL(wp) :: rn_bot_max ! max depth of ocean bottom (= ocean depth) (>0) (m)
+ REAL(wp) :: max_dep_env(max_nn_env) ! global maximum depth of envelopes
+ REAL(wp) :: min_dep_env(max_nn_env) ! global minimum depth of envelopes
+ REAL(wp) :: rn_e_hc(max_nn_env) ! Array specifing critical depth for transition to stretched
+ ! coordinates of each envelope
+ REAL(wp) :: rn_e_th(max_nn_env) ! Array specifing surface control parameter (0<=th<=20) of SH94
+ ! or rn_zs of SF12 for each vertical sub-zone
+ REAL(wp) :: rn_e_bb(max_nn_env) ! Array specifing bottom control parameter (0<=bb<=1) of SH94
+ ! or rn_zb_b of SF12 for each vertical sub-zone
+ REAL(wp) :: rn_e_ba(max_nn_env) ! Array specifing bottom control parameter rn_zb_a of SF12 for
+ ! each vertical sub-zone
+ REAL(wp) :: rn_e_al(max_nn_env) ! Array specifing stretching parameter rn_alpha of SF12 for
+ ! each vertical sub-zone
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: envl ! array for envelopes
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3 ! arrays for debugging
+ !! * Substitutions
+# include "do_loop_substitute.h90"
+!
+CONTAINS
+
+ SUBROUTINE zgr_mes
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE zgr_mes ***
+ !!
+ !! ** Purpose : Generate a gloabal or localised
+ !! ME s-coordinates systems
+ !!----------------------------------------------------------------------
+ INTEGER :: ji, jj, jk, je
+ !! Debugging variables
+ INTEGER :: eii, eij, irnk ! dummy loop argument
+ INTEGER, DIMENSION(2) :: e_ij
+ REAL(wp) :: min_dep
+ REAL(wp), DIMENSION(max_nn_env) :: rn_ebot_max
+ REAL(wp), DIMENSION(jpk) :: z_gsigw1, z_gsigt1 ! 1D arrays for debugging
+ REAL(wp), DIMENSION(jpk) :: gdepw1, gdept1 ! 1D arrays for debugging
+ REAL(wp), DIMENSION(jpk) :: e3t1, e3w1 ! 1D arrays for debugging
+ REAL(wp), DIMENSION(jpi,jpj) :: pmsk ! working array
+ !-----------------------------------------------------------------------
+ !
+ ! Initialise ME s-coordinates
+ ! -----------------------------------------------
+ CALL mes_ini
+ !
+ ! Generating a global MEs vertical grid
+ ! -----------------------------------------------
+ CALL mes_bld
+ !
+ ! Some global domain printings
+ ! -----------------------------------------------
+ IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
+ & ' MAX ', MAXVAL( mbathy(:,:) )
+
+ IF( lwp ) THEN ! min max values over the local domain
+ WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
+ WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
+ & ' w ', MINVAL( gdepw_0(:,:,:) )
+ WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
+ & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
+ & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
+ & ' w ', MINVAL( e3w_0 (:,:,:) )
+
+ WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
+ & ' w ', MAXVAL( gdepw_0(:,:,:) )
+ WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
+ & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
+ & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
+ & ' w ', MAXVAL( e3w_0 (:,:,:) )
+ ENDIF
+ !
+ ! Check coordinates makes sense
+ ! -----------------------------------------------
+ ! Extracting MEs depth profile in the shallowest point of the deepest
+ ! envelope for a first check of monotonicity of transformation
+ ! (also useful for debugging)
+ pmsk(:,:) = 0.0
+ WHERE ( bathy > 0.0 ) pmsk = 1.0
+ CALL mpp_minloc( 'zgr_mes', envl(:,:,nn_env), pmsk(:,:), min_dep, e_ij )
+ eii = e_ij(1)
+ eij = e_ij(2)
+ IF ((mi0(eii)>=1 .AND. mi0(eii)<=jpi) .AND. (mj0(eij)>=1 .AND. mj0(eij)<=jpj)) THEN
+ irnk = mpprank
+ DO je = 1, nn_env
+ rn_ebot_max(je) = envl(mi0(eii),mj0(eij),je)
+ END DO
+ z_gsigt1(:) = z_gsigt3(mi0(eii),mj0(eij),:)
+ z_gsigw1(:) = z_gsigw3(mi0(eii),mj0(eij),:)
+ gdept1(:) = gdept_0(mi0(eii),mj0(eij),:)
+ gdepw1(:) = gdepw_0(mi0(eii),mj0(eij),:)
+ e3t1(:) = e3t_0(mi0(eii),mj0(eij),:)
+ e3w1(:) = e3w_0(mi0(eii),mj0(eij),:)
+ ELSE
+ irnk = -1
+ END IF
+ IF( lk_mpp ) CALL mppsync
+ IF( lk_mpp ) CALL mpp_max('zgr_mes', irnk)
+ IF( lk_mpp ) CALL mpp_bcast_real(rn_ebot_max, max_nn_env, irnk )
+ IF( lk_mpp ) CALL mpp_bcast_real(z_gsigt1 , jpk , irnk )
+ IF( lk_mpp ) CALL mpp_bcast_real(z_gsigw1 , jpk , irnk )
+ IF( lk_mpp ) CALL mpp_bcast_real(gdept1 , jpk , irnk )
+ IF( lk_mpp ) CALL mpp_bcast_real(gdepw1 , jpk , irnk )
+ IF( lk_mpp ) CALL mpp_bcast_real(e3t1 , jpk , irnk )
+ IF( lk_mpp ) CALL mpp_bcast_real(e3w1 , jpk , irnk )
+
+ IF( lwp ) THEN
+ WRITE(numout,*) ""
+ WRITE(numout,*) "mes_build:"
+ WRITE(numout,*) "~~~~~~~~~"
+ WRITE(numout,*) ""
+ WRITE(numout,*) " FIRST CHECK: Checking MEs-levels profile in the shallowest point of the last envelope:"
+ WRITE(numout,*) " it is the most likely point where monotonicty of splines may be violeted. "
+ WRITE(numout,*) ""
+ DO je = 1, nn_env
+ WRITE(numout,*) ' * depth of envelope ', je, ' at point (',eii,',',eij,') is ', rn_ebot_max(je)
+ END DO
+ WRITE(numout,*) ""
+ WRITE(numout,*) " MEs-coordinates"
+ WRITE(numout,*) ""
+ WRITE(numout,*) " k z_gsigw1 z_gsigt1"
+ WRITE(numout,*) ""
+ DO jk = 1, jpk
+ WRITE(numout,*) ' ', jk, z_gsigw1(jk), z_gsigt1(jk)
+ ENDDO
+ WRITE(numout,*) ""
+ WRITE(numout,*) "-----------------------------------------------------------------"
+ WRITE(numout,*) ""
+ WRITE(numout,*) " MEs-levels depths and scale factors"
+ WRITE(numout,*) ""
+ WRITE(numout,*) " k gdepw1 e3w1 gdept1 e3t1"
+ WRITE(numout,*) ""
+ DO jk = 1, jpk
+ WRITE(numout,*) ' ', jk, gdepw1(jk), e3w1(jk), gdept1(jk), e3t1(jk)
+ ENDDO
+ WRITE(numout,*) "-----------------------------------------------------------------"
+ ENDIF
+
+ ! Checking monotonicity for cubic splines
+ DO jk = 1, jpk-1
+ IF ( gdept1(jk+1) < gdept1(jk) ) THEN
+ WRITE(ctlmes,*) 'NOT MONOTONIC gdept_0: change envelopes'
+ CALL ctl_stop( ctlmes )
+ END IF
+ IF ( gdepw1(jk+1) < gdepw1(jk) ) THEN
+ WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_0: change envelopes'
+ CALL ctl_stop( ctlmes )
+ END IF
+ IF ( e3t1(jk) < 0.0 ) THEN
+ WRITE(ctlmes,*) 'NEGATIVE e3t_0: change envelopes'
+ CALL ctl_stop( ctlmes )
+ END IF
+ IF ( e3w1(jk) < 0.0 ) THEN
+ WRITE(ctlmes,*) 'NEGATIVE e3w_0: change envelopes'
+ CALL ctl_stop( ctlmes )
+ END IF
+ END DO
+
+ ! CHECKING THE WHOLE DOMAIN
+ DO ji = 1, jpi
+ DO jj = 1, jpj
+ DO jk = 1, jpk !mbathy(ji,jj)
+ ! check coordinate is monotonically increasing
+ IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
+ WRITE(ctmp1,*) 'ERROR mes_build: e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
+ WRITE(numout,*) 'ERROR mes_build: e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
+ WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
+ WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
+ CALL ctl_stop( ctmp1 )
+ ENDIF
+ ! and check it has never gone negative
+ IF ( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
+ WRITE(ctmp1,*) 'ERROR mes_build: gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
+ WRITE(numout,*) 'ERROR mes_build: gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
+ WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
+ WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
+ CALL ctl_stop( ctmp1 )
+ ENDIF
+ END DO
+ END DO
+ END DO
+
+ END SUBROUTINE zgr_mes
+
+ SUBROUTINE mes_bld
+ !!-----------------------------------------------------------------------------
+ !! *** ROUTINE mes_bld ***
+ !!
+ !! ** Purpose : define the Multi-Envelope s-coordinate (MEs) system
+ !!
+ !! ** Method : Generalised s-coordinates are defined with respect to
+ !! multiple arbitrarily defined envelopes as detailed in
+ !!
+ !! Bruciaferri, Shapiro, Wobus, 2018. Oce. Dyn.
+ !! https://doi.org/10.1007/s10236-018-1189-x
+ !!
+ !! The ME method relies on n envelopes (i.e., arbitrarily
+ !! defined depth surfaces) dividing the ocean model vertical
+ !! domain into n subzones D_(i), with 1<=i<=n, each one bounded
+ !! by envelope He_(i-1) at the top and envelope He_(i) at the bottom,
+ !! with He_0 = eta (the free surface). Then:
+ !!
+ !! (a) FOR ODD i, the transormation from computational space
+ !! (MEs-space) to physical space (depth z-space) is
+ !!
+ !! z = He_0 + hc*s - C(s)*(He_1 - hc - He_0)
+ !!
+ !! where the depth z and envelopes are downward positive
+ !! defined, -1 <= s <= 0, with s(He_0)=0 and s(He_1)=-1
+ !! and C(s) is a stretching function.
+ !!
+ !! (b) FOR EVEN i: the transormation from MEs-space to z-space
+ !! is given by
+ !!
+ !! z = P3(C(s))
+ !!
+ !! where P3 is a 3rd order polynomial whose coefficients
+ !! are computed locally requiring monotonicity of the
+ !! transformation and continuity of its Jacobian and C(s)
+ !! is stretching function
+ !!
+ !! Three options for stretching are given:
+ !!
+ !! *) nn_strt(i) = 0 : Madec et al 1996 cosh/tanh function
+ !! *) nn_strt(i) = 1 : Song and Haidvogel 1994 sinh/tanh function
+ !! *) nn_strt(i) = 2 : Siddorn and Furner gamma function
+ !!
+ !!-----------------------------------------------------------------------------------
+ !! SKETCH of the GEOMETRY OF A MEs-COORDINATE SYSTEM
+ !!
+ !! 3 envelopes are used in this example, such that
+ !! 0 < He1 < He2 < He3 :
+ !!
+ !! === 1st envelope He1 (the shallowest)
+ !! ¬¬¬ 2nd envelope He2
+ !! ___ 3rd envelope he3 (the deepest)
+ !! --- W-levels
+ !!
+ !! D1: W-levels marked as D1 belong to the upper
+ !! sub-zone D1:
+ !! *) The number of discrete levels is controlled
+ !! by the nn_slev(1) namelist parameter.
+ !! *) The transormation from computational space
+ !! to physical space is (a)
+ !!
+ !! Depth first W-lev: 0 m (surface)
+ !! Depth last W-lev: depth of 1st envelope
+ !!
+ !! D2: W-levels marked as D2 belong to the second
+ !! sub-zone D2:
+ !! *) The number of discrete levels is controlled
+ !! by the nn_slev(2) namelist parameter.
+ !! *) The transormation from computational space
+ !! to physical space is (b)
+ !!
+ !! Depth last W-lev: depth of 2nd envelope
+ !!
+ !! D3: W-levels marked as D3 belong to the third
+ !! sub-zone D3:
+ !! *) The number of discrete levels is controlled
+ !! by the nn_slev(3) namelist parameter.
+ !! *) The transormation from computational space
+ !! to physical space is (a)
+ !!
+ !! Depth last W-lev: depth of 3rd envelope
+ !!
+ !! |~~~~~~~~~~~~~~~~~~~~D1~~~~~~~~~~~~~~~~~~~ SURFACE
+ !! |
+ !! |--------------------D1------------------- nn_slev(1)=3
+ !! |
+ !! |====================D1=================== ENVELOPE 1
+ !! |
+ !! |--------------------D2-------------------
+ !! | nn_slev(2)=3
+ !! |--------------------D2-------------------
+ !! |
+ !! |¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬D2¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ ENVELOPE 2
+ !! |
+ !! |--------------------D3------------------- nn_slev(3)=2
+ !! |
+ !! |____________________D3___________________ ENVELOPE 3
+ !! z \|/
+ !!----------------------------------------------------------------------
+ !! Coordinate transformation variables
+ INTEGER :: nn_s, jks, jks1 ! for loops over envelopes
+ INTEGER :: nn_s_top, nn_s_bot ! for loops over envelopes
+ REAL(wp) :: coefa, coefb ! for loops over envelopes
+ REAL(wp) :: alpha, env_hc ! for loops over envelopes
+ REAL(wp) :: zcoeft, zcoefw ! temporary scalars
+ REAL(wp) :: max_env_top, min_env_top ! to identify geopotential levels
+ REAL(wp) :: max_env_bot, min_env_bot ! to identify geopotential levels
+ REAL(wp), DIMENSION(jpi,jpj) :: env_top, env_bot ! for loops over envelopes
+ !! Cubic Splines variables
+ REAL(wp) :: dCds_top, dCds_bot ! for loops over envelopes
+ INTEGER :: ibcbeg, ibcend ! boundary condition for cubic splines
+ REAL(wp) :: c(4,2) ! matrix for cubic splines coefficents
+ REAL(wp) :: d(2,2) ! matrix for depth values and depth
+ ! derivative at intervals' boundaries
+ REAL(wp) :: tau(2) ! abscissas of intervals' boundaries
+ REAL(wp), DIMENSION(jpi,jpj) :: env0, env1, env2, env3 ! for loops over envelopes
+ ! for cubic splines
+ !
+ INTEGER :: ji, jj, jk, je ! dummy loop argument
+ !!------------------------------------------------------------------------
+
+ ! Initialise mbathy to the maximum ocean level available
+ mbathy(:,:) = jpkm1
+
+ ! Define scosrf
+ scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
+
+ !============================
+ ! Compute 3D MEs-levels depth
+ !============================
+ jks1 = 2
+ nn_s = 0
+ env_top = scosrf ! surface
+ !
+ DO je = 1, nn_env ! LOOP over all the requested envelopes
+ !
+ jks1 = jks1 + nn_s - 1
+ nn_s = nn_slev(je)
+ if ( je > 1 ) THEN
+ nn_s = nn_s + 1
+ env_top = envl(:,:,je-1)
+ END IF
+ env_bot = envl(:,:,je)
+ !
+ IF ( isodd(je) ) THEN
+ !
+ DO jk = 1, nn_s
+ jks = jks1 + jk - 1
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ !
+ ! -------------------------------------------
+ ! Stretching coefficients
+ !
+ env_hc = rn_e_hc(je)
+ coefa = rn_e_th(je)
+ coefb = rn_e_bb(je)
+ alpha = rn_e_al(je)
+ ! Coefficients for SF12 stretching
+ IF (nn_strt(je) == 2) THEN
+ coefa = rn_e_th(je) / (env_bot(ji,jj)-env_top(ji,jj))
+ coefb = rn_e_bb(je) + ((env_bot(ji,jj)-env_top(ji,jj))*rn_e_ba(je))
+ coefb = 1.0_wp-(coefb/(env_bot(ji,jj)-env_top(ji,jj)))
+ END IF
+ !
+ ! 1) Computing MEs-coordinates (-1 <= MEs <= 0)
+ ! ===============================================
+ !
+ ! SHALLOW WATER, uniform sigma
+ IF (env_bot(ji,jj) < env_hc) THEN
+ zcoefw = -sigma('W', jk, nn_s) ! W-GRID
+ zcoeft = -sigma('T', jk, nn_s) ! T-GRID
+ ! DEEP WATER, stretched s
+ ELSE
+ zcoefw = -stretch( sigma('W',jk,nn_s), coefa, &
+ & coefb, alpha, nn_s, nn_strt(je) ) ! W-GRID
+ zcoeft = -stretch( sigma('T',jk,nn_s), coefa, &
+ & coefb, alpha, nn_s, nn_strt(je) ) ! T-GRID
+ ENDIF
+ !
+ ! For debugging
+ z_gsigw3(ji,jj,jks) = zcoefw
+ z_gsigt3(ji,jj,jks) = zcoeft
+ !
+ ! 2) Computing model levels depths
+ ! ===============================================
+ !
+ IF (nn_strt(je) == 2) env_hc = 0._wp
+ gdepw_0(ji,jj,jks) = z_mes(env_top(ji,jj), env_bot(ji,jj), zcoefw, &
+ & -sigma('W', jk, nn_s), env_hc ) ! W-GRID
+ gdept_0(ji,jj,jks) = z_mes(env_top(ji,jj), env_bot(ji,jj), zcoeft, &
+ & -sigma('T', jk, nn_s), env_hc ) ! T-GRID
+ END DO
+ END DO
+ END DO
+ !
+ ELSE
+ !
+ IF ( je == 2 ) THEN
+ env0 = scosrf
+ ELSE
+ env0 = envl(:,:,je-2)
+ END IF
+ env1 = env_top
+ env2 = env_bot
+ IF ( je < nn_env ) env3 = envl(:,:,je+1)
+ !
+ DO jj = 1,jpj
+ DO ji = 1,jpi
+ d(1,1) = env_top(ji,jj)
+ d(1,2) = env_bot(ji,jj)
+ !
+ ! 1) Computing boundary conditions
+ ! for the 1st or 2nd derivative
+ ! to constrain splines
+ !===================================
+ !
+ ! SUB-ZONE ABOVE (je-1)
+ !
+ ibcbeg = 1 ! Bound. cond for the 1st derivative
+ ! -------------------------------------------
+ ! Stretching coefficients
+ !
+ env_hc = rn_e_hc(je-1)
+ alpha = rn_e_al(je-1)
+ coefa = rn_e_th(je-1)
+ coefb = rn_e_bb(je-1)
+ ! Coefficients for SF12 stretching
+ IF (nn_strt(je-1) == 2) THEN
+ env_hc = 0._wp
+ coefa = rn_e_th(je-1) / (env1(ji,jj)-env0(ji,jj))
+ coefb = rn_e_bb(je-1) + ((env1(ji,jj)-env0(ji,jj))*rn_e_ba(je-1))
+ coefb = 1.0_wp-(coefb/(env1(ji,jj)-env0(ji,jj)))
+ END IF
+ ! -------------------------------------------
+ ! Analytical derivative of coordinate
+ ! transformation at envelope (je-1)
+ !
+ zcoefw = -1._wp
+ dCds_top = dstretch( zcoefw, coefa, coefb, alpha, &
+ & nn_slev(je-1), nn_strt(je-1) )
+ d(2,1) = dzdmes(env0(ji,jj), env1(ji,jj), dCds_top, env_hc)
+ !
+ ! SUB-ZONE BELOW (je+1)
+ !
+ IF ( je < nn_env ) THEN
+ ibcend = 1 ! Bound. cond for the 1st derivative
+ ! -------------------------------------------
+ ! Stretching coefficients
+ !
+ env_hc = rn_e_hc(je+1)
+ alpha = rn_e_al(je+1)
+ coefa = rn_e_th(je+1)
+ coefb = rn_e_bb(je+1)
+ ! Coefficients for SF12 stretching
+ IF (nn_strt(je+1) == 2) THEN
+ env_hc = 0._wp
+ coefa = rn_e_th(je+1) / (env3(ji,jj)-env2(ji,jj))
+ coefb = rn_e_bb(je+1) + ((env3(ji,jj)-env2(ji,jj))*rn_e_ba(je+1))
+ coefb = 1.0_wp-(coefb/(env3(ji,jj)-env2(ji,jj)))
+ END IF
+ ! -------------------------------------------
+ ! Analytical derivative of coordinate
+ ! transformation at envelope (je+1)
+ !
+ zcoefw = 0._wp
+ dCds_bot = dstretch( zcoefw, coefa, coefb, alpha, &
+ & nn_slev(je+1), nn_strt(je+1))
+ d(2,2) = dzdmes(env2(ji,jj), env3(ji,jj), dCds_bot, env_hc)
+ ELSE
+ ! Splines are applied in the bottom most sub-zone:
+ ibcend = 2 ! Bound. cond for the 2nd derivative
+ d(2,2) = 0._wp ! 2nd der = 0
+ ENDIF
+ !
+ ! 2) Computing spline's coefficients
+ ! ===============================================
+ ! Interval's endpoints TAU are 0 and 1.
+ tau(1) = 0._wp
+ tau(2) = 1._wp
+ c = cub_spl(tau, d, 2, ibcbeg, ibcend )
+ ! TO ADD CONSTRAINED CUBIC SPLINES CASE
+ !
+ ! 3) Stretching coefficients for this sub-zone je
+ ! ===============================================
+ env_hc = rn_e_hc(je)
+ coefa = rn_e_th(je)
+ coefb = rn_e_bb(je)
+ alpha = rn_e_al(je)
+ !
+ DO jk = 1, nn_s
+ jks = jks1 + jk - 1
+ !
+ ! 4) Computing stretched distribution of
+ ! interpolation points -> they represent
+ ! MEs-coordinates (-1 <= MEs <= 0)
+ ! ===============================================
+ zcoefw = -stretch( sigma('W',jk,nn_s), coefa, &
+ & coefb, alpha, nn_s, nn_strt(je) ) ! W-GRID
+ zcoeft = -stretch( sigma('T',jk,nn_s), coefa, &
+ & coefb, alpha, nn_s, nn_strt(je) ) ! T-GRID
+ !
+ ! For debugging
+ z_gsigw3(ji,jj,jks) = zcoefw
+ z_gsigt3(ji,jj,jks) = zcoeft
+ !
+ ! 5) Computing model levels depths
+ ! ===============================================
+ gdept_0(ji,jj,jks) = z_cub_spl(zcoeft, tau, c)
+ gdepw_0(ji,jj,jks) = z_cub_spl(zcoefw, tau, c)
+ !
+ END DO ! jk
+ END DO ! ji
+ END DO ! jj
+ END IF
+ END DO ! je
+
+ !=================================
+ ! COMPUTE e3t, e3w for ALL levels
+ ! as finite differences
+ ! ================================
+ DO jk=1,jpkm1
+ e3t_0(:,:,jk) = gdepw_0(:,:,jk+1) - gdepw_0(:,:,jk)
+ e3w_0(:,:,jk+1) = gdept_0(:,:,jk+1) - gdept_0(:,:,jk)
+ ENDDO
+ ! Surface
+ jk = 1
+ e3w_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,1) - gdepw_0(:,:,1))
+ !
+ ! Bottom
+ jk = jpk
+ e3t_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,jk) - gdepw_0(:,:,jk))
+
+ ! Lateral B.C.
+ ! -----------------------------------------------
+ CALL lbc_lnk('dommes',e3t_0,'T', 1.)
+ CALL lbc_lnk('dommes',e3w_0,'W', 1.)
+
+ !==============================
+ ! Computing mbathy
+ !==============================
+ IF( lwp ) WRITE(numout,*) ''
+ IF( lwp ) WRITE(numout,*) 'MEs mbathy:'
+ IF( lwp ) WRITE(numout,*) ''
+
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ DO jk = 1, jpkm1
+ IF( bathy(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
+ END DO
+ END DO
+ END DO
+
+ WHERE (bathy(:,:)<=0) mbathy(:,:) = 0
+
+ !==============================================
+ ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0
+ !==============================================
+ DO jk = 1, jpk
+ DO jj = 1,jpjm1
+ DO ji = 1,jpim1
+ e3u_0(ji,jj,jk) = ( MIN(1,mbathy(ji ,jj)) * e3t_0(ji ,jj,jk) + &
+ & MIN(1,mbathy(ji+1,jj)) * e3t_0(ji+1,jj,jk) ) / &
+ & MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)))
+
+ e3v_0(ji,jj,jk) = ( MIN(1,mbathy(ji,jj )) * e3t_0(ji,jj ,jk) + &
+ & MIN(1,mbathy(ji,jj+1)) * e3t_0(ji,jj+1,jk) ) / &
+ & MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)))
+
+ e3uw_0(ji,jj,jk) = ( MIN(1,mbathy(ji ,jj)) * e3w_0(ji ,jj,jk) + &
+ & MIN(1,mbathy(ji+1,jj)) * e3w_0(ji+1,jj,jk) ) / &
+ & MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)))
+
+ e3vw_0(ji,jj,jk) = ( MIN(1,mbathy(ji,jj )) * e3w_0(ji,jj ,jk) + &
+ & MIN(1,mbathy(ji,jj+1)) * e3w_0(ji,jj+1,jk) ) / &
+ & MAX(1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)))
+
+ e3f_0(ji,jj,jk) = ( MIN(1,mbathy(ji ,jj )) * e3t_0(ji ,jj ,jk) + &
+ & MIN(1,mbathy(ji+1,jj )) * e3t_0(ji+1,jj ,jk) + &
+ & MIN(1,mbathy(ji+1,jj+1)) * e3t_0(ji+1,jj+1,jk) + &
+ & MIN(1,mbathy(ji ,jj+1)) * e3t_0(ji ,jj+1,jk) ) / &
+ & MAX(1, MIN(1,mbathy(ji ,jj))+MIN(1,mbathy(ji ,jj+1)) &
+ & + MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1)))
+
+ END DO
+ END DO
+ END DO
+ !==============================================
+ ! Adjusting for geopotential levels, if any
+ !==============================================
+ jks1 = 2
+ nn_s = 0
+ max_env_top = 0._wp ! surface
+ min_env_top = 0._wp ! surface
+ DO je = 1, nn_env ! LOOP over all the requested envelopes
+ jks1 = jks1 + nn_s - 1
+ nn_s = nn_slev(je)
+ if ( je > 1 ) THEN
+ nn_s = nn_s + 1
+ max_env_top = max_dep_env(je-1)
+ min_env_top = min_dep_env(je-1)
+ END IF
+ max_env_bot = max_dep_env(je)
+ min_env_bot = min_dep_env(je)
+
+ IF ( max_env_top == min_env_top .AND. max_env_bot == min_env_bot ) THEN
+ IF( lwp ) WRITE(numout,*) '... Adjusting for geopotential levels ...'
+ DO jj = 1,jpjm1
+ DO ji = 1,jpim1
+ DO jk = 1, nn_s
+ jks = jks1 + jk - 1
+ !
+ e3u_0 (ji,jj,jks) = MIN( e3t_0(ji,jj,jks), e3t_0(ji+1,jj ,jks))
+ e3uw_0(ji,jj,jks) = MIN( e3w_0(ji,jj,jks), e3w_0(ji+1,jj ,jks))
+ e3v_0 (ji,jj,jks) = MIN( e3t_0(ji,jj,jks), e3t_0(ji ,jj+1,jks))
+ e3vw_0(ji,jj,jks) = MIN( e3w_0(ji,jj,jks), e3w_0(ji ,jj+1,jks))
+ e3f_0 (ji,jj,jks) = MIN( e3t_0(ji,jj,jks), e3t_0(ji+1,jj ,jks))
+ !
+ END DO
+ END DO
+ END DO
+ END IF
+ END DO
+ !
+ ! Lateral B.C.
+ ! -----------------------------------------------
+ CALL lbc_lnk( 'dommes', e3u_0 , 'U', 1._wp )
+ CALL lbc_lnk( 'dommes', e3v_0 , 'V', 1._wp )
+ CALL lbc_lnk( 'dommes', e3f_0 , 'F', 1._wp )
+ CALL lbc_lnk( 'dommes', e3uw_0, 'U', 1._wp )
+ CALL lbc_lnk( 'dommes', e3vw_0, 'V', 1._wp )
+ !
+ WHERE (e3t_0 (:,:,:) == 0.0) e3t_0(:,:,:) = 1.0
+ WHERE (e3u_0 (:,:,:) == 0.0) e3u_0(:,:,:) = 1.0
+ WHERE (e3v_0 (:,:,:) == 0.0) e3v_0(:,:,:) = 1.0
+ WHERE (e3f_0 (:,:,:) == 0.0) e3f_0(:,:,:) = 1.0
+ WHERE (e3w_0 (:,:,:) == 0.0) e3w_0(:,:,:) = 1.0
+ WHERE (e3uw_0 (:,:,:) == 0.0) e3uw_0(:,:,:) = 1.0
+ WHERE (e3vw_0 (:,:,:) == 0.0) e3vw_0(:,:,:) = 1.0
+
+ !!----------------------------------------------------------------------
+ END SUBROUTINE mes_bld
+
+ SUBROUTINE mes_ini
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE mes_ini ***
+ !!
+ !! ** Purpose : Initialise a ME s-coordinates systems
+ !!
+ !!----------------------------------------------------------------------
+ CHARACTER(lc) :: env_name ! name of the externally defined envelopes
+ INTEGER :: ji, jj, je
+ INTEGER :: ierr, inum, ios
+
+ NAMELIST/namzgr_mes/rn_bot_min, rn_bot_max, nn_strt, &
+ & nn_slev , rn_e_hc , rn_e_th, &
+ & rn_e_bb , rn_e_ba , rn_e_al
+ !-----------------------------------------------------------------------
+ !
+ ! Initialising some variables and arrays
+ ALLOCATE( envl(jpi, jpj, max_nn_env), STAT = ierr )
+ ALLOCATE(z_gsigw3(jpi, jpj, jpk), z_gsigt3(jpi, jpj, jpk), STAT = ierr )
+ nn_env = 0
+ nn_strt(:) = 0
+ nn_slev(:) = 0
+ max_dep_env(:) = 0.0
+ min_dep_env(:) = 0.0
+
+ gdept_0(:,:,:) = 0._wp ; gdepw_0(:,:,:) = 0._wp ;
+ e3t_0 (:,:,:) = 0._wp ; e3w_0 (:,:,:) = 0._wp ;
+ e3u_0 (:,:,:) = 0._wp ; e3v_0 (:,:,:) = 0._wp ;
+ e3uw_0 (:,:,:) = 0._wp ; e3vw_0 (:,:,:) = 0._wp ;
+ !
+ IF(lwp) THEN ! control print
+ WRITE(numout,*) ''
+ WRITE(numout,*) 'mes_ini: Setting a Multi-Envelope s-coordinate system (Bruciaferri et al. 2018)'
+ WRITE(numout,*) '~~~~~~~~'
+ END IF
+ !
+ ! Namelist namzgr_mes in reference namelist:
+ !REWIND( numnam_ref )
+ READ ( numnam_ref, namzgr_mes, IOSTAT = ios, ERR = 901)
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in reference namelist')
+ !
+ ! Namelist namzgr_mes in configuration namelist:
+ !REWIND( numnam_cfg )
+ READ ( numnam_cfg, namzgr_mes, IOSTAT = ios, ERR = 902 )
+902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in configuration namelist')
+ IF(lwm) WRITE ( numond, namzgr_mes )
+
+ ! 1) Reading Bathymetry and envelopes
+ IF ( ntopo == 1 ) THEN
+ IF ( lwp ) THEN
+ WRITE(numout,*) ' Reading bathymetry and envelopes ...'
+ WRITE(numout,*) ''
+ END IF
+ !
+ CALL iom_open ( cn_topo, inum )
+ CALL iom_get ( inum, jpdom_auto, cn_bath, bathy )
+ !
+ DO je = 1, max_nn_env
+ WRITE(env_name, '(A6,I0)') 'hbatt_', je
+ IF ( iom_varid( inum, TRIM(env_name), ldstop = .FALSE. ) > 0 ) THEN
+ CALL iom_get ( inum, jpdom_auto, TRIM(env_name), envl(:,:,je) )
+ nn_env = nn_env + 1
+ END IF
+ END DO
+ !
+ CALL iom_close( inum )
+ ELSE
+ WRITE(ctlmes,*) 'parameter , ntopo = ', ntopo
+ CALL ctl_stop( ctlmes )
+ ENDIF
+ !
+ ! 2) Checking consistency of envelopes
+ DO je = 1, nn_env-1
+ WRITE(ctlmes,*) 'Envelope ', je+1, ' is shallower that Envelope ', je
+ IF (MAXVAL(envl(:,:,je+1)) < MAXVAL(envl(:,:,je))) CALL ctl_stop( ctlmes )
+ ENDDO
+ !
+ ! 3) Checking SF12 stretching function
+ ! is used only in the upper sub-zone
+ DO je = 2, nn_env
+ IF ( nn_strt(je) == 2 ) THEN
+ WRITE(ctlmes,*) 'SF12 stretching function MUST be used only in the upper sub-zone'
+ CALL ctl_stop( ctlmes )
+ END IF
+ END DO
+ !
+ ! 4) Computing max and min depths of envelopes
+ DO je = 1, nn_env
+ max_dep_env(je) = MAXVAL(envl(:,:,je))
+ min_dep_env(je) = MINVAL(envl(:,:,je))
+ IF( lk_mpp ) CALL mpp_max( 'mes_ini', max_dep_env(je) )
+ IF( lk_mpp ) CALL mpp_min( 'mes_ini', min_dep_env(je) )
+ END DO
+ IF( lk_mpp ) CALL mppsync
+ !
+ ! 5) Set maximum and minimum ocean depth
+ bathy(:,:) = MIN( rn_bot_max, bathy(:,:) )
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_bot_min, bathy(ji,jj) )
+ END DO
+ END DO
+ !
+ IF( lwp ) THEN ! control print
+ WRITE(numout,*)
+ WRITE(numout,*) ' GEOMETRICAL FEATURES OF THE REQUESTED MEs GRID:'
+ WRITE(numout,*) ' -----------------------------------------------'
+ WRITE(numout,*) ' Minimum depth of the ocean rn_bot_min = ', rn_bot_min
+ WRITE(numout,*) ' Maximum depth of the ocean rn_bot_max = ', rn_bot_max
+ WRITE(numout,*) ''
+ WRITE(numout,*) ' ENVELOPES:'
+ WRITE(numout,*) ' ========= '
+ WRITE(numout,*) ''
+ DO je = 1, nn_env
+ WRITE(numout,*) ' * envelope ', je, ': min depth = ', min_dep_env(je)
+ WRITE(numout,*) ' max depth = ', max_dep_env(je)
+ WRITE(numout,*) ''
+ END DO
+ WRITE(numout,*) ''
+ WRITE(numout,*) ' SUBDOMAINS:'
+ WRITE(numout,*) ' ========== '
+ WRITE(numout,*) ''
+ DO je = 1, nn_env
+ WRITE(numout,*) ' * subdomain ', je,':'
+ IF ( je == 1) THEN
+ WRITE(numout,*) ' Shallower envelope: free surface'
+ ELSE
+ WRITE(numout,*) ' Shallower envelope: envelope ',je-1
+ END IF
+ WRITE(numout,*) ' Deeper envelope: envelope ',je
+ WRITE(numout,*) ' Number of MEs-levs: nn_slev(',je,') = ',nn_slev(je)
+ IF( isodd(je) ) THEN
+ WRITE(numout,*) ' Stretched s-coordinates: '
+ ELSE
+ WRITE(numout,*) ' Stretched CUBIC SPLINES: '
+ END IF
+ IF ( nn_strt(je) == 0 ) WRITE(numout,*) ' M96 stretching function'
+ IF ( nn_strt(je) == 1 ) WRITE(numout,*) ' SH94 stretching function'
+ IF ( nn_strt(je) == 2 ) WRITE(numout,*) ' SF12 stretching function'
+ WRITE(numout,*) ' critical depth rn_e_hc(',je,') = ',rn_e_hc(je)
+ WRITE(numout,*) ' surface stretc. coef. rn_e_th(',je,') = ',rn_e_th(je)
+ IF( nn_strt(je) == 2 ) THEN
+ WRITE(numout,*) ' bottom stretc. coef. rn_e_ba(',je,') = ',rn_e_ba(je)
+ END IF
+ WRITE(numout,*) ' bottom stretc. coef. rn_e_bb(',je,') = ',rn_e_bb(je)
+ IF( nn_strt(je) == 2 ) THEN
+ WRITE(numout,*) ' bottom stretc. coef. rn_e_al(',je,') = ',rn_e_al(je)
+ END IF
+ WRITE(numout,*) ' ------------------------------------------------------------------'
+ ENDDO
+ ENDIF
+ END SUBROUTINE mes_ini
+ !!
+
+ FUNCTION isodd( a ) RESULT( odd )
+ !!---------------------------------------------------------------------
+ !! *** FUNCTION isodd ***
+ !!
+ !! ** Purpose : Determine whether an interger is odd or not
+ !!
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT (in) :: a
+ LOGICAL :: odd
+ odd = btest(a, 0)
+ END FUNCTION isodd
+
+ FUNCTION sech( x ) RESULT( seh )
+ !!---------------------------------------------------------------------
+ !! *** FUNCTION sech ***
+ !!
+ !! ** Purpose : Compute sech of a real number
+ !!
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in ) :: x
+ REAL(wp) :: seh
+ seh = 1._wp / COSH(x)
+ END FUNCTION sech
+
+ FUNCTION sigma( vgrid, pk, kmax ) RESULT( ps1 )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE sigma ***
+ !!
+ !! ** Purpose : provide the analytical function for sigma-coordinate
+ !! (not stretched s-coordinate).
+ !!
+ !! ** Method : the function provide the non-dimensional position of
+ !! T and W points (i.e. between 0 and 1).
+ !! T-points at integer values (between 1 and jpk)
+ !! W-points at integer values - 1/2 (between 0.5 and
+ !! jpk-0.5)
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT (in) :: pk ! continuous k coordinate
+ INTEGER, INTENT (in) :: kmax ! number of levels
+ CHARACTER(LEN=1), INTENT (in) :: vgrid ! type of vertical grid: T, W
+ REAL(wp) :: kindx ! index of T or W points
+ REAL(wp) :: ps1 ! value of sigma coordinate (-1 <= ps1 <=0)
+ !
+ SELECT CASE (vgrid)
+ CASE ('T') ! T-points at integer values (between 1 and jpk)
+ kindx = REAL(pk,wp)
+ CASE ('W') ! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+ kindx = REAL(pk,wp) - 0.5_wp
+ CASE DEFAULT
+ WRITE(ctlmes,*) 'No valid vertical grid option selected'
+ END SELECT
+ ps1 = -(kindx - 0.5) / REAL(kmax-1) ! The surface is at the first W level
+ ! while the bottom coincides with the
+ ! last W level
+ END FUNCTION sigma
+
+ FUNCTION stretch( s, ca, cb, alpha, kmax, ftype ) RESULT ( pf1 )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE stretch ***
+ !!
+ !! ** Purpose : provide the analytical stretching function
+ !! for s-coordinate.
+ !!
+ !! ** Method : if ftype = 2: Siddorn and Furner 2012
+ !! analytical function is used
+ !! if ftype = 1: Song and Haidvogel 1994
+ !! analytical function is used
+ !! if ftype = 0: Madec et al. 1996
+ !! analytical function is used
+ !! (pag 65 of NEMO Manual)
+ !! Reference : Madec, Lott, Delecluse and
+ !! Crepon, 1996. JPO, 26, 1393-1408
+ !!
+ !! s MUST be NEGATIVE: -1 <= s <= 0
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in) :: s ! not stretched sigma-coordinate
+ REAL(wp), INTENT(in) :: ca ! surface stretch. coeff
+ REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff
+ REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12
+ INTEGER, INTENT (in) :: kmax ! number of levels
+ INTEGER, INTENT(in) :: ftype ! type of stretching function
+ REAL(wp) :: pf1 ! stretched s-coordinate (-1 <= pf1 <=0)
+ ! SF12 stretch. funct. parameters
+ REAL(wp) :: psmth ! smoothing parameter for transition
+ ! between shallow and deep areas
+ REAL(wp) :: za1,za2,za3 ! local variables
+ REAL(wp) :: zn1,zn2,ps ! local variables
+ REAL(wp) :: za,zb,zx ! local variables
+ !!----------------------------------------------------------------------
+ SELECT CASE (ftype)
+ CASE (0) ! M96 stretching function
+ pf1 = ( TANH( ca * ( s + cb ) ) - TANH( cb * ca ) ) &
+ & * ( COSH( ca ) &
+ & + COSH( ca * ( 2.e0 * cb - 1.e0 ) ) ) &
+ & / ( 2._wp * SINH( ca ) )
+ CASE (1) ! SH94 stretching function
+ IF ( ca == 0 ) THEN ! uniform sigma
+ pf1 = s
+ ELSE ! stretched sigma
+ pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) + &
+ & cb * ( ( TANH(ca*(s + 0.5_wp)) - TANH(0.5_wp*ca) ) / &
+ & (2._wp*TANH(0.5_wp*ca)) )
+ END IF
+ CASE (2) ! SF12 stretching function
+ psmth = 1.0_wp ! We consider only the case for efold = 0
+ ps = -s
+ zn1 = 1. / REAL(kmax-1)
+ zn2 = 1. - zn1
+
+ za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp)
+ za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp)
+ za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
+
+ za = cb - za3*(ca-za1)-za2
+ za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) )
+ zb = (ca - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
+ zx = 1.0_wp-za/2.0_wp-zb
+
+ pf1 = za*(ps*(1.0_wp-ps/2.0_wp))+zb*ps**3.0_wp + &
+ & zx*( (alpha+2.0_wp)*ps**(alpha+1.0_wp) - &
+ & (alpha+1.0_wp)*ps**(alpha+2.0_wp) )
+ pf1 = -pf1*psmth+ps*(1.0_wp-psmth)
+ END SELECT
+
+ END FUNCTION stretch
+
+ FUNCTION dstretch( s, ca, cb, alpha, kmax, ftype) RESULT( pf1 )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE dstretch ***
+ !!
+ !! ** Purpose : provides the 1st derivative of the analytical
+ !! stretching function dC(sigma)/dsigma.
+ !!
+ !! ** Method : if ftype = 2: Siddorn and Furner 2012
+ !! analytical function is used
+ !! if ftype = 1: Song and Haidvogel 1994
+ !! analytical function is used
+ !! if ftype = 0: Madec et al. 1996
+ !! analytical function is used
+ !! (pag 65 of NEMO Manual)
+ !! Reference : Madec, Lott, Delecluse and
+ !! Crepon, 1996. JPO, 26, 1393-1408
+ !!
+ !! s MUST be NEGATIVE: -1 <= s <= 0
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in ) :: s ! not stretched sigma-coordinate
+ REAL(wp), INTENT(in) :: ca ! surface stretch. coeff
+ REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff
+ REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12
+ INTEGER, INTENT (in) :: kmax ! number of levels
+ INTEGER, INTENT(in ) :: ftype ! type of stretching function
+ REAL(wp) :: pf1 ! first derivative
+ ! SF12 stretch. funct. parameters
+ REAL(wp) :: psmth ! smoothing parameter for transition
+ ! between shallow and deep areas
+ REAL(wp) :: za1,za2,za3 ! local variables
+ REAL(wp) :: zn1,zn2,ps ! local variables
+ REAL(wp) :: za,zb,zx ! local variables
+ REAL(wp) :: zt1,zt2,zt3 ! local variables
+ !!----------------------------------------------------------------------
+
+ SELECT CASE (ftype)
+
+ CASE (0) ! M96 stretching function
+ pf1 = ( ca * ( COSH(ca*(2._wp * cb - 1._wp)) + COSH(ca)) * &
+ sech(ca * (s+cb))**2 ) / (2._wp * SINH(ca))
+ CASE (1) ! SH94 stretching function
+ IF ( ca == 0 ) then ! uniform sigma
+ pf1 = 1._wp / REAL(kmax,wp)
+ ELSE ! stretched sigma
+ !pf1 = (1._wp - cb) * ca * COSH(ca*s) / (SINH(ca) * REAL(kmax,wp)) + &
+ ! cb * ca * &
+ ! (sech((s+0.5_wp)*ca)**2) / (2._wp * TANH(0.5_wp*ca) * REAL(kmax,wp))
+ pf1 = (1._wp - cb) * ca * COSH(ca*s) / SINH(ca) + cb * ca * &
+ (sech((s+0.5_wp)*ca)**2) / (2._wp * TANH(0.5_wp*ca))
+ END IF
+ CASE (2) ! SF12 stretching function
+ ps = -s
+ zn1 = 1. / REAL(kmax-1)
+ zn2 = 1. - zn1
+
+ za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp)
+ za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp)
+ za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
+
+ za = cb - za3*(ca-za1)-za2
+ za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) )
+ zb = (ca - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
+ zx = (alpha+2.0_wp)*(alpha+1.0_wp)
+
+ zt1 = 0.5_wp*za*(ps-1._wp)*((ps**2.0_wp + 3._wp*ps + 2._wp)*ps**alpha - 2._wp)
+ zt2 = zb*(zx*ps**(alpha+1.0_wp) - zx*ps**(alpha) + 3._wp*ps**2._wp)
+ zt3 = zx*(ps-1._wp)*ps**(alpha)
+
+ pf1 = zt1 + zt2 - zt3
+
+ END SELECT
+
+ END FUNCTION dstretch
+
+ FUNCTION z_mes(dep_top, dep_bot, Cs, s, hc) RESULT( z )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE z_mes ***
+ !!
+ !! ** Purpose : provide the analytical trasformation from the
+ !! computational space (mes-coordinate) to the
+ !! physical space (depth z)
+ !!
+ !! N.B.: z is downward positive defined as well as envelope surfaces.
+ !! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1.
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in ) :: dep_top ! shallower envelope
+ REAL(wp), INTENT(in ) :: dep_bot ! deeper envelope
+ REAL(wp), INTENT(in ) :: Cs ! stretched s coordinate
+ REAL(wp), INTENT(in ) :: s ! not stretched s coordinate
+ REAL(wp), INTENT(in ) :: hc ! critical depth
+ REAL :: z ! downward positive depth
+ !!----------------------------------------------------------------------
+ !
+ z = dep_top + hc*s + Cs*(dep_bot - hc - dep_top)
+
+ END FUNCTION z_mes
+
+ FUNCTION dzdmes(dep_top, dep_bot, dCds, hc) RESULT( d1z )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE Dzdmes ***
+ !!
+ !! ** Purpose : provide the 1st derivative of the analytical
+ !! trasformation from the computational space
+ !! (mes-coordinate) to the physical space (depth z)
+ !!
+ !! N.B.: z is downward positive defined as well as envelope surfaces.
+ !! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1.
+ !!----------------------------------------------------------------------
+ REAL(wp), INTENT(in ) :: dep_top ! shallower envelope
+ REAL(wp), INTENT(in ) :: dep_bot ! deeper envelope
+ REAL(wp), INTENT(in ) :: dCds ! 1st derivative of the
+ ! stretched s coordinate
+ REAL(wp), INTENT(in ) :: hc ! critical depth
+ REAL(wp) :: d1z ! 1st derivative of downward
+ ! positive depth
+ !!----------------------------------------------------------------------
+ !
+ d1z = hc + dCds * (dep_bot - hc - dep_top)
+
+ END FUNCTION Dzdmes
+
+ FUNCTION cub_spl(tau, d, n, ibcbeg, ibcend) RESULT ( c )
+ !!----------------------------------------------------------------------
+ !!
+ !! CUBSPL defines an interpolatory cubic spline.
+ !!
+ !! Discussion:
+ !!
+ !! A tridiagonal linear system for the unknown slopes S(I) of
+ !! F at TAU(I), I=1,..., N, is generated and then solved by Gauss
+ !! elimination, with S(I) ending up in C(2,I), for all I.
+ !!
+ !! Author: Carl de Boor
+ !!
+ !! Reference: Carl de Boor, Practical Guide to Splines,
+ !! Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27.
+ !!
+ !! Parameters:
+ !!
+ !! Input:
+ !! TAU(N) : the abscissas or X values of the data points.
+ !! The entries of TAU are assumed to be strictly
+ !! increasing.
+ !!
+ !! N : the number of data points. N is assumed to be
+ !! at least 2.
+ !!
+ !! IBCBEG : boundary condition indicator at TAU(1).
+ !! = 0 no boundary condition at TAU(1) is given.
+ !! In this case, the "not-a-knot condition"
+ !! is used.
+ !! = 1 the 1st derivative at TAU(1) is equal to
+ !! the input value D(2,1).
+ !! = 2 the 2nd derivative at TAU(1) is equal to
+ !! the input value D(2,1).
+ !!
+ !! IBCEND : boundary condition indicator at TAU(N).
+ !! = 0 no boundary condition at TAU(N) is given.
+ !! In this case, the "not-a-knot condition"
+ !! is used.
+ !! = 1 the 1st derivative at TAU(N) is equal to
+ !! the input value D(2,2).
+ !! = 2 the 2nd derivative at TAU(N) is equal to
+ !! the input value D(2,2).
+ !!
+ !! D(1,1): value of the function at TAU(1)
+ !! D(1,1): value of the function at TAU(N)
+ !! D(2,1): if IBCBEG is 1 (2) it is the value of the
+ !! 1st (2nd) derivative at TAU(1)
+ !! D(2,2): if IBCBEG is 1 (2) it is the value of the
+ !! 1st (2nd) derivative at TAU(N)
+ !!
+ !! Output:
+ !! C(4,N) : contains the polynomial coefficients of the
+ !! cubic interpolating spline.
+ !! In the interval interval (TAU(I), TAU(I+1)),
+ !! the spline F is given by
+ !!
+ !! F(X) =
+ !! C(1,I) +
+ !! C(2,I) * H +
+ !! C(3,I) * H^2 / 2 +
+ !! C(4,I) * H^3 / 6.
+ !!
+ !! where H = X - TAU(I).
+ !!
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT(in ) :: n
+ REAL(wp), INTENT(in ) :: tau(n)
+ REAL(wp), INTENT(in ) :: d(2,2)
+ INTEGER, INTENT(in ) :: ibcbeg, ibcend
+ REAL(wp) :: c(4,n)
+ REAL(wp) :: divdf1
+ REAL(wp) :: divdf3
+ REAL(wp) :: dtau
+ REAL(wp) :: g
+ INTEGER(wp) :: i
+ !!----------------------------------------------------------------------
+ ! Initialise c and copy d values to c
+ c(:,:) = 0.0D+00
+ c(1,1) = d(1,1)
+ c(1,n) = d(1,2)
+ c(2,1) = d(2,1)
+ c(2,n) = d(2,2)
+ !
+ ! C(3,*) and C(4,*) are used initially for temporary storage.
+ ! Store first differences of the TAU sequence in C(3,*).
+ ! Store first divided difference of data in C(4,*).
+ DO i = 2, n
+ c(3,i) = tau(i) - tau(i-1)
+ c(4,i) = ( c(1,i) - c(1,i-1) ) / c(3,i)
+ END DO
+
+ ! Construct the first equation from the boundary condition
+ ! at the left endpoint, of the form:
+ !
+ ! C(4,1) * S(1) + C(3,1) * S(2) = C(2,1)
+ !
+ ! IBCBEG = 0: Not-a-knot
+ IF ( ibcbeg == 0 ) THEN
+
+ IF ( n <= 2 ) THEN
+ c(4,1) = 1._wp
+ c(3,1) = 1._wp
+ c(2,1) = 2._wp * c(4,2)
+ ELSE
+ c(4,1) = c(3,3)
+ c(3,1) = c(3,2) + c(3,3)
+ c(2,1) = ( ( c(3,2) + 2._wp * c(3,1) ) * c(4,2) &
+ * c(3,3) + c(3,2)**2 * c(4,3) ) / c(3,1)
+ END IF
+
+ ! IBCBEG = 1: derivative specified.
+ ELSE IF ( ibcbeg == 1 ) then
+
+ c(4,1) = 1._wp
+ c(3,1) = 0._wp
+
+ ! IBCBEG = 2: Second derivative prescribed at left end.
+ ELSE IF ( ibcbeg == 2 ) then
+
+ c(4,1) = 2._wp
+ c(3,1) = 1._wp
+ c(2,1) = 3._wp * c(4,2) - c(3,2) / 2._wp * c(2,1)
+
+ ELSE
+ WRITE(ctlmes,*) 'CUBSPL - Error, invalid IBCBEG input option!'
+ CALL ctl_stop( ctlmes )
+ END IF
+
+ ! If there are interior knots, generate the corresponding
+ ! equations and carry out the forward pass of Gauss,
+ ! elimination after which the I-th equation reads:
+ !
+ ! C(4,I) * S(I) + C(3,I) * S(I+1) = C(2,I).
+
+ IF ( n > 2 ) THEN
+
+ DO i = 2, n-1
+ g = -c(3,i+1) / c(4,i-1)
+ c(2,i) = g * c(2,i-1) + 3._wp * &
+ ( c(3,i) * c(4,i+1) + c(3,i+1) * c(4,i) )
+ c(4,i) = g * c(3,i-1) + 2._wp * ( c(3,i) + c(3,i+1))
+ END DO
+
+ ! Construct the last equation from the second boundary,
+ ! condition of the form
+ !
+ ! -G * C(4,N-1) * S(N-1) + C(4,N) * S(N) = C(2,N)
+ !
+ ! If 1st der. is prescribed at right end ( ibcend == 1 ),
+ ! one can go directly to back-substitution, since the C
+ ! array happens to be set up just right for it at this
+ ! point.
+
+ IF ( ibcend < 1 ) THEN
+ ! Not-a-knot and 3 <= N, and either 3 < N or also
+ ! not-a-knot at left end point.
+ IF ( n /= 3 .OR. ibcbeg /= 0 ) THEN
+ g = c(3,n-1) + c(3,n)
+ c(2,n) = ( ( c(3,n) + 2._wp * g ) * c(4,n) * c(3,n-1) &
+ + c(3,n)**2 * ( c(1,n-1) - c(1,n-2) ) / &
+ c(3,n-1) ) / g
+ g = - g / c(4,n-1)
+ c(4,n) = c(3,n-1)
+ c(4,n) = c(4,n) + g * c(3,n-1)
+ c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
+ ELSE
+ ! N = 3 and not-a-knot also at left.
+ c(2,n) = 2._wp * c(4,n)
+ c(4,n) = 1._wp
+ g = -1._wp / c(4,n-1)
+ c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
+ c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
+ END IF
+
+ ELSE IF ( ibcend == 2 ) THEN
+ ! IBCEND = 2: Second derivative prescribed at right endpoint.
+ c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n)
+ c(4,n) = 2._wp
+ g = -1._wp / c(4,n-1)
+ c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
+ c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
+ END IF
+
+ ELSE
+ ! N = 2 (assumed to be at least equal to 2!).
+
+ IF ( ibcend == 2 ) THEN
+
+ c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n)
+ c(4,n) = 2._wp
+ g = -1._wp / c(4,n-1)
+ c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
+ c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
+
+ ELSE IF ( ibcend == 0 .AND. ibcbeg /= 0 ) THEN
+
+ c(2,n) = 2._wp * c(4,n)
+ c(4,n) = 1._wp
+ g = -1._wp / c(4,n-1)
+ c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
+ c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
+
+ ELSE IF ( ibcend == 0 .AND. ibcbeg == 0 ) THEN
+
+ c(2,n) = c(4,n)
+
+ END IF
+
+ END IF
+
+ ! Back solve the upper triangular system
+ !
+ ! C(4,I) * S(I) + C(3,I) * S(I+1) = B(I)
+ !
+ ! for the slopes C(2,I), given that S(N) is already known.
+ !
+ DO i = n-1, 1, -1
+ c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i)
+ END DO
+
+ ! Generate cubic coefficients in each interval, that is, the
+ ! derivatives at its left endpoint, from value and slope at its
+ ! endpoints.
+ DO i = 2, n
+ dtau = c(3,i)
+ divdf1 = ( c(1,i) - c(1,i-1) ) / dtau
+ divdf3 = c(2,i-1) + c(2,i) - 2._wp * divdf1
+ c(3,i-1) = 2._wp * ( divdf1 - c(2,i-1) - divdf3 ) / dtau
+ c(4,i-1) = 6._wp * divdf3 / dtau**2
+ END DO
+
+ END FUNCTION cub_spl
+
+ FUNCTION z_cub_spl(s, tau, c) RESULT ( z_cs )
+ !----------------------------------------------------------------------
+ ! *** function z_cub_spl ***
+ !
+ ! ** Purpose : evaluate the cubic spline F in the interval
+ ! (TAU(I), TAU(I+1)). F is given by
+ !
+ !
+ ! F(X) = C1 + C2*H + (C3*H^2)/2 + (C4*H^3)/6
+ !
+ ! where H = S-TAU(I).
+ !
+ ! s MUST be positive: 0 <= s <= 1
+ !---------------------------------------------------------------------
+ !INTEGER, PARAMETER :: nn = 2
+ REAL(wp), INTENT(in ) :: s
+ REAL(wp), INTENT(in ) :: tau(2)
+ REAL(wp), INTENT(in ) :: c(4,2)
+ REAL(wp) :: z_cs
+ REAL(wp) :: c1, c2, c3, c4
+ REAL(wp) :: H
+ !---------------------------------------------------------------------
+ c1 = c(1,1)
+ c2 = c(2,1)
+ c3 = c(3,1)
+ c4 = c(4,1)
+
+ H = s - tau(1)
+
+ z_cs = c1 + c2 * H + (c3 * H**2._wp)/2._wp + (c4 * H**3._wp)/6._wp
+
+ END FUNCTION z_cub_spl
+
+END MODULE dommes
+
diff --git a/tools/DOMAINcfg/src/domwri.F90 b/tools/DOMAINcfg/src/domwri.F90
index bc5c3c6a..6ff6b50e 100644
--- a/tools/DOMAINcfg/src/domwri.F90
+++ b/tools/DOMAINcfg/src/domwri.F90
@@ -81,9 +81,9 @@ CONTAINS
! ! domain characteristics
CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
! ! type of vertical coordinate
- IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
- IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
- IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
+ IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF
+ IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF
+ IF( ln_sco .OR. ln_mes ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF
CALL iom_rstput( 0, 0, inum, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )
CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )
CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )
@@ -175,7 +175,7 @@ CONTAINS
CALL iom_rstput( 0, 0, inum, 'gdept_0' , gdept_0 , ktype = jp_r8 )
CALL iom_rstput( 0, 0, inum, 'gdepw_0' , gdepw_0 , ktype = jp_r8 )
!
- IF( ln_sco ) THEN ! s-coordinate stiffness
+ IF( ln_sco .OR. ln_mes ) THEN ! s-coordinate stiffness
CALL dom_stiff( zprt )
CALL iom_rstput( 0, 0, inum, 'stiffness', zprt ) ! Max. grid stiffness ratio
ENDIF
diff --git a/tools/DOMAINcfg/src/domzgr.F90 b/tools/DOMAINcfg/src/domzgr.F90
index f45cadbc..065ac376 100644
--- a/tools/DOMAINcfg/src/domzgr.F90
+++ b/tools/DOMAINcfg/src/domzgr.F90
@@ -44,6 +44,8 @@ MODULE domzgr
USE lib_fortran
USE dombat
USE domisf
+ USE dommes ! Multi-Envelope s-coordinates
+ USE domloc ! Localise the vertical coordinate system
#if defined key_agrif
USE agrif_connect
@@ -112,15 +114,15 @@ CONTAINS
INTEGER :: ji,jj,jk
REAL(wp) :: zrefdep ! depth of the reference level (~10m)
- NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
+ NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_mes, ln_isfcav, ln_loczgr
!!----------------------------------------------------------------------
!
!
- ! REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
+ ! REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist')
- ! REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
+ ! REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist')
IF(lwm) WRITE ( numond, namzgr )
@@ -154,15 +156,14 @@ CONTAINS
!
! ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
- CALL zgr_top_bot( k_top, k_bot ) ! with a minimum value set to 1
+ CALL zgr_top_bot( k_top, k_bot ) ! with a minimum value set to 1
! ! deepest/shallowest W level Above/Below ~10m
!!gm BUG in s-coordinate this does not work!
- zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
- nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
- nla10 = nlb10 - 1 ! deepest W level Above ~10m
+ zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
+ nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
+ nla10 = nlb10 - 1 ! deepest W level Above ~10m
!!gm end bug
-
ENDIF
IF(lwp) THEN ! Control print
@@ -173,19 +174,23 @@ CONTAINS
WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
+ WRITE(numout,*) ' ME s-coordinate ln_mes = ', ln_mes
WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav
+ IF(ln_loczgr) WRITE(numout,*) ' The chosen vertical coordinate system is implemented locally'
ENDIF
ioptio = 0 ! Check Vertical coordinate options
IF( ln_zco ) ioptio = ioptio + 1
IF( ln_zps ) ioptio = ioptio + 1
IF( ln_sco ) ioptio = ioptio + 1
+ IF( ln_mes ) ioptio = ioptio + 1
IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
!
IF ( ln_isfcav ) CALL zgr_isf_nam
ioptio = 0
IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
+ IF ( ln_mes .AND. ln_isfcav ) ioptio = ioptio + 1
IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
#if defined key_agrif
@@ -194,24 +199,34 @@ CONTAINS
! ENDIF
#endif
- IF(.NOT.ln_read_cfg) THEN
- !
- ! Build the vertical coordinate system
- ! ------------------------------------
- CALL zgr_z ! Reference z-coordinate system (always called)
- CALL zgr_bat ! Bathymetry fields (levels and meters)
- IF( ln_zco ) CALL zgr_zco ! z-coordinate
- IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
- IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
- CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
- !
- ! final adjustment of mbathy & check
- ! -----------------------------------
- CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
- k_bot = mbkt
- CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
- k_top = mikt
- WHERE( bathy(:,:) <= 0._wp ) k_top(:,:) = 0 ! set k_top to zero over land
+ IF( .NOT.ln_read_cfg ) THEN
+ !
+ ! Build the vertical coordinate system
+ ! ------------------------------------
+ CALL zgr_z ! Reference z-coordinate system (always called)
+ CALL zgr_bat ! Bathymetry fields (levels and meters)
+ IF( ln_zco ) CALL zgr_zco ! z-coordinate
+ IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
+ IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
+ IF( ln_mes ) CALL zgr_mes ! ME s-coordinate
+ !
+ ! Localising the vert. coord. system if requested
+ ! -----------------------------------------------
+ l2g_msk(:,:) = 2.0 ! Initialise the localisation
+ ! mask to global
+ IF( ln_loczgr ) CALL loc_zgr
+ !
+ ! Check bathymetry (mbathy) and suppress isolated ocean points
+ ! ------------------------------------------------------------
+ CALL zgr_bat_ctl
+ !
+ ! final adjustment of mbathy & check
+ ! -----------------------------------
+ CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
+ k_bot = mbkt
+ CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
+ k_top = mikt
+ WHERE( bathy(:,:) <= 0._wp ) k_top(:,:) = 0 ! set k_top to zero over land
ENDIF
!
#if defined key_agrif
@@ -264,102 +279,6 @@ CONTAINS
!
END SUBROUTINE dom_zgr
- SUBROUTINE zgr_read( 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 zgr_read ***
- !!
- !! ** Purpose : Read the vertical information in the domain configuration file
- !!
- !!----------------------------------------------------------------------
- LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags
- 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 vertical scale factors [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 ! - - -
- INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level
- !
- INTEGER :: jk ! dummy loop index
- INTEGER :: inum ! local logical unit
- REAL(WP) :: z_zco, z_zps, z_sco, z_cav
- REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
- !!----------------------------------------------------------------------
- !
- IF(lwp) THEN
- WRITE(numout,*)
- WRITE(numout,*) ' zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
- WRITE(numout,*) ' ~~~~~~~~'
- ENDIF
- !
- CALL iom_open( cn_domcfg, inum )
- !
- ! !* type of vertical coordinate
- CALL iom_get( inum, 'ln_zco' , z_zco )
- CALL iom_get( inum, 'ln_zps' , z_zps )
- CALL iom_get( inum, 'ln_sco' , z_sco )
- IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF
- IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF
- IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF
- !
- ! !* ocean cavities under iceshelves
- CALL iom_get( inum, 'ln_isfcav', z_cav )
- IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF
- !
- ! !* vertical scale factors
- CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate
- CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d )
- !
- CALL iom_get( inum, jpdom_global, 'e3t_0' , pe3t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy ) ! 3D coordinate
- CALL iom_get( inum, jpdom_global, 'e3u_0' , pe3u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
- CALL iom_get( inum, jpdom_global, 'e3v_0' , pe3v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
- CALL iom_get( inum, jpdom_global, 'e3f_0' , pe3f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
- CALL iom_get( inum, jpdom_global, 'e3w_0' , pe3w , cd_type = 'W', psgn = 1._wp, kfill = jpfillcopy )
- CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
- CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
- !
- ! !* depths
- ! !- old depth definition (obsolescent feature)
- IF( iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0 .AND. &
- & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0 .AND. &
- & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0 .AND. &
- & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0 ) THEN
- CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', &
- & ' depths at t- and w-points read in the domain configuration file')
- CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )
- CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
- CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept, kfill = jpfillcopy )
- CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw, kfill = jpfillcopy )
- !
- ELSE !- depths computed from e3. scale factors
- CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) ! 1D reference depth
- CALL e3_to_depth( pe3t , pe3w , pdept , pdepw ) ! 3D depths
- IF(lwp) THEN
- 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
- ENDIF
- !
- ! !* ocean top and bottom level
- CALL iom_get( inum, jpdom_global, 'top_level' , z2d ) ! 1st wet T-points (ISF)
- k_top(:,:) = NINT( z2d(:,:) )
- CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d ) ! last wet T-points
- k_bot(:,:) = NINT( z2d(:,:) )
- !
- ! reference depth for negative bathy (wetting and drying only)
- ! IF( ll_wd ) CALL iom_get( inum, 'rn_wd_ref_depth' , ssh_ref )
- !
- CALL iom_close( inum )
- !
- END SUBROUTINE zgr_read
-
SUBROUTINE zgr_top_bot( k_top, k_bot )
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_top_bot ***
@@ -797,7 +716,7 @@ CONTAINS
CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
ENDIF
!
- IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
+ IF (( .NOT. ln_sco).AND.(.NOT. ln_mes )) THEN !== set a minimum depth ==!
IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
ENDIF
@@ -1404,11 +1323,11 @@ CONTAINS
#endif
bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
- DO jj = 1, jpj
- DO ji = 1, jpi
- IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
- END DO
+ DO jj = 1, jpj
+ DO ji = 1, jpi
+ IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
END DO
+ END DO
! ! =============================
! ! Define the envelop bathymetry (hbatt)
! ! =============================
@@ -1494,10 +1413,10 @@ CONTAINS
END DO
!!bug: key_helsinki a verifer
- hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
- hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
- hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
- hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+ hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
+ hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
+ hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
+ hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
IF( lwp ) THEN
WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
@@ -1539,14 +1458,14 @@ CONTAINS
CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
!
- WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp
- WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp
- WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp
- WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp
- WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp
- WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp
- WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp
-!!
+ WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp
+ WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp
+ WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp
+ WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp
+ WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp
+ WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp
+ WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp
+
! HYBRID :
DO jj = 1, jpj
DO ji = 1, jpi
diff --git a/tools/DOMAINcfg/src/lib_mpp.F90 b/tools/DOMAINcfg/src/lib_mpp.F90
index e6a276b0..160c351a 100644
--- a/tools/DOMAINcfg/src/lib_mpp.F90
+++ b/tools/DOMAINcfg/src/lib_mpp.F90
@@ -51,6 +51,7 @@ MODULE lib_mpp
!! mpp_ini_north : initialisation of north fold
!! mpp_lbc_north_icb : alternative to mpp_nfd for extra outer halo with icebergs
!! mpp_bcast_nml : broadcast/receive namelist character buffer from reading process to all others
+ !! mpp_bcast_real : broadcast/receive double precision array buffer from reading process to all others
!!----------------------------------------------------------------------
USE dom_oce ! ocean space and time domain
USE in_out_manager ! I/O manager
@@ -70,6 +71,7 @@ MODULE lib_mpp
PUBLIC mppsend_dp, mpprecv_dp ! needed by TAM and ICB routines
PUBLIC mpp_report
PUBLIC mpp_bcast_nml
+ PUBLIC mpp_bcast_real
PUBLIC tic_tac
#if ! defined key_mpp_mpi
PUBLIC MPI_Wtime
@@ -665,7 +667,29 @@ CONTAINS
!
END SUBROUTINE mpp_bcast_nml
-
+ SUBROUTINE mpp_bcast_real( kvals, kno, kroot )
+ REAL(wp), DIMENSION(kno), INTENT(inout) :: kvals ! Array to send on kroot, receive for non-kroot
+ INTEGER , INTENT(in ) :: kno ! Number of elements in array
+ INTEGER , INTENT(in ) :: kroot ! Processor to send data
+ !!----------------------------------------------------------------------
+ !! *** routine mpp_bcast_real ***
+ !!
+ !! ** Purpose : Send array kvals to all processors
+ !!
+ !! ** Method : MPI broadcast
+ !!
+ !!-----------------------------------------------------------------------
+ !!
+ INTEGER :: iflag
+ !!-----------------------------------------------------------------------
+ !
+#if defined key_mpp_mpi
+ call MPI_BCAST( kvals, kno, mpi_double_precision, kroot, mpi_comm_oce, iflag )
+ call MPI_BARRIER(mpi_comm_oce, iflag)
+#endif
+ !
+ END SUBROUTINE mpp_bcast_real
+
!!----------------------------------------------------------------------
!! *** mppmax_a_int, mppmax_int, mppmax_a_real, mppmax_real ***
!!
--
GitLab