From 35a1b5016781215c5766c6a358cd481fc58b7bf4 Mon Sep 17 00:00:00 2001
From: Andrew Coward <acc@noc.ac.uk>
Date: Mon, 20 Mar 2023 18:41:41 +0000
Subject: [PATCH] Resolve "[WP2022 PHYS - NOC] Add GEOMETRIC parameterisation
(2)"
---
cfgs/SHARED/field_def_nemo-oce.xml | 16 +
cfgs/SHARED/namelist_ref | 31 ++
sette/prepare_job.sh | 3 +
src/OCE/DYN/dynldf_lev.F90 | 28 +-
src/OCE/DYN/dynldf_lev_rot_scheme.h90 | 11 +
src/OCE/LDF/ldfdyn.F90 | 8 +-
src/OCE/LDF/ldfeke.F90 | 610 ++++++++++++++++++++++++++
src/OCE/LDF/ldftra.F90 | 24 +-
src/OCE/nemogcm.F90 | 9 +-
src/OCE/step_oce.F90 | 1 +
src/OCE/stpmlf.F90 | 6 +
src/OFF/nemogcm.F90 | 16 +-
12 files changed, 744 insertions(+), 19 deletions(-)
create mode 100644 src/OCE/LDF/ldfeke.F90
diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml
index 6ad835d9..d66e26cc 100644
--- a/cfgs/SHARED/field_def_nemo-oce.xml
+++ b/cfgs/SHARED/field_def_nemo-oce.xml
@@ -229,6 +229,15 @@ that are available in the tidal-forcing implementation (see
<field id="cfl_cv" long_name="v-courant number" unit="#" grid_ref="grid_T_2D_inner" />
<field id="cfl_cw" long_name="w-courant number" unit="#" grid_ref="grid_T_2D_inner" />
+ <!-- GEOMETRIC fields (requires nn_aei_ijk_t = 32) -->
+ <field id="eke" long_name="total EKE (EKE+EPE)" unit="m3/s2" />
+ <field id="trd_eke_adv_ubt" long_name="ubt advective trend of EKE (LHS)" unit="m3/s3" grid_ref="grid_T_2D_inner" />
+ <field id="trd_eke_adv_wav" long_name="wav advective trend of EKE (LHS)" unit="m3/s3" grid_ref="grid_T_2D_inner" />
+ <field id="trd_eke_lap" long_name="diffusive trend of EKE (RHS)" unit="m3/s3" grid_ref="grid_T_2D_inner" />
+ <field id="trd_eke_peS" long_name="PE to EKE source trend (RHS)" unit="m3/s3" grid_ref="grid_T_2D_inner" />
+ <field id="trd_eke_keS" long_name="KE to EKE source trend (RHS)" unit="m3/s3" grid_ref="grid_T_2D_inner" />
+ <field id="trd_eke_dis" long_name="dissipation trend of EKE (RHS)" unit="m3/s3" grid_ref="grid_T_2D_inner" />
+
<!-- variables available with ln_zdfmfc=.true. -->
<field id="mf_Tp" long_name="plume_temperature" standard_name="plume_temperature" unit="degC" grid_ref="grid_T_3D_inner" />
<field id="mf_Sp" long_name="plume_salinity" standard_name="plume_salinity" unit="1e-3" grid_ref="grid_T_3D_inner" />
@@ -754,6 +763,13 @@ that are available in the tidal-forcing implementation (see
<!-- EOS -->
<field id="bn2" long_name="squared Brunt-Vaisala frequency" unit="s-2" />
+ <!-- GEOMETRIC fields (requires nn_aei_ijk_t = 32) -->
+ <field id="aeiv_geom" long_name="3D w-EIV coefficient from GEOMETRIC param." unit="m2/s" />
+ <field id="rossby_rad" long_name="internal Rossby deformation radius" unit="m" grid_ref="grid_W_2D_inner"/>
+ <field id="bn2" long_name="squared Brunt-Vaisala frequency" unit="s-1" />
+ <field id="c1_vert" long_name="1st baroclinic mode phase speed" unit="m/s" grid_ref="grid_W_2D_inner"/>
+ <field id="c_ros" long_name="long Rossby phase speed" unit="m/s" grid_ref="grid_W_2D"/>
+
<!-- dissipation diagnostics (note: ediss_k is only available with tke scheme) -->
<field id="avt_k" long_name="vertical eddy diffusivity from closure schemes" standard_name="ocean_vertical_eddy_diffusivity" unit="m2/s" grid_ref="grid_W_3D_inner" />
<field id="avm_k" long_name="vertical eddy viscosity from closure schemes" standard_name="ocean_vertical_eddy_viscosity" unit="m2/s" />
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index e52985c3..49000349 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -980,11 +980,42 @@
! ! = 20 F(i,j) =ldf_c2d
! ! = 21 F(i,j,t) =Treguier et al. JPO 1997 formulation
! ! = 30 F(i,j,k) =ldf_c2d * ldf_c1d
+ ! ! = 32 F(i,j,t) = GEOMETRIC parameterization (=> fill namldf_eke)
! ! time invariant coefficients: aei0 = 1/2 Ue*Le
rn_Ue = 0.02 ! lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30)
rn_Le = 200.e+3 ! lateral diffusive length [m] (nn_aht_ijk_t= 0, 10)
!
ln_ldfeiv_dia =.false. ! diagnose eiv stream function and velocities
+ ln_eke_equ =.false. ! switch on update of GEOMETRIC eddy energy equation (=> fill namldf_eke)
+ ! forced to be true if nn_aei_ijk_t = 32
+/
+!----------------------------------------------------------------------------------
+&namldf_eke ! GEOMETRIC param. (total EKE equation) (nn_aei_ijk_t = 32)
+!----------------------------------------------------------------------------------
+ rn_ekedis = 180. ! dissipation time scale of EKE [days]
+ nn_eke_dis = 0 ! dissipation option
+ ! ! = 0 constant in space
+ ! ! =-20 read in geom_diss_2D.nc file
+ rn_geom = 0.05 ! geometric parameterization master coefficient (>0 & <1)
+ rn_eke_init = 1.e-1 ! initial total EKE value
+ rn_eke_min = 1.e+0 ! background value of total EKE
+ rn_ross_min = 4.e+3 ! tapering of aeiv based on min Rossby radius [m]
+ ! ! set to zero to not taper it
+ rn_eke_lap = 500. ! Laplacian diffusion coefficient of EKE
+ ! ! this is in all options below, so set it to zero and nothing is done
+ rn_aeiv_min = 1.e+1 ! minimum bound of eiv coefficient
+ rn_aeiv_max = 1.5e+4 ! maximum bound of eiv coefficient
+ rn_SFmin = 1.0 ! minimum bound of Structure Function
+ rn_SFmax = 1.0 ! maximum bound of Structure Function
+ nn_eke_opt = 1 ! options for terms to include in EKE budget
+ ! ! = 0 PE->EKE conversion, dissipation only
+ ! ! = 1 as 0 but with advection
+ ! ! = 2 as 1 but with additional KE->EKE conversion
+ ! ! for testing purposes:
+ ! ! = 88 only advection by depth-averaged flow
+ ! ! = 99 only Laplacian diffusion
+ ln_adv_wav = .false. ! include advection at long Rossby speed
+ ln_beta_plane = .false. ! beta plane option for computing long Rossby speed (default: sphere option)
/
!-----------------------------------------------------------------------
&namtra_dmp ! tracer: T & S newtonian damping (default: OFF)
diff --git a/sette/prepare_job.sh b/sette/prepare_job.sh
index 2c76c99d..1a95f861 100755
--- a/sette/prepare_job.sh
+++ b/sette/prepare_job.sh
@@ -204,6 +204,9 @@ fi
X86_ARCHER2-Gnu*)
MK_TEMPLATE=$( /work/n01/shared/nemo/mkslurm_settejob_4.2_Gnu -S $NXIO_PROC -s 8 -m 4 -C $NB_PROC -g 2 -a n01-CLASS -j sette_job -t 20:00 > ${SETTE_DIR}/job_batch_template )
;;
+ ANEMONE-ifort*)
+ MK_TEMPLATE=$( /home/acc/mkslurm_settejob_4.2 -S $NXIO_PROC -s 8 -m 4 -C $NB_PROC -g 2 -j sette_job -t 20:00 > ${SETTE_DIR}/job_batch_template )
+ ;;
XC40_METO*) #Setup for Met Office XC40 with any compiler
# ocean cores are packed 32 to a node
# If we need more than one node then have to use parallel queue and XIOS must have a node to itself
diff --git a/src/OCE/DYN/dynldf_lev.F90 b/src/OCE/DYN/dynldf_lev.F90
index 863151cb..07180003 100644
--- a/src/OCE/DYN/dynldf_lev.F90
+++ b/src/OCE/DYN/dynldf_lev.F90
@@ -18,6 +18,8 @@ MODULE dynldf_lev
USE domutl, ONLY : is_tile
USE ldfdyn ! lateral diffusion: eddy viscosity coef.
USE ldfslp ! iso-neutral slopes
+ USE ldftra , ONLY : l_ldfeke ! GEOMETRIC param. activation
+ USE ldfeke , ONLY : eke_keS, nn_eke_opt ! GEOMETRIC source term of eke due to KE dissipation
USE zdf_oce ! ocean vertical physics
!
USE in_out_manager ! I/O manager
@@ -59,13 +61,14 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), DIMENSION(T2D(1)) :: zwf, zwt
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zah_cur2, zah_div2 ! temporaries used to calculate GEOMTERIC source term
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
- IF( kt == nit000 ) THEN
- IF(lwp) WRITE(numout,*)
- IF(lwp) WRITE(numout,*) 'dynldf_lev_lap : laplacian operator momentum '
- IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
+ IF( kt == nit000 .AND. lwp ) THEN
+ WRITE(numout,*)
+ WRITE(numout,*) 'dynldf_lev_lap : laplacian operator momentum '
+ WRITE(numout,*) '~~~~~~~~~~~~~~'
ENDIF
ENDIF
!
@@ -78,6 +81,12 @@ CONTAINS
!
CASE ( np_typ_rot ) !== Vorticity-Divergence operator ==!
!
+ IF( l_ldfeke .AND. nn_eke_opt == 2 ) THEN ! GEOMETRIC source term
+ ALLOCATE( zah_cur2(T2D(1)) , zah_div2(T2D(1)) )
+ zah_cur2(:,:) = 0._wp
+ zah_div2(:,:) = 0._wp
+ ENDIF
+ !
# define zcur zwf
# define zdiv zwt
!
@@ -88,6 +97,17 @@ CONTAINS
# undef zcur
# undef zdiv
!
+ IF( l_ldfeke .AND. nn_eke_opt == 2 ) THEN ! GEOMETRIC source term
+ zah_cur2(T2D(1)) = zah_cur2(T2D(1)) * e1e2f(T2D(1))
+ DO_2D( 0, 0, 0, 0 )
+ eke_keS(ji,jj) = zah_div2(ji,jj) + ( ( zah_cur2(ji-1,jj ) + zah_cur2(ji,jj ) ) & ! add () for NP repro
+ & + ( zah_cur2(ji-1,jj-1) + zah_cur2(ji,jj-1) ) ) & ! add () for NP repro
+ & / MAX( 1._wp , fmask (ji-1,jj ,1) + fmask (ji,jj ,1) &
+ & + fmask (ji-1,jj-1,1) + fmask (ji,jj-1,1) ) * r1_e1e2t(ji,jj)
+ END_2D
+ DEALLOCATE( zah_cur2 , zah_div2 )
+ ENDIF
+ !
CASE ( np_typ_sym ) !== Symmetric operator ==!
!
# define zshe zwf
diff --git a/src/OCE/DYN/dynldf_lev_rot_scheme.h90 b/src/OCE/DYN/dynldf_lev_rot_scheme.h90
index 81a3841e..f96cc2e7 100644
--- a/src/OCE/DYN/dynldf_lev_rot_scheme.h90
+++ b/src/OCE/DYN/dynldf_lev_rot_scheme.h90
@@ -51,3 +51,14 @@
& ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) &
& + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) )
END_2D
+ !
+#if defined lap
+ IF( l_ldfeke .AND. nn_eke_opt == 2 ) THEN ! GEOMETRIC source term
+ DO_2D( INN, INN+1, INN, INN+1 )
+ zah_cur2(ji-1,jj-1) = zah_cur2(ji-1,jj-1) + &
+ & zcur(ji-1,jj-1)**2 / MAX( 1._wp , ahmf(ji-1,jj-1,jk) ) * fmask(ji-1,jj-1,jk)
+ zah_div2(ji ,jj ) = zah_div2(ji ,jj ) + e3t(ji ,jj ,jk,Kbb) * &
+ & zdiv(ji ,jj )**2 / MAX( 1._wp , ahmt(ji ,jj ,jk) ) * tmask(ji,jj,jk)
+ END_2D
+ ENDIF
+#endif
diff --git a/src/OCE/LDF/ldfdyn.F90 b/src/OCE/LDF/ldfdyn.F90
index 9f9bf187..02056548 100644
--- a/src/OCE/LDF/ldfdyn.F90
+++ b/src/OCE/LDF/ldfdyn.F90
@@ -17,6 +17,7 @@ MODULE ldfdyn
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
USE ldfslp ! lateral diffusion: slopes of mixing orientation
+ USE ldftra , ONLY : ln_eke_equ ! for compatability check only (GEOMETRIC)
USE ldfc1d_c2d ! lateral diffusion: 1D and 2D cases
!
USE in_out_manager ! I/O manager
@@ -170,8 +171,11 @@ CONTAINS
IF(.NOT.ln_dynldf_OFF ) THEN !== direction ==>> type of operator ==!
!
SELECT CASE( nn_dynldf_typ ) ! div-rot or symmetric
- CASE( np_typ_rot ) ; IF(lwp) WRITE(numout,*) ' ==>>> use div-rot operator '
- CASE( np_typ_sym ) ; IF(lwp) WRITE(numout,*) ' ==>>> use symmetric operator '
+ CASE( np_typ_rot )
+ IF(lwp) WRITE(numout,*) ' ==>>> use div-rot operator '
+ CASE( np_typ_sym )
+ IF(lwp) WRITE(numout,*) ' ==>>> use symmetric operator '
+ IF( ln_eke_equ ) CALL ctl_stop('STOP', 'ldf_dyn_init : GEOMETRIC parameterisation is only available with the div-rot operator')
CASE DEFAULT ! error
CALL ctl_stop('ldf_dyn_init: wrong value for nn_dynldf_typ (0 or 1)' )
END SELECT
diff --git a/src/OCE/LDF/ldfeke.F90 b/src/OCE/LDF/ldfeke.F90
new file mode 100644
index 00000000..082e1111
--- /dev/null
+++ b/src/OCE/LDF/ldfeke.F90
@@ -0,0 +1,610 @@
+MODULE ldfeke
+ !!======================================================================
+ !! *** MODULE ldfeke ***
+ !! Ocean physics: Eddy induced velocity coefficient according to the
+ !! GEOMETRIC parameterization scheme
+ !!=====================================================================
+ !! History : 4.0 ! 2017-11 (J. Mak, G. Madec) original code
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! ldf_eke : time step depth-integrated EKE and update aeiw (and
+ !! from that aeiu and aeiv) according to the GEOMETRIC
+ !! parameterization scheme
+ !! ldf_eke_init : initialization, namelist read, and parameters control
+ !! ldf_eke_rst : read/write eke restart in ocean restart file
+ !!----------------------------------------------------------------------
+ USE oce ! ocean: dynamics and active tracers variables
+ USE phycst ! physical constants
+ USE dom_oce ! domain: ocean
+ USE ldfslp ! lateral physics: slope of iso-neutral surfaces
+ USE ldftra ! lateral physics: eddy coefficients
+ USE zdfmxl ! vertical physics: mixed layer
+
+ !
+ USE in_out_manager ! I/O manager
+ USE iom ! I/O manager library
+ USE lib_mpp ! MPP library
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
+ USE prtctl ! Print control
+ USE timing ! Timing
+ USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+
+ IMPLICIT NONE
+ PRIVATE
+
+ PUBLIC ldf_eke ! routine called in step module
+ PUBLIC ldf_eke_init ! routine called in opa module
+
+
+ ! !!** Namelist namldf_eke **
+ REAL(wp) :: rn_ekedis ! dissipation time scale of EKE [days]
+ REAL(wp) :: rn_geom ! geometric parameterization master coefficient [-]
+ REAL(wp) :: rn_eke_lap ! diffusion of EKE [m2/s]
+ REAL(wp) :: rn_eke_init ! initial value of total EKE [m2/s2]
+ REAL(wp) :: rn_eke_min ! background value of total EKE [m2/s2]
+ REAL(wp) :: rn_ross_min ! tapering based of minimum Rossby radius [m]
+ REAL(wp) :: rn_aeiv_min, rn_aeiv_max ! min and max bounds of aeiv coefficient [m2/s]
+ REAL(wp) :: rn_SFmin, rn_SFmax ! min and max bounds of Structure Function [-]
+ REAL(wp) :: zf0, zbeta ! f0 and beta for computing Rossby speed
+ !
+ INTEGER, PUBLIC :: nn_eke_opt ! option for what term to include in eddy energy budget
+ INTEGER :: nn_eke_dis ! option for taking constant or spatially varying linear dissipation
+ !
+ LOGICAL :: ln_adv_wav ! option for having advection by Rossby wave or not
+ LOGICAL :: ln_beta_plane ! option for computing long Rossby phase speed
+ !
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: eke_geom ! vertical sum of total Eddy Kinetic Energy [m3/s2]
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ekedis ! linear dissipation rate (= 1/rn_ekedis) [ /s ]
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zc1, zc_ros, zadv_wav ! 1st baroclinic mode and long Rossby speed [m /s ]
+ !
+ REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: eke_keS ! source term of eke due to KE dissipation
+
+ !! * Substitutions
+# include "do_loop_substitute.h90"
+# include "domzgr_substitute.h90"
+ !!----------------------------------------------------------------------
+ !! NEMO/OPA 4.0 , NEMO Consortium (2017)
+ !! $Id: ldfeke.F90 7813 2017-03-20 16:17:45Z clem $
+ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE ldf_eke( kt, Kbb, Kmm, Kaa )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE ldf_eke ***
+ !!
+ !! ** Purpose : implements the GEOMETRIC parameterization
+ !!
+ !! ** Notes : If nn_aei_ijk_t = 32 then eke and aeiw are BOTH updated
+ !! If ln_eke_equ = .true. in namtra_ldfeiv but nn_aei_ijk_t
+ !! is something else, then ONLY eddy equation is updated
+ !! (but the eddy energy is passive and doesn't do anything)
+ !!
+ !! ** Method : GEOMETRIC calculates the Gent-McWilliams / eddy induced
+ !! velocity coefficient according to
+ !!
+ !! aeiw = alpha * (\int E dz) / (\int S M^2 / N dz),
+ !!
+ !! where (\int E dz) is the depth-integrated eddy energy
+ !! (at the previous time level), informed by a parameterized
+ !! depth-integrated eddy energy, where
+ !!
+ !! nn_eke_opt = 0 => default: just PE->EKE growth and linear dissipation
+ !! !
+ !! = 1 => default + advection by depth-averaged flow
+ !! !
+ !! = 2 => default + advection + contribution to EKE from
+ !! momentum dissipation
+ !! !
+ !! = 88 => ONLY advection
+ !! !
+ !! = 99 => ONLY Laplacian diffusion
+ !!
+ !! S is a structure function, and M and N are horizontal and
+ !! vertical buoyancy frequencies
+ !!
+ !! linear dissipation may be specified by
+ !!
+ !! nn_eke_dus = 0 => constant
+ !! !
+ !! =-20 => read in a geom_diss_2D.nc field (give it in days)
+ !!
+ !! ** Action : * time step depth-integrated eddy energy budget
+ !! * calculate aeiw
+ !!
+ !! References : Marshall, Maddison, Berloff JPO 2012 ; Mak et al JPO 2018
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT(in) :: kt ! ocean time step
+ INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! time-level indicators
+ !
+ INTEGER :: ji, jj, jk ! dummy loop arguments
+ INTEGER :: ik ! local integer
+ REAL(wp) :: zn2_min = 1.e-8_wp ! minimum value of N2 used to compute structure function
+ REAL(wp) :: z1_f20, zfw, ze3w ! local scalar
+ REAL(wp) :: zfp_ui, zfp_vj ! - -
+ REAL(wp) :: zfm_ui, zfm_vj ! - -
+ REAL(wp) :: zn_slp2, zn2 ! - -
+ REAL(wp) :: zmsku, zaeiu_w ! - -
+ REAL(wp) :: zmskv, zaeiv_w ! - -
+ REAL(wp) :: zen, zed ! - -
+ REAL(wp) :: zeke_rhs ! - -
+ REAL(wp) :: zck, zwslpi, zwslpj ! - - tapering near coasts
+ REAL(wp) :: zc_rosu ! - -
+ REAL(wp), DIMENSION(A2D(0)) :: zeke_peS, zn_slp ! 2D workspace, PE-KE conversion
+ REAL(wp), DIMENSION(A2D(0)) :: zadv_ubt ! - - barotropic advection
+ REAL(wp), DIMENSION(A2D(0)) :: zlap ! - - diffusion
+ REAL(wp), DIMENSION(A2D(0)) :: zdis ! - - linear dissipation
+ REAL(wp), DIMENSION(A2D(0)) :: zn, zross ! - - tapering near coasts
+ REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy ! - - barotropic advection
+ REAL(wp), DIMENSION(A2D(nn_hls)) :: zaheeu, zaheev ! - - diffusion
+ REAL(wp), DIMENSION(A2D(nn_hls)) :: zaeiw ! 2D EIV coefficient
+ REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwrk3d ! 3D workspace (structure function and then 3D EIV coeff)
+ !!--------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('ldf_eke')
+ !
+ IF( kt == nit000 .AND. lwp) THEN !* Control print
+ WRITE(numout,*)
+ WRITE(numout,*) 'ldf_eke : GEOMETRIC parameterization (total EKE time evolution equation)'
+ WRITE(numout,*) '~~~~~~~'
+ ENDIF
+
+ ! !== EIV mean conversion to EKE (ah N^2 slp2) & N slp ==!
+ !
+ ! work out the 3D structure function (SF) here and
+ ! store in 3d workspace
+ !
+ ! current: Ferreria et al, SF = N^2 / N^2_ref, W-points
+ ! capped between rn_SFmax and rn_SFmin
+ zwrk3d(:,:,:) = 0._wp
+ DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+ IF( jk <= nmln(ji,jj) ) THEN ! above and at mixed layer
+ zwrk3d(ji,jj,jk) = rn_SFmax
+ ELSE
+ ik = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! one level below the mixed layer (MIN in case ML depth is the ocean depth)
+ zwrk3d(ji,jj,jk) = MAX( 0._wp , rn2b(ji,jj,jk) ) / MAX( zn2_min , rn2b(ji,jj,ik) ) ! Structure Function : N^2 / N^2_ref
+ zwrk3d(ji,jj,jk) = MAX( rn_SFmin , MIN( zwrk3d(ji,jj,jk) , rn_SFmax ) )
+ ENDIF
+ END_3D
+ !
+ ! !* parametrized PE_EKE conversion due to eddy induced velocity
+ zeke_peS(:,:) = 0._wp
+ zn_slp (:,:) = 0._wp
+ zn (:,:) = 0._wp
+ DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+ zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) &
+ & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp )
+ zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &
+ & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp )
+ !
+ zaeiu_w = ( ( aeiu(ji ,jj,jk-1) + aeiu(ji-1,jj,jk) ) & ! add () for NP repro
+ & + ( aeiu(ji-1,jj,jk-1) + aeiu(ji ,jj,jk) ) ) * zmsku ! add () for NP repro
+ zaeiv_w = ( ( aeiv(ji,jj ,jk-1) + aeiv(ji,jj-1,jk) ) & ! add () for NP repro
+ & + ( aeiv(ji,jj-1,jk-1) + aeiv(ji,jj ,jk) ) ) * zmskv ! add () for NP repro
+ !
+ zn_slp2 = ( ( zaeiu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) ) & ! (slope ** 2) * aeiv + add () for NP repro
+ & + ( zaeiv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) ) ! JM 28 Jun: undo slope reduction here too?
+ zn2 = MAX( 0._wp , rn2b(ji,jj,jk) )
+ !
+ ze3w = e3w(ji,jj,jk,Kbb) * tmask(ji,jj,jk)
+ zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * ze3w ! for working out taper at small rossby radius regions later
+ !
+ zeke_peS(ji,jj) = zeke_peS(ji,jj) + ze3w * zn2 * zn_slp2 ! note this is >=0
+ !
+ zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & ! taken from ldfslp, undo the slope reduction
+ & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25_wp ! near topographic features
+ zwslpi = wslpi(ji,jj,jk) / MAX( zck, 0.1_wp) ! (just to avoid dividing by zeros)
+ zwslpj = wslpj(ji,jj,jk) / MAX( zck, 0.1_wp)
+ !
+ zn_slp(ji,jj) = zn_slp(ji,jj) + ze3w * zwrk3d(ji,jj,jk) & ! note this >=0 and structure function weighted
+ & * SQRT( zn2 * ( zwslpi * zwslpi + zwslpj * zwslpj ) )
+ END_3D
+
+ ! !* upstream advection with initial mass fluxes & intermediate update
+ ! !* upstream tracer flux in the i and j direction
+ IF( .NOT. ln_linssh ) THEN
+ DO_2D( 1, 0, 1, 0 )
+ ! upstream scheme
+ zfp_ui = un_adv(ji,jj) + ABS( un_adv(ji,jj) )
+ zfm_ui = un_adv(ji,jj) - ABS( un_adv(ji,jj) )
+ zfp_vj = vn_adv(ji,jj) + ABS( vn_adv(ji,jj) )
+ zfm_vj = vn_adv(ji,jj) - ABS( vn_adv(ji,jj) )
+ zwx(ji,jj) = 0.5_wp * ( zfp_ui * eke_geom(ji,jj,Kbb) + zfm_ui * eke_geom(ji+1,jj ,Kbb) )
+ zwy(ji,jj) = 0.5_wp * ( zfp_vj * eke_geom(ji,jj,Kbb) + zfm_vj * eke_geom(ji ,jj+1,Kbb) )
+ END_2D
+ ! !* divergence of ubt advective fluxes
+ DO_2D( 0, 0, 0, 0 )
+ zadv_ubt(ji,jj) = - ( ( zwx(ji,jj) - zwx(ji-1,jj ) ) & ! add () for NP repro
+ & + ( zwy(ji,jj) - zwy(ji ,jj-1) ) ) * r1_e1e2t(ji,jj) ! add () for NP repro
+ END_2D
+ ELSE !* top value (linear free surf. only as zwz is multiplied by wmask)
+ DO_2D( 0, 0, 0, 0 )
+ zadv_ubt(ji,jj) = - ww(ji,jj,1) * eke_geom(ji,jj,Kbb) / ( ht_0(ji,jj) + 1._wp-tmask(ji,jj,1) ) ! jm (03 Mar 18): was eke_n
+ END_2D
+ ENDIF
+ !
+ ! !* same as above but for advection by Rossby waves
+ zadv_wav(:,:) = 0._wp
+ IF ( ln_adv_wav ) THEN
+ !
+ DO_2D( 0, 0, 0, 0 )
+ ! compute only for deep enough places
+ IF ( ht_0(ji,jj) > 300._wp ) THEN ! jm: should use something like ht_b really...
+ ! compute vertical mode phase speed on T point
+ zc1(ji,jj) = MIN( 10._wp, zn(ji,jj) / rpi )
+ ! compute long Rossby phase speed on T point (minus sign later)
+ IF ( ln_beta_plane ) THEN
+ zc_ros(ji,jj) = zc1(ji,jj) * zc1(ji,jj) * zbeta / (zf0 * zf0)
+ ELSE
+ zc_ros(ji,jj) = zc1(ji,jj) * zc1(ji,jj) * COS( rad * gphit(ji,jj) ) &
+ / ( ra * ff_t(ji,jj) * SIN( rad * gphit(ji,jj) ) &
+ + rsmall )
+ END IF
+ ! cap the Rossby phase speeds by the fastest equatorial Rossby wave speed
+ ! the minus sign for westward propagation goes here
+ zc_ros(ji,jj) = -MIN( zc1(ji,jj) / 3._wp, zc_ros(ji,jj) )
+ END IF
+ END_2D
+ !
+ CALL lbc_lnk( 'ldfeke', zc_ros, 'W', 1. )
+ zwx(:,:) = 0._wp ! wipe the advective contribution from above
+ !
+ DO_2D( 1, 0, 0, 0 )
+ zc_rosu = 0.5_wp * ( zc_ros(ji,jj) + zc_ros(ji+1,jj) ) * umask(ji,jj,1)
+ zfp_ui = zc_rosu + ABS( zc_rosu )
+ zfm_ui = zc_rosu - ABS( zc_rosu )
+ zwx(ji,jj) = 0.5_wp * ( zfp_ui * eke_geom(ji,jj,Kbb) + zfm_ui * eke_geom(ji+1,jj,Kbb) )
+ END_2D
+ ! !* divergence of wav advective fluxes (only e1 here)
+ z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) )
+ DO_2D( 0, 0, 0, 0 )
+ zadv_wav(ji,jj) = - ( zwx(ji,jj) - zwx(ji-1,jj) ) * r1_e1t(ji,jj)
+ zadv_wav(ji,jj) = zadv_wav(ji,jj) * MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) ! tropical decrease
+ END_2D
+ END IF
+ !
+ !* divergence of diffusive fluxes
+ ! code adapted from traldf_lap_blp.F90
+ IF ( rn_eke_lap >= 0._wp ) THEN
+ DO_2D( 1, 0, 1, 0 )
+ zaheeu(ji,jj) = rn_eke_lap * e2_e1u(ji,jj) * umask(ji,jj,1) ! rn_eke_lap is constant (for now) and NOT masked
+ zaheev(ji,jj) = rn_eke_lap * e1_e2v(ji,jj) * vmask(ji,jj,1) ! before it is pahu and pahv which IS masked
+ END_2D
+ DO_2D( 0, 0, 0, 0 )
+ zlap(ji,jj) = ( ( zaheeu(ji,jj) * eke_geom(ji+1,jj ,Kbb) + zaheeu(ji-1,jj ) * eke_geom(ji-1,jj ,Kbb) ) &
+ & - ( zaheeu(ji,jj) * eke_geom(ji ,jj ,Kbb) + zaheeu(ji-1,jj ) * eke_geom(ji ,jj ,Kbb) ) &
+ & + ( zaheev(ji,jj) * eke_geom(ji ,jj+1,Kbb) + zaheev(ji ,jj-1) * eke_geom(ji ,jj-1,Kbb) ) &
+ & - ( zaheev(ji,jj) * eke_geom(ji ,jj ,Kbb) + zaheev(ji ,jj-1) * eke_geom(ji ,jj ,Kbb) ) ) &
+ & / e1e2t(ji,jj)
+ END_2D
+ ELSE
+ zlap(:,:) = 0._wp
+ END IF
+ !* form the trend for linear dissipation
+ DO_2D( 0, 0, 0, 0 )
+ zdis(ji,jj) = - r1_ekedis(ji,jj) * (eke_geom(ji,jj,Kbb) - rn_eke_min) * tmask(ji,jj,1)
+ END_2D
+ ! !== time stepping of EKE Eq. ==!
+ !
+ ! note: the rn_eke_min term is a forcing in eddy equation, thus damping in mean equation
+ ! added to prevent overshoots and oscillations of energy from exponential growth/decay
+ ! maintains a background (depth-integrated) eddy energy level
+ !
+ zeke_rhs = 0._wp
+ DO_2D( 0, 0, 0, 0 )
+ !
+ SELECT CASE( nn_eke_opt ) ! Specification of what to include in EKE budget
+ !
+ CASE( 0 ) ! default: just PE->EKE growth and linear dissipation
+ zeke_rhs = zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj)
+ CASE( 1 ) ! as default but with full advection
+ zeke_rhs = - zadv_ubt(ji,jj) + zadv_wav(ji,jj) + zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj)
+ CASE( 2 ) ! full thing with additional KE->EKE growth
+ zeke_rhs = - zadv_ubt(ji,jj) + zadv_wav(ji,jj) + zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj) + eke_keS(ji,jj)
+ CASE( 88 ) ! ONLY advection by mean flow
+ zeke_rhs = - zadv_ubt(ji,jj)
+ CASE( 99 ) ! ONLY diffusion
+ zeke_rhs = zlap(ji,jj)
+ CASE DEFAULT
+ CALL ctl_stop('ldf_eke: wrong choice nn_eke_opt, set at least to 0 (default)')
+ END SELECT
+ !
+ ! leap-frog
+ eke_geom(ji,jj,Kaa) = eke_geom(ji,jj,Kbb) + rDt * zeke_rhs * ssmask(ji,jj)
+ !
+ END_2D
+ CALL lbc_lnk( 'ldfeke', eke_geom(:,:,Kaa), 'T', 1._wp ) ! Lateral boundary conditions on eke_geom (unchanged sign)
+
+ IF( .NOT. ( l_1st_euler .AND. kt == nit000 ) ) THEN ! Apply Asselin filter except if Euler time-stepping at first time-step
+ DO_2D( 1, 1, 1, 1 )
+ !
+ zen = eke_geom(ji,jj,Kmm)
+ zed = eke_geom(ji,jj,Kaa) - 2._wp * zen + eke_geom(ji,jj,Kbb) ! time laplacian on tracers
+ !
+ eke_geom(ji,jj,Kmm) = zen + rn_atfp * zed ! filtered eke_n (will be eke_b after level-indicators are shuffled in stpmlf.F90)
+ END_2D
+ ENDIF
+
+ ! initialise it here so XIOS stops complaining...
+ zross(:,:) = 0._wp
+ zaeiw(:,:) = 0._wp
+ IF( l_eke_eiv ) THEN
+ ! !== resulting EIV coefficient ==!
+ !
+ ! ! surface value
+ z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) )
+ DO_2D( 0, 0, 0, 0 )
+ ! Rossby radius at w-point, Ro = .5 * sum_jpk(N) / f
+ zfw = MAX( ABS( ff_t(ji,jj) ) , 1.e-10_wp )
+ zross(ji,jj) = 0.5_wp * zn(ji,jj) / zfw
+ !
+ zaeiw(ji,jj) = rn_geom * eke_geom(ji,jj,Kmm) / MAX( 1.e-10_wp , zn_slp(ji,jj) ) * tmask(ji,jj,1) ! zn_slp has SF multiplied to it
+ zaeiw(ji,jj) = MIN( rn_aeiv_max, zaeiw(ji,jj) ) ! bound aeiv from above
+ zaeiw(ji,jj) = zaeiw(ji,jj) &
+ ! tanh taper to deal with some some large values near coast
+ & * 0.5_wp * ( 1._wp - TANH( ( -ht_0(ji,jj) + 800._wp ) / 300._wp ) ) &
+ ! tanh taper of aeiv on internal Rossby radius
+ & * 0.5_wp * ( 1._wp + TANH( ( zross(ji,jj) - rn_ross_min ) * 0.5_wp ) )
+ zaeiw(ji,jj) = zaeiw(ji,jj) * MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) ! tropical decrease
+ zaeiw(ji,jj) = MAX( rn_aeiv_min, zaeiw(ji,jj) ) ! bound aeiv from below
+ END_2D
+ CALL lbc_lnk( 'ldfeke', zaeiw(:,:), 'W', 1. )
+ !
+ ! ! inner value
+ ! bottom value is already set to zero, use the un-masked zaeiw(ji,jj) with the structure function to
+ ! set interior values. Re-purpose 3D workspace, replacing structure function (SF) with 3D eiv coefficient
+ zwrk3d(:,:,1) = zaeiw(:,:)
+ DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+ zwrk3d(ji,jj,jk) = zaeiw(ji,jj) * zwrk3d(ji,jj,jk) * wmask(ji,jj,jk) ! SF has already been capped
+ END_3D
+ ! ! aei at u- and v-points
+ aeiu(:,:,jpk) = 0._wp
+ aeiv(:,:,jpk) = 0._wp
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+ aeiu(ji,jj,jk) = ( ( zwrk3d(ji,jj,jk ) + zwrk3d(ji+1,jj,jk ) ) & ! add () for NP repro
+ & + ( zwrk3d(ji,jj,jk+1) + zwrk3d(ji+1,jj,jk+1) ) ) & ! add () for NP repro
+ & / MAX( wmask(ji,jj,jk ) + wmask(ji+1,jj,jk ) &
+ & + wmask(ji,jj,jk+1) + wmask(ji+1,jj,jk+1) , 1._wp ) * umask(ji,jj,jk)
+ aeiv(ji,jj,jk) = ( ( zwrk3d(ji,jj,jk ) + zwrk3d(ji,jj+1,jk ) ) & ! add () for NP repro
+ & + ( zwrk3d(ji,jj,jk+1) + zwrk3d(ji,jj+1,jk+1) ) ) & ! add () for NP repro
+ & / MAX( wmask(ji,jj,jk ) + wmask(ji,jj+1,jk ) &
+ & + wmask(ji,jj,jk+1) + wmask(ji,jj+1,jk+1) , 1._wp ) * vmask(ji,jj,jk)
+ END_3D
+ ! !-- diagnostics --!
+ CALL lbc_lnk( 'ldfeke', aeiu , 'U', 1._wp , aeiv , 'V', 1._wp )
+ END IF
+ !
+ IF( lrst_oce .AND. kt == nitrst ) CALL ldf_eke_rst( kt, 'WRITE', Kmm, Kaa ) ! write Kmm and Kaa since time level
+ ! indicators will be shuffled before the
+ ! rest of the restart is written
+
+ ! !== output EKE related variables ==!
+ !
+ CALL iom_put( "eke" , eke_geom(:,:,Kaa) ) ! parameterized total EKE (EPE+ EKE)
+ CALL iom_put( "trd_eke_adv_ubt", zadv_ubt ) ! ubt advective trend of EKE(LHS)
+ CALL iom_put( "trd_eke_adv_wav", zadv_wav ) ! rossby advective trend of EKE(LHS)
+ CALL iom_put( "trd_eke_dis" , zdis ) ! dissipative trend of EKE (RHS)
+ CALL iom_put( "trd_eke_lap" , zlap ) ! diffusive trend of EKE (RHS)
+ CALL iom_put( "trd_eke_peS" , zeke_peS ) ! PE to EKE source trend (RHS)
+ CALL iom_put( "trd_eke_keS" , eke_keS ) ! KE to EKE source trend (RHS)
+ CALL iom_put( "aeiv_geom" , zwrk3d ) ! eddy induced coefficient from GEOMETRIC param
+ CALL iom_put( "rossby_rad" , zross ) ! internal Rossby deformation radius
+ CALL iom_put( "c1_vert" , zc1 ) ! 1st baroclinic mode phase speed
+ CALL iom_put( "c_ros" , zc_ros ) ! long Rossby phase speed
+!!
+!! jm: modified ldf_tra in ldftra.F90 to include the ln_eke_equ flag
+!! into consideration, otherwise iom_put is called twice aeiu and aeiv
+!!
+ CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff.
+ CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff.
+ CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff.
+ CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff.
+ !
+ IF( ln_timing ) CALL timing_stop('ldf_eke')
+ !
+ END SUBROUTINE ldf_eke
+
+ SUBROUTINE ldf_eke_init( Kbb, Kmm )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE ldf_eke_init ***
+ !!
+ !! ** Purpose : Initialization of the depth-integrated eke, vertical
+ !! structure function and max/min caps of aeiw when using the
+ !! GEOMETRIC parameterization
+ !!
+ !! ** Method : Read the namldf_eke namelist and check the parameters
+ !! called at the first timestep (nit000)
+ !!
+ !! ** input : Namlist namldf_eke
+ !!----------------------------------------------------------------------
+ INTEGER, INTENT( in ) :: Kbb, Kmm ! time level indicator
+ INTEGER :: ios, inum
+ INTEGER :: ji, jj
+ INTEGER :: ierr
+ !!
+ NAMELIST/namldf_eke/ rn_ekedis , nn_eke_dis , & ! GEOM master params (lambda and option, alpha, e0)
+ & rn_geom , rn_eke_init, &
+ & rn_eke_lap , & ! coeff of Laplacian diffusion of EKE
+ & rn_eke_min , & ! background value of (depth-integrated) EKE
+ & rn_ross_min, & ! taper aeiv based on Rossby internal radius
+ & rn_aeiv_max, rn_aeiv_min, & ! caps and options on result aeiv
+ & rn_SFmin, rn_SFmax, & ! caps on vertical structure function
+ & nn_eke_opt , & ! option for EKE budget terms
+ & ln_adv_wav , & ! option for advection by Rossby wave
+ & ln_beta_plane ! beta plane option for computing Rossby wave speed
+ !!----------------------------------------------------------------------
+ !
+ !
+ READ ( numnam_ref, namldf_eke, IOSTAT = ios, ERR = 901)
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namldf_eke in reference namelist' )
+
+ READ ( numnam_cfg, namldf_eke, IOSTAT = ios, ERR = 902 )
+902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namldf_eke in configuration namelist' )
+ IF(lwm) WRITE ( numond, namldf_eke )
+ !
+ IF(lwp) THEN !* Control print
+ WRITE(numout,*)
+ WRITE(numout,*) 'ldf_eke_init : GEOMETRIC parameterization (total EKE time evolution equation)'
+ WRITE(numout,*) '~~~~~~~~~~~~'
+ WRITE(numout,*) ' Namelist namldf_eke : set the GEOMETRIC parameters'
+ WRITE(numout,*) ' aeiw updated according to GEOMETRIC l_eke_eiv = ', l_eke_eiv
+ WRITE(numout,*) ' dissipation time scale of EKE in days rn_ekedis = ', rn_ekedis, ' days'
+ WRITE(numout,*) ' option for dissipation field nn_eke_dis = ', nn_eke_dis
+ SELECT CASE( nn_eke_dis )
+ !
+ CASE( 0 )
+ WRITE(numout,*) ' spatially constant '
+ CASE( 2 )
+ WRITE(numout,*) ' read in a 2d varying file geom_diss_2D.nc '
+ END SELECT
+ !
+ WRITE(numout,*) ' geometric parameterization master coefficient rn_geom = ', rn_geom
+ WRITE(numout,*) ' coeff of Lapican diffusion of EKE rn_eke_lap = ', rn_eke_lap
+ WRITE(numout,*) ' initial total EKE value rn_eke_init = ', rn_eke_init
+ WRITE(numout,*) ' background values of total EKE rn_eke_min = ', rn_eke_min
+ WRITE(numout,*) ' taper aeiv subject to min Rossby radius rn_ross_min = ', rn_ross_min
+ WRITE(numout,*) ' maximum bound of aeiv coeff rn_aeiv_max = ', rn_aeiv_max
+ WRITE(numout,*) ' minimum bound of aeiv coeff rn_aeiv_min = ', rn_aeiv_min
+ WRITE(numout,*) ' minimum bound of Structure Function rn_SFmin = ', rn_SFmin
+ WRITE(numout,*) ' maximum bound of Structure Function rn_SFmax = ', rn_SFmax
+ WRITE(numout,*) ' [set rnSFmin = rnSFmax = 1 to use aeiv = aeiv(x,y,t)] '
+ WRITE(numout,*) ' option for terms in eddy energy budget nn_eke_opt = ', nn_eke_opt
+ WRITE(numout,*) ' default: just PE->EKE growth and linear dissipation '
+ SELECT CASE( nn_eke_opt )
+ !
+ CASE( 0 )
+ WRITE(numout,*) ' default '
+ CASE( 1 )
+ WRITE(numout,*) ' default + advection '
+ CASE( 2 )
+ WRITE(numout,*) ' default + advection + KE->EKE '
+#if defined key_loop_fusion
+ CALL ctl_stop( 'STOP', 'ldf_eke_init : This option of the GEOMETRIC parameterisation is not available with loop fusion' , &
+ & ' The source term for KE (eke_keS) is not set in the loop-fused version of dynldf_lap_blp' )
+#endif
+ CASE( 88 )
+ WRITE(numout,*) ' ONLY advection by z-avg mean flow (no growth or dissipation)'
+ CASE( 99 )
+ WRITE(numout,*) ' ONLY Laplacian diffusion (no growth or dissipation)'
+ END SELECT
+ WRITE(numout,*) ' advection of energy by Rossby waves ln_adv_wav = ', ln_adv_wav
+ ENDIF
+ ! ! allocate eke arrays
+ ALLOCATE( eke_geom(A2D(nn_hls),jpt) , eke_keS(A2D(0)) , r1_ekedis(A2D(0)), &
+ & zc1(A2D(0)) , zc_ros(A2D(nn_hls)) , zadv_wav(A2D(0)), &
+ & STAT=ierr )
+ IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ldf_eke_init: failed to allocate arrays')
+ !
+ eke_keS(:,:) = 0._wp ! To avoid unset values in unused array if nn_eke_opt /= 2
+ ! ! (array is unused but may still be mistakeningly requested in XML output)
+ !
+ SELECT CASE( nn_eke_dis ) ! Specification of linear dissipation
+ !
+ CASE( 0 ) !== constant ==!
+ IF(lwp) WRITE(numout,*) ' linear EKE dissipation coef. = constant = ', rn_ekedis, ' days'
+ r1_ekedis(:,:) = 1._wp / (rn_ekedis * rday)
+ !
+ CASE( -20 ) !== fixed horizontal shape read in file ==!
+ IF(lwp) WRITE(numout,*) ' linear EKE dissipation coef. = F(i,j) read in geom_diss_2D.nc file'
+ CALL iom_open ( 'geom_diss_2D.nc', inum )
+ CALL iom_get ( inum, jpdom_auto, 'r1_ekedis', r1_ekedis, cd_type = 'T', psgn = 1._wp ) ! read in as time-scale in days...
+ CALL iom_close( inum )
+ DO_2D( 0, 0, 0, 0 )
+ r1_ekedis(ji,jj) = 1._wp / (r1_ekedis(ji,jj) * rday) ! ...convert rate in per second
+ END_2D
+ !
+ CASE DEFAULT
+ CALL ctl_stop('ldf_eke_init: wrong choice for nn_eke_dis, the option for linear dissipation in GEOMETRIC')
+ END SELECT
+ !
+ CALL ldf_eke_rst( nit000, 'READ', Kbb, Kmm ) !* read or initialize all required files
+ !
+ ! initialise some optional arrays (ln_adv_wave related)
+ zc1(:,:) = 0._wp
+ zc_ros(:,:) = 0._wp
+ !
+ ! if using beta_plane, compute beta and f0 using values from the 1st domain
+ IF ( ln_beta_plane ) THEN
+ IF ( narea .eq. 1 ) THEN
+ zf0 = ff_f(1,1)
+ zbeta = (ff_f(1,2) - zf0) / e2t(1,1)
+ ELSE
+ zf0 = HUGE(1._wp)
+ zbeta = HUGE(1._wp)
+ ENDIF
+ CALL mpp_min( 'ldfeke', zf0 )
+ CALL mpp_min( 'ldfeke', zbeta )
+ IF(lwp) WRITE(numout,*) ' beta plane option is ln_beta_plane =', ln_beta_plane
+ IF(lwp) WRITE(numout,*) ' f0 = ', zf0, ' s-1'
+ IF(lwp) WRITE(numout,*) ' beta = ', zbeta, 'm-1 s-1'
+ END IF
+ !
+ !
+ END SUBROUTINE ldf_eke_init
+
+
+ SUBROUTINE ldf_eke_rst( kt, cdrw, Kbb, Kmm )
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE eke_rst ***
+ !!
+ !! ** Purpose : Read or write EKE file (eke) in restart file
+ !!
+ !! ** Method : use of IOM library
+ !! if the restart does not contain EKE, eke is set
+ !! according to rn_eke_init, and aeiu = aeiv = 10 m s^-2
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in) :: kt ! ocean time-step
+ INTEGER , INTENT(in) :: Kbb, Kmm ! time level indicator (needed for 'ht')
+ CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
+ !
+ INTEGER :: jit, jk ! dummy loop indices
+ INTEGER :: id1, id2, id3, id4, id5 ! local integer
+ !!----------------------------------------------------------------------
+ !
+ IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
+ ! ! ---------------
+ IF( ln_rstart ) THEN !* Read the restart file
+ id1 = iom_varid( numror, 'eke_b' , ldstop = .FALSE. )
+ id2 = iom_varid( numror, 'eke_n' , ldstop = .FALSE. )
+ id3 = iom_varid( numror, 'aeiu' , ldstop = .FALSE. )
+ id4 = iom_varid( numror, 'aeiv' , ldstop = .FALSE. )
+ IF(lwp) write(numout,*) ' ldfeke init restart'
+ !
+ IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist
+ IF(lwp) write(numout,*) ' all files for ldfeke exist, loading ...'
+ CALL iom_get( numror, jpdom_auto, 'eke_b', eke_geom(:,:,Kbb), cd_type = 'T', psgn = 1._wp )
+ CALL iom_get( numror, jpdom_auto, 'eke_n', eke_geom(:,:,Kmm), cd_type = 'T', psgn = 1._wp )
+ CALL iom_get( numror, jpdom_auto, 'aeiu' , aeiu , cd_type = 'U', psgn = 1._wp )
+ CALL iom_get( numror, jpdom_auto, 'aeiv' , aeiv , cd_type = 'V', psgn = 1._wp )
+ ELSE ! one at least array is missing
+ IF(lwp) write(numout,*) ' not all files for ldfeke exist '
+ IF(lwp) write(numout,*) ' --- initialize from namelist'
+ eke_geom(:,:,Kbb) = rn_eke_init * ht(:,:,Kmm) * ssmask(:,:)
+ eke_geom(:,:,Kmm) = eke_geom(:,:,Kbb)
+ aeiu(:,:,:) = 10._wp * umask(:,:,:) ! bottom eddy coeff set to zero at last level
+ aeiv(:,:,:) = 10._wp * vmask(:,:,:)
+ ENDIF
+ ELSE !* Start from rest with a non zero value (required)
+ IF(lwp) write(numout,*) ' ldfeke init restart from namelist'
+ eke_geom(:,:,Kbb) = rn_eke_init * ht(:,:,Kmm) * ssmask(:,:)
+ eke_geom(:,:,Kmm) = eke_geom(:,:,Kbb)
+ aeiu(:,:,:) = 10._wp * umask(:,:,:) ! bottom eddy coeff set to zero at last level
+ aeiv(:,:,:) = 10._wp * vmask(:,:,:)
+ !
+ ENDIF
+ !
+ ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
+ ! ! -------------------
+ IF(lwp) WRITE(numout,*) '---- eke-rst ----'
+ CALL iom_rstput( kt, nitrst, numrow, 'eke_b', eke_geom(:,:,Kbb) )
+ CALL iom_rstput( kt, nitrst, numrow, 'eke_n', eke_geom(:,:,Kmm) )
+ CALL iom_rstput( kt, nitrst, numrow, 'aeiu' , aeiu )
+ CALL iom_rstput( kt, nitrst, numrow, 'aeiv' , aeiv )
+ !
+ ENDIF
+ !
+ END SUBROUTINE ldf_eke_rst
+
+ !!======================================================================
+END MODULE ldfeke
diff --git a/src/OCE/LDF/ldftra.F90 b/src/OCE/LDF/ldftra.F90
index 74c0c67c..3bfeb2db 100644
--- a/src/OCE/LDF/ldftra.F90
+++ b/src/OCE/LDF/ldftra.F90
@@ -8,6 +8,7 @@ MODULE ldftra
!! 2.0 ! 2005-11 (G. Madec)
!! 3.7 ! 2013-12 (F. Lemarie, G. Madec) restructuration/simplification of aht/aeiv specification,
!! ! add velocity dependent coefficient and optional read in file
+ !! 4.3 ! 2023-02 (J. Mak, A. C. Coward, G. Madec) added GEOMETRIC parameterization
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -85,6 +86,10 @@ MODULE ldftra
INTEGER , PUBLIC :: nldf_tra = 0 !: type of lateral diffusion used defined from ln_traldf_... (namlist logicals)
LOGICAL , PUBLIC :: l_ldftra_time = .FALSE. !: flag for time variation of the lateral eddy diffusivity coef.
LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. !: flag for time variation of the eiv coef.
+
+ LOGICAL , PUBLIC :: ln_eke_equ !: flag for having updates to eddy energy equation
+ LOGICAL , PUBLIC :: l_ldfeke = .FALSE. !: GEOMETRIC - total EKE flag
+ LOGICAL , PUBLIC :: l_eke_eiv = .FALSE. !: GEOMETRIC - aeiw flag
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahtu, ahtv !: eddy diffusivity coef. at U- and V-points [m2/s]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aeiu, aeiv !: eddy induced velocity coeff. [m2/s]
@@ -449,7 +454,7 @@ CONTAINS
CALL iom_put( "ahtu_3d", ahtu(:,:,:) ) ! 3D u-eddy diffusivity coeff.
CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff.
!
- IF( ln_ldfeiv ) THEN
+ IF( ln_ldfeiv .AND. (.NOT. ln_eke_equ) ) THEN
CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff.
CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff.
CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff.
@@ -485,7 +490,8 @@ CONTAINS
REAL(wp) :: zah_max, zUfac ! - scalar
!!
NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv)
- & nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient
+ & nn_aei_ijk_t, rn_Ue, rn_Le , & ! eiv coefficient
+ & ln_eke_equ ! GEOMETRIC eddy energy equation
!!----------------------------------------------------------------------
!
IF(lwp) THEN ! control print
@@ -593,6 +599,16 @@ CONTAINS
CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! surface value proportional to scale factor^inn
CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) ! reduction with depth
!
+ CASE( 32 ) !-- time varying 3D field --!
+ IF(lwp) WRITE(numout,*) ' eddy induced velocity coef. = F( latitude, longitude, depth, time )'
+ IF(lwp) WRITE(numout,*) ' = F( total EKE ) GEOMETRIC parameterization'
+ !
+ IF ( lk_RK3 ) CALL ctl_stop('ldf_tra_init: The GEOMETRIC parameterisation is not yet available with RK3 time-stepping')
+ IF(lwp .AND. .NOT. ln_eke_equ ) WRITE(numout,*) ' ln_eke_equ will be set to .true. '
+ ln_eke_equ = .TRUE. ! force the eddy energy equation to be updated
+ l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90
+ l_eke_eiv = .TRUE.
+ !
CASE DEFAULT
CALL ctl_stop('ldf_tra_init: wrong choice for nn_aei_ijk_t, the type of space-time variation of aei')
END SELECT
@@ -606,6 +622,10 @@ CONTAINS
!
ENDIF
!
+ IF( ln_eke_equ ) THEN
+ l_ldfeke = .TRUE. ! GEOMETRIC param initialization done in nemogcm_init
+ IF(lwp) WRITE(numout,*) ' GEOMETRIC eddy energy equation to be computed ln_eke_equ = ', ln_eke_equ
+ ENDIF
END SUBROUTINE ldf_eiv_init
diff --git a/src/OCE/nemogcm.F90 b/src/OCE/nemogcm.F90
index 20fdf05c..6ac75ede 100644
--- a/src/OCE/nemogcm.F90
+++ b/src/OCE/nemogcm.F90
@@ -464,10 +464,11 @@ CONTAINS
! ! Ocean physics
CALL zdf_phy_init( Nnn ) ! Vertical physics
- ! ! Lateral physics
- CALL ldf_tra_init ! Lateral ocean tracer physics
- CALL ldf_eiv_init ! eddy induced velocity param. must be done after ldf_tra_init
- CALL ldf_dyn_init ! Lateral ocean momentum physics
+ ! ! Lateral physics
+ CALL ldf_tra_init ! Lateral ocean tracer physics
+ CALL ldf_eiv_init ! eddy induced velocity param.
+ IF( l_ldfeke ) CALL ldf_eke_init( Nbb, Nnn ) ! GEOMETRIC param.
+ CALL ldf_dyn_init ! Lateral ocean momentum physics
! ! Active tracers
IF( ln_traqsr ) CALL tra_qsr_init ! penetrative solar radiation qsr
diff --git a/src/OCE/step_oce.F90 b/src/OCE/step_oce.F90
index 36d6d6ba..e9afda4a 100644
--- a/src/OCE/step_oce.F90
+++ b/src/OCE/step_oce.F90
@@ -63,6 +63,7 @@ MODULE step_oce
USE ldfslp ! iso-neutral slopes (ldf_slp routine)
USE ldfdyn ! lateral eddy viscosity coef. (ldf_dyn routine)
USE ldftra ! lateral eddy diffusive coef. (ldf_tra routine)
+ USE ldfeke ! GEOMETRIC parameterization (ldf_eke routine)
USE zdf_oce ! ocean vertical physics variables
USE zdfphy ! vertical physics manager (zdf_phy_init routine)
diff --git a/src/OCE/stpmlf.F90 b/src/OCE/stpmlf.F90
index 8bced301..92716992 100644
--- a/src/OCE/stpmlf.F90
+++ b/src/OCE/stpmlf.F90
@@ -279,6 +279,12 @@ CONTAINS
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF ( ln_diurnal ) CALL diurnal_layers( kstp )
+ !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ ! GEOMETRIC
+ !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ IF ( l_ldfeke ) CALL ldf_eke( kstp, Nbb, Nnn, Naa ) ! GEOMETRIC param. (time evolution of eiv coefficient)
+
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! diagnostics and outputs
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
diff --git a/src/OFF/nemogcm.F90 b/src/OFF/nemogcm.F90
index 536f3348..e1046b80 100644
--- a/src/OFF/nemogcm.F90
+++ b/src/OFF/nemogcm.F90
@@ -34,6 +34,7 @@ MODULE nemogcm
! ! ocean physics
USE ldftra ! lateral diffusivity setting (ldf_tra_init routine)
USE ldfslp ! slopes of neutral surfaces (ldf_slp_init routine)
+ USE ldfeke ! GEOMETRIC parameterisation (ldf_eke_init routine)
USE traqsr ! solar radiation penetration (tra_qsr_init routine)
USE trabbl ! bottom boundary layer (tra_bbl_init routine)
USE traldf ! lateral physics (tra_ldf_init routine)
@@ -353,15 +354,16 @@ CONTAINS
CALL sbc_init( Nbb, Nnn, Naa ) ! Forcings : surface module
CALL bdy_init ! Open boundaries initialisation
- CALL zdf_phy_init( Nnn ) ! Vertical physics
+ CALL zdf_phy_init( Nnn ) ! Vertical physics
! ! Tracer physics
- CALL ldf_tra_init ! Lateral ocean tracer physics
- CALL ldf_eiv_init ! Eddy induced velocity param. must be done after ldf_tra_init
- CALL tra_ldf_init ! lateral mixing
- IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing
- IF( ln_traqsr ) CALL tra_qsr_init ! penetrative solar radiation
- IF( ln_trabbl ) CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme
+ CALL ldf_tra_init ! Lateral ocean tracer physics
+ CALL ldf_eiv_init ! Eddy induced velocity param. must be done after ldf_tra_init
+ IF( l_ldfeke ) CALL ldf_eke_init( Nbb, Nnn ) ! GEOMETRIC param.
+ CALL tra_ldf_init ! lateral mixing
+ IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing
+ IF( ln_traqsr ) CALL tra_qsr_init ! penetrative solar radiation
+ IF( ln_trabbl ) CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme
! ! Passive tracers
CALL trc_nam_run ! Needed to get restart parameters for passive tracers
--
GitLab