diff --git a/.github/actions/testing-setup/action.yml b/.github/actions/testing-setup/action.yml index c47270af0d..1ab96aa3df 100644 --- a/.github/actions/testing-setup/action.yml +++ b/.github/actions/testing-setup/action.yml @@ -33,6 +33,7 @@ runs: echo "::group::Install linux packages" sudo apt-get update sudo apt-get install netcdf-bin libnetcdf-dev libnetcdff-dev mpich libmpich-dev + sudo apt-get install linux-tools-common echo "::endgroup::" - name: Compile FMS library diff --git a/.github/workflows/perfmon.yml b/.github/workflows/perfmon.yml new file mode 100644 index 0000000000..00e645c4fd --- /dev/null +++ b/.github/workflows/perfmon.yml @@ -0,0 +1,36 @@ +name: Performance Monitor + +on: [pull_request] + +jobs: + build-test-perfmon: + + runs-on: ubuntu-latest + defaults: + run: + working-directory: .testing + + steps: + - uses: actions/checkout@v2 + with: + submodules: recursive + + - uses: ./.github/actions/testing-setup + + - name: Compile optimized models + run: >- + make -j build.prof + MOM_TARGET_SLUG=$GITHUB_REPOSITORY + MOM_TARGET_LOCAL_BRANCH=$GITHUB_BASE_REF + DO_REGRESSION_TESTS=true + + - name: Generate profile data + run: >- + pip install f90nml && + make profile + DO_REGRESSION_TESTS=true + + - name: Generate perf data + run: | + sudo sysctl -w kernel.perf_event_paranoid=2 + make perf DO_REGRESSION_TESTS=true diff --git a/.testing/Makefile b/.testing/Makefile index 06b29dc690..1946b133d0 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -34,6 +34,7 @@ # Build configuration: # FCFLAGS_DEBUG Testing ("debug") compiler flags # FCFLAGS_REPRO Production ("repro") compiler flags +# FCFLAGS_OPT Aggressive optimization compiler flags # FCFLAGS_INIT Variable initialization flags # FCFLAGS_COVERAGE Code coverage flags # @@ -72,6 +73,7 @@ export MPIFC # NOTE: FMS will be built using FCFLAGS_DEBUG FCFLAGS_DEBUG ?= -g -O0 FCFLAGS_REPRO ?= -g -O2 +FCFLAGS_OPT ?= -g -O3 -mavx -fno-omit-frame-pointer FCFLAGS_INIT ?= FCFLAGS_COVERAGE ?= # Additional notes: @@ -96,6 +98,7 @@ DO_REPRO_TESTS ?= # Time measurement (configurable by the CI) TIME ?= time + #--- # Dependencies DEPS = deps @@ -114,6 +117,9 @@ CONFIGS := $(wildcard tc*) TESTS = grids layouts restarts nans dims openmps rotations DIMS = t l h z q r +# Set the framework +FRAMEWORK ?= fms1 + # REPRO tests enable reproducibility with optimization, and often do not match # the DEBUG results in older GCCs and vendor compilers, so we can optionally # disable them. @@ -122,6 +128,11 @@ ifeq ($(DO_REPRO_TESTS), true) TESTS += repros endif +# Profiling +ifeq ($(DO_PROFILE), false) + BUILDS += opt opt_target +endif + # The following variables are configured by Travis: # DO_REGRESSION_TESTS: true if $(TRAVIS_PULL_REQUEST) is a PR number # MOM_TARGET_SLUG: TRAVIS_REPO_SLUG @@ -195,6 +206,7 @@ endif .PHONY: all build.regressions all: $(foreach b,$(BUILDS),build/$(b)/MOM6) $(VENV_PATH) build.regressions: $(foreach b,symmetric target,build/$(b)/MOM6) +build.prof: $(foreach b,opt opt_target,build/$(b)/MOM6) # Executable BUILD_TARGETS = MOM6 Makefile path_names @@ -217,6 +229,7 @@ PATH_FMS = PATH="${PATH}:../../$(DEPS)/bin" SYMMETRIC_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(COVERAGE) $(FCFLAGS_FMS)" ASYMMETRIC_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" REPRO_FCFLAGS := FCFLAGS="$(FCFLAGS_REPRO) $(FCFLAGS_FMS)" +OPT_FCFLAGS := FCFLAGS="$(FCFLAGS_OPT) $(FCFLAGS_FMS)" OPENMP_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" TARGET_FCFLAGS := FCFLAGS="$(FCFLAGS_DEBUG) $(FCFLAGS_INIT) $(FCFLAGS_FMS)" @@ -230,6 +243,8 @@ build/asymmetric/Makefile: MOM_ENV=$(PATH_FMS) $(ASYMMETRIC_FCFLAGS) $(MOM_LDFLA build/repro/Makefile: MOM_ENV=$(PATH_FMS) $(REPRO_FCFLAGS) $(MOM_LDFLAGS) build/openmp/Makefile: MOM_ENV=$(PATH_FMS) $(OPENMP_FCFLAGS) $(MOM_LDFLAGS) build/target/Makefile: MOM_ENV=$(PATH_FMS) $(TARGET_FCFLAGS) $(MOM_LDFLAGS) +build/opt/Makefile: MOM_ENV=$(PATH_FMS) $(OPT_FCFLAGS) $(MOM_LDFLAGS) +build/opt_target/Makefile: MOM_ENV=$(PATH_FMS) $(OPT_FCFLAGS) $(MOM_LDFLAGS) build/coupled/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) build/nuopc/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) build/mct/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) @@ -240,12 +255,15 @@ build/asymmetric/Makefile: MOM_ACFLAGS=--enable-asymmetric build/repro/Makefile: MOM_ACFLAGS= build/openmp/Makefile: MOM_ACFLAGS=--enable-openmp build/target/Makefile: MOM_ACFLAGS= +build/opt/Makefile: MOM_ACFLAGS= +build/opt_target/Makefile: MOM_ACFLAGS= build/coupled/Makefile: MOM_ACFLAGS=--with-driver=coupled_driver build/nuopc/Makefile: MOM_ACFLAGS=--with-driver=nuopc_driver build/mct/Makefile: MOM_ACFLAGS=--with-driver=mct_driver # Fetch regression target source code build/target/Makefile: | $(TARGET_CODEBASE) +build/opt_target/Makefile: | $(TARGET_CODEBASE) # Define source code dependencies @@ -267,7 +285,7 @@ build/%/MOM6: build/%/Makefile build/%/Makefile: ../ac/configure ../ac/Makefile.in $(DEPS)/lib/libFMS.a $(MKMF) $(LIST_PATHS) mkdir -p $(@D) cd $(@D) \ - && $(MOM_ENV) ../../../ac/configure $(MOM_ACFLAGS) \ + && $(MOM_ENV) ../../../ac/configure $(MOM_ACFLAGS) --with-framework=$(FRAMEWORK) \ || (cat config.log && false) @@ -276,7 +294,8 @@ build/%/Makefile: ../ac/configure ../ac/Makefile.in $(DEPS)/lib/libFMS.a $(MKMF) # Fetch the regression target codebase -build/target/Makefile: $(TARGET_CODEBASE)/ac/configure $(DEPS)/lib/libFMS.a $(MKMF) $(LIST_PATHS) +build/target/Makefile build/opt_target/Makefile: \ + $(TARGET_CODEBASE)/ac/configure $(DEPS)/lib/libFMS.a $(MKMF) $(LIST_PATHS) mkdir -p $(@D) cd $(@D) \ && $(MOM_ENV) ../../$(TARGET_CODEBASE)/ac/configure $(MOM_ACFLAGS) \ @@ -549,7 +568,11 @@ $(eval $(call STAT_RULE,dim.q,symmetric,,Q_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.r,symmetric,,R_RESCALE_POWER=11,,1)) -# Restart tests require significant preprocessing, and are handled separately. +# Generate the half-period input namelist as follows: +# 1. Fetch DAYMAX and TIMEUNIT from MOM_input +# 2. Convert DAYMAX from TIMEUNIT to seconds +# 3. Apply seconds to `ocean_solo_nml` inside input.nml. +# NOTE: Assumes that runtime set by DAYMAX, will fail if set by input.nml work/%/restart/ocean.stats: build/symmetric/MOM6 $(VENV_PATH) rm -rf $(@D) mkdir -p $(@D) @@ -562,14 +585,13 @@ work/%/restart/ocean.stats: build/symmetric/MOM6 $(VENV_PATH) cd work/$*/restart; \ fi mkdir -p $(@D)/RESTART - # Generate the half-period input namelist - # TODO: Assumes that runtime set by DAYMAX, will fail if set by input.nml + # Set the half-period cd $(@D) \ && daymax=$$(grep DAYMAX MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ && timeunit=$$(grep TIMEUNIT MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ && if [ -z "$${timeunit}" ]; then timeunit="8.64e4"; fi \ && printf -v timeunit_int "%.f" "$${timeunit}" \ - && halfperiod=$$(printf "%.f" $$(bc <<< "scale=10; 0.5 * $${daymax} * $${timeunit_int}")) \ + && halfperiod=$$(awk -v t=$${daymax} -v dt=$${timeunit} 'BEGIN {printf "%.f", 0.5*t*dt}') \ && printf "\n&ocean_solo_nml\n seconds = $${halfperiod}\n/\n" >> input.nml # Remove any previous archived output rm -f results/$*/std.restart{1,2}.{out,err} @@ -621,6 +643,64 @@ test.summary: fi +#--- +# Profiling +# XXX: This is experimental work to track, log, and report changes in runtime +PCONFIGS = p0 + +.PHONY: profile +profile: $(foreach p,$(PCONFIGS), prof.$(p)) + +.PHONY: prof.p0 +prof.p0: work/p0/opt/clocks.json work/p0/opt_target/clocks.json + python tools/compare_clocks.py $^ + +work/p0/%/clocks.json: work/p0/%/std.out + python tools/parse_fms_clocks.py -d $(@D) $^ > $@ + +work/p0/opt/std.out: build/opt/MOM6 +work/p0/opt_target/std.out: build/opt_target/MOM6 + +work/p0/%/std.out: + mkdir -p $(@D) + cp -RL p0/* $(@D) + mkdir -p $(@D)/RESTART + echo -e "" > $(@D)/MOM_override + cd $(@D) \ + && $(MPIRUN) -n 1 ../../../$< 2> std.err > std.out + +#--- +# Same but with perf + +# TODO: This expects the -e flag, can I handle it in the command? +PERF_EVENTS ?= + +.PHONY: perf +perf: $(foreach p,$(PCONFIGS), perf.$(p)) + +.PHONY: prof.p0 +perf.p0: work/p0/opt/profile.json work/p0/opt_target/profile.json + python tools/compare_perf.py $^ + +work/p0/%/profile.json: work/p0/%/perf.data + python tools/parse_perf.py -f $< > $@ + +work/p0/opt/perf.data: build/opt/MOM6 +work/p0/opt_target/perf.data: build/opt_target/MOM6 + +work/p0/%/perf.data: + mkdir -p $(@D) + cp -RL p0/* $(@D) + mkdir -p $(@D)/RESTART + echo -e "" > $(@D)/MOM_override + cd $(@D) \ + && perf record \ + -F 3999 \ + ${PERF_EVENTS} \ + ../../../$< 2> std.perf.err > std.perf.out \ + || cat std.perf.err + + #---- # NOTE: These tests assert that we are in the .testing directory. diff --git a/.testing/p0/MOM_input b/.testing/p0/MOM_input new file mode 100644 index 0000000000..8f751d7bf1 --- /dev/null +++ b/.testing/p0/MOM_input @@ -0,0 +1,505 @@ +! This input file provides the adjustable run-time parameters for version 6 of the Modular Ocean Model (MOM6). +! Where appropriate, parameters use usually given in MKS units. + +! This particular file is for the example in benchmark. + +! This MOM_input file typically contains only the non-default values that are needed to reproduce this example. +! A full list of parameters for this example can be found in the corresponding MOM_parameter_doc.all file +! which is generated by the model at run-time. + +! === module MOM_domains === +NIGLOBAL = 32 ! + ! The total number of thickness grid points in the x-direction in the physical + ! domain. With STATIC_MEMORY_ this is set in MOM_memory.h at compile time. +NJGLOBAL = 32 ! + ! The total number of thickness grid points in the y-direction in the physical + ! domain. With STATIC_MEMORY_ this is set in MOM_memory.h at compile time. +LAYOUT = 1, 1 + ! The processor layout that was actually used. + +! === module MOM === +THICKNESSDIFFUSE = True ! [Boolean] default = False + ! If true, interface heights are diffused with a coefficient of KHTH. +THICKNESSDIFFUSE_FIRST = True ! [Boolean] default = False + ! If true, do thickness diffusion before dynamics. This is only used if + ! THICKNESSDIFFUSE is true. +DT = 900.0 ! [s] + ! The (baroclinic) dynamics time step. The time-step that is actually used will + ! be an integer fraction of the forcing time-step (DT_FORCING in ocean-only mode + ! or the coupling timestep in coupled mode.) +DT_THERM = 3600.0 ! [s] default = 900.0 + ! The thermodynamic and tracer advection time step. Ideally DT_THERM should be + ! an integer multiple of DT and less than the forcing or coupling time-step, + ! unless THERMO_SPANS_COUPLING is true, in which case DT_THERM can be an integer + ! multiple of the coupling timestep. By default DT_THERM is set to DT. +DTBT_RESET_PERIOD = 0.0 ! [s] default = 3600.0 + ! The period between recalculations of DTBT (if DTBT <= 0). If DTBT_RESET_PERIOD + ! is negative, DTBT is set based only on information available at + ! initialization. If 0, DTBT will be set every dynamics time step. The default + ! is set by DT_THERM. This is only used if SPLIT is true. +FRAZIL = True ! [Boolean] default = False + ! If true, water freezes if it gets too cold, and the accumulated heat deficit + ! is returned in the surface state. FRAZIL is only used if + ! ENABLE_THERMODYNAMICS is true. +C_P = 3925.0 ! [J kg-1 K-1] default = 3991.86795711963 + ! The heat capacity of sea water, approximated as a constant. This is only used + ! if ENABLE_THERMODYNAMICS is true. The default value is from the TEOS-10 + ! definition of conservative temperature. +SAVE_INITIAL_CONDS = True ! [Boolean] default = False + ! If true, write the initial conditions to a file given by IC_OUTPUT_FILE. +IC_OUTPUT_FILE = "GOLD_IC" ! default = "MOM_IC" + ! The file into which to write the initial conditions. + +! === module MOM_fixed_initialization === +INPUTDIR = "INPUT" ! default = "." + ! The directory in which input files are found. + +! === module MOM_grid_init === +GRID_CONFIG = "cartesian" ! + ! A character string that determines the method for defining the horizontal + ! grid. Current options are: + ! mosaic - read the grid from a mosaic (supergrid) + ! file set by GRID_FILE. + ! cartesian - use a (flat) Cartesian grid. + ! spherical - use a simple spherical grid. + ! mercator - use a Mercator spherical grid. +SOUTHLAT = -41.0 ! [degrees] + ! The southern latitude of the domain. +LENLAT = 41.0 ! [degrees] + ! The latitudinal length of the domain. +LENLON = 90.0 ! [degrees] + ! The longitudinal length of the domain. +ISOTROPIC = True ! [Boolean] default = False + ! If true, an isotropic grid on a sphere (also known as a Mercator grid) is + ! used. With an isotropic grid, the meridional extent of the domain (LENLAT), + ! the zonal extent (LENLON), and the number of grid points in each direction are + ! _not_ independent. In MOM the meridional extent is determined to fit the zonal + ! extent and the number of grid points, while grid is perfectly isotropic. +TOPO_CONFIG = "benchmark" ! + ! This specifies how bathymetry is specified: + ! file - read bathymetric information from the file + ! specified by (TOPO_FILE). + ! flat - flat bottom set to MAXIMUM_DEPTH. + ! bowl - an analytically specified bowl-shaped basin + ! ranging between MAXIMUM_DEPTH and MINIMUM_DEPTH. + ! spoon - a similar shape to 'bowl', but with an vertical + ! wall at the southern face. + ! halfpipe - a zonally uniform channel with a half-sine + ! profile in the meridional direction. + ! bbuilder - build topography from list of functions. + ! benchmark - use the benchmark test case topography. + ! Neverworld - use the Neverworld test case topography. + ! DOME - use a slope and channel configuration for the + ! DOME sill-overflow test case. + ! ISOMIP - use a slope and channel configuration for the + ! ISOMIP test case. + ! DOME2D - use a shelf and slope configuration for the + ! DOME2D gravity current/overflow test case. + ! Kelvin - flat but with rotated land mask. + ! seamount - Gaussian bump for spontaneous motion test case. + ! dumbbell - Sloshing channel with reservoirs on both ends. + ! shelfwave - exponential slope for shelfwave test case. + ! Phillips - ACC-like idealized topography used in the Phillips config. + ! dense - Denmark Strait-like dense water formation and overflow. + ! USER - call a user modified routine. + +! === module benchmark_initialize_topography === +MINIMUM_DEPTH = 1.0 ! [m] default = 0.0 + ! The minimum depth of the ocean. +MAXIMUM_DEPTH = 5500.0 ! [m] + ! The maximum depth of the ocean. + +! === module MOM_verticalGrid === +! Parameters providing information about the vertical grid. +NK = 75 ! [nondim] + ! The number of model layers. + +! === module MOM_EOS === + +! === module MOM_tracer_flow_control === +USE_IDEAL_AGE_TRACER = True ! [Boolean] default = False + ! If true, use the ideal_age_example tracer package. + +! === module ideal_age_example === + +! === module MOM_coord_initialization === +COORD_CONFIG = "ts_range" ! default = "none" + ! This specifies how layers are to be defined: + ! ALE or none - used to avoid defining layers in ALE mode + ! file - read coordinate information from the file + ! specified by (COORD_FILE). + ! BFB - Custom coords for buoyancy-forced basin case + ! based on SST_S, T_BOT and DRHO_DT. + ! linear - linear based on interfaces not layers + ! layer_ref - linear based on layer densities + ! ts_ref - use reference temperature and salinity + ! ts_range - use range of temperature and salinity + ! (T_REF and S_REF) to determine surface density + ! and GINT calculate internal densities. + ! gprime - use reference density (RHO_0) for surface + ! density and GINT calculate internal densities. + ! ts_profile - use temperature and salinity profiles + ! (read from COORD_FILE) to set layer densities. + ! USER - call a user modified routine. +TS_RANGE_T_LIGHT = 25.0 ! [degC] default = 10.0 + ! The initial temperature of the lightest layer when COORD_CONFIG is set to + ! ts_range. +TS_RANGE_T_DENSE = 3.0 ! [degC] default = 10.0 + ! The initial temperature of the densest layer when COORD_CONFIG is set to + ! ts_range. +TS_RANGE_RESOLN_RATIO = 5.0 ! [nondim] default = 1.0 + ! The ratio of density space resolution in the densest part of the range to that + ! in the lightest part of the range when COORD_CONFIG is set to ts_range. Values + ! greater than 1 increase the resolution of the denser water. + +! === module MOM_state_initialization === +THICKNESS_CONFIG = "benchmark" ! default = "uniform" + ! A string that determines how the initial layer thicknesses are specified for a + ! new run: + ! file - read interface heights from the file specified + ! thickness_file - read thicknesses from the file specified + ! by (THICKNESS_FILE). + ! coord - determined by ALE coordinate. + ! uniform - uniform thickness layers evenly distributed + ! between the surface and MAXIMUM_DEPTH. + ! list - read a list of positive interface depths. + ! DOME - use a slope and channel configuration for the + ! DOME sill-overflow test case. + ! ISOMIP - use a configuration for the + ! ISOMIP test case. + ! benchmark - use the benchmark test case thicknesses. + ! Neverworld - use the Neverworld test case thicknesses. + ! search - search a density profile for the interface + ! densities. This is not yet implemented. + ! circle_obcs - the circle_obcs test case is used. + ! DOME2D - 2D version of DOME initialization. + ! adjustment2d - 2D lock exchange thickness ICs. + ! sloshing - sloshing gravity thickness ICs. + ! seamount - no motion test with seamount ICs. + ! dumbbell - sloshing channel ICs. + ! soliton - Equatorial Rossby soliton. + ! rossby_front - a mixed layer front in thermal wind balance. + ! USER - call a user modified routine. + +! === module benchmark_initialize_thickness === +TS_CONFIG = "benchmark" ! + ! A string that determines how the initial tempertures and salinities are + ! specified for a new run: + ! file - read velocities from the file specified + ! by (TS_FILE). + ! fit - find the temperatures that are consistent with + ! the layer densities and salinity S_REF. + ! TS_profile - use temperature and salinity profiles + ! (read from TS_FILE) to set layer densities. + ! benchmark - use the benchmark test case T & S. + ! linear - linear in logical layer space. + ! DOME2D - 2D DOME initialization. + ! ISOMIP - ISOMIP initialization. + ! adjustment2d - 2d lock exchange T/S ICs. + ! sloshing - sloshing mode T/S ICs. + ! seamount - no motion test with seamount ICs. + ! dumbbell - sloshing channel ICs. + ! rossby_front - a mixed layer front in thermal wind balance. + ! SCM_CVMix_tests - used in the SCM CVMix tests. + ! USER - call a user modified routine. + +! === module MOM_diag_mediator === + +! === module MOM_lateral_mixing_coeffs === +USE_VARIABLE_MIXING = True ! [Boolean] default = False + ! If true, the variable mixing code will be called. This allows diagnostics to + ! be created even if the scheme is not used. If KHTR_SLOPE_CFF>0 or + ! KhTh_Slope_Cff>0, this is set to true regardless of what is in the parameter + ! file. +USE_VISBECK = True ! [Boolean] default = False + ! If true, use the Visbeck et al. (1997) formulation for + ! thickness diffusivity. +RESOLN_SCALED_KH = True ! [Boolean] default = False + ! If true, the Laplacian lateral viscosity is scaled away when the first + ! baroclinic deformation radius is well resolved. +RESOLN_SCALED_KHTH = True ! [Boolean] default = False + ! If true, the interface depth diffusivity is scaled away when the first + ! baroclinic deformation radius is well resolved. +RESOLN_SCALED_KHTR = True ! [Boolean] default = False + ! If true, the epipycnal tracer diffusivity is scaled away when the first + ! baroclinic deformation radius is well resolved. +KHTH_SLOPE_CFF = 0.1 ! [nondim] default = 0.0 + ! The nondimensional coefficient in the Visbeck formula for the interface depth + ! diffusivity +KHTR_SLOPE_CFF = 0.1 ! [nondim] default = 0.0 + ! The nondimensional coefficient in the Visbeck formula for the epipycnal tracer + ! diffusivity +VARMIX_KTOP = 6 ! [nondim] default = 2 + ! The layer number at which to start vertical integration of S*N for purposes of + ! finding the Eady growth rate. +VISBECK_L_SCALE = 3.0E+04 ! [m] default = 0.0 + ! The fixed length scale in the Visbeck formula. + +! === module MOM_set_visc === +PRANDTL_TURB = 0.0 ! [nondim] default = 1.0 + ! The turbulent Prandtl number applied to shear instability. +DYNAMIC_VISCOUS_ML = True ! [Boolean] default = False + ! If true, use a bulk Richardson number criterion to determine the mixed layer + ! thickness for viscosity. +ML_OMEGA_FRAC = 1.0 ! [nondim] default = 0.0 + ! When setting the decay scale for turbulence, use this fraction of the absolute + ! rotation rate blended with the local value of f, as sqrt((1-of)*f^2 + + ! of*4*omega^2). +HBBL = 10.0 ! [m] + ! The thickness of a bottom boundary layer with a viscosity of KVBBL if + ! BOTTOMDRAGLAW is not defined, or the thickness over which near-bottom + ! velocities are averaged for the drag law if BOTTOMDRAGLAW is defined but + ! LINEAR_DRAG is not. +DRAG_BG_VEL = 0.1 ! [m s-1] default = 0.0 + ! DRAG_BG_VEL is either the assumed bottom velocity (with LINEAR_DRAG) or an + ! unresolved velocity that is combined with the resolved velocity to estimate + ! the velocity magnitude. DRAG_BG_VEL is only used when BOTTOMDRAGLAW is + ! defined. +BBL_THICK_MIN = 0.1 ! [m] default = 0.0 + ! The minimum bottom boundary layer thickness that can be used with + ! BOTTOMDRAGLAW. This might be Kv/(cdrag*drag_bg_vel) to give Kv as the minimum + ! near-bottom viscosity. +KV = 1.0E-04 ! [m2 s-1] + ! The background kinematic viscosity in the interior. The molecular value, ~1e-6 + ! m2 s-1, may be used. + +! === module MOM_thickness_diffuse === +KHTH = 1.0 ! [m2 s-1] default = 0.0 + ! The background horizontal thickness diffusivity. +KHTH_MAX = 900.0 ! [m2 s-1] default = 0.0 + ! The maximum horizontal thickness diffusivity. + +! === module MOM_dynamics_split_RK2 === + +! === module MOM_continuity === + +! === module MOM_continuity_PPM === +ETA_TOLERANCE = 1.0E-06 ! [m] default = 1.1E-09 + ! The tolerance for the differences between the barotropic and baroclinic + ! estimates of the sea surface height due to the fluxes through each face. The + ! total tolerance for SSH is 4 times this value. The default is + ! 0.5*NK*ANGSTROM, and this should not be set less than about + ! 10^-15*MAXIMUM_DEPTH. +VELOCITY_TOLERANCE = 0.001 ! [m s-1] default = 3.0E+08 + ! The tolerance for barotropic velocity discrepancies between the barotropic + ! solution and the sum of the layer thicknesses. + +! === module MOM_CoriolisAdv === +BOUND_CORIOLIS = True ! [Boolean] default = False + ! If true, the Coriolis terms at u-points are bounded by the four estimates of + ! (f+rv)v from the four neighboring v-points, and similarly at v-points. This + ! option would have no effect on the SADOURNY Coriolis scheme if it were + ! possible to use centered difference thickness fluxes. + +! === module MOM_PressureForce === + +! === module MOM_PressureForce_AFV === + +! === module MOM_hor_visc === +AH_VEL_SCALE = 0.05 ! [m s-1] default = 0.0 + ! The velocity scale which is multiplied by the cube of the grid spacing to + ! calculate the biharmonic viscosity. The final viscosity is the largest of this + ! scaled viscosity, the Smagorinsky and Leith viscosities, and AH. +SMAGORINSKY_AH = True ! [Boolean] default = False + ! If true, use a biharmonic Smagorinsky nonlinear eddy viscosity. +SMAG_BI_CONST = 0.06 ! [nondim] default = 0.0 + ! The nondimensional biharmonic Smagorinsky constant, typically 0.015 - 0.06. + +! === module MOM_vert_friction === +MAXVEL = 10.0 ! [m s-1] default = 3.0E+08 + ! The maximum velocity allowed before the velocity components are truncated. + +! === module MOM_barotropic === +BOUND_BT_CORRECTION = True ! [Boolean] default = False + ! If true, the corrective pseudo mass-fluxes into the barotropic solver are + ! limited to values that require less than maxCFL_BT_cont to be accommodated. +NONLINEAR_BT_CONTINUITY = True ! [Boolean] default = False + ! If true, use nonlinear transports in the barotropic continuity equation. This + ! does not apply if USE_BT_CONT_TYPE is true. +BT_PROJECT_VELOCITY = True ! [Boolean] default = False + ! If true, step the barotropic velocity first and project out the velocity + ! tendency by 1+BEBT when calculating the transport. The default (false) is to + ! use a predictor continuity step to find the pressure field, and then to do a + ! corrector continuity step using a weighted average of the old and new + ! velocities, with weights of (1-BEBT) and BEBT. +BEBT = 0.2 ! [nondim] default = 0.1 + ! BEBT determines whether the barotropic time stepping uses the forward-backward + ! time-stepping scheme or a backward Euler scheme. BEBT is valid in the range + ! from 0 (for a forward-backward treatment of nonrotating gravity waves) to 1 + ! (for a backward Euler treatment). In practice, BEBT must be greater than about + ! 0.05. +DTBT = -0.95 ! [s or nondim] default = -0.98 + ! The barotropic time step, in s. DTBT is only used with the split explicit time + ! stepping. To set the time step automatically based the maximum stable value + ! use 0, or a negative value gives the fraction of the stable value. Setting + ! DTBT to 0 is the same as setting it to -0.98. The value of DTBT that will + ! actually be used is an integer fraction of DT, rounding down. + +! === module MOM_mixed_layer_restrat === +MIXEDLAYER_RESTRAT = True ! [Boolean] default = False + ! If true, a density-gradient dependent re-stratifying flow is imposed in the + ! mixed layer. Can be used in ALE mode without restriction but in layer mode can + ! only be used if BULKMIXEDLAYER is true. +FOX_KEMPER_ML_RESTRAT_COEF = 5.0 ! [nondim] default = 0.0 + ! A nondimensional coefficient that is proportional to the ratio of the + ! deformation radius to the dominant lengthscale of the submesoscale mixed layer + ! instabilities, times the minimum of the ratio of the mesoscale eddy kinetic + ! energy to the large-scale geostrophic kinetic energy or 1 plus the square of + ! the grid spacing over the deformation radius, as detailed by Fox-Kemper et al. + ! (2010) + +! === module MOM_diagnostics === + +! === module MOM_diabatic_driver === +! The following parameters are used for diabatic processes. + +! === module MOM_entrain_diffusive === +MAX_ENT_IT = 20 ! default = 5 + ! The maximum number of iterations that may be used to calculate the interior + ! diapycnal entrainment. +TOLERANCE_ENT = 1.0E-05 ! [m] default = 1.341640786499874E-05 + ! The tolerance with which to solve for entrainment values. + +! === module MOM_set_diffusivity === + +! === module MOM_bkgnd_mixing === +! Adding static vertical background mixing coefficients +KD = 2.0E-05 ! [m2 s-1] default = 0.0 + ! The background diapycnal diffusivity of density in the interior. Zero or the + ! molecular value, ~1e-7 m2 s-1, may be used. + +! === module MOM_kappa_shear === +! Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008 +USE_JACKSON_PARAM = True ! [Boolean] default = False + ! If true, use the Jackson-Hallberg-Legg (JPO 2008) shear mixing + ! parameterization. +MAX_RINO_IT = 25 ! [nondim] default = 50 + ! The maximum number of iterations that may be used to estimate the Richardson + ! number driven mixing. + +! === module MOM_diabatic_aux === +! The following parameters are used for auxiliary diabatic processes. +RECLAIM_FRAZIL = False ! [Boolean] default = True + ! If true, try to use any frazil heat deficit to cool any overlying layers down + ! to the freezing point, thereby avoiding the creation of thin ice when the SST + ! is above the freezing point. + +! === module MOM_mixed_layer === +MSTAR = 0.3 ! [nondim] default = 1.2 + ! The ratio of the friction velocity cubed to the TKE input to the mixed layer. +BULK_RI_ML = 0.05 ! [nondim] + ! The efficiency with which mean kinetic energy released by mechanically forced + ! entrainment of the mixed layer is converted to turbulent kinetic energy. +ABSORB_ALL_SW = True ! [Boolean] default = False + ! If true, all shortwave radiation is absorbed by the ocean, instead of passing + ! through to the bottom mud. +TKE_DECAY = 10.0 ! [nondim] default = 2.5 + ! TKE_DECAY relates the vertical rate of decay of the TKE available for + ! mechanical entrainment to the natural Ekman depth. +HMIX_MIN = 2.0 ! [m] default = 0.0 + ! The minimum mixed layer depth if the mixed layer depth is determined + ! dynamically. +LIMIT_BUFFER_DETRAIN = True ! [Boolean] default = False + ! If true, limit the detrainment from the buffer layers to not be too different + ! from the neighbors. +DEPTH_LIMIT_FLUXES = 0.1 ! [m] default = 0.2 + ! The surface fluxes are scaled away when the total ocean depth is less than + ! DEPTH_LIMIT_FLUXES. +CORRECT_ABSORPTION_DEPTH = True ! [Boolean] default = False + ! If true, the average depth at which penetrating shortwave radiation is + ! absorbed is adjusted to match the average heating depth of an exponential + ! profile by moving some of the heating upward in the water column. + +! === module MOM_opacity === +PEN_SW_SCALE = 15.0 ! [m] default = 0.0 + ! The vertical absorption e-folding depth of the penetrating shortwave + ! radiation. +PEN_SW_FRAC = 0.42 ! [nondim] default = 0.0 + ! The fraction of the shortwave radiation that penetrates below the surface. + +! === module MOM_tracer_advect === + +! === module MOM_tracer_hor_diff === +KHTR = 1.0 ! [m2 s-1] default = 0.0 + ! The background along-isopycnal tracer diffusivity. +KHTR_MAX = 900.0 ! [m2 s-1] default = 0.0 + ! The maximum along-isopycnal tracer diffusivity. +DIFFUSE_ML_TO_INTERIOR = True ! [Boolean] default = False + ! If true, enable epipycnal mixing between the surface boundary layer and the + ! interior. +ML_KHTR_SCALE = 0.0 ! [nondim] default = 1.0 + ! With Diffuse_ML_interior, the ratio of the truly horizontal diffusivity in the + ! mixed layer to the epipycnal diffusivity. The valid range is 0 to 1. + +! === module MOM_sum_output === +MAXTRUNC = 5000 ! [truncations save_interval-1] default = 0 + ! The run will be stopped, and the day set to a very large value if the velocity + ! is truncated more than MAXTRUNC times between energy saves. Set MAXTRUNC to 0 + ! to stop if there is any truncation of velocities. +ENERGYSAVEDAYS = 0.25 ! [days] default = 1.0 + ! The interval in units of TIMEUNIT between saves of the energies of the run and + ! other globally summed diagnostics. + +! === module MOM_surface_forcing === +BUOY_CONFIG = "linear" ! default = "zero" + ! The character string that indicates how buoyancy forcing is specified. Valid + ! options include (file), (zero), (linear), (USER), (BFB) and (NONE). +WIND_CONFIG = "gyres" ! default = "zero" + ! The character string that indicates how wind forcing is specified. Valid + ! options include (file), (2gyre), (1gyre), (gyres), (zero), and (USER). +TAUX_SIN_AMP = 0.1 ! [Pa] default = 0.0 + ! With the gyres wind_config, the sine amplitude in the zonal wind stress + ! profile: B in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L). +TAUX_N_PIS = 1.0 ! [nondim] default = 0.0 + ! With the gyres wind_config, the number of gyres in the zonal wind stress + ! profile: n in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L). +RESTOREBUOY = True ! [Boolean] default = False + ! If true, the buoyancy fluxes drive the model back toward some specified + ! surface state with a rate given by FLUXCONST. +FLUXCONST = 0.5 ! [m day-1] default = 0.0 + ! The constant that relates the restoring surface fluxes to the relative surface + ! anomalies (akin to a piston velocity). Note the non-MKS units. +SST_NORTH = 27.0 ! [deg C] default = 0.0 + ! With buoy_config linear, the sea surface temperature at the northern end of + ! the domain toward which to to restore. +SST_SOUTH = 3.0 ! [deg C] default = 0.0 + ! With buoy_config linear, the sea surface temperature at the southern end of + ! the domain toward which to to restore. +GUST_CONST = 0.02 ! [Pa] default = 0.0 + ! The background gustiness in the winds. + +! === module MOM_main (MOM_driver) === +DT_FORCING = 3600.0 ! [s] default = 900.0 + ! The time step for changing forcing, coupling with other components, or + ! potentially writing certain diagnostics. The default value is given by DT. +DAYMAX = 3.0 ! [days] + ! The final time of the whole simulation, in units of TIMEUNIT seconds. This + ! also sets the potential end time of the present run segment if the end time is + ! not set via ocean_solo_nml in input.nml. +RESTART_CONTROL = 3 ! default = 1 + ! An integer whose bits encode which restart files are written. Add 2 (bit 1) + ! for a time-stamped file, and odd (bit 0) for a non-time-stamped file. A + ! non-time-stamped restart file is saved at the end of the run segment for any + ! non-negative value. +RESTINT = 365.0 ! [days] default = 0.0 + ! The interval between saves of the restart file in units of TIMEUNIT. Use 0 + ! (the default) to not save incremental restart files at all. + +! === module MOM_write_cputime === +MAXCPU = 2.88E+04 ! [wall-clock seconds] default = -1.0 + ! The maximum amount of cpu time per processor for which MOM should run before + ! saving a restart file and quitting with a return value that indicates that a + ! further run is required to complete the simulation. If automatic restarts are + ! not desired, use a negative value for MAXCPU. MAXCPU has units of wall-clock + ! seconds, so the actual CPU time used is larger by a factor of the number of + ! processors used. + +! Debugging parameters set to non-default values +U_TRUNC_FILE = "U_velocity_truncations" ! default = "" + ! The absolute path to a file into which the accelerations leading to zonal + ! velocity truncations are written. Undefine this for efficiency if this + ! diagnostic is not needed. +V_TRUNC_FILE = "V_velocity_truncations" ! default = "" + ! The absolute path to a file into which the accelerations leading to meridional + ! velocity truncations are written. Undefine this for efficiency if this + ! diagnostic is not needed. diff --git a/.testing/p0/MOM_override b/.testing/p0/MOM_override new file mode 100644 index 0000000000..e69de29bb2 diff --git a/.testing/p0/diag_table b/.testing/p0/diag_table new file mode 100644 index 0000000000..68c71dd2c4 --- /dev/null +++ b/.testing/p0/diag_table @@ -0,0 +1,91 @@ +"MOM benchmark Experiment" +1 1 1 0 0 0 +"prog", 1,"days",1,"days","time", +#"ave_prog", 5,"days",1,"days","Time",365,"days" +#"cont", 5,"days",1,"days","Time",365,"days" + +#This is the field section of the diag_table. + +# Prognostic Ocean fields: +#========================= + +"ocean_model","u","u","prog","all",.false.,"none",2 +"ocean_model","v","v","prog","all",.false.,"none",2 +"ocean_model","h","h","prog","all",.false.,"none",1 +"ocean_model","e","e","prog","all",.false.,"none",2 +"ocean_model","temp","temp","prog","all",.false.,"none",2 +#"ocean_model","salt","salt","prog","all",.false.,"none",2 + +#"ocean_model","u","u","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","v","v","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","h","h","ave_prog_%4yr_%3dy","all",.true.,"none",1 +#"ocean_model","e","e","ave_prog_%4yr_%3dy","all",.true.,"none",2 + +# Auxilary Tracers: +#================== +#"ocean_model","vintage","vintage","prog_%4yr_%3dy","all",.false.,"none",2 +#"ocean_model","age","age","prog_%4yr_%3dy","all",.false.,"none",2 + +# Continuity Equation Terms: +#=========================== +#"ocean_model","dhdt","dhdt","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","wd","wd","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","uh","uh","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","vh","vh","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","h_rho","h_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","uh_rho","uh_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","vh_rho","vh_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","uhGM_rho","uhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","vhGM_rho","vhGM_rho","cont_%4yr_%3dy","all",.true.,"none",2 + +# +# Tracer Fluxes: +#================== +#"ocean_model","T_adx", "T_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","T_ady", "T_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","T_diffx","T_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","T_diffy","T_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_adx", "S_adx", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_ady", "S_ady", "ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_diffx","S_diffx","ave_prog_%4yr_%3dy","all",.true.,"none",2 +#"ocean_model","S_diffy","S_diffy","ave_prog_%4yr_%3dy","all",.true.,"none",2 + +# testing +# ======= +#"ocean_model","Kv_u","Kv_u","prog","all",.false.,"none",2 +#"ocean_model","Kv_v","Kv_v","prog","all",.false.,"none",2 + +#============================================================================================= +# +#===- This file can be used with diag_manager/v2.0a (or higher) ==== +# +# +# FORMATS FOR FILE ENTRIES (not all input values are used) +# ------------------------ +# +#"file_name", output_freq, "output_units", format, "time_units", "time_long_name", ... +# (opt) new_file_frequecy, (opt) "new_file_freq_units", "new_file_start_date" +# +# +#output_freq: > 0 output frequency in "output_units" +# = 0 output frequency every time step +# =-1 output frequency at end of run +# +#output_units = units used for output frequency +# (years, months, days, minutes, hours, seconds) +# +#time_units = units used to label the time axis +# (days, minutes, hours, seconds) +# +# +# FORMAT FOR FIELD ENTRIES (not all input values are used) +# ------------------------ +# +#"module_name", "field_name", "output_name", "file_name" "time_sampling", time_avg, "other_opts", packing +# +#time_avg = .true. or .false. +# +#packing = 1 double precision +# = 2 float +# = 4 packed 16-bit integers +# = 8 packed 1-byte (not tested?) diff --git a/.testing/p0/input.nml b/.testing/p0/input.nml new file mode 100644 index 0000000000..41555b8822 --- /dev/null +++ b/.testing/p0/input.nml @@ -0,0 +1,22 @@ + &MOM_input_nml + output_directory = './', + input_filename = 'n' + restart_input_dir = 'INPUT/', + restart_output_dir = 'RESTART/', + parameter_filename = 'MOM_input', + 'MOM_override' / + + &diag_manager_nml + / + + &fms_nml + clock_grain='ROUTINE' + clock_flags='SYNC' + !domains_stack_size = 955296 + domains_stack_size = 14256000 + stack_size =0 / + +!&ocean_solo_nml +! hours = 1 +! !days = 1 +!/ diff --git a/.testing/tc0/MOM_input b/.testing/tc0/MOM_input index ff64c55803..e4d1694e72 100644 --- a/.testing/tc0/MOM_input +++ b/.testing/tc0/MOM_input @@ -230,7 +230,6 @@ ENERGYSAVEDAYS = 1.0 DIAG_AS_CHKSUM = True DEBUG = True -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True USE_GM_WORK_BUG = True ! [Boolean] default = True FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False diff --git a/.testing/tc1/MOM_input b/.testing/tc1/MOM_input index 68674f7a86..151c093ff9 100644 --- a/.testing/tc1/MOM_input +++ b/.testing/tc1/MOM_input @@ -575,7 +575,6 @@ ENERGYSAVEDAYS = 0.125 ! [days] default = 3600.0 DIAG_AS_CHKSUM = True DEBUG = True USE_PSURF_IN_EOS = False ! [Boolean] default = False -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True INTERPOLATE_RES_FN = True ! [Boolean] default = True GILL_EQUATORIAL_LD = False ! [Boolean] default = False diff --git a/.testing/tc2/MOM_input b/.testing/tc2/MOM_input index 1818390192..ca84d1c382 100644 --- a/.testing/tc2/MOM_input +++ b/.testing/tc2/MOM_input @@ -610,7 +610,6 @@ DIAG_AS_CHKSUM = True DEBUG = True USE_GM_WORK_BUG = False USE_PSURF_IN_EOS = False ! [Boolean] default = False -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True REMAP_UV_USING_OLD_ALG = True ! [Boolean] default = True USE_LAND_MASK_FOR_HVISC = False ! [Boolean] default = False diff --git a/.testing/tc3/MOM_input b/.testing/tc3/MOM_input index 9112898b4c..a034960d1e 100644 --- a/.testing/tc3/MOM_input +++ b/.testing/tc3/MOM_input @@ -469,7 +469,6 @@ ENERGYSAVEDAYS = 3.0 ! [hours] default = 1.44E+04 ! energies of the run and other globally summed diagnostics. DIAG_AS_CHKSUM = True DEBUG = True -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True OBC_RADIATION_MAX = 10.0 ! [nondim] default = 10.0 GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True USE_GM_WORK_BUG = True ! [Boolean] default = True diff --git a/.testing/tc4/MOM_input b/.testing/tc4/MOM_input index 04598a9dc9..e33bf40bf6 100644 --- a/.testing/tc4/MOM_input +++ b/.testing/tc4/MOM_input @@ -409,7 +409,6 @@ DIAG_AS_CHKSUM = True DEBUG = True USE_PSURF_IN_EOS = False ! [Boolean] default = False -DEFAULT_2018_ANSWERS = True ! [Boolean] default = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True INTERPOLATE_RES_FN = True ! [Boolean] default = True GILL_EQUATORIAL_LD = False ! [Boolean] default = False diff --git a/.testing/tools/compare_clocks.py b/.testing/tools/compare_clocks.py new file mode 100755 index 0000000000..77198fda6a --- /dev/null +++ b/.testing/tools/compare_clocks.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python +import argparse +import json + +# Ignore timers below this threshold (in seconds) +DEFAULT_THRESHOLD = 0.05 + +# Thresholds for reporting +DT_WARN = 0.10 # Slowdown warning +DT_FAIL = 0.25 # Slowdown abort + +ANSI_RED = '\033[31m' +ANSI_GREEN = '\033[32m' +ANSI_YELLOW = '\033[33m' +ANSI_RESET = '\033[0m' + + +def main(): + desc = ( + 'Compare two FMS clock output files and report any differences within ' + 'a defined threshold.' + ) + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('expt') + parser.add_argument('ref') + parser.add_argument('--threshold') + parser.add_argument('--verbose', action='store_true') + args = parser.parse_args() + + threshold = float(args.threshold) if args.threshold else DEFAULT_THRESHOLD + verbose = args.verbose + + clock_cmp = {} + + print('{:35s}{:8s} {:8s}'.format('', 'Profile', 'Reference')) + print() + + with open(args.expt) as log_expt, open(args.ref) as log_ref: + clocks_expt = json.load(log_expt)['clocks'] + clocks_ref = json.load(log_ref)['clocks'] + + # Gather timers which appear in both clocks + clock_tags = [clk for clk in clocks_expt if clk in clocks_ref] + + for clk in clock_tags: + clock_cmp[clk] = {} + + # For now, we only comparge tavg, the rank-averaged timing + rec = 'tavg' + + t_expt = clocks_expt[clk][rec] + t_ref = clocks_ref[clk][rec] + + # Compare the relative runtimes + if all(t > threshold for t in (t_expt, t_ref)): + dclk = (t_expt - t_ref) / t_ref + else: + dclk = 0. + clock_cmp[clk][rec] = dclk + + # Skip trivially low clocks + if all(t < threshold for t in (t_expt, t_ref)) and not verbose: + continue + + # Report the time differences + ansi_color = ANSI_RESET + + if abs(t_expt - t_ref) > threshold: + if dclk > DT_FAIL: + ansi_color = ANSI_RED + elif dclk > DT_WARN: + ansi_color = ANSI_YELLOW + elif dclk < -DT_WARN: + ansi_color = ANSI_GREEN + + print('{}{}: {:7.3f}s, {:7.3f}s ({:.1f}%){}'.format( + ansi_color, + ' ' * (32 - len(clk)) + clk, + t_expt, + t_ref, + 100. * dclk, + ANSI_RESET, + )) + + +if __name__ == '__main__': + main() diff --git a/.testing/tools/compare_perf.py b/.testing/tools/compare_perf.py new file mode 100755 index 0000000000..e4a651c709 --- /dev/null +++ b/.testing/tools/compare_perf.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python +import argparse +import json + +# Ignore timers below this threshold (in seconds) +DEFAULT_THRESHOLD = 0.05 + +# Thresholds for reporting +DT_WARN = 0.10 # Slowdown warning +DT_FAIL = 0.25 # Slowdown abort + +ANSI_RED = '\033[31m' +ANSI_GREEN = '\033[32m' +ANSI_YELLOW = '\033[33m' +ANSI_RESET = '\033[0m' + + +def main(): + desc = ( + 'Compare two FMS clock output files and report any differences within ' + 'a defined threshold.' + ) + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('expt') + parser.add_argument('ref') + parser.add_argument('--threshold') + parser.add_argument('--verbose', action='store_true') + args = parser.parse_args() + + threshold = float(args.threshold) if args.threshold else DEFAULT_THRESHOLD + verbose = args.verbose + + clock_cmp = {} + + print('{:35s}{:8s} {:8s}'.format('', 'Profile', 'Reference')) + print() + + with open(args.expt) as profile_expt, open(args.ref) as profile_ref: + perf_expt = json.load(profile_expt) + perf_ref = json.load(profile_ref) + + events = [ev for ev in perf_expt if ev in perf_ref] + + for event in events: + # For now, only report the times + if event not in ('task-clock', 'cpu-clock'): + continue + + count_expt = perf_expt[event]['count'] + count_ref = perf_ref[event]['count'] + + symbols_expt = perf_expt[event]['symbol'] + symbols_ref = perf_ref[event]['symbol'] + + symbols = [ + s for s in symbols_expt + if s in symbols_ref + and not s.startswith('0x') + ] + + for symbol in symbols: + t_expt = float(symbols_expt[symbol]) / 1e9 + t_ref = float(symbols_ref[symbol]) / 1e9 + + # Compare the relative runtimes + if all(t > threshold for t in (t_expt, t_ref)): + dclk = (t_expt - t_ref) / t_ref + else: + dclk = 0. + + # Skip trivially low clocks + if all(t < threshold for t in (t_expt, t_ref)) and not verbose: + continue + + # Report the time differences + ansi_color = ANSI_RESET + + if abs(t_expt - t_ref) > threshold: + if dclk > DT_FAIL: + ansi_color = ANSI_RED + elif dclk > DT_WARN: + ansi_color = ANSI_YELLOW + elif dclk < -DT_WARN: + ansi_color = ANSI_GREEN + + # Remove module name + sname = symbol.split('_MOD_', 1)[-1] + + # Strip version from glibc calls + sname = sname.split('@')[0] + + # Remove GCC optimization renaming + sname = sname.replace('.constprop.0', '') + + if len(sname) > 32: + sname = sname[:29] + '...' + + print('{}{}: {:7.3f}s, {:7.3f}s ({:.1f}%){}'.format( + ansi_color, + ' ' * (32 - len(sname)) + sname, + t_expt, + t_ref, + 100. * dclk, + ANSI_RESET, + )) + + +if __name__ == '__main__': + main() diff --git a/.testing/tools/parse_fms_clocks.py b/.testing/tools/parse_fms_clocks.py new file mode 100755 index 0000000000..b57fc481ab --- /dev/null +++ b/.testing/tools/parse_fms_clocks.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python +import argparse +import collections +import json +import os +import sys + +import f90nml + +record_type = collections.defaultdict(lambda: float) +for rec in ('grain', 'pemin', 'pemax',): + record_type[rec] = int + + +def main(): + desc = 'Parse MOM6 model stdout and return clock data in JSON format.' + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('--format', '-f', action='store_true') + parser.add_argument('--dir', '-d') + parser.add_argument('log') + args = parser.parse_args() + + config = {} + + if args.dir: + # Gather model configuration + input_nml = os.path.join(args.dir, 'input.nml') + nml = f90nml.read(input_nml) + config['input.nml'] = nml.todict() + + parameter_filenames = [ + ('params', 'MOM_parameter_doc.all'), + ('layout', 'MOM_parameter_doc.layout'), + ('debug', 'MOM_parameter_doc.debugging'), + ] + for key, fname in parameter_filenames: + config[key] = {} + with open(os.path.join(args.dir, fname)) as param_file: + params = parse_mom6_param(param_file) + config[key].update(params) + + # Get log path + if os.path.isfile(args.log): + log_path = args.log + elif os.path.isfile(os.path.join(args.dir, args.log)): + log_path = os.path.join(args.dir, args.log) + else: + sys.exit('stdout log not found.') + + # Parse timings + with open(log_path) as log: + clocks = parse_clocks(log) + + config['clocks'] = clocks + + if args.format: + print(json.dumps(config, indent=4)) + else: + print(json.dumps(config)) + + +def parse_mom6_param(param_file): + params = {} + for line in param_file: + param_stmt = line.split('!')[0].strip() + if param_stmt: + key, val = [s.strip() for s in param_stmt.split('=')] + + # TODO: Convert to equivalent Python types + if val in ('True', 'False'): + params[key] = bool(val) + else: + params[key] = val + + return params + + +def parse_clocks(log): + clock_start_msg = 'Tabulating mpp_clock statistics across' + clock_end_msg = 'MPP_STACK high water mark=' + + fields = [] + for line in log: + if line.startswith(clock_start_msg): + npes = line.lstrip(clock_start_msg).split()[0] + + # Get records + fields = [] + line = next(log) + + # Skip blank lines + while line.isspace(): + line = next(log) + + fields = line.split() + + # Exit this loop, begin clock parsing + break + + clocks = {} + for line in log: + # Treat MPP_STACK usage as end of clock reports + if line.lstrip().startswith(clock_end_msg): + break + + record = line.split()[-len(fields):] + + clk = line.split(record[0])[0].strip() + clocks[clk] = {} + + for fld, rec in zip(fields, record): + rtype = record_type[fld] + clocks[clk][fld] = rtype(rec) + + return clocks + + +if __name__ == '__main__': + main() diff --git a/.testing/tools/parse_perf.py b/.testing/tools/parse_perf.py new file mode 100755 index 0000000000..b86b1cc106 --- /dev/null +++ b/.testing/tools/parse_perf.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python +import argparse +import collections +import json +import os +import shlex +import subprocess +import sys + + +def main(): + desc = 'Parse perf.data and return in JSON format.' + + parser = argparse.ArgumentParser(description=desc) + parser.add_argument('--format', '-f', action='store_true') + parser.add_argument('data') + args = parser.parse_args() + + profile = parse_perf_report(args.data) + + if args.format: + print(json.dumps(profile, indent=4)) + else: + print(json.dumps(profile)) + + +def parse_perf_report(perf_data_path): + profile = {} + + cmd = shlex.split( + 'perf report -s symbol,period -i {}'.format(perf_data_path) + ) + with subprocess.Popen(cmd, stdout=subprocess.PIPE, text=True) as proc: + event_name = None + for line in proc.stdout: + # Skip blank lines: + if not line or line.isspace(): + continue + + # Set the current event + if line.startswith('# Samples: '): + event_name = line.split()[-1].strip("'") + + # Remove perf modifiers for now + event_name = event_name.rsplit(':', 1)[0] + + profile[event_name] = {} + profile[event_name]['symbol'] = {} + + # Get total count + elif line.startswith('# Event count '): + event_count = int(line.split()[-1]) + profile[event_name]['count'] = event_count + + # skip all other 'comment' lines + elif line.startswith('#'): + continue + + # get per-symbol count + else: + tokens = line.split() + symbol = tokens[2] + period = int(tokens[3]) + + profile[event_name]['symbol'][symbol] = period + + return profile + + +if __name__ == '__main__': + main() diff --git a/ac/configure.ac b/ac/configure.ac index 9cb7147846..3d1af81b05 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -65,7 +65,7 @@ AS_IF([test "x$with_driver" != "x"], MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1 AC_ARG_WITH([framework], AS_HELP_STRING([--with-framework=fms1|fms2], [Select the model framework])) -AS_CASE([with_framework], +AS_CASE(["$with_framework"], [fms1], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1], [fms2], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS2], [MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1] diff --git a/config_src/infra/FMS1/MOM_interp_infra.F90 b/config_src/infra/FMS1/MOM_interp_infra.F90 index 170573f7ec..774f6a67d2 100644 --- a/config_src/infra/FMS1/MOM_interp_infra.F90 +++ b/config_src/infra/FMS1/MOM_interp_infra.F90 @@ -231,7 +231,7 @@ end subroutine time_interp_extern_3d !> initialize an external field integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & - threading, ierr, ignore_axis_atts ) + threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency ) character(len=*), intent(in) :: file !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -246,13 +246,20 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, logical, optional, intent(in) :: ignore_axis_atts !< If present and true, do not issue a !! fatal error if the axis Cartesian attribute is !! not set to a recognized value. + logical, optional, intent(in) :: correct_leap_year_inconsistency !< If present and true, + !! then if, (1) a calendar containing leap years + !! is in use, and (2) the modulo time period of the + !! data is an integer number of years, then map + !! a model date of Feb 29. onto a common year on Feb. 28. if (present(MOM_Domain)) then init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts) + verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency) else init_extern_field = init_external_field(file, fieldname, domain=domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts) + verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency) endif end function init_extern_field diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index 170573f7ec..b02beca313 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -231,7 +231,7 @@ end subroutine time_interp_extern_3d !> initialize an external field integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & - threading, ierr, ignore_axis_atts ) + threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency ) character(len=*), intent(in) :: file !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -246,13 +246,22 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, logical, optional, intent(in) :: ignore_axis_atts !< If present and true, do not issue a !! fatal error if the axis Cartesian attribute is !! not set to a recognized value. + logical, optional, intent(in) :: correct_leap_year_inconsistency !< If present and true, + !! then if, (1) a calendar containing leap years + !! is in use, and (2) the modulo time period of the + !! data is an integer number of years, then map + !! a model date of Feb 29. onto a common year on Feb. 28. + + if (present(MOM_Domain)) then init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts) + verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency) else init_extern_field = init_external_field(file, fieldname, domain=domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts) + verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency) endif end function init_extern_field diff --git a/docs/discrete_space.rst b/docs/discrete_space.rst index b954915256..08a41a5f2d 100644 --- a/docs/discrete_space.rst +++ b/docs/discrete_space.rst @@ -13,7 +13,6 @@ algorithm. :maxdepth: 2 api/generated/pages/Discrete_Grids - api/generated/pages/Finite_Difference_Operators api/generated/pages/PPM api/generated/pages/Discrete_Coriolis api/generated/pages/Discrete_PG diff --git a/docs/parameterizations_vertical.rst b/docs/parameterizations_vertical.rst index 27285034d7..0d22787294 100644 --- a/docs/parameterizations_vertical.rst +++ b/docs/parameterizations_vertical.rst @@ -26,6 +26,8 @@ Kappa-shear Internal-tide driven mixing The schemes of :cite:`st_laurent2002`, :cite:`polzin2009`, and :cite:`melet2012`, are all implemented through MOM_set_diffusivity and MOM_diabatic_driver. + :ref:`Internal_Tidal_Mixing` + Vertical friction ----------------- diff --git a/docs/tracers.rst b/docs/tracers.rst index 6190fe096d..8b5a21ee12 100644 --- a/docs/tracers.rst +++ b/docs/tracers.rst @@ -9,4 +9,5 @@ Tracers in MOM6 api/generated/pages/Horizontal_Diffusion.rst api/generated/pages/Vertical_Diffusion.rst api/generated/pages/Passive_Tracers + api/generated/pages/Frazil_Ice diff --git a/docs/zotero.bib b/docs/zotero.bib index f0e1a3b44d..957097f217 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -847,6 +847,16 @@ @article{melet2012 pages = {602--615} } +@article{simmons2004, + title = {Tidally driven mixing in a numerical model of the ocean general circulation}, + volume = {6}, + author = {Simmons, H. L. and S. R. Jayne and L. C. St. Laurent and A. J. Weaver}, + year = {2004}, + journal = {Ocean Modell.}, + pages = {245--263}, + doi = {10.1016/S1463-5003(03)00011-8} +} + @article{polzin2009, title = {An abyssal recipe}, volume = {30}, @@ -863,6 +873,15 @@ @article{polzin2009 pages = {298--309} } +@article{polzin2004, + title = {Idealized solutions for the energy balance of the finescale internal wave field}, + volume = {34}, + journal = {J. Phys. Oceanogr.}, + author = {Polzin, Kurt L.}, + year = {2004}, + pages = {231--246} +} + @article{white2009, title = {High-order regridding-remapping schemes for continuous isopycnal and generalized coordinates in ocean models}, volume = {228}, @@ -2524,3 +2543,13 @@ @article{hallberg2005 volume = {8}, doi = {10.1016/j.ocemod.2004.01.001} } + +@article{bell1975, + author = {T. H. Bell}, + year = {1975}, + title = {Lee wavews in stratified flows with simple harmonic time dependence"}, + journal = {J. Fluid Mech.}, + volume = {67}, + pages = {705--722} +} + diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 06850dca97..0681cb9e74 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1276,7 +1276,13 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call enable_averages(dtdia, Time_end_thermo, CS%diag) if (associated(CS%odaCS)) then - call apply_oda_tracer_increments(US%T_to_s*dtdia, G, GV, tv, h, CS%odaCS) + if (CS%debug) then + call MOM_thermo_chksum("Pre-oda ", tv, G, US, haloshift=0) + endif + call apply_oda_tracer_increments(US%T_to_s*dtdia, Time_end_thermo, G, GV, tv, h, CS%odaCS) + if (CS%debug) then + call MOM_thermo_chksum("Post-oda ", tv, G, US, haloshift=0) + endif endif if (associated(fluxes%p_surf) .or. associated(fluxes%p_surf_full)) then @@ -2807,7 +2813,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & (LEN_TRIM(dirs%input_filename) == 1)) if (CS%ensemble_ocean) then - call init_oda(Time, G, GV, CS%odaCS) + call init_oda(Time, G, GV, CS%diag, CS%odaCS) endif !### This could perhaps go here instead of in finish_MOM_initialization? diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 51f9a5cb85..4d19459bc7 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2685,6 +2685,16 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k) enddo ; enddo ; enddo endif + if (present(ADp) .and. (associated(ADp%visc_rem_u))) then + do k=1,nz ; do j=js,je ; do I=is-1,ie + ADp%visc_rem_u(I,j,k) = visc_rem_u(I,j,k) + enddo ; enddo ; enddo + endif + if (present(ADp) .and. (associated(ADp%visc_rem_u))) then + do k=1,nz ; do J=js-1,je ; do i=is,ie + ADp%visc_rem_v(i,J,k) = visc_rem_v(i,J,k) + enddo ; enddo ; enddo + endif if (G%nonblocking_updates) then if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 291703a242..a168fe1319 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -169,10 +169,12 @@ module MOM_dynamics_split_RK2 integer :: id_h_PFu = -1, id_h_PFv = -1 integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1 integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1 + integer :: id_PFu_visc_rem = -1, id_PFv_visc_rem = -1 ! integer :: id_hf_CAu = -1, id_hf_CAv = -1 integer :: id_h_CAu = -1, id_h_CAv = -1 integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1 integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1 + integer :: id_CAu_visc_rem = -1, id_CAv_visc_rem = -1 ! Split scheme only. integer :: id_uav = -1, id_vav = -1 @@ -181,6 +183,7 @@ module MOM_dynamics_split_RK2 integer :: id_h_u_BT_accel = -1, id_h_v_BT_accel = -1 integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1 integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1 + integer :: id_u_BT_accel_visc_rem = -1, id_v_BT_accel_visc_rem = -1 !>@} type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the @@ -360,6 +363,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s real, dimension(SZI_(G),SZJB_(G)) :: & intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2]. + ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv], + real, allocatable, dimension(:,:,:) :: & + PFu_visc_rem, PFv_visc_rem, & ! Pressure force accel. x visc_rem_[uv] [L T-2 ~> m s-2]. + CAu_visc_rem, CAv_visc_rem, & ! Coriolis force accel. x visc_rem_[uv] [L T-2 ~> m s-2]. + u_BT_accel_visc_rem, v_BT_accel_visc_rem ! barotropic correction accel. x visc_rem_[uv] [L T-2 ~> m s-2]. + real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. logical :: dyn_p_surf @@ -1108,6 +1117,61 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s deallocate(h_v_BT_accel) endif + if (CS%id_PFu_visc_rem > 0) then + allocate(PFu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + PFu_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%ADp%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_PFu_visc_rem, PFu_visc_rem, CS%diag) + deallocate(PFu_visc_rem) + endif + if (CS%id_PFv_visc_rem > 0) then + allocate(PFv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + PFv_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%ADp%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_PFv_visc_rem, PFv_visc_rem, CS%diag) + deallocate(PFv_visc_rem) + endif + if (CS%id_CAu_visc_rem > 0) then + allocate(CAu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + CAu_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%ADp%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_CAu_visc_rem, CAu_visc_rem, CS%diag) + deallocate(CAu_visc_rem) + endif + if (CS%id_CAv_visc_rem > 0) then + allocate(CAv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + CAv_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%ADp%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_CAv_visc_rem, CAv_visc_rem, CS%diag) + deallocate(CAv_visc_rem) + endif + if (CS%id_u_BT_accel_visc_rem > 0) then + allocate(u_BT_accel_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + u_BT_accel_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_u_BT_accel_visc_rem, u_BT_accel_visc_rem, CS%diag) + deallocate(u_BT_accel_visc_rem) + endif + if (CS%id_v_BT_accel_visc_rem > 0) then + allocate(v_BT_accel_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + v_BT_accel_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_v_BT_accel_visc_rem, v_BT_accel_visc_rem, CS%diag) + deallocate(v_BT_accel_visc_rem) + endif + if (CS%debug) then call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) @@ -1612,6 +1676,33 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) + CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, & + 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if (CS%id_PFu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) + CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, & + 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) + + CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, & + 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if (CS%id_CAu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) + CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, & + 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) + + CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, & + 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if (CS%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz) + CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, & + 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz) + id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 3f98a97052..ab85db8baf 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -195,6 +195,9 @@ module MOM_variables real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points + real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points + real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points + end type accel_diag_ptrs !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances. diff --git a/src/core/_Finite_difference.dox b/src/core/_Finite_difference.dox deleted file mode 100644 index ecbd37d8b7..0000000000 --- a/src/core/_Finite_difference.dox +++ /dev/null @@ -1,5 +0,0 @@ -/*! \page Finite_Difference_Operators Finite Difference Operators - -\brief Finite Difference Operators - -*/ diff --git a/src/core/_Sea_ice.dox b/src/core/_Sea_ice.dox index bec05af17c..232bac1bb8 100644 --- a/src/core/_Sea_ice.dox +++ b/src/core/_Sea_ice.dox @@ -1,5 +1,11 @@ /*! \page Sea_Ice Sea Ice Considerations -\section Frazil Ice Formation +\section section_seaice Sea Ice Considerations + +For realistic domains, it is assumed that MOM6 will be run in a coupled mode, such that either the +sea-ice model or the coupler will be computing atmospheric bulk fluxes and passing them to the ocean. +Likewise, MOM6 can compute the frazil ice formation as described in \ref section_frazil, which it +then passes to the sea-ice model, expecting to get back the rejected brine or melted fresh water in +return. */ diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 44b05cc081..cc16a25fc3 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -2341,7 +2341,6 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS) call safe_alloc_ptr(ADp%gradKEu,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%gradKEv,isd,ied,JsdB,JedB,nz) endif - if (associated(CS%KE_visc)) then call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) @@ -2395,7 +2394,7 @@ subroutine MOM_diagnostics_end(CS, ADp, CDp) if (associated(CS%vhGM_Rlay)) deallocate(CS%vhGM_Rlay) if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu) - if (associated(ADp%gradKEu)) deallocate(ADp%gradKEu) + if (associated(ADp%gradKEv)) deallocate(ADp%gradKEv) if (associated(ADp%du_dt_visc)) deallocate(ADp%du_dt_visc) if (associated(ADp%dv_dt_visc)) deallocate(ADp%dv_dt_visc) if (associated(ADp%du_dt_str)) deallocate(ADp%du_dt_str) diff --git a/src/equation_of_state/_Equation_of_State.dox b/src/equation_of_state/_Equation_of_State.dox index 2e18b49f54..791c7001b1 100644 --- a/src/equation_of_state/_Equation_of_State.dox +++ b/src/equation_of_state/_Equation_of_State.dox @@ -36,7 +36,7 @@ Compute the required quantities using the equation of state from \cite jackett19 Compute the required quantities using the equation of state from [TEOS-10](http://www.teos-10.org/). -\section TFREEZE Freezing Temperature of Sea Water +\section section_TFREEZE Freezing Temperature of Sea Water There are three choices for computing the freezing point of sea water: diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 1f8d45e88d..c24bd82f73 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -277,24 +277,24 @@ subroutine register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS) "ice sheet/shelf u-velocity", "m s-1", hor_grid='Bu') call register_restart_field(CS%v_shelf, "v_shelf", .false., restart_CS, & "ice sheet/shelf v-velocity", "m s-1", hor_grid='Bu') - call register_restart_field(CS%u_bdry_val, "u_bdry", .false., restart_CS, & + call register_restart_field(CS%u_bdry_val, "u_bdry_val", .false., restart_CS, & "ice sheet/shelf boundary u-velocity", "m s-1", hor_grid='Bu') - call register_restart_field(CS%v_bdry_val, "v_bdry", .false., restart_CS, & + call register_restart_field(CS%v_bdry_val, "v_bdry_val", .false., restart_CS, & "ice sheet/shelf boundary v-velocity", "m s-1", hor_grid='Bu') - call register_restart_field(CS%u_face_mask_bdry, "u_bdry_mask", .false., restart_CS, & + call register_restart_field(CS%u_face_mask_bdry, "u_face_mask_bdry", .false., restart_CS, & "ice sheet/shelf boundary u-mask", "nondim", hor_grid='Bu') - call register_restart_field(CS%v_face_mask_bdry, "v_bdry_mask", .false., restart_CS, & + call register_restart_field(CS%v_face_mask_bdry, "v_face_mask_bdry", .false., restart_CS, & "ice sheet/shelf boundary v-mask", "nondim", hor_grid='Bu') call register_restart_field(CS%OD_av, "OD_av", .true., restart_CS, & "Average open ocean depth in a cell","m") call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, & "fractional degree of grounding", "nondim") - call register_restart_field(CS%C_basal_friction, "tau_b_beta", .true., restart_CS, & + call register_restart_field(CS%C_basal_friction, "C_basal_friction", .true., restart_CS, & "basal sliding coefficients", "Pa (m s-1)^n_sliding") - call register_restart_field(CS%AGlen_visc, "A_Glen", .true., restart_CS, & + call register_restart_field(CS%AGlen_visc, "AGlen_visc", .true., restart_CS, & "ice-stiffness parameter", "Pa-3 s-1") - call register_restart_field(CS%h_bdry_val, "h_bdry", .false., restart_CS, & + call register_restart_field(CS%h_bdry_val, "h_bdry_val", .false., restart_CS, & "ice thickness at the boundary","m") endif @@ -503,6 +503,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ call pass_var(CS%ground_frac,G%domain) call pass_var(CS%ice_visc,G%domain) call pass_var(CS%basal_traction, G%domain) + call pass_var(CS%AGlen_visc, G%domain) call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) endif @@ -533,30 +534,31 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif ! initialize basal friction coefficients - call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file) - call pass_var(CS%C_basal_friction, G%domain) + if (new_sim) then + call initialize_ice_C_basal_friction(CS%C_basal_friction, G, US, param_file) + call pass_var(CS%C_basal_friction, G%domain) - ! initialize ice-stiffness AGlen - call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file) - call pass_var(CS%AGlen_visc, G%domain) + ! initialize ice-stiffness AGlen + call initialize_ice_AGlen(CS%AGlen_visc, G, US, param_file) + call pass_var(CS%AGlen_visc, G%domain) - !initialize boundary conditions - call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & + !initialize boundary conditions + call initialize_ice_shelf_boundary_from_file(CS%u_face_mask_bdry, CS%v_face_mask_bdry, & CS%u_bdry_val, CS%v_bdry_val, CS%umask, CS%vmask, CS%h_bdry_val, & CS%thickness_bdry_val, ISS%hmask, ISS%h_shelf, G, US, param_file ) - call pass_var(ISS%hmask, G%domain) - call pass_var(CS%h_bdry_val, G%domain) - call pass_var(CS%thickness_bdry_val, G%domain) - call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE) - call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE) - - !initialize ice flow velocities from file - call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, & + call pass_var(ISS%hmask, G%domain) + call pass_var(CS%h_bdry_val, G%domain) + call pass_var(CS%thickness_bdry_val, G%domain) + call pass_vector(CS%u_bdry_val, CS%v_bdry_val, G%domain, TO_ALL, BGRID_NE) + call pass_vector(CS%u_face_mask_bdry, CS%v_face_mask_bdry, G%domain, TO_ALL, BGRID_NE) + + !initialize ice flow velocities from file + call initialize_ice_flow_from_file(CS%bed_elev,CS%u_shelf, CS%v_shelf,CS%ground_frac, ISS%hmask,ISS%h_shelf, & G, US, param_file) - call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) - call pass_var(CS%bed_elev, G%domain,CENTER) - call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) - + call pass_vector(CS%u_shelf, CS%v_shelf, G%domain, TO_ALL, BGRID_NE) + call pass_var(CS%bed_elev, G%domain,CENTER) + call update_velocity_masks(CS, G, ISS%hmask, CS%umask, CS%vmask, CS%u_face_mask, CS%v_face_mask) + endif ! Register diagnostics. CS%id_u_shelf = register_diag_field('ice_shelf_model','u_shelf',CS%diag%axesB1, Time, & 'x-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) @@ -564,16 +566,15 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'y-velocity of ice', 'm yr-1', conversion=365.0*86400.0*US%L_T_to_m_s) ! I think that the conversion factors for the next two diagnostics are wrong. - RWH CS%id_taudx_shelf = register_diag_field('ice_shelf_model','taudx_shelf',CS%diag%axesB1, Time, & - 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) + 'x-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa) CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesB1, Time, & - 'y-driving stress of ice', 'kPa', conversion=1.e-9*US%L_T_to_m_s) + 'y-driving stress of ice', 'kPa', conversion=1.e-9*US%RL2_T2_to_Pa) CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesB1, Time, & 'mask for u-nodes', 'none') CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesB1, Time, & 'mask for v-nodes', 'none') CS%id_ground_frac = register_diag_field('ice_shelf_model','ice_ground_frac',CS%diag%axesT1, Time, & 'fraction of cell that is grounded', 'none') -! CS%id_col_thick = register_diag_field('ice_shelf_model','col_thick',CS%diag%axesT1, Time, & 'ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) CS%id_visc_shelf = register_diag_field('ice_shelf_model','ice_visc',CS%diag%axesT1, Time, & @@ -582,12 +583,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'taub', 'Pa yr m-1', conversion=1e-6*US%Z_to_m) CS%id_OD_av = register_diag_field('ice_shelf_model','OD_av',CS%diag%axesT1, Time, & 'intermediate ocean column thickness passed to ice model', 'm', conversion=US%Z_to_m) - if (new_sim) then - call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") - call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) - endif endif + call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: initialize ice velocity.") + call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) end subroutine initialize_ice_shelf_dyn @@ -683,46 +682,46 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled coupled_GL = .false. if (present(ocean_mass) .and. present(coupled_grounding)) coupled_GL = coupled_grounding ! - call ice_shelf_advect(CS, ISS, G, time_step, Time) - CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step - if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. - - if (coupled_GL) then - call update_OD_ffrac(CS, G, US, ocean_mass, update_ice_vel) - elseif (update_ice_vel) then - call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) - endif + call ice_shelf_advect(CS, ISS, G, time_step, Time) + CS%elapsed_velocity_time = CS%elapsed_velocity_time + time_step + if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) update_ice_vel = .true. + if (coupled_GL) then + call update_OD_ffrac(CS, G, US, ocean_mass, update_ice_vel) + elseif (update_ice_vel) then + call update_OD_ffrac_uncoupled(CS, G, ISS%h_shelf(:,:)) + endif - if (update_ice_vel) then - call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) - endif + + if (update_ice_vel) then + call ice_shelf_solve_outer(CS, ISS, G, US, CS%u_shelf, CS%v_shelf,CS%taudx_shelf,CS%taudy_shelf, iters, Time) + endif ! call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) - if (update_ice_vel) then - call enable_averages(CS%elapsed_velocity_time, Time, CS%diag) - if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) - if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) - if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) + if (update_ice_vel) then + call enable_averages(CS%elapsed_velocity_time, Time, CS%diag) + if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) + if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) + if (CS%id_v_shelf > 0) call post_data(CS%id_v_shelf, CS%v_shelf, CS%diag) ! if (CS%id_t_shelf > 0) call post_data(CS%id_t_shelf,CS%t_shelf,CS%diag) - if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) - if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) - if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) - if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) - if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) - if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag) + if (CS%id_taudx_shelf > 0) call post_data(CS%id_taudx_shelf, CS%taudx_shelf, CS%diag) + if (CS%id_taudy_shelf > 0) call post_data(CS%id_taudy_shelf, CS%taudy_shelf, CS%diag) + if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac,CS%diag) + if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) + if (CS%id_visc_shelf > 0) call post_data(CS%id_visc_shelf, CS%ice_visc,CS%diag) + if (CS%id_taub > 0) call post_data(CS%id_taub, CS%basal_traction,CS%diag) !! - if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) - if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) - if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask,CS%u_face_mask_bdry,CS%diag) - if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask,CS%v_face_mask_bdry,CS%diag) + if (CS%id_u_mask > 0) call post_data(CS%id_u_mask,CS%umask,CS%diag) + if (CS%id_v_mask > 0) call post_data(CS%id_v_mask,CS%vmask,CS%diag) + if (CS%id_ufb_mask > 0) call post_data(CS%id_ufb_mask,CS%u_face_mask_bdry,CS%diag) + if (CS%id_vfb_mask > 0) call post_data(CS%id_vfb_mask,CS%v_face_mask_bdry,CS%diag) ! if (CS%id_t_mask > 0) call post_data(CS%id_t_mask,CS%tmask,CS%diag) - call disable_averaging(CS%diag) + call disable_averaging(CS%diag) - CS%elapsed_velocity_time = 0.0 - endif + CS%elapsed_velocity_time = 0.0 + endif end subroutine update_ice_shelf @@ -869,13 +868,13 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite CS%ground_frac(:,:) = 0.0 allocate(Phisub(nsub,nsub,2,2,2,2)) ; Phisub(:,:,:,:,:,:) = 0.0 - do j=G%jsc,G%jec - do i=G%isc,G%iec - if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then + do j=G%jsc,G%jec + do i=G%isc,G%iec + if (rhoi_rhow * ISS%h_shelf(i,j) - G%bathyT(i,j) > 0) then float_cond(i,j) = 1.0 CS%ground_frac(i,j) = 1.0 - endif - enddo + endif + enddo enddo call calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, CS%OD_av) @@ -1043,16 +1042,15 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf,taudx,taudy, ite call MOM_mesg(mesg, 5) if (err_max <= CS%nonlinear_tolerance * err_init) then - write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init - call MOM_mesg(mesg) - write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" -! call MOM_mesg(mesg, 5) - call MOM_mesg(mesg) exit endif enddo + write(mesg,*) "ice_shelf_solve_outer: nonlinear fractional residual = ", err_max/err_init + call MOM_mesg(mesg) + write(mesg,*) "ice_shelf_solve_outer: exiting nonlinear solve after ",iter," iterations" + call MOM_mesg(mesg) deallocate(Phi) deallocate(Phisub) @@ -1086,6 +1084,7 @@ subroutine ice_shelf_solve_inner(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, H !! iterations have converged to the specified tolerance integer, intent(out) :: iters !< The number of iterations used in the solver. type(time_type), intent(in) :: Time !< The current model time + character(len=160) :: mesg ! The text of an error message real, dimension(8,4,SZDI_(G),SZDJ_(G)), & intent(in) :: Phi !< The gradients of bilinear basis elements at Gaussian !! quadrature points surrounding the cell vertices [L-1 ~> m-1]. @@ -2586,7 +2585,7 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (u_shlf(I,J-1) + (u_shlf(I-1,J-1) + u_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) vy = ((v_shlf(I,J) + (v_shlf(I-1,J)+ v_shlf(I+1,J))) - & (v_shlf(I,J-1) + (v_shlf(I-1,J-1)+ v_shlf(I+1,J-1)))) / (3*G%dyT(i,j)) -! CS%ice_visc(i,j) =1e14*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging +! CS%ice_visc(i,j) =1e15*(G%areaT(i,j) * ISS%h_shelf(i,j)) ! constant viscocity for debugging CS%ice_visc(i,j) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & (US%s_to_T**2 * (ux**2 + vy**2 + ux*vy + 0.25*(uy+vx)**2 + eps_min**2))**((1.-n_g)/(2.*n_g)) endif @@ -2938,7 +2937,7 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face do j=js,G%jed do i=is,G%ied - if (hmask(i,j) == 1) then + if ((hmask(i,j) == 1) .OR. (hmask(i,j) == 3)) then umask(I,j) = 1. vmask(I,j) = 1. @@ -2947,10 +2946,10 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%u_face_mask_bdry(I-1+k,j))) case (3) - ! vmask(I-1+k,J-1)=0. + vmask(I-1+k,J-1)=3. u_face_mask(I-1+k,j)=3. umask(I-1+k,J)=3. - !vmask(I-1+k,J)=0. + vmask(I-1+k,J)=3. vmask(I-1+k,J)=3. case (2) u_face_mask(I-1+k,j)=2. @@ -2973,9 +2972,9 @@ subroutine update_velocity_masks(CS, G, hmask, umask, vmask, u_face_mask, v_face select case (int(CS%v_face_mask_bdry(i,J-1+k))) case (3) vmask(I-1,J-1+k)=3. - umask(I-1,J-1+k)=0. + umask(I-1,J-1+k)=3. vmask(I,J-1+k)=3. - umask(I,J-1+k)=0. + umask(I,J-1+k)=3. v_face_mask(i,J-1+k)=3. case (2) v_face_mask(i,J-1+k)=2. diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index f3a5f210fc..73db36596e 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -101,7 +101,7 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U ! This subroutine reads ice thickness and area from a file and puts it into ! h_shelf [Z ~> m] and area_shelf_h [L2 ~> m2] (and dimensionless) and updates hmask character(len=200) :: filename,thickness_file,inputdir ! Strings for file/path - character(len=200) :: thickness_varname, area_varname ! Variable name in file + character(len=200) :: thickness_varname, area_varname, hmask_varname ! Variable name in file character(len=40) :: mdl = "initialize_ice_thickness_from_file" ! This subroutine's name. integer :: i, j, isc, jsc, iec, jec real :: len_sidestress, mask, udh @@ -125,45 +125,46 @@ subroutine initialize_ice_thickness_from_file(h_shelf, area_shelf_h, hmask, G, U call get_param(PF, mdl, "ICE_AREA_VARNAME", area_varname, & "The name of the area variable in ICE_THICKNESS_FILE.", & default="area_shelf_h") - + hmask_varname="h_mask" if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & " initialize_topography_from_file: Unable to open "//trim(filename)) call MOM_read_data(filename, trim(thickness_varname), h_shelf, G%Domain, scale=US%m_to_Z) call MOM_read_data(filename,trim(area_varname), area_shelf_h, G%Domain, scale=US%m_to_L**2) - + call MOM_read_data(filename, trim(hmask_varname), hmask, G%Domain) isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec - do j=jsc,jec - do i=isc,iec + if (len_sidestress > 0.) then + do j=jsc,jec + do i=isc,iec ! taper ice shelf in area where there is no sidestress - ! but do not interfere with hmask - if ((G%geoLonCv(i,j) > len_sidestress).and. & - (len_sidestress > 0.)) then - udh = exp(-(G%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j) - if (udh <= 25.0) then - h_shelf(i,j) = 0.0 - area_shelf_h(i,j) = 0.0 - else + if (G%geoLonCv(i,j) > len_sidestress) then + udh = exp(-(G%geoLonCv(i,j)-len_sidestress)/5.0) * h_shelf(i,j) + if (udh <= 25.0) then + h_shelf(i,j) = 0.0 + area_shelf_h(i,j) = 0.0 + else h_shelf(i,j) = udh + endif endif - endif ! update thickness mask - if (area_shelf_h (i,j) >= G%areaT(i,j)) then - hmask(i,j) = 1. - elseif (area_shelf_h (i,j) == 0.0) then - hmask(i,j) = 0. - elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= G%areaT(i,j))) then - hmask(i,j) = 2. - else - call MOM_error(FATAL,mdl// " AREA IN CELL OUT OF RANGE") - endif + if (area_shelf_h (i,j) >= G%areaT(i,j)) then + hmask(i,j) = 1. + area_shelf_h(i,j)=G%areaT(i,j) + elseif (area_shelf_h (i,j) == 0.0) then + hmask(i,j) = 0. + elseif ((area_shelf_h(i,j) > 0) .and. (area_shelf_h(i,j) <= G%areaT(i,j))) then + hmask(i,j) = 2. + else + call MOM_error(FATAL,mdl// " AREA IN CELL OUT OF RANGE") + endif + enddo enddo - enddo - + endif end subroutine initialize_ice_thickness_from_file !> Initialize ice shelf thickness for a channel configuration diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 8057234cdc..dd9c46ff90 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -9,23 +9,32 @@ module MOM_oda_driver_mod use MOM_domains, only : domain2d, global_field, get_domain_extent use MOM_domains, only : pass_var, redistribute_array, broadcast_domain use MOM_diag_mediator, only : register_diag_field, diag_axis_init, post_data +use MOM_diag_mediator, only : enable_averaging, disable_averaging +use MOM_diag_mediator, only : diag_update_remap_grids use MOM_ensemble_manager, only : get_ensemble_id, get_ensemble_size use MOM_ensemble_manager, only : get_ensemble_pelist, get_ensemble_filter_pelist use MOM_error_handler, only : stdout, stdlog, MOM_error use MOM_io, only : SINGLE_FILE +use MOM_interp_infra, only : init_extern_field, get_external_field_info +use MOM_interp_infra, only : time_interp_extern use MOM_time_manager, only : time_type, real_to_time, get_date use MOM_time_manager, only : operator(+), operator(>=), operator(/=) use MOM_time_manager, only : operator(==), operator(<) - +use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end, cpu_clock_id +use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer ! ODA Modules use ocean_da_types_mod, only : grid_type, ocean_profile_type, ocean_control_struct use ocean_da_core_mod, only : ocean_da_core_init, get_profiles +!This preprocessing directive enables the SPEAR online ensemble data assimilation +!configuration. Existing community based APIs for data assimilation are currently +!called offline for forecast applications using information read from a MOM6 state file. +!The SPEAR configuration (https://doi.org/10.1029/2020MS002149) calculated increments +!efficiently online. A community-based set of APIs should be implemented in place +!of the CPP directive when this is available. #ifdef ENABLE_ECDA use eakf_oda_mod, only : ensemble_filter #endif -use write_ocean_obs_mod, only : open_profile_file -use write_ocean_obs_mod, only : write_profile,close_profile_file -use kdtree, only : kd_root !# JEDI +use kdtree, only : kd_root !# A kd-tree object using JEDI APIs ! MOM Modules use MOM_io, only : slasher, MOM_read_data use MOM_diag_mediator, only : diag_ctrl, set_axes_info @@ -52,7 +61,16 @@ module MOM_oda_driver_mod implicit none ; private public :: init_oda, oda_end, set_prior_tracer, get_posterior_tracer -public :: set_analysis_time, oda, save_obs_diff, apply_oda_tracer_increments +public :: set_analysis_time, oda, apply_oda_tracer_increments + +!>@{ CPU time clock ID +integer :: id_clock_oda_init +integer :: id_clock_oda_filter +integer :: id_clock_bias_adjustment +integer :: id_clock_apply_increments +integer :: id_clock_oda_prior +integer :: id_clock_oda_posterior +!>@} #include @@ -61,13 +79,23 @@ module MOM_oda_driver_mod type(domain2d), pointer :: mpp_domain => NULL() !< pointer to a domain2d end type ptr_mpp_domain +!> A structure containing integer handles for bias adjustment of tracers +type :: INC_CS + integer :: fldno = 0 !< The number of tracers + integer :: T_id !< The integer handle for the temperature file + integer :: S_id !< The integer handle for the salinity file +end type INC_CS + !> Control structure that contains a transpose of the ocean state across ensemble members. type, public :: ODA_CS ; private type(ocean_control_struct), pointer :: Ocean_prior=> NULL() !< ensemble ocean prior states in DA space type(ocean_control_struct), pointer :: Ocean_posterior=> NULL() !< ensemble ocean posterior states !! or increments to prior in DA space + type(ocean_control_struct), pointer :: Ocean_increment=> NULL() !< A separate structure for + !! increment diagnostics integer :: nk !< number of vertical layers used for DA type(ocean_grid_type), pointer :: Grid => NULL() !< MOM6 grid type and decomposition for the DA + type(ocean_grid_type), pointer :: G => NULL() !< MOM6 grid type and decomposition for the model type(MOM_domain_type), pointer, dimension(:) :: domains => NULL() !< Pointer to mpp_domain objects !! for ensemble members type(verticalGrid_type), pointer :: GV => NULL() !< vertical grid for DA @@ -78,12 +106,17 @@ module MOM_oda_driver_mod type(grid_type), pointer :: oda_grid !< local tracer grid real, pointer, dimension(:,:,:) :: h => NULL() ! m or kg m-2] for DA type(thermo_var_ptrs), pointer :: tv => NULL() !< pointer to thermodynamic variables + type(thermo_var_ptrs), pointer :: tv_bc => NULL() !< pointer to thermodynamic bias correction integer :: ni !< global i-direction grid size integer :: nj !< global j-direction grid size logical :: reentrant_x !< grid is reentrant in the x direction logical :: reentrant_y !< grid is reentrant in the y direction logical :: tripolar_N !< grid is folded at its north edge logical :: symmetric !< Values at C-grid locations are symmetric + logical :: use_basin_mask !< If true, use a basin file to delineate weakly coupled ocean basins + logical :: do_bias_adjustment !< If true, use spatio-temporally varying climatological tendency + !! adjustment for Temperature and Salinity + real :: bias_adjustment_multiplier !< A scaling for the bias adjustment integer :: assim_method !< Method: NO_ASSIM,EAKF_ASSIM or OI_ASSIM integer :: ensemble_size !< Size of the ensemble integer :: ensemble_id = 0 !< id of the current ensemble member @@ -99,7 +132,10 @@ module MOM_oda_driver_mod type(regridding_CS) :: regridCS !< ALE control structure for regridding type(remapping_CS) :: remapCS !< ALE control structure for remapping type(time_type) :: Time !< Current Analysis time - type(diag_ctrl) :: diag_cs ! NULL() ! initialize First_guess (prior) and Analysis grid !! information for all ensemble members -subroutine init_oda(Time, G, GV, CS) +subroutine init_oda(Time, G, GV, diag_CS, CS) type(time_type), intent(in) :: Time !< The current model time. type(ocean_grid_type), pointer :: G !< domain and grid information for ocean model type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(diag_ctrl), target, intent(inout) :: diag_CS !< A pointer to a diagnostic control structure type(ODA_CS), pointer, intent(inout) :: CS !< The DA control structure ! Local variables @@ -133,6 +170,7 @@ subroutine init_oda(Time, G, GV, CS) integer :: isg,ieg,jsg,jeg integer :: idg_offset, jdg_offset integer :: stdout_unit + integer, dimension(4) :: fld_sz character(len=32) :: assim_method integer :: npes_pm, ens_info(6), ni, nj character(len=128) :: mesg @@ -140,14 +178,21 @@ subroutine init_oda(Time, G, GV, CS) character(len=30) :: coord_mode character(len=200) :: inputdir, basin_file logical :: reentrant_x, reentrant_y, tripolar_N, symmetric + character(len=80) :: bias_correction_file, inc_file if (associated(CS)) call MOM_error(FATAL, 'Calling oda_init with associated control structure') allocate(CS) + + id_clock_oda_init=cpu_clock_id('(ODA initialization)') + id_clock_oda_prior=cpu_clock_id('(ODA setting prior)') + id_clock_oda_filter=cpu_clock_id('(ODA filter computation)') + id_clock_oda_posterior=cpu_clock_id('(ODA getting posterior)') + call cpu_clock_begin(id_clock_oda_init) + ! Use ens1 parameters , this could be changed at a later time ! if it were desirable to have alternate parameters, e.g. for the grid ! for the analysis call get_MOM_input(PF,dirs,ensemble_num=0) - call unit_scaling_init(PF, CS%US) call get_param(PF, "MOM", "ASSIM_METHOD", assim_method, & @@ -167,6 +212,20 @@ subroutine init_oda(Time, G, GV, CS) "Use tripolar connectivity at the northern edge of the "//& "domain. With TRIPOLAR_N, NIGLOBAL must be even.", & default=.false.) + call get_param(PF,"MOM", "APPLY_TRACER_TENDENCY_ADJUSTMENT", CS%do_bias_adjustment, & + "If true, add a spatio-temporally varying climatological adjustment "//& + "to temperature and salinity.", & + default=.false.) + if (CS%do_bias_adjustment) then + call get_param(PF,"MOM", "TRACER_ADJUSTMENT_FACTOR", CS%bias_adjustment_multiplier, & + "A multiplicative scaling factor for the climatological tracer tendency adjustment ", & + default=1.0) + endif + call get_param(PF,"MOM", "USE_BASIN_MASK", CS%use_basin_mask, & + "If true, add a basin mask to delineate weakly connected "//& + "ocean basins for the purpose of data assimilation.", & + default=.false.) + call get_param(PF,"MOM", "NIGLOBAL", CS%ni, & "The total number of thickness grid points in the "//& "x-direction in the physical domain.") @@ -200,19 +259,19 @@ subroutine init_oda(Time, G, GV, CS) call set_PElist(CS%filter_pelist) allocate(CS%domains(CS%ensemble_size)) - CS%domains(CS%ensemble_id)%mpp_domain => G%Domain%mpp_domain + CS%domains(CS%ensemble_id)%mpp_domain => G%Domain%mpp_domain ! this should go away do n=1,CS%ensemble_size if (.not. associated(CS%domains(n)%mpp_domain)) allocate(CS%domains(n)%mpp_domain) - call set_rootPE(CS%ensemble_pelist(n,1)) + call set_rootPE(CS%ensemble_pelist(n,1)) ! this line is not in Feiyu's version (needed?) call broadcast_domain(CS%domains(n)%mpp_domain) enddo - call set_rootPE(CS%filter_pelist(1)) + call set_rootPE(CS%filter_pelist(1)) ! this line is not in Feiyu's version (needed?) + CS%G => G allocate(CS%Grid) ! params NIHALO_ODA, NJHALO_ODA set the DA halo size call MOM_domains_init(CS%Grid%Domain,PF,param_suffix='_ODA') allocate(HI) - call hor_index_init(CS%Grid%Domain, HI, PF, & - local_indexing=.false.) ! Use global indexing for DA + call hor_index_init(CS%Grid%Domain, HI, PF) call verticalGridInit( PF, CS%GV, CS%US ) allocate(dG) call create_dyn_horgrid(dG, HI) @@ -222,7 +281,7 @@ subroutine init_oda(Time, G, GV, CS) call MOM_initialize_coord(CS%GV, CS%US, PF, .false., & dirs%output_directory, tv_dummy, dG%max_depth) call ALE_init(PF, CS%GV, CS%US, dG%max_depth, CS%ALE_CS) - call MOM_grid_init(CS%Grid, PF, global_indexing=.true.) + call MOM_grid_init(CS%Grid, PF) call ALE_updateVerticalGridType(CS%ALE_CS, CS%GV) call copy_dyngrid_to_MOM_grid(dG, CS%Grid, CS%US) CS%mpp_domain => CS%Grid%Domain%mpp_domain @@ -233,7 +292,9 @@ subroutine init_oda(Time, G, GV, CS) call init_ocean_ensemble(CS%Ocean_prior,CS%Grid,CS%GV,CS%ensemble_size) allocate(CS%Ocean_posterior) call init_ocean_ensemble(CS%Ocean_posterior,CS%Grid,CS%GV,CS%ensemble_size) - allocate(CS%tv) + allocate(CS%Ocean_increment) + call init_ocean_ensemble(CS%Ocean_increment,CS%Grid,CS%GV,CS%ensemble_size) + call get_param(PF, 'oda_driver', "REGRIDDING_COORDINATE_MODE", coord_mode, & "Coordinate mode for vertical regridding.", & @@ -241,76 +302,80 @@ subroutine init_oda(Time, G, GV, CS) call initialize_regridding(CS%regridCS, CS%GV, CS%US, dG%max_depth,PF,'oda_driver',coord_mode,'','') call initialize_remapping(CS%remapCS,'PLM') call set_regrid_params(CS%regridCS, min_thickness=0.) + isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed + ! breaking with the MOM6 convention and using global indices - call get_domain_extent(G%Domain,is,ie,js,je,isd,ied,jsd,jed,& - isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric) - isd=isd+idg_offset; ied=ied+idg_offset - jsd=jsd+jdg_offset; jed=jed+jdg_offset - !call mpp_get_data_domain(G%Domain%mpp_domain,isd,ied,jsd,jed) + !call get_domain_extent(G%Domain,is,ie,js,je,isd,ied,jsd,jed,& + ! isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric) + !isd=isd+idg_offset; ied=ied+idg_offset ! using global indexing within the DA module + !jsd=jsd+jdg_offset; jed=jed+jdg_offset ! TODO: switch to local indexing? (mjh) + if (.not. associated(CS%h)) then - allocate(CS%h(isd:ied,jsd:jed,CS%GV%ke)); CS%h(:,:,:)=0.0 + allocate(CS%h(isd:ied,jsd:jed,CS%GV%ke)); CS%h(:,:,:)=CS%GV%Angstrom_m*CS%GV%H_to_m ! assign thicknesses call ALE_initThicknessToCoord(CS%ALE_CS,G,CS%GV,CS%h) endif + allocate(CS%tv) allocate(CS%tv%T(isd:ied,jsd:jed,CS%GV%ke)); CS%tv%T(:,:,:)=0.0 allocate(CS%tv%S(isd:ied,jsd:jed,CS%GV%ke)); CS%tv%S(:,:,:)=0.0 - call set_axes_info(CS%Grid, CS%GV, CS%US, PF, CS%diag_cs, set_vertical=.true.) - ! get domain extents for the analysis grid and use global indexing - !call get_domain_extent(CS%Grid%Domain,is,ie,js,je,isd,ied,jsd,jed,& - ! isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric) - !isd=isd+idg_offset; ied=ied+idg_offset - !jsd=jsd+jdg_offset; jed=jed+jdg_offset - !call mpp_get_data_domain(CS%mpp_domain,isd,ied,jsd,jed) +! call set_axes_info(CS%Grid, CS%GV, CS%US, PF, CS%diag_cs, set_vertical=.true.) ! missing in Feiyu's fork allocate(CS%oda_grid) CS%oda_grid%x => CS%Grid%geolonT CS%oda_grid%y => CS%Grid%geolatT - call get_param(PF, 'oda_driver', "BASIN_FILE", basin_file, & + + if (CS%use_basin_mask) then + call get_param(PF, 'oda_driver', "BASIN_FILE", basin_file, & "A file in which to find the basin masks, in variable 'basin'.", & default="basin.nc") - basin_file = trim(inputdir) // trim(basin_file) - allocate(CS%oda_grid%basin_mask(isd:ied,jsd:jed)) - CS%oda_grid%basin_mask(:,:) = 0.0 - call MOM_read_data(basin_file,'basin',CS%oda_grid%basin_mask,CS%Grid%domain, timelevel=1) - -! get global grid information from ocean_model - allocate(T_grid) - allocate(T_grid%x(CS%ni,CS%nj)) - allocate(T_grid%y(CS%ni,CS%nj)) - allocate(T_grid%basin_mask(CS%ni,CS%nj)) - call global_field(CS%mpp_domain, CS%Grid%geolonT, T_grid%x) - call global_field(CS%mpp_domain, CS%Grid%geolatT, T_grid%y) - call global_field(CS%mpp_domain, CS%oda_grid%basin_mask, T_grid%basin_mask) - T_grid%ni = CS%ni - T_grid%nj = CS%nj - T_grid%nk = CS%nk - allocate(T_grid%mask(CS%ni,CS%nj,CS%nk)) - allocate(T_grid%z(CS%ni,CS%nj,CS%nk)) - allocate(global2D(CS%ni,CS%nj)) - allocate(global2D_old(CS%ni,CS%nj)) - T_grid%mask(:,:,:) = 0.0 - T_grid%z(:,:,:) = 0.0 - - do k = 1, CS%nk - call global_field(G%Domain%mpp_domain, CS%h(:,:,k), global2D) - do i=1,CS%ni ; do j=1,CS%nj - if ( global2D(i,j) > 1 ) then - T_grid%mask(i,j,k) = 1.0 - endif - enddo ; enddo - if (k == 1) then - T_grid%z(:,:,k) = global2D/2 - else - T_grid%z(:,:,k) = T_grid%z(:,:,k-1) + (global2D + global2D_old)/2 - endif - global2D_old = global2D - enddo + basin_file = trim(inputdir) // trim(basin_file) + allocate(CS%oda_grid%basin_mask(isd:ied,jsd:jed)) + CS%oda_grid%basin_mask(:,:) = 0.0 + call MOM_read_data(basin_file,'basin',CS%oda_grid%basin_mask,CS%Grid%domain, timelevel=1) + endif - call ocean_da_core_init(CS%mpp_domain, T_grid, CS%Profiles, Time) + ! set up diag variables for analysis increments + CS%diag_CS => diag_CS + CS%id_inc_t=register_diag_field('ocean_model','temp_increment',diag_CS%axesTL,& + Time,'ocean potential temperature increments','degC') + CS%id_inc_s=register_diag_field('ocean_model','salt_increment',diag_CS%axesTL,& + Time,'ocean salinity increments','psu') + + !! get global grid information from ocean model needed for ODA initialization + T_grid=>NULL() + call set_up_global_tgrid(T_grid, CS, G) + call ocean_da_core_init(CS%mpp_domain, T_grid, CS%Profiles, Time) + deallocate(T_grid) CS%Time=Time !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) + + if (CS%do_bias_adjustment) then + call get_param(PF, "MOM", "TEMP_SALT_ADJUSTMENT_FILE", bias_correction_file, & + "The name of the file containing temperature and salinity "//& + "tendency adjustments", default='temp_salt_adjustment.nc') + + inc_file = trim(inputdir) // trim(bias_correction_file) + CS%INC_CS%T_id = init_extern_field(inc_file, "temp_increment", & + correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) + CS%INC_CS%S_id = init_extern_field(inc_file, "salt_increment", & + correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) + call get_external_field_info(CS%INC_CS%T_id,size=fld_sz) + CS%INC_CS%fldno = 2 + if (CS%nk .ne. fld_sz(3)) call MOM_error(FATAL,'Increment levels /= ODA levels') + allocate(CS%tv_bc) ! storage for increment + allocate(CS%tv_bc%T(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke)); CS%tv_bc%T(:,:,:)=0.0 + allocate(CS%tv_bc%S(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke)); CS%tv_bc%S(:,:,:)=0.0 + endif + + call cpu_clock_end(id_clock_oda_init) + +! if (CS%write_obs) then +! temp_fid = open_profile_file("temp_"//trim(obs_file)) +! salt_fid = open_profile_file("salt_"//trim(obs_file)) +! end if + end subroutine init_oda !> Copy ensemble member tracers to ensemble vector. @@ -340,7 +405,8 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) !! switch to global pelist call set_PElist(CS%filter_pelist) - call MOM_mesg('Setting prior') + !call MOM_mesg('Setting prior') + call cpu_clock_begin(id_clock_oda_prior) ! computational domain for the analysis grid isc=CS%Grid%isc;iec=CS%Grid%iec;jsc=CS%Grid%jsc;jec=CS%Grid%jec @@ -367,6 +433,7 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) call pass_var(CS%Ocean_prior%S(:,:,:,m),CS%Grid%domain) enddo + call cpu_clock_end(id_clock_oda_prior) !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) @@ -379,28 +446,31 @@ end subroutine set_prior_tracer subroutine get_posterior_tracer(Time, CS, h, tv, increment) type(time_type), intent(in) :: Time !< the current model time type(ODA_CS), pointer :: CS !< ocean DA control structure - real, dimension(:,:,:), pointer :: h !< Layer thicknesses [H ~> m or kg m-2] - type(thermo_var_ptrs), pointer :: tv !< A structure pointing to various thermodynamic variables + real, dimension(:,:,:), pointer, optional :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), pointer, optional :: tv !< A structure pointing to various thermodynamic variables logical, optional, intent(in) :: increment !< True if returning increment only type(ocean_control_struct), pointer :: Ocean_increment=>NULL() integer :: i, j, m logical :: used, get_inc + integer :: seconds_per_hour = 3600. ! return if not analysis time (retain pointers for h and tv) - if (Time < CS%Time) return + if (Time < CS%Time .or. CS%assim_method .eq. NO_ASSIM) return !! switch to global pelist call set_PElist(CS%filter_pelist) call MOM_mesg('Getting posterior') - + call cpu_clock_begin(id_clock_oda_posterior) + if (present(h)) h => CS%h ! get analysis thickness + !! Calculate and redistribute increments to CS%tv right after assimilation + !! Retain CS%tv to calculate increments for IAU updates CS%tv_inc otherwise get_inc = .true. if (present(increment)) get_inc = increment if (get_inc) then allocate(Ocean_increment) - call init_ocean_ensemble(Ocean_increment,CS%Grid,CS%GV,CS%ensemble_size) Ocean_increment%T = CS%Ocean_posterior%T - CS%Ocean_prior%T Ocean_increment%S = CS%Ocean_posterior%S - CS%Ocean_prior%S endif @@ -418,17 +488,28 @@ subroutine get_posterior_tracer(Time, CS, h, tv, increment) endif enddo - tv => CS%tv - h => CS%h + if (present(tv)) tv => CS%tv + if (present(h)) h => CS%h + + call cpu_clock_end(id_clock_oda_posterior) + !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) + call pass_var(CS%tv%T,CS%domains(CS%ensemble_id)) + call pass_var(CS%tv%S,CS%domains(CS%ensemble_id)) + + !convert to a tendency (degC or PSU per second) + CS%tv%T = CS%tv%T / (CS%assim_frequency * seconds_per_hour) + CS%tv%S = CS%tv%S / (CS%assim_frequency * seconds_per_hour) + + end subroutine get_posterior_tracer !> Gather observations and call ODA routines subroutine oda(Time, CS) type(time_type), intent(in) :: Time !< the current model time - type(oda_CS), intent(inout) :: CS !< the ocean DA control structure + type(oda_CS), pointer :: CS !< A pointer the ocean DA control structure integer :: i, j integer :: m @@ -438,20 +519,61 @@ subroutine oda(Time, CS) !! switch to global pelist call set_PElist(CS%filter_pelist) - + call cpu_clock_begin(id_clock_oda_filter) call get_profiles(Time, CS%Profiles, CS%CProfiles) #ifdef ENABLE_ECDA call ensemble_filter(CS%Ocean_prior, CS%Ocean_posterior, CS%CProfiles, CS%kdroot, CS%mpp_domain, CS%oda_grid) #endif - + call cpu_clock_end(id_clock_oda_filter) !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) + call get_posterior_tracer(Time, CS, increment=.true.) + if (CS%do_bias_adjustment) call get_bias_correction_tracer(Time, CS) endif return end subroutine oda +subroutine get_bias_correction_tracer(Time, CS) + type(time_type), intent(in) :: Time !< the current model time + type(ODA_CS), pointer :: CS !< ocean DA control structure + + integer :: i,j,k + real, allocatable, dimension(:,:,:) :: T_bias, S_bias + real, allocatable, dimension(:,:,:) :: mask_z + real, allocatable, dimension(:), target :: z_in, z_edges_in + real :: missing_value + integer,dimension(3) :: fld_sz + + call cpu_clock_begin(id_clock_bias_adjustment) + call horiz_interp_and_extrap_tracer(CS%INC_CS%T_id,Time,1.0,CS%G,T_bias,& + mask_z,z_in,z_edges_in,missing_value,.true.,.false.,.false.,.true.) + call horiz_interp_and_extrap_tracer(CS%INC_CS%S_id,Time,1.0,CS%G,S_bias,& + mask_z,z_in,z_edges_in,missing_value,.true.,.false.,.false.,.true.) + + ! This should be replaced to use mask_z instead of the following lines + ! which are intended to zero land values using an arbitrary limit. + fld_sz=shape(T_bias) + do i=1,fld_sz(1) + do j=1,fld_sz(2) + do k=1,fld_sz(3) + if (T_bias(i,j,k) .gt. 1.0E-3) T_bias(i,j,k) = 0.0 + if (S_bias(i,j,k) .gt. 1.0E-3) S_bias(i,j,k) = 0.0 + enddo + enddo + enddo + + CS%tv_bc%T = T_bias * CS%bias_adjustment_multiplier + CS%tv_bc%S = S_bias * CS%bias_adjustment_multiplier + + call pass_var(CS%tv_bc%T, CS%domains(CS%ensemble_id)) + call pass_var(CS%tv_bc%S, CS%domains(CS%ensemble_id)) + + call cpu_clock_end(id_clock_bias_adjustment) + + end subroutine get_bias_correction_tracer + !> Finalize DA module subroutine oda_end(CS) type(ODA_CS), intent(inout) :: CS !< the ocean DA control structure @@ -513,45 +635,127 @@ subroutine set_analysis_time(Time,CS) end subroutine set_analysis_time -!> Write observation differences to a file -subroutine save_obs_diff(filename,CS) - character(len=*), intent(in) :: filename !< name of output file - type(ODA_CS), pointer, intent(in) :: CS !< pointer to DA control structure - - integer :: fid ! profile file handle - type(ocean_profile_type), pointer :: Prof=>NULL() - - fid = open_profile_file(trim(filename), nvar=2, thread=SINGLE_FILE, fset=SINGLE_FILE) - Prof=>CS%CProfiles - - !! switch to global pelist - !call set_PElist(CS%filter_pelist) - - do while (associated(Prof)) - call write_profile(fid,Prof) - Prof=>Prof%cnext - enddo - call close_profile_file(fid) - - !! switch back to ensemble member pelist - !call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) - - return -end subroutine save_obs_diff - !> Apply increments to tracers -subroutine apply_oda_tracer_increments(dt, G, GV, tv, h, CS) +subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) real, intent(in) :: dt !< The tracer timestep [s] + type(time_type), intent(in) :: Time_end !< Time at the end of the interval type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< layer thickness [H ~> m or kg m-2] - type(ODA_CS), intent(inout) :: CS !< the data assimilation structure + type(ODA_CS), pointer :: CS !< the data assimilation structure + + !! local variables + integer :: yr, mon, day, hr, min, sec + integer :: i, j, k + integer :: isc, iec, jsc, jec + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_inc !< an adjustment to the temperature + !! tendency [degC T-1 -> degC s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_inc !< an adjustment to the salinity + !! tendency [g kg-1 T-1 -> g kg-1 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: T !< The updated temperature [degC] + real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: S !< The updated salinity [g kg-1] + real :: missing_value + + if (.not. associated(CS)) return + if (CS%assim_method .eq. NO_ASSIM .and. (.not. CS%do_bias_adjustment)) return + + call cpu_clock_begin(id_clock_apply_increments) + + T_inc(:,:,:) = 0.0; S_inc(:,:,:) = 0.0; T(:,:,:) = 0.0; S(:,:,:) = 0.0 + if (CS%assim_method > 0 ) then + T = T + CS%tv%T + S = S + CS%tv%S + endif + if (CS%do_bias_adjustment ) then + T = T + CS%tv_bc%T + S = S + CS%tv_bc%S + endif + + isc=G%isc; iec=G%iec; jsc=G%jsc; jec=G%jec + do j=jsc,jec; do i=isc,iec + call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), T(i,j,:), & + G%ke, h(i,j,:), T_inc(i,j,:)) + call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), S(i,j,:), & + G%ke, h(i,j,:), S_inc(i,j,:)) + enddo; enddo + + + call pass_var(T_inc, G%Domain) + call pass_var(S_inc, G%Domain) + + tv%T(isc:iec,jsc:jec,:)=tv%T(isc:iec,jsc:jec,:)+T_inc(isc:iec,jsc:jec,:)*dt + tv%S(isc:iec,jsc:jec,:)=tv%S(isc:iec,jsc:jec,:)+S_inc(isc:iec,jsc:jec,:)*dt + + call pass_var(tv%T, G%Domain) + call pass_var(tv%S, G%Domain) + + call enable_averaging(dt, Time_end, CS%diag_CS) + if (CS%id_inc_t > 0) call post_data(CS%id_inc_t, T_inc, CS%diag_CS) + if (CS%id_inc_s > 0) call post_data(CS%id_inc_s, S_inc, CS%diag_CS) + call disable_averaging(CS%diag_CS) + + call diag_update_remap_grids(CS%diag_CS) + call cpu_clock_end(id_clock_apply_increments) + end subroutine apply_oda_tracer_increments + subroutine set_up_global_tgrid(T_grid, CS, G) + type(grid_type), pointer :: T_grid !< global tracer grid + type(ODA_CS), pointer, intent(in) :: CS !< A pointer to DA control structure. + type(ocean_grid_type), pointer :: G !< domain and grid information for ocean model + + ! local variables + real, dimension(:,:), allocatable :: global2D, global2D_old + integer :: i, j, k + + ! get global grid information from ocean_model + T_grid=>NULL() + !if (associated(T_grid)) call MOM_error(FATAL,'MOM_oda_driver:set_up_global_tgrid called with associated T_grid') + + allocate(T_grid) + T_grid%ni = CS%ni + T_grid%nj = CS%nj + T_grid%nk = CS%nk + allocate(T_grid%x(CS%ni,CS%nj)) + allocate(T_grid%y(CS%ni,CS%nj)) + allocate(T_grid%bathyT(CS%ni,CS%nj)) + call global_field(CS%mpp_domain, CS%Grid%geolonT, T_grid%x) + call global_field(CS%mpp_domain, CS%Grid%geolatT, T_grid%y) + call global_field(CS%domains(CS%ensemble_id)%mpp_domain, G%bathyT, T_grid%bathyT) + if (CS%use_basin_mask) then + allocate(T_grid%basin_mask(CS%ni,CS%nj)) + call global_field(CS%mpp_domain, CS%oda_grid%basin_mask, T_grid%basin_mask) + endif + allocate(T_grid%mask(CS%ni,CS%nj,CS%nk)) + allocate(T_grid%z(CS%ni,CS%nj,CS%nk)) + allocate(global2D(CS%ni,CS%nj)) + allocate(global2D_old(CS%ni,CS%nj)) + T_grid%mask(:,:,:) = 0.0 + T_grid%z(:,:,:) = 0.0 + + do k = 1, CS%nk + call global_field(G%Domain%mpp_domain, CS%h(:,:,k), global2D) + do i=1,CS%ni ; do j=1,CS%nj + if ( global2D(i,j) > 1 ) then + T_grid%mask(i,j,k) = 1.0 + endif + enddo; enddo + if (k == 1) then + T_grid%z(:,:,k) = global2D/2 + else + T_grid%z(:,:,k) = T_grid%z(:,:,k-1) + (global2D + global2D_old)/2 + endif + global2D_old = global2D + enddo + + deallocate(global2D) + deallocate(global2D_old) + end subroutine set_up_global_tgrid + !> \namespace MOM_oda_driver_mod !! !! \section section_ODA The Ocean data assimilation (DA) and Ensemble Framework diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 8779968bcd..6497890040 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -82,7 +82,9 @@ module MOM_MEKE real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. - + logical :: fixed_total_depth !< If true, use the nominal bathymetric depth as the estimate of + !! the time-varying ocean depth. Otherwise base the depth on the total + !! ocean mass per unit area. logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled. logical :: initialize !< If True, invokes a steady state solver to calculate MEKE. logical :: debug !< If true, write out checksums of data for debugging @@ -131,8 +133,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h MEKE_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1]. drag_rate_visc, & ! Near-bottom velocity contribution to bottom dratg [L T-1 ~> m s-1] drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1]. - drag_rate_J15, & ! The MEKE spindown timescale due to bottom drag with the Jansen 2015 scheme. - ! Unfortunately, as written the units seem inconsistent. [T-1 ~> s-1]. del2MEKE, & ! Laplacian of MEKE, used for bi-harmonic diffusion [T-2 ~> s-2]. del4MEKE, & ! Time-integrated MEKE tendency arising from the biharmonic of MEKE [L2 T-2 ~> m2 s-2]. LmixScale, & ! Eddy mixing length [L ~> m]. @@ -177,8 +177,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h if (.not.associated(MEKE)) call MOM_error(FATAL, & "MOM_MEKE: MEKE must be initialized before it is used.") - if ((CS%MEKE_damping > 0.0) .or. (CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) & - .or. CS%visc_drag) then + if ((CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) .or. CS%visc_drag) then use_drag_rate = .true. else use_drag_rate = .false. @@ -234,12 +233,6 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo endif - if (CS%MEKE_Cd_scale == 0.0 .and. .not. CS%visc_drag) then - !$OMP parallel do default(shared) private(ldamping) - do j=js,je ; do i=is,ie - drag_rate(i,j) = 0. ; drag_rate_J15(i,j) = 0. - enddo ; enddo - endif ! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow if (CS%visc_drag) then @@ -283,12 +276,17 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo enddo - !$OMP parallel do default(shared) - do j=js-1,je+1 ; do i=is-1,ie+1 - depth_tot(i,j) = G%bathyT(i,j) - !### Try changing this to: - ! depth_tot(i,j) = mass(i,j) * I_Rho0 - enddo ; enddo + if (CS%fixed_total_depth) then + !$OMP parallel do default(shared) + do j=js-1,je+1 ; do i=is-1,ie+1 + depth_tot(i,j) = G%bathyT(i,j) + enddo ; enddo + else + !$OMP parallel do default(shared) + do j=js-1,je+1 ; do i=is-1,ie+1 + depth_tot(i,j) = mass(i,j) * I_Rho0 + enddo ; enddo + endif if (CS%initialize) then call MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot) @@ -363,6 +361,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h drag_rate(i,j) = (US%L_to_Z*Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) enddo ; enddo + else + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + drag_rate(i,j) = 0. + enddo ; enddo endif ! First stage of Strang splitting @@ -531,23 +534,19 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h drag_rate(i,j) = (US%L_to_Z*Rho0 * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + & cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) ) enddo ; enddo - !$OMP parallel do default(shared) - do j=js,je ; do i=is,ie - ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) - if (MEKE%MEKE(i,j) < 0.) ldamping = 0. - ! notice that the above line ensures a damping only if MEKE is positive, - ! while leaving MEKE unchanged if it is negative - MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) - MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) - enddo ; enddo endif + !$OMP parallel do default(shared) + do j=js,je ; do i=is,ie + ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j) + if (MEKE%MEKE(i,j) < 0.) ldamping = 0. + ! notice that the above line ensures a damping only if MEKE is positive, + ! while leaving MEKE unchanged if it is negative + MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping) + MEKE_decay(i,j) = ldamping*G%mask2dT(i,j) + enddo ; enddo endif endif ! MEKE_KH>=0 - ! do j=js,je ; do i=is,ie - ! MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0) - ! enddo ; enddo - call cpu_clock_begin(CS%id_clock_pass) call do_group_pass(CS%pass_MEKE, G%Domain) call cpu_clock_end(CS%id_clock_pass) @@ -711,8 +710,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then beta_topo_x = 0. ; beta_topo_y = 0. else - !### Consider different combinations of these estimates of topographic beta, and the use - ! of the water column thickness instead of the bathymetric depth. + !### Consider different combinations of these estimates of topographic beta. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) & / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) & @@ -887,14 +885,13 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, EKE, depth_tot, & FatH = 0.25* ( ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1) ) + & ( G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1) ) ) ! Coriolis parameter at h points - ! If bathyT is zero, then a division by zero FPE will be raised. In this + ! If depth_tot is zero, then a division by zero FPE will be raised. In this ! case, we apply Adcroft's rule of reciprocals and set the term to zero. ! Since zero-bathymetry cells are masked, this should not affect values. if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then beta_topo_x = 0. ; beta_topo_y = 0. else - !### Consider different combinations of these estimates of topographic beta, and the use - ! of the water column thickness instead of the bathymetric depth. + !### Consider different combinations of these estimates of topographic beta. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) & / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) & @@ -1167,6 +1164,11 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) "If positive, is a fixed length contribution to the expression "//& "for mixing length used in MEKE-derived diffusivity.", & units="m", default=0.0, scale=US%m_to_L) + call get_param(param_file, mdl, "MEKE_FIXED_TOTAL_DEPTH", CS%fixed_total_depth, & + "If true, use the nominal bathymetric depth as the estimate of the "//& + "time-varying ocean depth. Otherwise base the depth on the total ocean mass"//& + "per unit area.", default=.true.) + call get_param(param_file, mdl, "MEKE_ALPHA_DEFORM", CS%aDeform, & "If positive, is a coefficient weighting the deformation scale "//& "in the expression for mixing length used in MEKE-derived diffusivity.", & diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index b4f857dec4..84eabc9317 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -190,6 +190,7 @@ module MOM_hor_visc integer :: id_h_diffu = -1, id_h_diffv = -1 integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1 integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1 + integer :: id_diffu_visc_rem = -1, id_diffv_visc_rem = -1 integer :: id_Ah_h = -1, id_Ah_q = -1 integer :: id_Kh_h = -1, id_Kh_q = -1 integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1 @@ -276,10 +277,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_vel_mag_h, & ! Magnitude of the velocity gradient tensor squared at h-points [T-2 ~> s-2] grad_vel_mag_bt_h, & ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2] grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2] + GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim] + htot, & ! The total thickness of all layers [Z ~> m] boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] real, allocatable, dimension(:,:,:) :: h_diffu ! h x diffu [H L T-2 ~> m2 s-2] real, allocatable, dimension(:,:,:) :: h_diffv ! h x diffv [H L T-2 ~> m2 s-2] + real, allocatable, dimension(:,:,:) :: diffu_visc_rem ! diffu x visc_rem_u [L T-2 ~> m s-2] + real, allocatable, dimension(:,:,:) :: diffv_visc_rem ! diffv x visc_rem_v [L T-2 ~> m s-2] real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1] @@ -301,6 +306,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] ! This form guarantees that hq/hu < 4. grad_vel_mag_bt_q, & ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] + GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim] boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: & @@ -338,6 +344,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity ! points where masks are applied [H ~> m or kg m-2]. + real :: h_arith_q ! The arithmetic mean total thickness at q points [Z ~> m] + real :: h_harm_q ! The harmonic mean total thickness at q points [Z ~> m] + real :: I_hq ! The inverse of the arithmetic mean total thickness at q points [Z-1 ~> m-1] + real :: I_GME_h0 ! The inverse of GME tapering scale [Z-1 ~> m-1] real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2] real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6] real :: h_min ! Minimum h at the 4 neighboring velocity points [H ~> m] @@ -501,6 +511,35 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2) enddo ; enddo + do j=js-1,je+1 ; do i=is-1,ie+1 + htot(i,j) = 0.0 + enddo ; enddo + do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + htot(i,j) = htot(i,j) + GV%H_to_Z*h(i,j,k) + enddo ; enddo ; enddo + + I_GME_h0 = 1.0 / CS%GME_h0 + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (grad_vel_mag_bt_h(i,j)>0) then + GME_effic_h(i,j) = CS%GME_efficiency * boundary_mask_h(i,j) * & + (MIN(htot(i,j) * I_GME_h0, 1.0)**2) + else + GME_effic_h(i,j) = 0.0 + endif + enddo ; enddo + + do J=js-1,Jeq ; do I=is-1,Ieq + if (grad_vel_mag_bt_q(I,J)>0) then + h_arith_q = 0.25 * ((htot(i,j) + htot(i+1,j+1)) + (htot(i+1,j) + htot(i,j+1))) + I_hq = 1.0 / h_arith_q + h_harm_q = 0.25 * h_arith_q * ((htot(i,j)*I_hq + htot(i+1,j+1)*I_hq) + & + (htot(i+1,j)*I_hq + htot(i,j+1)*I_hq)) + GME_effic_q(I,J) = CS%GME_efficiency * boundary_mask_q(I,J) * (MIN(h_harm_q * I_GME_h0, 1.0)**2) + else + GME_effic_q(I,J) = 0.0 + endif + enddo ; enddo + endif ! use_GME !$OMP parallel do default(none) & @@ -509,7 +548,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, & !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & !$OMP use_MEKE_Ku, use_MEKE_Au, boundary_mask_h, boundary_mask_q, & - !$OMP backscat_subround, GME_coeff_limiter, & + !$OMP backscat_subround, GME_coeff_limiter, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, FWfrac, inv_PI3, inv_PI6, H0_GME, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, & @@ -1388,15 +1427,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call pass_vector(KH_u_GME, KH_v_GME, G%Domain) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if (grad_vel_mag_bt_h(i,j)>0) then - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & - (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)+KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) - else - GME_coeff = 0.0 - endif - - ! apply mask - GME_coeff = GME_coeff * boundary_mask_h(i,j) + GME_coeff = GME_effic_h(i,j) * 0.25 * & + ((KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) GME_coeff = MIN(GME_coeff, CS%GME_limiter) if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff @@ -1405,17 +1437,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq - if (grad_vel_mag_bt_q(I,J)>0) then - !### This expression is not rotationally invariant - bathyT is to the SW of q points, - ! and it needs parentheses in the sum of the 4 diffusivities. - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & - (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)+KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) - else - GME_coeff = 0.0 - endif - - ! apply mask - GME_coeff = GME_coeff * boundary_mask_q(I,J) + GME_coeff = GME_effic_q(I,J) * 0.25 * & + ((KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) GME_coeff = MIN(GME_coeff, CS%GME_limiter) if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff @@ -1424,8 +1447,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo ! Applying GME diagonal term. This is linear and the arguments can be rescaled. - call smooth_GME(CS,G,GME_flux_h=str_xx_GME) - call smooth_GME(CS,G,GME_flux_q=str_xy_GME) + call smooth_GME(CS, G, GME_flux_h=str_xx_GME) + call smooth_GME(CS, G, GME_flux_q=str_xy_GME) do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = (str_xx(i,j) + str_xx_GME(i,j)) * (h(i,j,k) * CS%reduction_xx(i,j)) @@ -1446,7 +1469,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo endif - else ! use_GME + else ! .not. use_GME do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo @@ -1706,6 +1729,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, deallocate(h_diffv) endif + if (present(ADp) .and. (CS%id_diffu_visc_rem > 0)) then + allocate(diffu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) + diffu_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + diffu_visc_rem(I,j,k) = diffu(I,j,k) * ADp%visc_rem_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_diffu_visc_rem, diffu_visc_rem, CS%diag) + deallocate(diffu_visc_rem) + endif + if (present(ADp) .and. (CS%id_diffv_visc_rem > 0)) then + allocate(diffv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke)) + diffv_visc_rem(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + diffv_visc_rem(i,J,k) = diffv(i,J,k) * ADp%visc_rem_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_diffv_visc_rem, diffv_visc_rem, CS%diag) + deallocate(diffv_visc_rem) + endif + end subroutine horizontal_viscosity !> Allocates space for and calculates static variables used by horizontal_viscosity(). @@ -2450,6 +2492,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) endif + CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, & + 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if ((CS%id_diffu_visc_rem > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%visc_rem_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke) + endif + + CS%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, Time, & + 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if ((CS%id_diffv_visc_rem > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%visc_rem_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) + endif + if (CS%biharmonic) then CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time, & 'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, & diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 4d66471408..0452a83a23 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -96,6 +96,9 @@ module MOM_internal_tides real :: decay_rate !< A constant rate at which internal tide energy is !! lost to the interior ocean internal wave field. real :: cdrag !< The bottom drag coefficient [nondim]. + real :: drag_min_depth !< The minimum total ocean thickness that will be used in the denominator + !! of the quadratic drag terms for internal tides when + !! INTERNAL_TIDE_QUAD_DRAG is true [Z ~> m] logical :: apply_background_drag !< If true, apply a drag due to background processes as a sink. logical :: apply_bottom_drag @@ -185,6 +188,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & tot_En, & ! energy summed over angles, modes, frequencies [R Z3 T-2 ~> J m-2] tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_Froude_loss, tot_allprocesses_loss, & ! energy loss rates summed over angle, freq, and mode [R Z3 T-3 ~> W m-2] + htot, & ! The vertical sum of the layer thicknesses [H ~> m or kg m-2] drag_scale, & ! bottom drag scale [T-1 ~> s-1] itidal_loss_mode, allprocesses_loss_mode ! energy loss rates for a given mode and frequency (summed over angles) [R Z3 T-3 ~> W m-2] @@ -200,7 +204,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & real :: En_initial, Delta_E_check ! Energies for debugging [R Z3 T-2 ~> J m-2] real :: TKE_Froude_loss_check, TKE_Froude_loss_tot ! Energy losses for debugging [R Z3 T-3 ~> W m-2] character(len=160) :: mesg ! The text of an error message - integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm + integer :: a, m, fr, i, j, k, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging) type(group_pass_type), save :: pass_test, pass_En type(time_type) :: time_end @@ -212,6 +216,10 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & I_rho0 = 1.0 / GV%Rho0 cn_subRO = 1e-100*US%m_s_to_L_T ! The hard-coded value here might need to increase. + ! init local arrays + drag_scale(:,:) = 0. + Ub(:,:,:,:) = 0. + ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.********************** ! This is wrong, of course, but it works reasonably in some cases. ! Uncomment if wave_speed is not used to calculate the true values (BDM). @@ -360,9 +368,12 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! Extract the energy for mixing due to bottom drag------------------------------- if (CS%apply_bottom_drag) then + do j=jsd,jed ; do i=isd,ied ; htot(i,j) = 0.0 ; enddo ; enddo + do k=1,GV%ke ; do j=jsd,jed ; do i=isd,ied + htot(i,j) = htot(i,j) + h(i,j,k) + enddo ; enddo ; enddo do j=jsd,jed ; do i=isd,ied - ! Note the 1 m dimensional scale here. Should this be a parameter? - I_D_here = 1.0 / (max(G%bathyT(i,j), 1.0*US%m_to_Z)) + I_D_here = 1.0 / (max(GV%H_to_Z*htot(i,j), CS%drag_min_depth)) drag_scale(i,j) = CS%cdrag * sqrt(max(0.0, US%L_to_Z**2*vel_btTide(i,j)**2 + & tot_En(i,j) * I_rho0 * I_D_here)) * I_D_here enddo ; enddo @@ -2143,20 +2154,20 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) ! Local variables real :: Angle_size ! size of wedges, rad real, allocatable :: angles(:) ! orientations of wedge centers, rad - real, allocatable, dimension(:,:) :: h2 ! topographic roughness scale, m^2 - real :: kappa_itides, kappa_h2_factor - ! characteristic topographic wave number - ! and a scaling factor - real, allocatable :: ridge_temp(:,:) - ! array for temporary storage of flags + real, dimension(:,:), allocatable :: h2 ! topographic roughness scale squared [Z2 ~> m2] + real :: kappa_itides ! characteristic topographic wave number [L-1 ~> m-1] + real, dimension(:,:), allocatable :: ridge_temp ! array for temporary storage of flags ! of cells with double-reflecting ridges - logical :: use_int_tides, use_temperature - real :: period_1 ! The period of the gravest modeled mode [T ~> s] + logical :: use_int_tides, use_temperature + real :: kappa_h2_factor ! A roughness scaling factor [nondim] + real :: RMS_roughness_frac ! The maximum RMS topographic roughness as a fraction of the + ! nominal ocean depth, or a negative value for no limit [nondim] + real :: period_1 ! The period of the gravest modeled mode [T ~> s] integer :: num_angle, num_freq, num_mode, m, fr integer :: isd, ied, jsd, jed, a, id_ang, i, j type(axes_grp) :: axes_ang ! This include declares and sets the variable "version". -#include "version_variable.h" +# include "version_variable.h" character(len=40) :: mdl = "MOM_internal_tides" ! This module's name. character(len=16), dimension(8) :: freq_name character(len=40) :: var_name @@ -2206,7 +2217,10 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) ! Allocate and populate frequency array (each a multiple of first for now) allocate(CS%frequency(num_freq)) - call get_param(param_file, mdl, "FIRST_MODE_PERIOD", period_1, units="s", scale=US%s_to_T) + call get_param(param_file, mdl, "FIRST_MODE_PERIOD", period_1, & + "The period of the first mode for internal tides", default=44567., & + units="s", scale=US%s_to_T) + do fr=1,num_freq CS%frequency(fr) = (8.0*atan(1.0) * (real(fr)) / period_1) ! ADDED BDM enddo @@ -2280,16 +2294,20 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "1st-order upwind advection. This scheme is highly "//& "continuity solver. This scheme is highly "//& "diffusive but may be useful for debugging.", default=.false.) - call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", & - CS%apply_background_drag, "If true, the internal tide "//& - "ray-tracing advection uses a background drag term as a sink.",& - default=.false.) + call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", CS%apply_background_drag, & + "If true, the internal tide ray-tracing advection uses a background drag "//& + "term as a sink.", default=.false.) call get_param(param_file, mdl, "INTERNAL_TIDE_QUAD_DRAG", CS%apply_bottom_drag, & "If true, the internal tide ray-tracing advection uses "//& "a quadratic bottom drag term as a sink.", default=.false.) call get_param(param_file, mdl, "INTERNAL_TIDE_WAVE_DRAG", CS%apply_wave_drag, & "If true, apply scattering due to small-scale roughness as a sink.", & default=.false.) + call get_param(param_file, mdl, "INTERNAL_TIDE_DRAG_MIN_DEPTH", CS%drag_min_depth, & + "The minimum total ocean thickness that will be used in the denominator "//& + "of the quadratic drag terms for internal tides.", & + units="m", default=1.0, scale=US%m_to_Z, do_not_log=.not.CS%apply_bottom_drag) + CS%drag_min_depth = MAX(CS%drag_min_depth, GV%H_subroundoff * GV%H_to_Z) call get_param(param_file, mdl, "INTERNAL_TIDE_FROUDE_DRAG", CS%apply_Froude_drag, & "If true, apply wave breaking as a sink.", & default=.false.) @@ -2313,7 +2331,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "The default is 2pi/10 km, as in St.Laurent et al. 2002.", & units="m-1", default=8.e-4*atan(1.0), scale=US%L_to_m) call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, & - "A scaling factor for the roughness amplitude with n"//& + "A scaling factor for the roughness amplitude with "//& "INT_TIDE_DISSIPATION.", units="nondim", default=1.0) ! Allocate various arrays needed for loss rates @@ -2340,10 +2358,18 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) fail_if_missing=.true.) filename = trim(CS%inputdir) // trim(h2_file) call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename) - call MOM_read_data(filename, 'h2', h2, G%domain, scale=US%m_to_Z) + call get_param(param_file, mdl, "INTERNAL_TIDE_ROUGHNESS_FRAC", RMS_roughness_frac, & + "The maximum RMS topographic roughness as a fraction of the nominal ocean depth, "//& + "or a negative value for no limit.", units="nondim", default=0.1) + + call MOM_read_data(filename, 'h2', h2, G%domain, scale=US%m_to_Z**2) do j=G%jsc,G%jec ; do i=G%isc,G%iec - ! Restrict rms topo to 10 percent of column depth. - h2(i,j) = min(0.01*(G%bathyT(i,j))**2, h2(i,j)) + ! Restrict RMS topographic roughness to a fraction (10 percent by default) of the column depth. + if (RMS_roughness_frac >= 0.0) then + h2(i,j) = max(min((RMS_roughness_frac*G%bathyT(i,j))**2, h2(i,j)), 0.0) + else + h2(i,j) = max(h2(i,j), 0.0) + endif ! Compute the fixed part; units are [R L-2 Z3 ~> kg m-2] here ! will be multiplied by N and the squared near-bottom velocity to get into [R Z3 T-3 ~> W m-2] CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*GV%Rho0 * US%L_to_Z*kappa_itides * h2(i,j) @@ -2513,8 +2539,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) do a=1,num_angle ; angles(a) = (real(a) - 1) * Angle_size ; enddo id_ang = diag_axis_init("angle", angles, "Radians", "N", "Angular Orienation of Fluxes") - call define_axes_group(diag, (/ diag%axesT1%handles(1), diag%axesT1%handles(2), id_ang /), axes_ang) - + call define_axes_group(diag, (/ diag%axesT1%handles(1), diag%axesT1%handles(2), id_ang /), & + axes_ang, is_h_point=.true.) do fr=1,CS%nFreq ; write(freq_name(fr), '("freq",i1)') fr ; enddo do m=1,CS%nMode ; do fr=1,CS%nFreq ! Register 2-D energy density (summed over angles) for each freq and mode diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 3b3d72576c..0c8f25b42e 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -145,6 +145,8 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp KH_u_CFL ! The maximum stable interface height diffusivity at u grid points [L2 T-1 ~> m2 s-1] real, dimension(SZI_(G), SZJB_(G)) :: & KH_v_CFL ! The maximum stable interface height diffusivity at v grid points [L2 T-1 ~> m2 s-1] + real, dimension(SZI_(G), SZJ_(G)) :: & + htot ! The sum of the total layer thicknesses [H ~> m or kg m-2] real :: Khth_Loc_u(SZIB_(G), SZJ_(G)) real :: Khth_Loc_v(SZI_(G), SZJB_(G)) real :: Khth_Loc(SZIB_(G), SZJB_(G)) ! locally calculated thickness diffusivity [L2 T-1 ~> m2 s-1] @@ -479,38 +481,41 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp ! thicknesses across u and v faces, converted to 0/1 mask ! layer average of the interface diffusivities KH_u and KH_v do j=js,je ; do I=is-1,ie - hu(I,j) = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect) - if (hu(I,j) /= 0.0) hu(I,j) = 1.0 - !### The same result would be accomplished with the following without a division: - ! hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0 + ! This expression uses harmonic mean thicknesses: + ! hu(I,j) = 2.0*h(i,j,k)*h(i+1,j,k) / (h(i,j,k)+h(i+1,j,k)+h_neglect) + ! This expression is a 0/1 mask based on depths where there are thick layers: + hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0 KH_u_lay(I,j) = 0.5*(KH_u(I,j,k)+KH_u(I,j,k+1)) enddo ; enddo do J=js-1,je ; do i=is,ie - hv(i,J) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect) - if (hv(i,J) /= 0.0) hv(i,J) = 1.0 - !### The same result would be accomplished with the following without a division: - ! hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0 + ! This expression uses harmonic mean thicknesses: + ! hv(i,J) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect) + ! This expression is a 0/1 mask based on depths where there are thick layers: + hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0 KH_v_lay(i,J) = 0.5*(KH_v(i,J,k)+KH_v(i,J,k+1)) enddo ; enddo - ! diagnose diffusivity at T-point - !### Because hu and hv are nondimensional here, the denominator is dimensionally inconsistent. + ! diagnose diffusivity at T-points do j=js,je ; do i=is,ie - Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j)+hu(I,j)*KH_u_lay(I,j)) & - +(hv(i,J-1)*KH_v_lay(i,J-1)+hv(i,J)*KH_v_lay(i,J))) & - / (hu(I-1,j)+hu(I,j)+hv(i,J-1)+hv(i,J)+h_neglect) + Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j) + hu(I,j)*KH_u_lay(I,j)) + & + (hv(i,J-1)*KH_v_lay(i,J-1) + hv(i,J)*KH_v_lay(i,J))) / & + ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + 1.0e-20) + ! Use this denominator instead if hu and hv are actual thicknesses rather than a 0/1 mask: + ! ((hu(I-1,j)+hu(I,j)) + (hv(i,J-1)+hv(i,J)) + h_neglect) enddo ; enddo enddo if (CS%Use_KH_in_MEKE) then MEKE%Kh_diff(:,:) = 0.0 + htot(:,:) = 0.0 do k=1,nz do j=js,je ; do i=is,ie MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) + Kh_t(i,j,k) * h(i,j,k) + htot(i,j) = htot(i,j) + h(i,j,k) enddo ; enddo enddo do j=js,je ; do i=is,ie - MEKE%Kh_diff(i,j) = GV%H_to_Z * MEKE%Kh_diff(i,j) / MAX(1.0*US%m_to_Z, G%bathyT(i,j)) + MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0*GV%m_to_H, htot(i,j)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 9eeb8867db..604ac03626 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -218,7 +218,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) call get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", CS%hor_regrid_answers_2018, & - "If true, use the order of arithmetic for horizonal regridding that recovers "//& + "If true, use the order of arithmetic for horizontal regridding that recovers "//& "the answers from the end of 2018. Otherwise, use rotationally symmetric "//& "forms of the same expressions.", default=default_2018_answers) call get_param(param_file, mdl, "REENTRANT_X", CS%reentrant_x, & @@ -478,7 +478,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) call get_param(param_file, mdl, "HOR_REGRID_2018_ANSWERS", CS%hor_regrid_answers_2018, & - "If true, use the order of arithmetic for horizonal regridding that recovers "//& + "If true, use the order of arithmetic for horizontal regridding that recovers "//& "the answers from the end of 2018 and retain a bug in the 3-dimensional mask "//& "returned in certain cases. Otherwise, use rotationally symmetric "//& "forms of the same expressions and initialize the mask properly.", & diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index ee27c6c5df..072bc1445e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -45,7 +45,7 @@ module MOM_diabatic_aux real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m]. logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to !! to cool the topmost layer down to the freezing - !! point. The default is false. + !! point. The default is true. logical :: pressure_dependent_frazil !< If true, use a pressure dependent !! freezing temperature when making frazil. The !! default is false, which will be faster but is diff --git a/src/parameterizations/vertical/_Frazil.dox b/src/parameterizations/vertical/_Frazil.dox new file mode 100644 index 0000000000..06321231a3 --- /dev/null +++ b/src/parameterizations/vertical/_Frazil.dox @@ -0,0 +1,33 @@ +/*! \page Frazil_Ice Frazil Ice Formation + +\section section_frazil Frazil Ice Formation + +Frazil ice forms in the model when the in situ temperature drops below +the local freezing point, taking into account the in situ salinity and +pressure. Starting at the bottom and working up through the water column, +if the water is below freezing, set it to freezing and add the heat +required to the heat deficit. If the water above is warmer than freezing, +use that heat to take away the heat deficit and to cool the water. If +you get all the way to the surface with a heat deficit, that quantity +is passed to the ice model as a heat flux it will need to provide to +the ocean. + +The local freezing point code is provided by the equation of state being +used by MOM6. See \ref section_TFREEZE for the MOM6 options. + +The salinity is adjusted only at the surface when frazil ice is +formed. This happens when the ice model creates ice with the heat deficit, +taking salt out of the surface waters. We inherit this behavior from +older versions of MOM, but the effect of not adjusting the in situ +salinity is thought to be small. + +Note that versions simply whisking all the heat deficit to the surface +without checking for warm water above tended to produce rapidly-melting +ice floes in warm waters. This was deemed unphysical and was corrected. + +A similar process that we are also omitting is the formation of salt +crystals when the salinity becomes too high. The salt crystals should +form and sink, leaving a layer on the bed that will be diluted when the +salinity drops again. This process can be seen in a lake in Death Valley. + +*/ diff --git a/src/parameterizations/vertical/_Internal_tides.dox b/src/parameterizations/vertical/_Internal_tides.dox new file mode 100644 index 0000000000..882b73dd1b --- /dev/null +++ b/src/parameterizations/vertical/_Internal_tides.dox @@ -0,0 +1,216 @@ +/*! \page Internal_Tidal_Mixing Internal Tidal Mixing + +Two parameterizations of vertical mixing due to internal tides are +available with the option INT_TIDE_DISSIPATION. The first is that of +\cite st_laurent2002 while the second is that of \cite polzin2009. Choose +between them with the INT_TIDE_PROFILE option. There are other relevant +paramters which can be seen in MOM_parameter_doc.all once the main tidal +dissipation switch is turned on. + +\section section_st_laurent St Laurent et al. + +The estimated turbulent dissipation rate of +internal tide energy \f$\epsilon\f$ is: + +\f[ + \epsilon = \frac{q E(x,y)}{\rho} F(z). +\f] + +where \f$\rho\f$ is the reference density of seawater, \f$E(x,y)\f$ is +the energy flux per unit area transferred from barotropic to baroclinic +tides, \f$q\f$ is the fraction of the internal-tide energy dissipated +locally, and \f$F(z)\f$ is the vertical structure of the dissipation. +This \f$q\f$ is estimated to be roughly 0.3 based on observations. The +term \f$E(x,y)\f$ is given by \cite st_laurent2002 as: + +\f[ + E(x,y) \simeq \frac{1}{2} \rho N_b \kappa h^2 \langle U^2 \rangle +\f] + +where \f$\rho\f$ is the reference density of seawater, \f$N_b\f$ is +the buoyancy frequency along the seafloor, and \f$(\kappa, h)\f$ are +the wavenumber and amplitude scales for the topographic roughness, and +\f$\langle U^2 \rangle\f$ is the barotropic tide variance. It is assumed +that the model will read in topographic roughness squared \f$h^2\f$ +from a file (the variable must be named "h2"). + +To convert from energy dissipation to vertical diffusion \f$K_d\f$, +the simple estimate is: + +\f[ + K_d \approx \frac{\Gamma q E(x,y) F(z)}{\rho N^2} +\f] + +where \f$\Gamma\f$ is the mixing efficiency, generally set to 0.2 +and \f$F(z)\f$ is a vertical structure function with exponential decay +from the bottom: + +\f[ + F(z) = \frac{e^{-(H+z)/\zeta}}{\zeta (1 - e^{H/\zeta}}. +\f] + +Here, \f$\zeta\f$ is a vertical decay scale with a default of 500 meters. +One change in MOM6 from the St. Laurent scheme is to use this form of \f$\Gamma\f$: + +\f[ + \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} +\f] + +instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity +of the Earth. This allows the buoyancy fluxes to tend to zero in regions +of very weak stratification, allowing a no-flux bottom boundary condition +to be satisfied. + +This \f$K_d\f$ gets added to the diffusivity due to the background and +other contributions unless you set BBL_MIXING_AS_MAX to True, in which +case the maximum of all the contributions is used. + +\section section_polzin Polzin + +The vertical diffusion profile of \cite polzin2009 is a WKB-stretched +algebraic decay profile. It is based on a radiation balance equation, +which links the dissipation profile associtated with internal breaking to +the finescale internal wave shear producing that dissipation. The vertical +profile of internal-tide driven energy dissipation can then vary in time +and space, and evolve in a changing climate (\cite melet2012). \cite melet2012 +describes how the Polzin scheme is implemented in MOM6, +copied here. + +The parameterization of \cite polzin2009 links the energy dissipation +profile to the finescale internal wave shear producing that +dissipation, using an idealized vertical wavenumber energy spectrum +to identify analytic solutions to a radiation balance equation +(\cite polzin2004). These solutions yield a dissipation profile +\f$\epsilon(z)\f$: + +\f[ + \epsilon = \frac{\epsilon_0}{[1 + (z/z_p)]^2}, +\f] + +where the magnitude \f$\epsilon_0\f$ and scale height \f$z_p\f$ can be expressed in terms of the +spectral amplitude and bandwidth of the idealized vertical wavenumber energy spectrum in uniform +stratification (\cite polzin2009). + +To take into account the nonuniform stratification, \cite polzin2009 applied a buoyancy scaling +using the Wentzel-Kramers-Brillouin (WKB) approximation. As a result, the vertical wavenumber of a +wave packet varies in proportion to the buoyancy frequency \f$N\f$, which in turn implies an +additional transport of energy to smaller scales, and thus a possible enhanced mixing in regions of +strong stratification. Such effects can be described by buoyancy scaling the vertical coordinate +\f$z\f$ as + +\f[ + z^{\ast}(z) = \int_{0}^{z} \left[ \frac{N^2 (z^\prime )}{N_b^2} \right] dz^{\prime} , +\f] + +with \f$z^\prime\f$ being positive upward relative to the bottom of the ocean. The turbulent +dissipation rate then becomes + +\f[ + \epsilon = \frac{\epsilon_0}{[1 + (z^{\ast} /z_p)]^2} \frac{N^2(z)}{N_b^2} . +\f] + +The spectral amplitude and bandwidth of the idealized vertical wavenumber +energy spectrum are identified after WKB scaling using a quasi-linear +spectral model of internal-tide generation that incorporates horizontal +advection of the barotropic tide into the momentum equation (\cite bell1975). +As a result, Polzin's formulation leads to an expression for +the spatially and temporally varying dissipation of internal tide energy +at the bottom \f$\epsilon_0\f$, and the vertical scale of decay for the +dissipation of internal tide energy \f$z_p\f$. + +\subsection subsection_energy_conserving Energy-conserving form + +To satisfy energy conservation (the integral of the vertical structure for the turbulent dissipation +over depth should be unity), the dissipation is rewritten as + +\f[ + \epsilon = \frac{\epsilon_0 z_p}{1 + (z^\ast/z_p)]^2} \frac{N^2(z)}{N^2_b} \left[ + \frac{1}{z^{\ast(z=H)}} + \frac{1}{z_p} \right] . +\f] + +In the MOM6 implementation, we use the \cite st_laurent2002 template for the vertical flux of energy +at the ocean floor, so that in both formulations: + +\f[ + \int_{0}^{H} \epsilon (z) dz = \frac{qE}{\rho} . +\f] + +Whereas \cite polzin2009 assumed tthat the total dissipation was locally in balance with the +barotropic to baroclinic energy conversion rate \f$(q=1)\f$, here we use the \cite simmons2004 value +of \f$q=1/3\f$ to retain as much consistency as passible between both parameterizations. + +\subsection subsection_vertical_decay_scale Vertical decay-scale reformulation + +We follow the \cite polzin2009 prescription for the vertical scale of +decay for the dissipation of internal-tide energy. However, we assume +that the topographic power law, denoted as \f$\nu\f$ in \cite polzin2009, +is equal to 1 (instead of 0.9) and we reformulated the expression of +\f$z_p\f$ to put it in a more readable form: + +\f[ + z_p = \frac{z_p^\mbox{ref} (\kappa^\mbox{ref})^2 (h^\mbox{ref})^2 (N_b^\mbox{ref})^3} {U^\mbox{ref}} + \frac{U}{h^2 \kappa^2 N_b^3} . +\f] + +The superscript ref refers to reference values of the various parameters, as given by +observations from the Brazil basin. Therefore, the above can be rewritten as + +\f[ + z_p = \mu (N_b^\mbox{ref} )^2 + \frac{U}{h^2 \kappa^2 N_b^3} . +\f] + +where \f$\mu\f$ is a nondimensional constant \f$(\mu = 0.06970)\f$ and \f$N_b^\mbox{ref} = 9.6 \times +10^{-4} s^{-1}\f$. Finally, a minimum decay scale of \f$z_p = 100 m\f$ is imposed in our +implementation. + +\subsection subsection_reformulation_WKB Reformulation of the WKB scaling + +Since the dissipation is expressed as a function of the ratio \f$z^\ast / z_p\f$, a different WKB +scaling can be used so long as we modify \f$z_p\f$ accordingly. In the implemented parameterization, +we define the scaled height coordinate \f$z^\ast\f$ by + +\f[ + z^\ast (z) = \frac{1}{\overline{N^2 (z)}^z} \int_{0}^{z} N^2(z^\prime ) dz ^\prime , +\f] + +with \f$z^\prime\f$ defined to be the height above the ocean bottom. By normalizing \f$N^2\f$ by its +vertical mean \f$\overline{N^2}^z\f$, \f$z^\ast\f$ ranges from \f$0\f$ to \f$H\f$, the depth of the +ocean. + +The WKB-scaled vertical decay scale for the Polzin formulation becomes + +\f[ + z^\ast_p = \mu(N_b^\mbox{ref})^2 \frac{U}{h^2 \kappa^2 N_b \overline{N^2}^z} . +\f] + +Unlike the \cite st_laurent2002 parameterization, the vertical decay scale now depends on physical +variables and can evolve with a changing climate. + +Finally, the Polzin vertical profile of dissipation implemented in the model is given by + +\f[ + \epsilon = \frac{qE(x,y)}{\rho [1 + (z^\ast/z_p^\ast)]^2} \frac{N^2(z)}{\overline{N^2}^z} + \left( \frac{1}{H} + \frac{1}{z_p^\ast} \right) . +\f] + +In both parameterizations, turbulent diapycnal diffusivities are inferred from the dissipation +\f$\epsilon\f$ by: + +\f[ + K_d = \frac{\Gamma \epsilon}{N^2} +\f] + +and using this form of \f$\Gamma\f$: + +\f[ + \Gamma = 0.2 \frac{N^2}{N^2 + \Omega^2} +\f] + +instead of \f$\Gamma = 0.2\f$, where \f$\Omega\f$ is the angular velocity +of the Earth. This allows the buoyancy fluxes to tend to zero in regions +of very weak stratification, allowing a no-flux bottom boundary condition +to be satisfied. + +*/ + diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index c20eda7745..44421c7387 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -75,7 +75,7 @@ function register_DOME_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) character(len=48) :: flux_units ! The units for tracer fluxes, usually ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1. character(len=200) :: inputdir - real, pointer :: tr_ptr(:,:,:) => NULL() + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to the tracer field logical :: register_DOME_tracer integer :: isd, ied, jsd, jed, nz, m isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke @@ -166,13 +166,15 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & character(len=48) :: units ! The dimensions of the variable. character(len=48) :: flux_units ! The units for tracer fluxes, usually ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1. - real, pointer :: tr_ptr(:,:,:) => NULL() + real, pointer :: tr_ptr(:,:,:) => NULL() ! A pointer to the tracer field real :: PI ! 3.1415926... calculated as 4*atan(1) real :: tr_y ! Initial zonally uniform tracer concentrations. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. - real :: e(SZK_(GV)+1), e_top, e_bot ! Heights [Z ~> m]. - real :: d_tr ! A change in tracer concentraions, in tracer units. + real :: e(SZK_(GV)+1) ! Interface heights relative to the sea surface (negative down) [Z ~> m] + real :: e_top ! Height of the top of the tracer band relative to the sea surface [Z ~> m] + real :: e_bot ! Height of the bottom of the tracer band relative to the sea surface [Z ~> m] + real :: d_tr ! A change in tracer concentrations, in tracer units. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m integer :: IsdB, IedB, JsdB, JedB @@ -215,9 +217,9 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, US, h, diag, OBC, CS, & if (NTR > 7) then do j=js,je ; do i=is,ie - e(nz+1) = -G%bathyT(i,j) - do k=nz,1,-1 - e(K) = e(K+1) + h(i,j,k)*GV%H_to_Z + e(1) = 0.0 + do k=1,nz + e(K+1) = e(K) - h(i,j,k)*GV%H_to_Z do m=7,NTR e_top = (-600.0*real(m-1) + 3000.0) * US%m_to_Z e_bot = (-600.0*real(m-1) + 2700.0) * US%m_to_Z diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 26ef197ae2..60ff9c981c 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -62,6 +62,8 @@ module MOM_tracer_flow_control use boundary_impulse_tracer, only : boundary_impulse_tracer_column_physics, boundary_impulse_tracer_surface_state use boundary_impulse_tracer, only : boundary_impulse_stock, boundary_impulse_tracer_end use boundary_impulse_tracer, only : boundary_impulse_tracer_CS +use nw2_tracers, only : nw2_tracers_CS, register_nw2_tracers, nw2_tracer_column_physics +use nw2_tracers, only : initialize_nw2_tracers, nw2_tracers_end implicit none ; private @@ -84,6 +86,7 @@ module MOM_tracer_flow_control logical :: use_pseudo_salt_tracer = .false. !< If true, use the psuedo_salt tracer package logical :: use_boundary_impulse_tracer = .false. !< If true, use the boundary impulse tracer package logical :: use_dyed_obc_tracer = .false. !< If true, use the dyed OBC tracer package + logical :: use_nw2_tracers = .false. !< If true, use the ideal age tracer package !>@{ Pointers to the control strucures for the tracer packages type(USER_tracer_example_CS), pointer :: USER_tracer_example_CSp => NULL() type(DOME_tracer_CS), pointer :: DOME_tracer_CSp => NULL() @@ -98,6 +101,7 @@ module MOM_tracer_flow_control type(pseudo_salt_tracer_CS), pointer :: pseudo_salt_tracer_CSp => NULL() type(boundary_impulse_tracer_CS), pointer :: boundary_impulse_tracer_CSp => NULL() type(dyed_obc_tracer_CS), pointer :: dyed_obc_tracer_CSp => NULL() + type(nw2_tracers_CS), pointer :: nw2_tracers_CSp => NULL() !>@} end type tracer_flow_control_CS @@ -206,6 +210,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) call get_param(param_file, mdl, "USE_DYED_OBC_TRACER", CS%use_dyed_obc_tracer, & "If true, use the dyed_obc_tracer tracer package.", & default=.false.) + call get_param(param_file, mdl, "USE_NW2_TRACERS", CS%use_nw2_tracers, & + "If true, use the NeverWorld2 tracers.", & + default=.false.) ! Add other user-provided calls to register tracers for restarting here. Each ! tracer package registration call returns a logical false if it cannot be run @@ -249,7 +256,8 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) if (CS%use_dyed_obc_tracer) CS%use_dyed_obc_tracer = & register_dyed_obc_tracer(HI, GV, param_file, CS%dyed_obc_tracer_CSp, & tr_Reg, restart_CS) - + if (CS%use_nw2_tracers) CS%use_ideal_age = & + register_nw2_tracers(HI, GV, param_file, CS%nw2_tracers_CSp, tr_Reg, restart_CS) end subroutine call_tracer_register @@ -328,6 +336,8 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag sponge_CSp, tv) if (CS%use_dyed_obc_tracer) & call initialize_dyed_obc_tracer(restart, day, G, GV, h, diag, OBC, CS%dyed_obc_tracer_CSp) + if (CS%use_nw2_tracers) & + call initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS%nw2_tracers_CSp) end subroutine tracer_flow_control_init @@ -486,8 +496,11 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, G, GV, US, CS%dyed_obc_tracer_CSp, & evap_CFL_limit=evap_CFL_limit, & minimum_forcing_depth=minimum_forcing_depth) - - + if (CS%use_nw2_tracers) & + call nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, tv, CS%nw2_tracers_CSp, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) else ! Apply tracer surface fluxes using ea on the first layer if (CS%use_USER_tracer_example) & call tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & @@ -532,10 +545,10 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, if (CS%use_dyed_obc_tracer) & call dyed_obc_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, CS%dyed_obc_tracer_CSp) - + if (CS%use_nw2_tracers) call nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, tv, CS%nw2_tracers_CSp) endif - end subroutine call_tracer_column_fns !> This subroutine calls all registered tracer packages to enable them to @@ -776,6 +789,7 @@ subroutine tracer_flow_control_end(CS) if (CS%use_pseudo_salt_tracer) call pseudo_salt_tracer_end(CS%pseudo_salt_tracer_CSp) if (CS%use_boundary_impulse_tracer) call boundary_impulse_tracer_end(CS%boundary_impulse_tracer_CSp) if (CS%use_dyed_obc_tracer) call dyed_obc_tracer_end(CS%dyed_obc_tracer_CSp) + if (CS%use_nw2_tracers) call nw2_tracers_end(CS%nw2_tracers_CSp) if (associated(CS)) deallocate(CS) end subroutine tracer_flow_control_end diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index 2919f2d95f..2bf3cd94ed 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -203,9 +203,10 @@ subroutine initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS, sponge_C character(len=48) :: units ! The dimensions of the variable. character(len=48) :: flux_units ! The units for age tracer fluxes, either ! years m3 s-1 or years kg s-1. + real :: z_bot ! Height of the bottom of the layer relative to the sea surface [Z ~> m] + real :: z_center ! Height of the center of the layer relative to the sea surface [Z ~> m] logical :: OK integer :: i, j, k, m - real :: z_bot, z_center if (.not.associated(CS)) return if (CS%ntr < 1) return @@ -221,14 +222,14 @@ subroutine initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS, sponge_C CS%dye_source_minlat(m)=G%geoLatT(i,j) .and. & G%mask2dT(i,j) > 0.0 ) then - z_bot = -G%bathyT(i,j) - do k = GV%ke, 1, -1 + z_bot = 0.0 + do k = 1, GV%ke + z_bot = z_bot - h(i,j,k)*GV%H_to_Z z_center = z_bot + 0.5*h(i,j,k)*GV%H_to_Z if ( z_center > -CS%dye_source_maxdepth(m) .and. & z_center < -CS%dye_source_mindepth(m) ) then CS%tr(i,j,k,m) = 1.0 endif - z_bot = z_bot + h(i,j,k)*GV%H_to_Z enddo endif enddo ; enddo @@ -273,9 +274,10 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US real :: sfc_val ! The surface value for the tracers. real :: Isecs_per_year ! The number of seconds in a year. real :: year ! The time in years. + real :: z_bot ! Height of the bottom of the layer relative to the sea surface [Z ~> m] + real :: z_center ! Height of the center of the layer relative to the sea surface [Z ~> m] integer :: secs, days ! Integer components of the time type. integer :: i, j, k, is, ie, js, je, nz, m - real :: z_bot, z_center is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -305,14 +307,14 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US CS%dye_source_minlat(m)=G%geoLatT(i,j) .and. & G%mask2dT(i,j) > 0.0 ) then - z_bot = -G%bathyT(i,j) - do k=nz,1,-1 + z_bot = 0.0 + do k=1,nz + z_bot = z_bot - h_new(i,j,k)*GV%H_to_Z z_center = z_bot + 0.5*h_new(i,j,k)*GV%H_to_Z if ( z_center > -CS%dye_source_maxdepth(m) .and. & z_center < -CS%dye_source_mindepth(m) ) then CS%tr(i,j,k,m) = 1.0 endif - z_bot = z_bot + h_new(i,j,k)*GV%H_to_Z enddo endif enddo ; enddo diff --git a/src/tracer/nw2_tracers.F90 b/src/tracer/nw2_tracers.F90 new file mode 100644 index 0000000000..e75c5c5d38 --- /dev/null +++ b/src/tracer/nw2_tracers.F90 @@ -0,0 +1,314 @@ +!> Ideal tracers designed to help diagnose a tracer diffusivity tensor in NeverWorld2 +module nw2_tracers + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing +use MOM_grid, only : ocean_grid_type +use MOM_hor_index, only : hor_index_type +use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc +use MOM_restart, only : query_initialized, MOM_restart_CS +use MOM_time_manager, only : time_type, time_type_to_real +use MOM_tracer_registry, only : register_tracer, tracer_registry_type +use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface, thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +public register_nw2_tracers +public initialize_nw2_tracers +public nw2_tracer_column_physics +public nw2_tracers_end + +!> The control structure for the nw2_tracers package +type, public :: nw2_tracers_CS ; private + integer :: ntr = 0 !< The number of tracers that are actually used. + type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. + type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the tracer registry + real, pointer :: tr(:,:,:,:) => NULL() !< The array of tracers used in this package, in g m-3? + real, allocatable , dimension(:) :: restore_rate !< The exponential growth rate for restoration value [year-1]. + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart controls structure +end type nw2_tracers_CS + +contains + +!> Register the NW2 tracer fields to be used with MOM. +logical function register_nw2_tracers(HI, GV, param_file, CS, tr_Reg, restart_CS) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracer. + type(tracer_registry_type), pointer :: tr_Reg !< A pointer that is set to point to the control + !! structure for the tracer advection and + !! diffusion module + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure + +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mdl = "nw2_tracers" ! This module's name. + character(len=200) :: inputdir ! The directory where the input files are. + character(len=8) :: var_name ! The variable's name. + real, pointer :: tr_ptr(:,:,:) => NULL() + logical :: do_nw2 + integer :: isd, ied, jsd, jed, nz, m, ig + integer :: n_groups ! Number of groups of three tracers (i.e. # tracers/3) + real, allocatable, dimension(:) :: timescale_in_days + type(vardesc) :: tr_desc ! Descriptions and metadata for the tracers + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + + if (associated(CS)) then + call MOM_error(FATAL, "register_nw2_tracer called with an "// & + "associated control structure.") + return + endif + allocate(CS) + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "NW2_TRACER_GROUPS", n_groups, & + "The number of tracer groups where a group is of three tracers "//& + "initialized and restored to sin(x), y and z, respectively. Each "//& + "group is restored with an independent restoration rate.", & + default=3) + allocate(timescale_in_days(n_groups)) + timescale_in_days = (/365., 730., 1460./) + call get_param(param_file, mdl, "NW2_TRACER_RESTORE_TIMESCALE", timescale_in_days, & + "A list of timescales, one for each tracer group.", & + units="days") + + CS%ntr = 3 * n_groups + allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0 + allocate(CS%restore_rate(CS%ntr)) + + do m=1,CS%ntr + write(var_name(1:8),'(a6,i2.2)') 'tracer',m + tr_desc = var_desc(var_name, "1", "Ideal Tracer", caller=mdl) + ! This is needed to force the compiler not to do a copy in the registration + ! calls. Curses on the designers and implementers of Fortran90. + tr_ptr => CS%tr(:,:,:,m) + ! Register the tracer for horizontal advection, diffusion, and restarts. + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=tr_desc, & + registry_diags=.true., restart_CS=restart_CS, mandatory=.false.) + ig = int( (m+2)/3 ) ! maps (1,2,3)->1, (4,5,6)->2, ... + CS%restore_rate(m) = 1.0 / ( timescale_in_days(ig) * 86400.0 ) + enddo + + CS%tr_Reg => tr_Reg + CS%restart_CSp => restart_CS + register_nw2_tracers = .true. +end function register_nw2_tracers + +!> Sets the NW2 traces to their initial values and sets up the tracer output +subroutine initialize_nw2_tracers(restart, day, G, GV, US, h, tv, diag, CS) + logical, intent(in) :: restart !< .true. if the fields have already + !! been read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(in) :: 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 + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables + type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate + !! diagnostic output. + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracer. + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights + real :: rscl ! z* scaling factor + character(len=8) :: var_name ! The variable's name. + integer :: i, j, k, m + + if (.not.associated(CS)) return + + CS%Time => day + CS%diag => diag + + ! Calculate z* interface positions + if (GV%Boussinesq) then + ! First calculate interface positions in z-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,GV%ke+1) = - G%mask2dT(i,j) * G%bathyT(i,j) + enddo ; enddo + do k=GV%ke,1,-1 ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h(i,j,k) * GV%H_to_Z + enddo ; enddo ; enddo + ! Re-calculate for interface positions in z*-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (G%bathyT(i,j)>0.) then + rscl = G%bathyT(i,j) / ( eta(i,j,1) + G%bathyT(i,j) ) + do K=GV%ke, 1, -1 + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h(i,j,k) * GV%H_to_Z * rscl + enddo + endif + enddo ; enddo + else + call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode") + endif + + do m=1,CS%ntr + ! Initialize only if this is not a restart or we are using a restart + ! in which the tracers were not present + write(var_name(1:8),'(a6,i2.2)') 'tracer',m + if ((.not.restart) .or. & + (.not. query_initialized(CS%tr(:,:,:,m),var_name,CS%restart_CSp))) then + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + CS%tr(i,j,k,m) = nw2_tracer_dist(m, G, GV, eta, i, j, k) + enddo ; enddo ; enddo + endif ! restart + enddo ! Tracer loop + +end subroutine initialize_nw2_tracers + +!> Applies diapycnal diffusion, aging and regeneration at the surface to the NW2 tracers +subroutine nw2_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, tv, CS, & + evap_CFL_limit, minimum_forcing_depth) + 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(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: ea !< an array to which the amount of fluid entrained + !! from the layer above during this call will be + !! added [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: eb !< an array to which the amount of fluid entrained + !! from the layer below during this call will be + !! added [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic + !! and tracer forcing fields. Unused fields have NULL ptrs. + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracer. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which + !! fluxes can be applied [H ~> m or kg m-2] +! This subroutine applies diapycnal diffusion and any other column +! tracer physics or chemistry to the tracers from this file. +! This is a simple example of a set of advected passive tracers. + +! The arguments to this subroutine are redundant in that +! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta ! Interface heights + integer :: i, j, k, m + real :: dt_x_rate ! dt * restoring rate + real :: rscl ! z* scaling factor + real :: target_value ! tracer value + +! if (.not.associated(CS)) return + + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + do m=1,CS%ntr + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CS%tr(:,:,:,m), dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV) + enddo + else + do m=1,CS%ntr + call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV) + enddo + endif + + ! Calculate z* interface positions + if (GV%Boussinesq) then + ! First calculate interface positions in z-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,GV%ke+1) = - G%mask2dT(i,j) * G%bathyT(i,j) + enddo ; enddo + do k=GV%ke,1,-1 ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h_new(i,j,k) * GV%H_to_Z + enddo ; enddo ; enddo + ! Re-calculate for interface positions in z*-space (m) + do j=G%jsc,G%jec ; do i=G%isc,G%iec + if (G%bathyT(i,j)>0.) then + rscl = G%bathyT(i,j) / ( eta(i,j,1) + G%bathyT(i,j) ) + do K=GV%ke, 1, -1 + eta(i,j,K) = eta(i,j,K+1) + G%mask2dT(i,j) * h_new(i,j,k) * GV%H_to_Z * rscl + enddo + endif + enddo ; enddo + else + call MOM_error(FATAL, "NW2 tracers assume Boussinesq mode") + endif + + do m=1,CS%ntr + dt_x_rate = ( dt * CS%restore_rate(m) ) * US%T_to_s +!$OMP parallel do default(private) shared(CS,G,dt,dt_x_rate) + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + target_value = nw2_tracer_dist(m, G, GV, eta, i, j, k) + CS%tr(i,j,k,m) = CS%tr(i,j,k,m) + G%mask2dT(i,j) * dt_x_rate * ( target_value - CS%tr(i,j,k,m) ) + enddo ; enddo ; enddo + enddo + +end subroutine nw2_tracer_column_physics + +!> The target value of a NeverWorld2 tracer label m at non-dimensional +!! position x=lon/Lx, y=lat/Ly, z=eta/H +real function nw2_tracer_dist(m, G, GV, eta, i, j, k) + integer, intent(in) :: m !< Indicates the NW2 tracer + 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(SZI_(G),SZJ_(G),0:SZK_(G)), & + intent(in) :: eta !< Interface position [m] + integer, intent(in) :: i !< Cell index i + integer, intent(in) :: j !< Cell index j + integer, intent(in) :: k !< Layer index k + ! Local variables + real :: pi ! 3.1415... + real :: x, y, z ! non-dimensional positions + pi = 2.*acos(0.) + x = ( G%geolonT(i,j) - G%west_lon ) / G%len_lon ! 0 ... 1 + y = -G%geolatT(i,j) / G%south_lat ! -1 ... 1 + z = - 0.5 * ( eta(i,j,K-1) + eta(i,j,K) ) / GV%max_depth ! 0 ... 1 + select case ( mod(m-1,3) ) + case (0) ! sin(2 pi x/L) + nw2_tracer_dist = sin( 2.0 * pi * x ) + case (1) ! y/L + nw2_tracer_dist = y + case (2) ! -z/L + nw2_tracer_dist = -z + case default + stop 'This should not happen. Died in nw2_tracer_dist()!' + end select + nw2_tracer_dist = nw2_tracer_dist * G%mask2dT(i,j) +end function nw2_tracer_dist + +!> Deallocate any memory associated with this tracer package +subroutine nw2_tracers_end(CS) + type(nw2_tracers_CS), pointer :: CS !< The control structure returned by a previous + !! call to register_nw2_tracers. + + integer :: m + + if (associated(CS)) then + if (associated(CS%tr)) deallocate(CS%tr) + deallocate(CS) + endif +end subroutine nw2_tracers_end + +!> \namespace nw2_tracers +!! +!! TBD + +end module nw2_tracers