diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index cdacb620b0..7ad78049f3 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -32,6 +32,7 @@ setup: - git clone --recursive http://gitlab.gfdl.noaa.gov/ogrp/Gaea-stats-MOM6-examples.git tests && cd tests # Install / update testing scripts - git clone https://github.com/adcroft/MRS.git MRS + - (cd MRS ; git checkout xanadu-fms) # Update MOM6-examples and submodules - (cd MOM6-examples && git checkout . && git checkout dev/gfdl && git pull && git submodule init && git submodule update) - (cd MOM6-examples/src/MOM6 && git submodule update) @@ -73,7 +74,7 @@ gnu:ice-ocean-nolibs: - time tar zxf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz && cd tests - make -f MRS/Makefile.build build/gnu/env && cd build/gnu # mkdir -p build/gnu/repro/symmetric_dynamic/ocean_only && cd build/gnu/repro/symmetric_dynamic/ocean_only - - ../../MOM6-examples/src/mkmf/bin/list_paths -l ../../../config_src/{coupled_driver,dynamic} ../../../src ../../MOM6-examples/src/{FMS,coupler,SIS2,icebergs,ice_ocean_extras,land_null,atmos_null} + - ../../MOM6-examples/src/mkmf/bin/list_paths -l ../../../config_src/{coupled_driver,dynamic} ../../../src ../../MOM6-examples/src/{FMS,coupler,SIS2,icebergs,ice_param,land_null,atmos_null} - ../../MOM6-examples/src/mkmf/bin/mkmf -t ../../MOM6-examples/src/mkmf/templates/ncrc-gnu.mk -p MOM6 -c"-Duse_libMPI -Duse_netCDF -D_USE_LEGACY_LAND_ -Duse_AM3_physics" path_names - time (source ./env ; make NETCDF=3 REPRO=1 MOM6 -s -j) diff --git a/.testing/Makefile b/.testing/Makefile index d0f098e411..ee561375a3 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -18,7 +18,7 @@ else EXPERIMENTS ?= unit_tests double_gyre flow_downslope/z CVmix_SCM_tests/cooling_only/EPBL endif -FMS_PACKAGES ?= platform,include,memutils,constants,mpp,fms,time_manager,diag_manager,data_override,coupler/ensemble_manager.F90,axis_utils,horiz_interp,time_interp,astronomy,mosaic,random_numbers +FMS_PACKAGES ?= platform,include,memutils,constants,mpp,fms,time_manager,diag_manager,data_override,coupler/coupler_types.F90,coupler/ensemble_manager.F90,axis_utils,horiz_interp,time_interp,astronomy,mosaic,random_numbers TEMPLATE ?= .testing/linux-ubuntu-xenial-gnu.mk MPIRUN ?= mpirun diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 75ad9c427f..3364943222 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -999,13 +999,11 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if (CS%use_RK2) then call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, & - CS%Barotropic_CSp, CS%thickness_diffuse_CSp) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE) else call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, & - CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, & - CS%Barotropic_CSp, CS%thickness_diffuse_CSp, Waves=Waves) + CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, Waves=Waves) endif if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)") diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 31cb35b54d..2422ac7028 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -689,8 +689,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & ! diffu = horizontal viscosity terms (u_av) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & - MEKE, Varmix, CS%barotropic_CSp, thickness_diffuse_CSp, & - G, GV, US, CS%hor_visc_CSp, OBC=CS%OBC) + MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, & + OBC=CS%OBC, BT=CS%barotropic_CSp) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)") @@ -1147,14 +1147,11 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, & CS%tides_CSp) -! CS%barotropic_CSp => Barotropic_CSp -! CS%thickness_diffuse_CSp => thickness_diffuse_CSp - if (.not. query_initialized(CS%diffu,"diffu",restart_CS) .or. & .not. query_initialized(CS%diffv,"diffv",restart_CS)) & call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, & - CS%barotropic_CSp, thickness_diffuse_CSp, & - G, GV, US, CS%hor_visc_CSp, OBC=CS%OBC) + G, GV, US, CS%hor_visc_CSp, & + OBC=CS%OBC, BT=CS%barotropic_CSp) if (.not. query_initialized(CS%u_av,"u2", restart_CS) .or. & .not. query_initialized(CS%u_av,"v2", restart_CS)) then CS%u_av(:,:,:) = u(:,:,:) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index afdc99f2d7..dd03e11f42 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -183,7 +183,7 @@ module MOM_dynamics_unsplit !! 3rd order (for the inviscid momentum equations) order scheme subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE, Barotropic, thickness_diffuse, Waves) + VarMix, MEKE, Waves) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -218,10 +218,6 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & !! that specify the spatially variable viscosities. type(MEKE_type), pointer :: MEKE !< A pointer to a structure containing !! fields related to the Mesoscale Eddy Kinetic Energy. - type(barotropic_CS), pointer :: Barotropic!< Pointer to a structure containing - !! barotropic velocities - type(thickness_diffuse_CS), pointer :: thickness_diffuse!< Pointer to a structure containing - !! interface height diffusivities type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing !! fields related to the surface wave conditions @@ -260,7 +256,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & call enable_averaging(dt,Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, & - Barotropic, thickness_diffuse, G, GV, US, CS%hor_visc_CSp) + G, GV, US, CS%hor_visc_CSp) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index f88c08ac45..2cb22b12fe 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -186,7 +186,7 @@ module MOM_dynamics_unsplit_RK2 !> Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, & p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, & - VarMix, MEKE, Barotropic, thickness_diffuse) + VarMix, MEKE) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -232,10 +232,6 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, type(MEKE_type), pointer :: MEKE !< A pointer to a structure containing !! fields related to the Mesoscale !! Eddy Kinetic Energy. - type(barotropic_CS), pointer :: Barotropic!< Pointer to a structure containing - !! barotropic velocities - type(thickness_diffuse_CS), pointer :: thickness_diffuse!< Pointer to a structure containing - !! interface height diffusivities ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up @@ -271,7 +267,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call enable_averaging(dt,Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_in, v_in, h_in, CS%diffu, CS%diffv, MEKE, VarMix, & - Barotropic, thickness_diffuse, G, GV, US, CS%hor_visc_CSp) + G, GV, US, CS%hor_visc_CSp) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) call pass_vector(CS%diffu, CS%diffv, G%Domain, clock=id_clock_pass) diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index 7c1ee90f12..23d6f19e4e 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -69,6 +69,9 @@ subroutine find_obsolete_params(param_file) hint="Instead use OBC_SEGMENT_XXX_DATA.") call obsolete_char(param_file, "EXTEND_OBC_SEGMENTS", & hint="This option is no longer needed, nor supported.") + call obsolete_char(param_file, "MEKE_VISCOSITY_COEFF", & + hint="This option has been replaced by MEKE_VISCOSITY_COEFF_KU and \n" //& + " MEKE_VISCOSITY_COEFF_AU. Please set these parameters instead.") nseg = 0 call read_param(param_file, "OBC_NUMBER_OF_SEGMENTS", nseg) do l_seg = 1,nseg diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index da0310c012..21ad5a9800 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -1420,9 +1420,13 @@ end subroutine MEKE_end !! \subsection section_MEKE_viscosity Viscosity derived from MEKE !! !! As for \f$ \kappa_M \f$, the predicted eddy velocity scale can be -!! used to form an eddy viscosity: +!! used to form a harmonic eddy viscosity, !! -!! \f[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } . \f] +!! \f[ \kappa_u = \gamma_u \sqrt{ U_e^2 A_\Delta } \f] +!! +!! as well as a biharmonic eddy viscosity, +!! +!! \f[ \kappa_4 = \gamma_4 \sqrt{ U_e^2 A_\Delta^3 } \f] !! !! \subsection section_MEKE_limit_case Limit cases for local source-dissipative balance !! @@ -1459,7 +1463,8 @@ end subroutine MEKE_end !! | \f$ \kappa_4 \f$ | MEKE_K4 | !! | \f$ \gamma_\kappa \f$ | MEKE_KHCOEFF | !! | \f$ \gamma_M \f$ | MEKE_KHMEKE_FAC | -!! | \f$ \gamma_u \f$ | MEKE_VISCOSITY_COEFF | +!! | \f$ \gamma_u \f$ | MEKE_VISCOSITY_COEFF_KU | +!! | \f$ \gamma_4 \f$ | MEKE_VISCOSITY_COEFF_AU | !! | \f$ \gamma_{min}^2 \f$| MEKE_MIN_GAMMA2 | !! | \f$ \alpha_d \f$ | MEKE_ALPHA_DEFORM | !! | \f$ \alpha_f \f$ | MEKE_ALPHA_FRICT | diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 7d0b5a8dca..977d9b9228 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -190,8 +190,8 @@ module MOM_hor_visc !! u[is-2:ie+2,js-2:je+2] !! v[is-2:ie+2,js-2:je+2] !! h[is-1:ie+1,js-1:je+1] -subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, & - thickness_diffuse, G, GV, US, CS, OBC) +subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & + CS, OBC, BT) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & @@ -199,7 +199,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & intent(in) :: v !< The meridional velocity [m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(inout) :: h !< Layer thicknesses, in H + intent(inout) :: h !< Layer thicknesses, in H !! (usually m or kg m-2). real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & intent(out) :: diffu !< Zonal acceleration due to convergence of @@ -211,14 +211,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, !! related to Mesoscale Eddy Kinetic Energy. type(VarMix_CS), pointer :: VarMix !< Pointer to a structure with fields that !! specify the spatially variable viscosities - type(barotropic_CS), pointer :: Barotropic !< Pointer to a structure containing - !! barotropic velocities - type(thickness_diffuse_CS), pointer :: thickness_diffuse !< Pointer to a structure containing - !! interface height diffusivities type(hor_visc_CS), pointer :: CS !< Control structure returned by a previous type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type !! call to hor_visc_init. type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type + type(barotropic_CS), optional, pointer :: BT !< Pointer to a structure containing + !! barotropic velocities. ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & @@ -386,20 +384,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, use_MEKE_Ku = associated(MEKE%Ku) use_MEKE_Au = associated(MEKE%Au) -!$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h, & -!$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, & -!$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, & -!$OMP find_FrictWork,FrictWork,use_MEKE_Ku, & -!$OMP use_MEKE_Au, MEKE, hq, & -!$OMP mod_Leith, legacy_bound, div_xx_h, vort_xy_q) & -!$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, & -!$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, & -!$OMP sh_xx_bt, sh_xy_bt, dvdx_bt, dudy_bt, & -!$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, & -!$OMP vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, & -!$OMP div_xx, div_xx_dx, div_xx_dy,local_strain, & -!$OMP Shear_mag, h2uq, h2vq, Kh_scale, hrat_min) - do j=js,je ; do i=is,ie boundary_mask(i,j) = (G%mask2dCu(I,j) * G%mask2dCv(i,J) * G%mask2dCu(I-1,j) * G%mask2dCv(i,J-1)) enddo ; enddo @@ -419,7 +403,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, ! call pass_var(boundary_mask, G%Domain, complete=.true.) ! Get barotropic velocities and their gradients - call barotropic_get_tav(Barotropic, ubtav, vbtav, G) + call barotropic_get_tav(BT, ubtav, vbtav, G) call pass_vector(ubtav, vbtav, G%Domain) do j=js,je ; do i=is,ie @@ -480,6 +464,19 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, Barotropic, endif ! use_GME + !$OMP parallel do default(none) shared(Isq,Ieq,Jsq,Jeq,nz,CS,G,GV,u,v,is,js,ie,je,h, & + !$OMP rescale_Kh,VarMix,h_neglect,h_neglect3, & + !$OMP Kh_h,Ah_h,Kh_q,Ah_q,diffu,apply_OBC,OBC,diffv, & + !$OMP find_FrictWork,FrictWork,use_MEKE_Ku, & + !$OMP use_MEKE_Au, MEKE, hq, & + !$OMP mod_Leith, legacy_bound, div_xx_h, vort_xy_q) & + !$OMP private(u0, v0, sh_xx, str_xx, visc_bound_rem, & + !$OMP sh_xy, str_xy, Ah, Kh, AhSm, KhSm, dvdx, dudy, & + !$OMP sh_xx_bt, sh_xy_bt, dvdx_bt, dudy_bt, & + !$OMP bhstr_xx, bhstr_xy,FatH,RoScl, hu, hv, h_u, h_v, & + !$OMP vort_xy,vort_xy_dx,vort_xy_dy,Vort_mag,AhLth,KhLth, & + !$OMP div_xx, div_xx_dx, div_xx_dy,local_strain, & + !$OMP Shear_mag, h2uq, h2vq, Kh_scale, hrat_min) do k=1,nz ! The following are the forms of the horizontal tension and horizontal @@ -1410,6 +1407,8 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS) logical :: get_all ! If true, read and log all parameters, regardless of ! whether they are used, to enable spell-checking of ! valid parameters. + logical :: split ! If true, use the split time stepping scheme. + ! If false and USE_GME = True, issue a FATAL error. character(len=64) :: inputdir, filename real :: deg2rad ! Converts degrees to radians real :: slat_fn ! sin(lat)**Kh_pwr_of_sine @@ -1563,7 +1562,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS) "the grid spacing to calculate the biharmonic viscosity. \n"//& "The final viscosity is the largest of this scaled \n"//& "viscosity, the Smagorinsky and Leith viscosities, and AH.", & - units="m s-1", default=0.1) + units="m s-1", default=0.0) call get_param(param_file, mdl, "SMAGORINSKY_AH", CS%Smagorinsky_Ah, & "If true, use a biharmonic Smagorinsky nonlinear eddy \n"//& "viscosity.", default=.false.) @@ -1642,6 +1641,14 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS) "If true, use the GM+E backscatter scheme in association \n"//& "with the Gent and McWilliams parameterization.", default=.false.) + if (CS%use_GME) then + call get_param(param_file, mdl, "SPLIT", split, & + "Use the split time stepping if true.", default=.true., & + do_not_log=.true.) + if (.not. split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// & + "cannot be used with SPLIT=False.") + endif + if (CS%bound_Kh .or. CS%bound_Ah .or. CS%better_bound_Kh .or. CS%better_bound_Ah) & call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step.", units = "s", & diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 7a85292e54..68b747182d 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -63,8 +63,10 @@ module MOM_thickness_diffuse !! with GME closure. logical :: MEKE_GEOMETRIC !< If true, uses the GM coefficient formulation from the GEOMETRIC !! framework (Marshall et al., 2012) - real :: MEKE_GEOMETRIC_alpha !< The nondimensional coefficient governing the efficiency of + real :: MEKE_GEOMETRIC_alpha!< The nondimensional coefficient governing the efficiency of !! the GEOMETRIC thickness difussion [nondim] + real :: MEKE_GEOMETRIC_epsilon !< Minimum Eady growth rate for the GEOMETRIC thickness + !! diffusivity [s-1]. logical :: Use_KH_in_MEKE !< If true, uses the thickness diffusivity calculated here to diffuse MEKE. logical :: GM_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather !! than the streamfunction for the GM source term. @@ -144,7 +146,6 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp real :: hv(SZI_(G), SZJ_(G)) ! v-thickness [H ~> m or kg m-2] real :: KH_u_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [m2 s-1] real :: KH_v_lay(SZI_(G), SZJ_(G)) ! layer ave thickness diffusivities [m2 s-1] - real :: epsilon if (.not. associated(CS)) call MOM_error(FATAL, "MOM_thickness_diffuse:"// & "Module must be initialized before it is used.") @@ -176,7 +177,6 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp cg1 => null() endif - epsilon = 1.e-6 !$OMP parallel do default(none) shared(is,ie,js,je,KH_u_CFL,dt,G,CS) do j=js,je ; do I=is-1,ie @@ -216,7 +216,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (CS%MEKE_GEOMETRIC) then do j=js,je ; do I=is-1,ie Khth_Loc_u(I,j) = Khth_Loc_u(I,j) + G%mask2dCu(I,j) * CS%MEKE_GEOMETRIC_alpha * 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i+1,j)) / & - (VarMix%SN_u(I,j) + epsilon) + (VarMix%SN_u(I,j) + CS%MEKE_GEOMETRIC_epsilon) enddo ; enddo else do j=js,je ; do I=is-1,ie @@ -294,7 +294,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (CS%MEKE_GEOMETRIC) then do j=js-1,je ; do I=is,ie Khth_Loc(I,j) = Khth_Loc(I,j) + G%mask2dCv(i,J) * CS%MEKE_GEOMETRIC_alpha * 0.5*(MEKE%MEKE(i,j)+MEKE%MEKE(i,j+1)) / & - (VarMix%SN_v(i,J) + epsilon) + (VarMix%SN_v(i,J) + CS%MEKE_GEOMETRIC_epsilon) enddo ; enddo else do J=js-1,je ; do i=is,ie @@ -362,7 +362,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp if (CS%MEKE_GEOMETRIC) then do j=js,je ; do I=is,ie MEKE%Kh(i,j) = CS%MEKE_GEOMETRIC_alpha * MEKE%MEKE(i,j) / & - (0.25*(VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)) + epsilon) + (0.25*(VarMix%SN_u(I,j)+VarMix%SN_u(I-1,j)+VarMix%SN_v(i,J)+VarMix%SN_v(i,J-1)) + & + CS%MEKE_GEOMETRIC_epsilon) enddo ; enddo endif endif ; endif @@ -1923,9 +1924,17 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & "If true, uses the GM coefficient formulation \n"//& "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.) - call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, & + if (CS%MEKE_GEOMETRIC) then + + call get_param(param_file, mdl, "MEKE_GEOMETRIC_EPSILON", CS%MEKE_GEOMETRIC_epsilon, & + "Minimum Eady growth rate used in the calculation of \n"//& + "GEOMETRIC thickness diffusivity.", units="s-1", default=1.0e-7) + + call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", CS%MEKE_GEOMETRIC_alpha, & "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//& "thickness diffusion.", units="nondim", default=0.05) + endif + call get_param(param_file, mdl, "USE_KH_IN_MEKE", CS%Use_KH_in_MEKE, & "If true, uses the thickness diffusivity calculated here to diffuse \n"//& "MEKE.", default=.false.)