diff --git a/.gitmodules b/.gitmodules index e69de29bb..f9a8115c8 100644 --- a/.gitmodules +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "src/GOTM5.2/code"] + path = src/GOTM5.2/code + url = https://github.com/gotm-model/code.git diff --git a/cmake/SCHISM.local.build b/cmake/SCHISM.local.build index 858576d04..72c6d66eb 100644 --- a/cmake/SCHISM.local.build +++ b/cmake/SCHISM.local.build @@ -27,7 +27,8 @@ set (PREC_EVAP OFF CACHE BOOLEAN "Include precipitation and evaporation calculat set (USE_BULK_FAIRALL OFF CACHE BOOLEAN "Enable Fairall bulk scheme for air-sea exchange") #IMPOSE_NET_FLUX will be removed eventually set (IMPOSE_NET_FLUX OFF CACHE BOOLEAN "Specify net heat and salt fluxes in sflux") -##Older versions of GOTM (3.*) have issues with netcdf v4, so are not maintained +## GOTM5.2 is installed by default, and has been tested successfully, defined GOTM_BASE to use a non-standard GOTM source +##set( GOTM_BASE /work2/03473/seatonc/frontera/GOTM5.2/code CACHE STRING "Path to non-standard GOTM base directory") set (USE_GOTM OFF CACHE BOOLEAN "Use GOTM turbulence model. This just enables the build -- GOTM must still be selected in param.nml") set (USE_HA OFF CACHE BOOLEAN "Enable harmonic analysis output modules") set( USE_MARSH OFF CACHE BOOLEAN "Use marsh module") diff --git a/docs/input-output/param.md b/docs/input-output/param.md index be4c2b6c4..bd7107b8d 100644 --- a/docs/input-output/param.md +++ b/docs/input-output/param.md @@ -258,15 +258,16 @@ If `itur=2`, the zero-equation Pacanowski and Philander closure is used. In this If `itur=3`, then the two-equation closure schemes from the GLS model of Umlauf and Burchard (2003) are used. In this case, 2 additional parameters are required: `mid`, `stab`, which specify the closure scheme and stability function used: `mid=` `MY` is Mellor & Yamada; `KL` is GLS as k-kl; `KE` is GLS as $k-\varepsilon$; `KW` is GLS as $k-\omega$; `UB` is Umlauf & Burchard's optimal. `stab=GA` is Galperin's clipping (only for MY); `KC` is Kantha & Clayson's stability function). Also the user needs to specify max/min diffusivity/viscosity in `diffmax.gr3` and `diffmin.gr3`, as well as a surface mixing length scale constant `xlsc0`. -If `itur=4`, GOTM turbulence model is invoked; the user needs to compile the GOTM libraries first -(see README inside GOTM/ for instructions), and turn on pre-processing flag `USE_GOTM` in makefile -and recompile. In this case, the minimum and maximum viscosity/diffusivity are still specified in `diffmin.gr3` -and `diffmax.gr3`. In addition, GOTM also requires an input called `gotmturb.inp`. There are some -ready-made samples for this input in the source code bundle. If you wish to tune some parameters +If `itur=4`, GOTM turbulence model is invoked; the user needs to turn on pre-processing flag `USE_GOTM` in makefile +and recompile (GOTM5.2 uses cmake, so does not need to be pre-compiled). In this case, the minimum and maximum +viscosity/diffusivity are still specified in `diffmin.gr3` +and `diffmax.gr3`. In addition, GOTM also requires an input file called `gotmturb.inp` (the file is in NML format. If it is absent in a GOTM run, GOTM will error out with a notice that gotmturb.nml is missing, but SCHISM expects the file to be named .inp). There are some +ready-made samples for this input in the source code bundle. An example for the Columbia River +(with steady state Richardson number set to ri_st= 0.15) is included in 'sample_inputs/'. If you wish to tune some parameters inside, you may consult [gotm.net](https://gotm.net) for more details. !!! note - 1. We only tested GOTM v3, not newer versions of GOTM. Using `itur=3` generally gave similar results. + 1. GOTM has only been tested up to v5.2, not newer versions of GOTM. Using `itur=3` generally gave similar results, except in the Columbia River under high flow conditions, where GOTM produces substantially better stratification. ### meth_sink=1 (int) Option for sinks. If `meth_sink =1`, the sink value is reset to `0` if an element is dry with diff --git a/sample_inputs/gotmturb.inp b/sample_inputs/gotmturb.inp new file mode 100644 index 000000000..4e6f6d155 --- /dev/null +++ b/sample_inputs/gotmturb.inp @@ -0,0 +1,304 @@ +!$Id: gotmturb.proto,v 1.1.1.1 2003/03/11 13:38:58 kbk Exp $ +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! What type of equations are solved in the turbulence model? +! +! turb_method -> type of turbulence closure +! +! 0: convective adjustment +! 1: analytical eddy visc. and diff. profiles, not coded yet +! 2: turbulence Model calculating TKE and length scale +! (specify stability function below) +! 3: second-order model (see "scnd" namelist below) +! 99: KPP model (requires "kpp.inp" with specifications) +! +! +! tke_method -> type of equation for TKE +! +! 1: algebraic equation +! 2: dynamic equation (k-epsilon style) +! 3: dynamic equation (Mellor-Yamada style) +! +! +! len_scale_method -> type of model for dissipative length scale +! +! 1: parabolic shape +! 2: triangle shape +! 3: Xing and Davies [1995] +! 4: Robert and Ouellet [1987] +! 5: Blackadar (two boundaries) [1962] +! 6: Bougeault and Andre [1986] +! 7: Eifler and Schrimpf (ISPRAMIX) [1992] +! 8: dynamic dissipation rate equation +! 9: dynamic Mellor-Yamada q^2l-equation +! 10: generic length scale (GLS) +! +! +! stab_method -> type of stability function +! +! 1: constant stability functions +! 2: Munk and Anderson [1954] +! 3: Schumann and Gerz [1995] +! 4: Eifler and Schrimpf [1992] +! +!------------------------------------------------------------------------------- + &turbulence + turb_method= 3 + tke_method= 2 + len_scale_method=8 + stab_method= 1 + / + +!------------------------------------------------------------------------------- +! What boundary conditions are used? +! +! k_ubc, k_lbc -> upper and lower boundary conditions +! for k-equation +! 0: prescribed BC +! 1: flux BC +! +! psi_ubc, psi_lbc -> upper and lower boundary conditions +! for the length-scale equation (e.g. +! epsilon, kl, omega, GLS) +! 0: prescribed BC +! 1: flux BC +! +! +! ubc_type -> type of upper boundary layer +! 0: viscous sublayer (not yet impl.) +! 1: logarithmic law of the wall +! 2: tke-injection (breaking waves) +! +! lbc_type -> type of lower boundary layer +! 0: viscous sublayer (not yet impl.) +! 1: logarithmic law of the wall +! +!------------------------------------------------------------------------------- + &bc + k_ubc= 1 + k_lbc= 1 + psi_ubc= 1 + psi_lbc= 1 + ubc_type= 1 + lbc_type= 1 + / + +!------------------------------------------------------------------------------- +! What turbulence parameters have been described? +! +! cm0_fix -> value of cm0 for turb_method=2 +! Prandtl0_fix -> value of the turbulent Prandtl-number for stab_method=1-4 +! cw -> constant of the wave-breaking model +! (Craig & Banner (1994) use cw=100) +! compute_kappa -> compute von Karman constant from model parameters +! kappa -> the desired von Karman constant (if compute_kappa=.true.) +! compute_c3 -> compute c3 (E3 for Mellor-Yamada) for given Ri_st +! Ri_st -> the desired steady-state Richardson number (if compute_c3=.true.) +! length_lim -> apply length scale limitation (see Galperin et al. 1988) +! galp -> coef. for length scale limitation +! const_num -> minimum eddy diffusivity (only with turb_method=0) +! const_nuh -> minimum heat diffusivity (only with turb_method=0) +! k_min -> minimun TKE +! eps_min -> minimum dissipation rate +! kb_min -> minimun buoyancy variance +! epsb_min -> minimum buoyancy variance destruction rate +! +!------------------------------------------------------------------------------- + &turb_param + cm0_fix= 0.5477 + Prandtl0_fix= 0.74 + cw= 100. + compute_kappa= .true. + kappa= 0.4 + compute_c3= .true. + ri_st= 0.15 + length_lim= .false. + galp= 0.53 + const_num= 5.e-4 + const_nuh= 5.e-4 + k_min= 1.e-10 + eps_min= 1.e-14 + kb_min= 1.e-10 + epsb_min= 1.e-14 + / + +!------------------------------------------------------------------------------- +! The generic model (Umlauf & Burchard, J. Mar. Res., 2003) +! +! This part is active only, when len_scale_method=10 has been set. +! +! compute_param -> compute the model parameters: +! if this is .false., you have to set all +! model parameters (m,n,cpsi1,...) explicitly +! if this is .true., all model parameters +! set by you (except m) will be ignored and +! re-computed from kappa, d, alpha, etc. +! (see Umlauf&Burchard 2002) +! +! m: -> exponent for k +! n: -> exponent for l +! p: -> exponent for cm0 +! +! Examples: +! +! k-epsilon (Rodi 1987) : m=3/2, n=-1, p=3 +! k-omega (Umlauf et al. 2003) : m=1/2, n=-1, p=-1 +! +! cpsi1 -> emp. coef. in psi equation +! cpsi2 -> emp. coef. in psi equation +! cpsi3minus -> cpsi3 for stable stratification +! cpsi3plus -> cpsi3 for unstable stratification +! sig_kpsi -> Schmidt number for TKE diffusivity +! sig_psi -> Schmidt number for psi diffusivity +! +!------------------------------------------------------------------------------- + &generic + compute_param= .false. + gen_m= 1.5 + gen_n= -1.0 + gen_p= 3.0 + cpsi1= 1.44 + cpsi2= 1.92 + cpsi3minus= 0.0 + cpsi3plus = 0.0 + sig_kpsi= 1.0 + sig_psi= 1.3 + gen_d= -1.087 + gen_alpha= -4.97 + gen_l= 0.09 + / + +!------------------------------------------------------------------------------- +! The k-epsilon model (Rodi 1987) +! +! This part is active only, when len_scale_method=8 has been set. +! +! ce1 -> emp. coef. in diss. eq. +! ce2 -> emp. coef. in diss. eq. +! ce3minus -> ce3 for stable stratification, overwritten if compute_c3=.true. +! ce3plus -> ce3 for unstable stratification (Rodi 1987: ce3plus=1.0) +! sig_k -> Schmidt number for TKE diffusivity +! sig_e -> Schmidt number for diss. diffusivity +! sig_peps -> if .true. -> the wave breaking parameterisation suggested +! by Burchard (JPO 31, 2001, 3133-3145) will be used. +!------------------------------------------------------------------------------- + &keps + ce1= 1.44 + ce2= 1.92 + ce3minus= 0.0 + ce3plus= 1.0 + sig_k= 1.0 + sig_e= 1.3 + sig_peps= .false. + / + +!------------------------------------------------------------------------------- +! The Mellor-Yamada model (Mellor & Yamada 1982) +! +! This part is active only, when len_scale_method=9 has been set! +! +! e1 -> coef. in MY q**2 l equation +! e2 -> coef. in MY q**2 l equation +! e3 -> coef. in MY q**2 l equation, overwritten if compute_c3=.true. +! sq -> turbulent diffusivities of q**2 (= 2k) +! sl -> turbulent diffusivities of q**2 l +! my_length -> prescribed barotropic lengthscale in q**2 l equation of MY +! 1: parabolic +! 2: triangular +! 3: lin. from surface +! new_constr -> stabilisation of Mellor-Yamada stability functions +! according to Burchard & Deleersnijder (2001) +! (if .true.) +! +!------------------------------------------------------------------------------- + &my + e1= 1.8 + e2= 1.33 + e3= 1.8 + sq= 0.2 + sl= 0.2 + my_length= 1 + new_constr= .false. + / + +!------------------------------------------------------------------------------- +! The second-order model +! +! scnd_method -> type of second-order model +! 1: EASM with quasi-equilibrium +! 2: EASM with weak equilibrium, buoy.-variance algebraic +! 3: EASM with weak equilibrium, buoy.-variance from PDE +! +! kb_method -> type of equation for buoyancy variance +! +! 1: algebraic equation for buoyancy variance +! 2: PDE for buoyancy variance +! +! +! epsb_method -> type of equation for variance destruction +! +! 1: algebraic equation for variance destruction +! 2: PDE for variance destruction +! +! +! scnd_coeff -> coefficients of second-order model +! +! 0: read the coefficients from this file +! 1: coefficients of Gibson and Launder (1978) +! 2: coefficients of Mellor and Yamada (1982) +! 3: coefficients of Kantha and Clayson (1994) +! 4: coefficients of Luyten et al. (1996) +! 5: coefficients of Canuto et al. (2001) (version A) +! 6: coefficients of Canuto et al. (2001) (version B) +! 7: coefficients of Cheng et al. (2002) +! +!------------------------------------------------------------------------------- + &scnd + scnd_method= 2 + kb_method= 1 + epsb_method= 1 + scnd_coeff= 5 + cc1= 5.0 + cc2= 0.8000 + cc3= 1.9680 + cc4= 1.1360 + cc5= 0.0000 + cc6= 0.4000 + ct1= 5.9500 + ct2= 0.6000 + ct3= 1.0000 + ct4= 0.0000 + ct5= 0.3333 + ctt= 0.7200 + / + +!------------------------------------------------------------------------------- +! The internal wave model +! +! iw_model -> method to compute internal wave mixing +! 0: no internal waves mixing parameterisation +! 1: Mellor 1989 internal wave mixing +! 2: Large et al. 1994 internal wave mixing +! +! alpha -> coeff. for Mellor IWmodel (0: no IW, 0.7 Mellor 1989) +! +! The following six empirical parameters are used for the +! Large et al. 1994 shear instability and internal wave breaking +! parameterisations (iw_model = 2, all viscosities are in m**2/s): +! +! klimiw -> critcal value of TKE +! rich_cr -> critical Richardson number for shear instability +! numshear -> background diffusivity for shear instability +! numiw -> background viscosity for internal wave breaking +! nuhiw -> background diffusivity for internal wave breaking +!------------------------------------------------------------------------------- + &iw + iw_model= 0 + alpha= 0.0 + klimiw= 1e-6 + rich_cr= 0.7 + numshear= 5.e-3 + numiw= 1.e-4 + nuhiw= 1.e-5 + / diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index badedcffa..dee654465 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -328,8 +328,9 @@ string(REPLACE _HYDRO_SCHISM "" mod_tag_rev ${mod_tag_rev}) # call cmake with -DUSE_GOTM=ON -DGOTM_BASE=[gotm code directory] # if(USE_GOTM) -find_path(GOTM_BASE src/gotm/gotm.F90 DOC "GOTM source directory") +find_path(GOTM_BASE src/gotm/gotm.F90 $GOTM_BASE ${CMAKE_CURRENT_SOURCE_DIR}/GOTM5.2/code DOC "GOTM source directory") if(GOTM_BASE) + message(STATUS "using GOTM_BASE ${GOTM_BASE}") set(GOTM_BUILD_LIBRARIES_ONLY ON) set(GOTM_USE_FABM OFF CACHE BOOL "use GOTM without FABM") set(GOTM_USE_FLEXIBLE_OUTPUT OFF CACHE BOOL "use GOTM without flexible output") diff --git a/src/GOTM5.2/README.md b/src/GOTM5.2/README.md new file mode 100644 index 000000000..3e88627cc --- /dev/null +++ b/src/GOTM5.2/README.md @@ -0,0 +1,4 @@ +### Source information for GOTM5.2 + +GOTM5.2 directory copied from https://github.com/gotm-model/code/tree/v5.2 +License information for GOTM is included in the code/COPYING file. diff --git a/src/GOTM5.2/code b/src/GOTM5.2/code new file mode 160000 index 000000000..24749f915 --- /dev/null +++ b/src/GOTM5.2/code @@ -0,0 +1 @@ +Subproject commit 24749f91505b43001e532d3f1ac24ebc2b576d7e