diff --git a/.github/workflows/macos-regression.yml b/.github/workflows/macos-regression.yml index 422c50b68a..d769e15131 100644 --- a/.github/workflows/macos-regression.yml +++ b/.github/workflows/macos-regression.yml @@ -10,6 +10,7 @@ jobs: env: CC: gcc FC: gfortran + FMS_COMMIT: 2019.01.03 defaults: run: diff --git a/.github/workflows/macos-stencil.yml b/.github/workflows/macos-stencil.yml index 36a5841bb2..6e77a5c4a6 100644 --- a/.github/workflows/macos-stencil.yml +++ b/.github/workflows/macos-stencil.yml @@ -10,6 +10,7 @@ jobs: env: CC: gcc FC: gfortran + FMS_COMMIT: 2019.01.03 defaults: run: diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 2c4d92c424..55494696ae 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -32,6 +32,8 @@ p:clone: tags: - ncrc5 script: + # NOTE: We could sweep any builds older than 3 days here if needed + #- find $HOME/ci/[0-9]* -mtime +3 -delete 2> /dev/null || true - .gitlab/pipeline-ci-tool.sh create-job-dir #.gitlab/pipeline-ci-tool.sh clean-job-dir @@ -353,4 +355,5 @@ cleanup: before_script: - echo Skipping usual preamble script: + - rm -rf $HOME/ci/$CI_PIPELINE_ID - rm -rf $JOB_DIR diff --git a/.gitlab/pipeline-ci-tool.sh b/.gitlab/pipeline-ci-tool.sh index 77409d29ef..d948b72008 100755 --- a/.gitlab/pipeline-ci-tool.sh +++ b/.gitlab/pipeline-ci-tool.sh @@ -94,6 +94,13 @@ create-job-dir () { make -f tools/MRS/Makefile.clone clone_gfdl -j # Extras and link to datasets bash tools/MRS/generate_manifest.sh . tools/MRS/excluded-expts.txt > manifest.mk mkdir -p results + # Temporarily move build directory to $HOME to circumvent poor F5 performance + mkdir -p $HOME/ci/$CI_PIPELINE_ID/build + ln -s $HOME/ci/$CI_PIPELINE_ID/build build + # Builds need non-mangled access to src/. + ln -s "$(pwd)"/src $HOME/ci/$CI_PIPELINE_ID/src + # Static builds need access to ocean_only/ + ln -s "$(pwd)"/ocean_only $HOME/ci/$CI_PIPELINE_ID/ocean_only fi section-end create-job-dir } diff --git a/.testing/Makefile b/.testing/Makefile index a1461fe7ab..f5da44342d 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -20,7 +20,6 @@ # # General test configuration: # MPIRUN MPI job launcher (mpirun, srun, etc) -# FRAMEWORK Model framework (fms1 or fms2) # DO_REPRO_TESTS Enable production ("repro") testing equivalence # DO_REGRESSION_TESTS Enable regression tests (usually dev/gfdl) # DO_COVERAGE Enable code coverage and generate .gcov reports @@ -74,8 +73,11 @@ AC_SRCDIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST))))../ac # User-defined configuration -include config.mk -# Set the infra framework -FRAMEWORK ?= fms2 +# Set the FMS library +FMS_COMMIT ?= 2023.03 +FMS_URL ?= https://github.com/NOAA-GFDL/FMS.git +export FMS_COMMIT +export FMS_URL # Set the MPI launcher here # TODO: This needs more automated configuration @@ -133,9 +135,6 @@ TIME ?= time WORKSPACE ?= . # Set directories for build/ and work/ -#BUILD ?= $(WORKSPACE)build -#DEPS ?= $(BUILD)/deps -#WORK ?= $(WORKSPACE)work BUILD ?= $(WORKSPACE)/build DEPS ?= $(BUILD)/deps WORK ?= $(WORKSPACE)/work @@ -207,34 +206,6 @@ else endif -# List of source files to link this Makefile's dependencies to model Makefiles -# Assumes a depth of two, and the following extensions: F90 inc c h -# (1): Root directory -# NOTE: extensions could be a second variable -SOURCE = \ - $(foreach ext,F90 inc c h,$(wildcard $(1)/*/*.$(ext) $(1)/*/*/*.$(ext))) - -MOM_SOURCE = \ - $(call SOURCE,../src) \ - $(wildcard ../config_src/drivers/solo_driver/*.F90) \ - $(wildcard ../config_src/ext*/*/*.F90) - -TARGET_SOURCE = \ - $(call SOURCE,$(BUILD)/target_codebase/src) \ - $(wildcard $(BUILD)/target_codebase/config_src/drivers/solo_driver/*.F90) \ - $(wildcard $(BUILD)target_codebase/config_src/ext*/*.F90) - -ifeq ($(FRAMEWORK), fms1) - MOM_SOURCE += $(wildcard ../config_src/infra/FMS1/*.F90) - TARGET_SOURCE += $(wildcard $(BUILD)/target_codebase/config_src/infra/FMS1/*.F90) -else - MOM_SOURCE +=$(wildcard ../config_src/infra/FMS2/*.F90) - TARGET_SOURCE += $(wildcard $(BUILD)/target_codebase/config_src/infra/FMS2/*.F90) -endif - -FMS_SOURCE = $(call SOURCE,$(DEPS)/fms/src) - - ## Rules .PHONY: all build.regressions build.prof @@ -286,7 +257,6 @@ $(BUILD)/unit/Makefile: MOM_ENV += $(COV_FCFLAGS) $(COV_LDFLAGS) $(BUILD)/timing/Makefile: MOM_ENV += $(OPT_FCFLAGS) $(MOM_LDFLAGS) # Configure script flags -MOM_ACFLAGS := --with-framework=$(FRAMEWORK) $(BUILD)/openmp/Makefile: MOM_ACFLAGS += --enable-openmp $(BUILD)/coupled/Makefile: MOM_ACFLAGS += --with-driver=FMS_cap $(BUILD)/nuopc/Makefile: MOM_ACFLAGS += --with-driver=nuopc_cap @@ -298,11 +268,21 @@ $(BUILD)/timing/Makefile: MOM_ACFLAGS += --with-driver=timing_tests $(BUILD)/unit/test_%: $(BUILD)/unit/Makefile FORCE cd $(@D) && $(TIME) $(MAKE) $(@F) -j $(BUILD)/unit/Makefile: $(foreach e,$(UNIT_EXECS),../config_src/drivers/unit_tests/$(e).F90) + $(BUILD)/timing/time_%: $(BUILD)/timing/Makefile FORCE cd $(@D) && $(TIME) $(MAKE) $(@F) -j $(BUILD)/timing/Makefile: $(foreach e,$(TIMING_EXECS),../config_src/drivers/timing_tests/$(e).F90) + $(BUILD)/%/MOM6: $(BUILD)/%/Makefile FORCE cd $(@D) && $(TIME) $(MAKE) $(@F) -j + +# Target codebase should use its own build system +$(BUILD)/target/MOM6: $(BUILD)/target FORCE | $(TARGET_CODEBASE) + $(MAKE) -C $(TARGET_CODEBASE)/.testing build/symmetric/MOM6 + +$(BUILD)/target: | $(TARGET_CODEBASE) + ln -s $(abspath $(TARGET_CODEBASE))/.testing/build/symmetric $@ + FORCE: @@ -334,27 +314,12 @@ $(BUILD)/%/configure.ac: ../ac/configure.ac | $(BUILD)/%/ $(BUILD)/%/m4/: ../ac/m4/ | $(BUILD)/%/ cp -r ../ac/m4 $(@D) -ALL_EXECS = symmetric asymmetric repro openmp target opt opt_target coupled \ - nuopc cov unit timing +ALL_EXECS = symmetric asymmetric repro openmp opt opt_target coupled nuopc \ + cov unit timing $(foreach b,$(ALL_EXECS),$(BUILD)/$(b)/): mkdir -p $@ # Fetch the regression target codebase - -$(BUILD)/target/config.status: $(BUILD)/target/configure $(DEPS)/lib/libFMS.a - cd $(@D) && $(MOM_ENV) ./configure -n \ - --srcdir=$(abspath $(BUILD))/target_codebase/ac $(MOM_ACFLAGS) \ - || (cat config.log && false) - -$(BUILD)/target/Makefile.in: | $(TARGET_CODEBASE) $(BUILD)/target/ - cp $(TARGET_CODEBASE)/ac/Makefile.in $(@D) - -$(BUILD)/target/configure.ac: | $(TARGET_CODEBASE) $(BUILD)/target/ - cp $(TARGET_CODEBASE)/ac/configure.ac $(@D) - -$(BUILD)/target/m4/: | $(TARGET_CODEBASE) $(BUILD)/target/ - cp -r $(TARGET_CODEBASE)/ac/m4 $(@D) - $(TARGET_CODEBASE): git clone --recursive $(MOM_TARGET_URL) $@ cd $@ && git checkout --recurse-submodules $(MOM_TARGET_BRANCH) diff --git a/.testing/README.rst b/.testing/README.rst index 49103da718..a84eeea80e 100644 --- a/.testing/README.rst +++ b/.testing/README.rst @@ -47,9 +47,11 @@ Several of the following may require configuration for particular systems. Name of the MPI launcher. Often this is ``mpirun`` or ``mpiexec`` but may all need to run through a scheduler, e.g. ``srun`` if using Slurm. -``FRAMEWORK`` (*default:* ``fms1``) - Select either the legacy FMS framework (``fms1``) or an FMS2 I/O compatible - version (``fms2``). +``FMS_COMMIT`` (*default:* ``2023.03``) + Set the FMS version, either by tag or commit (as defined in ``FMS_URL``). + +``FMS_URL`` (*default*: ``https://github.com/NOAA-GFDL/FMS.git``) + Set the URL of the FMS repository. ``DO_REPRO_TESTS`` (*default:* *none*) Set to ``true`` to test the REPRO build and confirm equivalence of DEBUG and diff --git a/.testing/tools/parse_perf.py b/.testing/tools/parse_perf.py index 7cbffd995d..76c6be5bcb 100755 --- a/.testing/tools/parse_perf.py +++ b/.testing/tools/parse_perf.py @@ -3,10 +3,20 @@ import collections import json import os +import re import shlex import subprocess import sys +perf_scanner = re.Scanner([ + (r'<', lambda scanner, token: token), + (r'>', lambda scanner, token: token), + (r'\(', lambda scanner, token: token), + (r'\)', lambda scanner, token: token), + (r'[ \t]+', lambda scanner, token: token), + (r'[^<>() \t]+', lambda scanner, token: token), +]) + def main(): desc = 'Parse perf.data and return in JSON format.' @@ -58,15 +68,55 @@ def parse_perf_report(perf_data_path): # get per-symbol count else: + tokens, remainder = perf_scanner.scan(line) + if remainder: + print('Line could not be tokenized', file=sys.stderr) + print(' line:', repr(line), file=sys.stderr) + print(' tokens:', tokens, file=sys.stderr) + print(' remainder:', remainder, file=sys.stderr) + sys.exit(os.EX_DATAERR) + + # Construct record from tokens + # (NOTE: Not a proper grammar, just dumb bracket counting) + record = [] + bracks = 0 + parens = 0 + + for tok in tokens: + if tok == '<': + bracks += 1 + + if tok == '(': + parens += 1 + + rec = record[-1] if record else None + + inside_bracket = rec and (bracks > 0 or parens > 0) + lead_rec = tok in '<(' and rec and not rec.isspace() + tail_rec = not tok.isspace() and rec and rec[-1] in '>)' + + if inside_bracket or lead_rec or tail_rec: + record[-1] += tok + else: + record.append(tok) + + if tok == '>': + bracks -= 1 + if tok == '(': + parens -= 1 + + # Strip any whitespace tokens + record = [rec for rec in record if not rec.isspace()] + try: - tokens = line.split() - symbol = tokens[2] - period = int(tokens[3]) - except ValueError: + symbol = record[2] + period = int(record[3]) + except: print("parse_perf.py: Error extracting symbol count", - file=sys.stderr) + file=sys.stderr) print("line:", repr(line), file=sys.stderr) print("tokens:", tokens, file=sys.stderr) + print("record:", record, file=sys.stderr) raise profile[event_name]['symbol'][symbol] = period diff --git a/ac/configure.ac b/ac/configure.ac index c774c39129..8196e2eb01 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -79,18 +79,6 @@ AS_IF([test "x$with_driver" != "x"], # used to configure a header based on a template. #AC_CONFIG_HEADERS(["$MEM_LAYOUT/MOM_memory.h"]) -# Select the model framework (default: FMS1) -# NOTE: We can phase this out after the FMS1 I/O has been removed from FMS and -# replace with a detection test. For now, it is a user-defined switch. -MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS2 -AC_ARG_WITH([framework], - AS_HELP_STRING([--with-framework=fms1|fms2], [Select the model 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/FMS2] -) - # Explicitly assume free-form Fortran AC_LANG(Fortran) @@ -220,7 +208,6 @@ AX_FC_CHECK_LIB([FMS], [fms_init], [fms_mod], ] ) - # Verify that FMS is at least 2019.01.02 # NOTE: 2019.01.02 introduced two changes: # - diag_axis_init supports an optional domain_position argument @@ -236,6 +223,14 @@ AC_COMPILE_IFELSE( ] ) +# Determine the FMS IO implementation. +AX_FC_CHECK_MODULE([fms2_io_mod], [ + MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS2 +],[ + MODEL_FRAMEWORK=${srcdir}/config_src/infra/FMS1 +]) + + # Python interpreter test # Declare the Python interpreter variable diff --git a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 index c4be8c769d..787fa7e7d1 100644 --- a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 @@ -26,6 +26,7 @@ program Shelf_main use MOM_debugging, only : MOM_debugging_init use MOM_diag_mediator, only : diag_mediator_init, diag_mediator_infrastructure_init, set_axes_info use MOM_diag_mediator, only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration + use MOM_diag_manager_infra, only : diag_manager_set_time_end_infra use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_domains, only : MOM_domains_init, clone_MOM_domain, pass_var use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid @@ -325,6 +326,8 @@ program Shelf_main Time_end = daymax endif + call diag_manager_set_time_end_infra (Time_end) + if (Time >= Time_end) call MOM_error(FATAL, & "Shelf_driver: The run has been started at or after the end time of the run.") diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 3574943918..3e3abba674 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -136,6 +136,7 @@ module MOM_cap_mod logical :: grid_attach_area = .false. logical :: use_coldstart = .true. logical :: use_mommesh = .true. +logical :: restart_eor = .false. character(len=128) :: scalar_field_name = '' integer :: scalar_field_count = 0 integer :: scalar_field_idx_grid_nx = 0 @@ -381,6 +382,13 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) geomtype = ESMF_GEOMTYPE_GRID endif + ! Read end of run restart config option + call NUOPC_CompAttributeGet(gcomp, name="write_restart_at_endofrun", value=value, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + if (trim(value) .eq. '.true.') restart_eor = .true. + end if end subroutine @@ -1637,6 +1645,8 @@ subroutine ModelAdvance(gcomp, rc) real(8) :: MPI_Wtime, timers logical :: write_restart logical :: write_restartfh + logical :: write_restart_eor + rc = ESMF_SUCCESS if(profile_memory) call ESMF_VMLogMemInfo("Entering MOM Model_ADVANCE: ") @@ -1776,7 +1786,6 @@ subroutine ModelAdvance(gcomp, rc) !--------------- ! Get the stop alarm !--------------- - call ESMF_ClockGetAlarm(clock, alarmname='stop_alarm', alarm=stop_alarm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -1807,7 +1816,18 @@ subroutine ModelAdvance(gcomp, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if - if (write_restart .or. write_restartfh) then + write_restart_eor = .false. + if (restart_eor) then + if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then + if (ChkErr(rc,__LINE__,u_FILE_u)) return + write_restart_eor = .true. + ! turn off the alarm + call ESMF_AlarmRingerOff(stop_alarm, rc=rc ) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end if + end if + + if (write_restart .or. write_restartfh .or. write_restart_eor) then ! determine restart filename call ESMF_ClockGetNextTime(clock, MyTime, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/config_src/drivers/solo_driver/MOM_driver.F90 b/config_src/drivers/solo_driver/MOM_driver.F90 index 0e355f8638..9b85fafb8d 100644 --- a/config_src/drivers/solo_driver/MOM_driver.F90 +++ b/config_src/drivers/solo_driver/MOM_driver.F90 @@ -28,6 +28,7 @@ program MOM6 use MOM_cpu_clock, only : CLOCK_COMPONENT use MOM_data_override, only : data_override_init use MOM_diag_mediator, only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration + use MOM_diag_manager_infra, only : diag_manager_set_time_end_infra use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only : extract_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized @@ -375,6 +376,8 @@ program MOM6 Time_end = daymax endif + call diag_manager_set_time_end_infra(Time_end) + call get_param(param_file, mod_name, "SINGLE_STEPPING_CALL", single_step_call, & "If true, advance the state of MOM with a single step "//& "including both dynamics and thermodynamics. If false "//& diff --git a/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 b/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 index b7ee7de684..5d78e0d501 100644 --- a/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 +++ b/config_src/external/GFDL_ocean_BGC/FMS_coupler_util.F90 @@ -22,6 +22,8 @@ subroutine extract_coupler_values(BC_struc, BC_index, BC_element, array_out, ilb integer, optional, intent(in) :: js !< The j- limits of array_out to be filled integer, optional, intent(in) :: je !< The j- limits of array_out to be filled real, optional, intent(in) :: conversion !< A number that every element is multiplied by + + array_out(:,:) = -1. end subroutine extract_coupler_values !> Set element and index of a boundary condition diff --git a/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 b/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 index ea9f225a27..fec9c80461 100644 --- a/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 +++ b/config_src/external/GFDL_ocean_BGC/generic_tracer_utils.F90 @@ -143,6 +143,17 @@ subroutine g_tracer_get_common(isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau,& integer, optional, dimension(:,:), pointer :: grid_mask_coast !< Unknown integer, optional, dimension(:,:), pointer :: grid_kmt !< Unknown type(g_diag_ctrl), optional, pointer :: diag_CS !< Unknown + + isc = -1 + iec = -1 + jsc = -1 + jec = -1 + isd = -1 + ied = -1 + jsd = -1 + jed = -1 + nk = -1 + ntau = -1 end subroutine g_tracer_get_common !> Unknown @@ -177,6 +188,8 @@ subroutine g_tracer_get_4D_val(g_tracer_list,name,member,array,isd,jsd) integer, intent(in) :: isd !< Unknown integer, intent(in) :: jsd !< Unknown real, dimension(isd:,jsd:,:,:), intent(out):: array !< Unknown + + array(:,:,:,:) = -1. end subroutine g_tracer_get_4D_val !> Unknown @@ -190,6 +203,8 @@ subroutine g_tracer_get_3D_val(g_tracer_list,name,member,array,isd,jsd,ntau,posi logical, optional, intent(in) :: positive !< Unknown real, dimension(isd:,jsd:,:), intent(out):: array !< Unknown character(len=fm_string_len), parameter :: sub_name = 'g_tracer_get_3D_val' + + array(:,:,:) = -1. end subroutine g_tracer_get_3D_val !> Unknown @@ -200,6 +215,8 @@ subroutine g_tracer_get_2D_val(g_tracer_list,name,member,array,isd,jsd) integer, intent(in) :: isd !< Unknown integer, intent(in) :: jsd !< Unknown real, dimension(isd:,jsd:), intent(out):: array !< Unknown + + array(:,:) = -1. end subroutine g_tracer_get_2D_val !> Unknown @@ -208,6 +225,8 @@ subroutine g_tracer_get_real(g_tracer_list,name,member,value) character(len=*), intent(in) :: member !< Unknown type(g_tracer_type), pointer :: g_tracer_list !< Unknown real, intent(out):: value !< Unknown + + value = -1 end subroutine g_tracer_get_real !> Unknown @@ -216,6 +235,8 @@ subroutine g_tracer_get_string(g_tracer_list,name,member,string) character(len=*), intent(in) :: member !< Unknown type(g_tracer_type), pointer :: g_tracer_list !< Unknown character(len=fm_string_len), intent(out) :: string !< Unknown + + string = "" end subroutine g_tracer_get_string !> Unknown @@ -268,18 +289,24 @@ end subroutine g_tracer_send_diag subroutine g_tracer_get_name(g_tracer,string) type(g_tracer_type), pointer :: g_tracer !< Unknown character(len=*), intent(out) :: string !< Unknown + + string = "" end subroutine g_tracer_get_name !> Unknown subroutine g_tracer_get_alias(g_tracer,string) type(g_tracer_type), pointer :: g_tracer !< Unknown character(len=*), intent(out) :: string !< Unknown + + string = "" end subroutine g_tracer_get_alias !> Is the tracer prognostic? function g_tracer_is_prog(g_tracer) logical :: g_tracer_is_prog type(g_tracer_type), pointer :: g_tracer !< Pointer to tracer node + + g_tracer_is_prog = .false. end function g_tracer_is_prog !> get the next tracer in the list @@ -297,6 +324,8 @@ subroutine g_tracer_get_obc_segment_props(g_tracer_list, name, obc_has, src_file real, optional,intent(out):: lfac_out !< OBC reservoir inverse lengthscale factor character(len=*),optional,intent(out):: src_file !< OBC source file character(len=*),optional,intent(out):: src_var_name !< OBC source variable in file + + obc_has = .false. end subroutine g_tracer_get_obc_segment_props !>Vertical Diffusion of a tracer node diff --git a/config_src/external/database_comms/database_client_interface.F90 b/config_src/external/database_comms/database_client_interface.F90 index 9b57628921..8b05b83daf 100644 --- a/config_src/external/database_comms/database_client_interface.F90 +++ b/config_src/external/database_comms/database_client_interface.F90 @@ -317,6 +317,7 @@ function unpack_tensor_float_1d(self, name, data, dims) result(code) integer :: code code = -1 + data(:) = -1. end function unpack_tensor_float_1d !> Unpack a 32-bit real 2d tensor from the database @@ -328,6 +329,7 @@ function unpack_tensor_float_2d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:) = -1. end function unpack_tensor_float_2d !> Unpack a 32-bit real 3d tensor from the database @@ -339,6 +341,7 @@ function unpack_tensor_float_3d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:,:) = -1. end function unpack_tensor_float_3d !> Unpack a 32-bit real 4d tensor from the database @@ -350,6 +353,7 @@ function unpack_tensor_float_4d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:,:,:) = -1. end function unpack_tensor_float_4d !> Unpack a 64-bit real 1d tensor from the database @@ -361,6 +365,7 @@ function unpack_tensor_double_1d(self, name, data, dims) result(code) integer :: code code = -1 + data(:) = -1. end function unpack_tensor_double_1d !> Unpack a 64-bit real 2d tensor from the database @@ -372,6 +377,7 @@ function unpack_tensor_double_2d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:) = -1. end function unpack_tensor_double_2d !> Unpack a 64-bit real 3d tensor from the database @@ -383,6 +389,7 @@ function unpack_tensor_double_3d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:,:) = -1. end function unpack_tensor_double_3d !> Unpack a 64-bit real 4d tensor from the database @@ -394,6 +401,7 @@ function unpack_tensor_double_4d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:,:,:) = -1. end function unpack_tensor_double_4d !> Unpack a 32-bit integer 1d tensor from the database @@ -405,6 +413,7 @@ function unpack_tensor_int32_1d(self, name, data, dims) result(code) integer :: code code = -1 + data(:) = -1_int32 end function unpack_tensor_int32_1d !> Unpack a 32-bit integer 2d tensor from the database @@ -416,6 +425,7 @@ function unpack_tensor_int32_2d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:) = -1_int32 end function unpack_tensor_int32_2d !> Unpack a 32-bit integer 3d tensor from the database @@ -427,6 +437,7 @@ function unpack_tensor_int32_3d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:,:) = -1_int32 end function unpack_tensor_int32_3d !> Unpack a 32-bit integer 4d tensor from the database @@ -438,6 +449,7 @@ function unpack_tensor_int32_4d(self, name, data, dims) result(code) integer :: code code = -1 + data(:,:,:,:) = -1_int32 end function unpack_tensor_int32_4d !> Move a tensor to a new name @@ -479,6 +491,7 @@ function get_model(self, name, model) result(code) integer :: code code = -1 + model = "" end function get_model !> Load the machine learning model from a file and set the configuration @@ -621,6 +634,7 @@ function get_script(self, name, script) result(code) integer :: code code = -1 + script = "" end function get_script !> Set a script (from file) in the database for future execution @@ -735,7 +749,12 @@ function get_dataset(self, name, dataset) result(code) type(dataset_type), intent( out) :: dataset !< receives the dataset integer :: code + type(dataset_type) :: dataset_out + ! Placeholder dataset to prevent compiler warnings + ! Since dataset_type contains no data, any declared instance should work. + code = -1 + dataset = dataset_out end function get_dataset !> Rename a dataset stored in the database diff --git a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 index 18ccdaae67..232986f480 100644 --- a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 @@ -57,6 +57,8 @@ module MOM_diag_manager_infra public MOM_diag_manager_init public MOM_diag_manager_end public send_data_infra +public diag_send_complete_infra +public diag_manager_set_time_end_infra public MOM_diag_field_add_attribute public register_diag_field_infra public register_static_field_infra @@ -451,4 +453,13 @@ subroutine MOM_diag_field_add_attribute_i1d(diag_field_id, att_name, att_value) end subroutine MOM_diag_field_add_attribute_i1d +!> Needed for backwards compatibility, does nothing +subroutine diag_send_complete_infra () +end subroutine diag_send_complete_infra + +!> Needed for backwards compatibility, does nothing +subroutine diag_manager_set_time_end_infra(time) + type(time_type), intent(in) :: time !< The model time that simulation ends +end subroutine diag_manager_set_time_end_infra + end module MOM_diag_manager_infra diff --git a/config_src/infra/FMS2/MOM_diag_manager_infra.F90 b/config_src/infra/FMS2/MOM_diag_manager_infra.F90 index 18ccdaae67..f05baa4474 100644 --- a/config_src/infra/FMS2/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS2/MOM_diag_manager_infra.F90 @@ -14,13 +14,15 @@ module MOM_diag_manager_infra use diag_data_mod, only : null_axis_id use diag_manager_mod, only : fms_diag_manager_init => diag_manager_init use diag_manager_mod, only : fms_diag_manager_end => diag_manager_end +use diag_manager_mod, only : diag_send_complete +use diag_manager_mod, only : diag_manager_set_time_end use diag_manager_mod, only : send_data_fms => send_data use diag_manager_mod, only : fms_diag_field_add_attribute => diag_field_add_attribute use diag_manager_mod, only : DIAG_FIELD_NOT_FOUND use diag_manager_mod, only : register_diag_field_fms => register_diag_field use diag_manager_mod, only : register_static_field_fms => register_static_field use diag_manager_mod, only : get_diag_field_id_fms => get_diag_field_id -use MOM_time_manager, only : time_type +use MOM_time_manager, only : time_type, set_time use MOM_domain_infra, only : MOM_domain_type use MOM_error_infra, only : MOM_error => MOM_err, FATAL, WARNING @@ -57,6 +59,8 @@ module MOM_diag_manager_infra public MOM_diag_manager_init public MOM_diag_manager_end public send_data_infra +public diag_send_complete_infra +public diag_manager_set_time_end_infra public MOM_diag_field_add_attribute public register_diag_field_infra public register_static_field_infra @@ -451,4 +455,19 @@ subroutine MOM_diag_field_add_attribute_i1d(diag_field_id, att_name, att_value) end subroutine MOM_diag_field_add_attribute_i1d +!> Finishes the diag manager reduction methods as needed for the time_step +subroutine diag_send_complete_infra () + !! The time_step in the diag_send_complete call is a dummy argument, needed for backwards compatibility + !! It won't be used at all when diag_manager_nml::use_modern_diag=.true. + !! It won't have any impact when diag_manager_nml::use_modern_diag=.false. + call diag_send_complete (set_time(0)) +end subroutine diag_send_complete_infra + +!> Sets the time that the simulation ends in the diag manager +subroutine diag_manager_set_time_end_infra(time) + type(time_type), optional, intent(in) :: time !< The time the simulation ends + + call diag_manager_set_time_end(time) +end subroutine diag_manager_set_time_end_infra + end module MOM_diag_manager_infra diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index abcd821790..0fdf80bf52 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -393,9 +393,9 @@ subroutine check_reconstructions_1d(n0, h0, u0, deg, boundary_extrapolation, & real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A] integer, intent(in) :: deg !< Degree of polynomial reconstruction logical, intent(in) :: boundary_extrapolation !< Extrapolate at boundaries if true - real, dimension(n0,deg+1),intent(out) :: ppoly_r_coefs !< Coefficients of polynomial [A] - real, dimension(n0,2), intent(out) :: ppoly_r_E !< Edge value of polynomial [A] - real, dimension(n0,2), intent(out) :: ppoly_r_S !< Edge slope of polynomial [A H-1] + real, dimension(n0,deg+1),intent(in) :: ppoly_r_coefs !< Coefficients of polynomial [A] + real, dimension(n0,2), intent(in) :: ppoly_r_E !< Edge value of polynomial [A] + real, dimension(n0,2), intent(in) :: ppoly_r_S !< Edge slope of polynomial [A H-1] ! Local variables integer :: i0, n real :: u_l, u_c, u_r ! Cell averages [A] diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index fb53c638db..0b602be944 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3114,10 +3114,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & if (CS%debug) then call uvchksum("Post ALE adjust init cond [uv]", CS%u, CS%v, G%HI, haloshift=1) - call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=1, scale=GV%H_to_MKS) + call hchksum(CS%h, "Post ALE adjust init cond h", G%HI, haloshift=2, scale=GV%H_to_MKS) if (use_temperature) then - call hchksum(CS%tv%T, "Post ALE adjust init cond T", G%HI, haloshift=1, scale=US%C_to_degC) - call hchksum(CS%tv%S, "Post ALE adjust init cond S", G%HI, haloshift=1, scale=US%S_to_ppt) + call hchksum(CS%tv%T, "Post ALE adjust init cond T", G%HI, haloshift=2, scale=US%C_to_degC) + call hchksum(CS%tv%S, "Post ALE adjust init cond S", G%HI, haloshift=2, scale=US%S_to_ppt) endif endif endif @@ -3221,13 +3221,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & if (CS%split) then allocate(eta(SZI_(G),SZJ_(G)), source=0.0) if (CS%use_alt_split) then - call initialize_dyn_split_RK2b(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, & + call initialize_dyn_split_RK2b(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, & G, GV, US, param_file, diag, CS%dyn_split_RK2b_CSp, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil) else - call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, & + call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%tv, CS%uh, CS%vh, eta, Time, & G, GV, US, param_file, diag, CS%dyn_split_RK2_CSp, restart_CSp, & CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, & CS%thickness_diffuse_CSp, CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index b73524e362..6d80d4dd55 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -1753,6 +1753,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t real :: F_guess, F_l, F_r ! Fractional positions [nondim] real :: GxRho ! The product of the gravitational acceleration and reference density [R L2 Z-1 T-2 ~> Pa m-1] real :: Pa, Pa_left, Pa_right, Pa_tol ! Pressure anomalies, P = integral of g*(rho-rho_ref) dz [R L2 T-2 ~> Pa] + integer :: m ! A counter for how many iterations have been done in the while loop character(len=240) :: msg GxRho = G_e * rho_ref @@ -1780,8 +1781,14 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t F_guess = F_l - Pa_left / (Pa_right - Pa_left) * (F_r - F_l) Pa = Pa_right - Pa_left ! To get into iterative loop + m = 0 ! Reset the counter for the loop to be zero do while ( abs(Pa) > Pa_tol ) + m = m + 1 + if (m > 30) then ! Call an error, because convergence to the tolerance has not been achieved + write(msg,*) Pa_left,Pa,Pa_right,P_t-P_tgt,P_b-P_tgt + call MOM_error(FATAL, 'find_depth_of_pressure_in_cell completes too many iterations: '//msg) + endif z_out = z_t + ( z_b - z_t ) * F_guess Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS) - ( P_tgt - P_t ) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 0557ec7cd5..a184cf473f 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -840,7 +840,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f if (CS%debug) then call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) - call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_MKS) + call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, scale=GV%H_to_MKS) ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US) call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s) call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) @@ -849,7 +849,7 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f ! diffu = horizontal viscosity terms (u_av) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & - MEKE, Varmix, G, GV, US, CS%hor_visc, & + MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt, & OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v) call cpu_clock_end(id_clock_horvisc) @@ -1296,7 +1296,7 @@ end subroutine remap_dyn_split_RK2_aux_vars !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. -subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, & +subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, & diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & @@ -1310,6 +1310,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param intent(inout) :: v !< merid velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< layer thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -1518,7 +1519,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. & .not. query_initialized(CS%diffv, "diffv", restart_CS)) then call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, VarMix, G, GV, US, CS%hor_visc, & - OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & + tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v) call set_initialized(CS%diffu, "diffu", restart_CS) call set_initialized(CS%diffv, "diffv", restart_CS) diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index 44a0b0bf5c..fa78938477 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -558,7 +558,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc ! diffu = horizontal viscosity terms (u_av) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, & - OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%AD_pred) + tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%AD_pred) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with predictor horizontal_viscosity (step_MOM_dyn_split_RK2b)") @@ -837,7 +837,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc if (CS%debug) then call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) - call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_MKS) + call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, scale=GV%H_to_MKS) ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US) call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s) call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T) @@ -845,7 +845,7 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc ! diffu = horizontal viscosity terms (u_av) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, & + call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt, & OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, ADp=CS%ADp) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2b)") @@ -1222,7 +1222,7 @@ end subroutine remap_dyn_split_RK2b_aux_vars !> This subroutine initializes all of the variables that are used by this !! dynamic core, including diagnostics and the cpu clocks. -subroutine initialize_dyn_split_RK2b(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, & +subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, & diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, & VarMix, MEKE, thickness_diffuse_CSp, & OBC, update_OBC_CSp, ALE_CSp, set_visc, & @@ -1236,6 +1236,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, uh, vh, eta, Time, G, GV, US, para intent(inout) :: v !< merid velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< layer thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index f9e4aa0efe..579ddead2d 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -263,7 +263,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! diffu = horizontal viscosity terms (u,h) call enable_averages(dt, Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc) + call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc, tv, dt) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 1c589f509c..65b3bdf50e 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -276,7 +276,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call enable_averages(dt,Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_in, v_in, h_in, CS%diffu, CS%diffv, MEKE, VarMix, & - G, GV, US, CS%hor_visc) + G, GV, US, CS%hor_visc, tv, dt) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) call pass_vector(CS%diffu, CS%diffv, G%Domain, clock=id_clock_pass) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 167b2e81f3..fd8057c38f 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -136,6 +136,7 @@ module MOM_diagnostics integer :: id_sst = -1, id_sst_sq = -1, id_sstcon = -1 integer :: id_sss = -1, id_sss_sq = -1, id_sssabs = -1 integer :: id_ssu = -1, id_ssv = -1 + integer :: id_ssu_east = -1, id_ssv_north = -1 ! Diagnostic IDs for heat and salt flux fields integer :: id_fraz = -1 @@ -1283,6 +1284,8 @@ subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh) ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: speed ! The surface speed [L T-1 ~> m s-1] + real :: ssu_east(SZI_(G),SZJ_(G)) ! Surface velocity due east component [L T-1 ~> m s-1] + real :: ssv_north(SZI_(G),SZJ_(G)) ! Surface velocity due north component [L T-1 ~> m s-1] integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -1304,6 +1307,17 @@ subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh) call post_data(IDs%id_speed, speed, diag, mask=G%mask2dT) endif + if (IDs%id_ssu_east > 0 .or. IDs%id_ssv_north > 0) then + do j=js,je ; do i=is,ie + ssu_east(i,j) = ((0.5*(sfc_state%u(I-1,j) + sfc_state%u(I,j))) * G%cos_rot(i,j)) + & + ((0.5*(sfc_state%v(i,J-1) + sfc_state%v(i,J))) * G%sin_rot(i,j)) + ssv_north(i,j) = ((0.5*(sfc_state%v(i,J-1) + sfc_state%v(i,J))) * G%cos_rot(i,j)) - & + ((0.5*(sfc_state%u(I-1,j) + sfc_state%u(I,j))) * G%sin_rot(i,j)) + enddo ; enddo + if (IDs%id_ssu_east > 0 ) call post_data(IDs%id_ssu_east, ssu_east, diag, mask=G%mask2dT) + if (IDs%id_ssv_north > 0 ) call post_data(IDs%id_ssv_north, ssv_north, diag, mask=G%mask2dT) + endif + end subroutine post_surface_dyn_diags @@ -1912,6 +1926,10 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv) 'Sea Surface Meridional Velocity', 'm s-1', conversion=US%L_T_to_m_s) IDs%id_speed = register_diag_field('ocean_model', 'speed', diag%axesT1, Time, & 'Sea Surface Speed', 'm s-1', conversion=US%L_T_to_m_s) + IDs%id_ssu_east = register_diag_field('ocean_model', 'ssu_east', diag%axesT1, Time, & + 'Eastward velocity', 'm s-1', conversion=US%L_T_to_m_s) + IDs%id_ssv_north = register_diag_field('ocean_model', 'ssv_north', diag%axesT1, Time, & + 'Northward velocity', 'm s-1', conversion=US%L_T_to_m_s) if (associated(tv%T)) then IDs%id_sst = register_diag_field('ocean_model', 'SST', diag%axesT1, Time, & diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index 2567e7591b..a590ae3893 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -152,6 +152,9 @@ subroutine find_obsolete_params(param_file) call obsolete_logical(param_file, "VERT_FRICTION_2018_ANSWERS", & hint="Instead use VERT_FRICTION_ANSWER_DATE.") + call obsolete_logical(param_file, "USE_GRID_SPACE_DIAGNOSTIC_AXES", & + hint="Instead use USE_INDEX_DIAGNOSTIC_AXIS.") + ! Write the file version number to the model log. call log_version(param_file, mdl, version) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 18dcd9a825..b3194af3d8 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -14,6 +14,7 @@ module MOM_diag_mediator use MOM_diag_manager_infra, only : send_data_infra, MOM_diag_field_add_attribute, EAST, NORTH use MOM_diag_manager_infra, only : register_diag_field_infra, register_static_field_infra use MOM_diag_manager_infra, only : get_MOM_diag_field_id, DIAG_FIELD_NOT_FOUND +use MOM_diag_manager_infra, only : diag_send_complete_infra use MOM_diag_remap, only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field @@ -243,7 +244,7 @@ module MOM_diag_mediator !! This file is open if available_diag_doc_unit is > 0. logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics logical :: show_call_tree !< Display the call tree while running. Set by VERBOSITY level. - logical :: grid_space_axes !< If true, diagnostic horizontal coordinates axes are in grid space. + logical :: index_space_axes !< If true, diagnostic horizontal coordinates axes are in index space. ! The following fields are used for the output of the data. integer :: is !< The start i-index of cell centers within the computational domain integer :: ie !< The end i-index of cell centers within the computational domain @@ -371,7 +372,7 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical - if (diag_cs%grid_space_axes) then + if (diag_cs%index_space_axes) then allocate(IaxB(G%IsgB:G%IegB)) do i=G%IsgB, G%IegB Iaxb(i)=real(i) @@ -392,7 +393,7 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) ! Horizontal axes for the native grids if (G%symmetric) then - if (diag_cs%grid_space_axes) then + if (diag_cs%index_space_axes) then id_xq = diag_axis_init('iq', IaxB(G%isgB:G%iegB), 'none', 'x', & 'q point grid-space longitude', G%Domain, position=EAST) id_yq = diag_axis_init('jq', JaxB(G%jsgB:G%jegB), 'none', 'y', & @@ -404,7 +405,7 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) 'q point nominal latitude', G%Domain, position=NORTH) endif else - if (diag_cs%grid_space_axes) then + if (diag_cs%index_space_axes) then id_xq = diag_axis_init('Iq', IaxB(G%isg:G%ieg), 'none', 'x', & 'q point grid-space longitude', G%Domain, position=EAST) id_yq = diag_axis_init('Jq', JaxB(G%jsg:G%jeg), 'none', 'y', & @@ -417,11 +418,11 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) endif endif - if (diag_cs%grid_space_axes) then + if (diag_cs%index_space_axes) then id_xh = diag_axis_init('ih', iax(G%isg:G%ieg), 'none', 'x', & - 'h point grid-space longitude', G%Domain, position=EAST) + 'h point grid-space longitude', G%Domain) id_yh = diag_axis_init('jh', jax(G%jsg:G%jeg), 'none', 'y', & - 'h point grid space latitude', G%Domain, position=NORTH) + 'h point grid space latitude', G%Domain) else id_xh = diag_axis_init('xh', G%gridLonT(G%isg:G%ieg), G%x_axis_units, 'x', & 'h point nominal longitude', G%Domain) @@ -579,7 +580,7 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) endif enddo - if (diag_cs%grid_space_axes) then + if (diag_cs%index_space_axes) then deallocate(IaxB, iax, JaxB, jax) endif !Define the downsampled axes @@ -2078,6 +2079,7 @@ end subroutine enable_averages subroutine disable_averaging(diag_cs) type(diag_ctrl), intent(inout) :: diag_CS !< Structure used to regulate diagnostic output + call diag_send_complete_infra() diag_cs%time_int = 0.0 diag_cs%ave_enabled = .false. @@ -3287,7 +3289,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) - call get_param(param_file, mdl, 'USE_GRID_SPACE_DIAGNOSTIC_AXES', diag_cs%grid_space_axes, & + call get_param(param_file, mdl, 'USE_INDEX_DIAGNOSTIC_AXES', diag_cs%index_space_axes, & 'If true, use a grid index coordinate convention for diagnostic axes. ',& default=.false.) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index b718f65b5e..ce739cba55 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -125,7 +125,7 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, real, dimension(SZI_(G),SZJ_(G)) :: good_new ! The values of good_ to use for the next iteration [nondim] real :: east, west, north, south ! Valid neighboring values or 0 for invalid values [arbitrary] - real :: ge, gw, gn, gs ! Flags indicating which neighbors have valid values [nondim] + real :: ge, gw, gn, gs ! Flags set to 0 or 1 indicating which neighbors have valid values [nondim] real :: ngood ! The number of valid values in neighboring points [nondim] real :: nfill ! The remaining number of points to fill [nondim] real :: nfill_prev ! The previous value of nfill [nondim] @@ -227,23 +227,30 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, ! Do Laplacian smoothing for the points that have been filled in. do k=1,npass call pass_var(aout,G%Domain) - do j=js,je ; do i=is,ie - if (fill(i,j) == 1) then - east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j)) - north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1)) - if (ans_2018) then + + a_chg(:,:) = 0.0 + if (ans_2018) then + do j=js,je ; do i=is,ie + if (fill(i,j) == 1) then + east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j)) + north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1)) a_chg(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + & west*aout(i-1,j)+east*aout(i+1,j) - & (south+north+west+east)*aout(i,j)) - else - a_chg(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + & - (west*aout(i-1,j)+east*aout(i+1,j))) - & - ((south+north)+(west+east))*aout(i,j) ) endif - else - a_chg(i,j) = 0. - endif - enddo ; enddo + enddo ; enddo + else + do j=js,je ; do i=is,ie + if (fill(i,j) == 1) then + ge = max(good(i+1,j),fill(i+1,j)) ; gw = max(good(i-1,j),fill(i-1,j)) + gn = max(good(i,j+1),fill(i,j+1)) ; gs = max(good(i,j-1),fill(i,j-1)) + a_chg(i,j) = relax_coeff*( ((gs*aout(i,j-1) + gn*aout(i,j+1)) + & + (gw*aout(i-1,j) + ge*aout(i+1,j))) - & + ((gs + gn) + (gw + ge))*aout(i,j) ) + endif + enddo ; enddo + endif + ares = 0.0 do j=js,je ; do i=is,ie aout(i,j) = a_chg(i,j) + aout(i,j) diff --git a/src/framework/MOM_io_file.F90 b/src/framework/MOM_io_file.F90 index 6697e56f68..261d4b628d 100644 --- a/src/framework/MOM_io_file.F90 +++ b/src/framework/MOM_io_file.F90 @@ -1709,6 +1709,8 @@ subroutine read_field_chksum_nc(handle, field, chksum, valid_chksum) !< If true, chksum has been successfully read call MOM_error(FATAL, 'read_field_chksum over netCDF is not yet implemented.') + chksum = -1_int64 + valid_chksum = .false. end subroutine read_field_chksum_nc diff --git a/src/framework/posix.F90 b/src/framework/posix.F90 index fffb619cba..1087958939 100644 --- a/src/framework/posix.F90 +++ b/src/framework/posix.F90 @@ -84,7 +84,7 @@ function stat_posix(path, buf) result(rc) bind(c, name="stat") character(kind=c_char), dimension(*), intent(in) :: path !< Pathname of a POSIX file - type(stat_buf), intent(in) :: buf + type(stat_buf), intent(inout) :: buf !< Information describing the file if it exists integer(kind=c_int) :: rc !< Function return code diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index dabb075cf3..42fa63fd95 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -514,9 +514,9 @@ function register_MOM_IS_diag_field(module_name, field_name, axes, init_time, & end function register_MOM_IS_diag_field !> Registers a static diagnostic, returning an integer handle -integer function register_MOM_IS_static_field(module_name, field_name, axes, & - long_name, units, missing_value, range, mask_variant, standard_name, & - do_not_log, interp_method, tile_count) +function register_MOM_IS_static_field(module_name, field_name, axes, & + long_name, units, missing_value, range, mask_variant, standard_name, & + do_not_log, interp_method, tile_count) result(register_static_field) integer :: register_static_field !< The returned diagnostic handle character(len=*), intent(in) :: module_name !< Name of this module, usually "ice_model" character(len=*), intent(in) :: field_name !< Name of the diagnostic field diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 4f11233c93..c18752c83d 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -343,7 +343,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & just_read=just_read) case ("dumbbell"); call dumbbell_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("soliton"); call soliton_initialize_thickness(dz, depth_tot, G, GV, US) + case ("soliton"); call soliton_initialize_thickness(dz, depth_tot, G, GV, US, PF, & + just_read=just_read) case ("phillips"); call Phillips_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) case ("rossby_front") @@ -508,7 +509,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & case ("phillips"); call Phillips_initialize_velocity(u, v, G, GV, US, PF, just_read) case ("rossby_front"); call Rossby_front_initialize_velocity(u, v, h, & G, GV, US, PF, just_read) - case ("soliton"); call soliton_initialize_velocity(u, v, G, GV, US) + case ("soliton"); call soliton_initialize_velocity(u, v, G, GV, US, PF, just_read) case ("USER"); call user_initialize_velocity(u, v, G, GV, US, PF, just_read) case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& "Unrecognized velocity configuration "//trim(config)) diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index d28a925c03..cac8a5cd6c 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -200,8 +200,10 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ enddo ; enddo if (h_is_in_Z_units) then + ! Because h is in units of [Z ~> m], dzSrc is already in the right units, but we need to + ! specify negligible thickness values with the right units. dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) - call ALE_remap_scalar(remapCS, G, GV, kd, hSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date, & + call ALE_remap_scalar(remapCS, G, GV, kd, dzSrc, tr_z, h, tr, all_cells=.false., answer_date=remap_answer_date, & H_neglect=dz_neglect, H_neglect_edge=dz_neglect_edge) else ! Equation of state data is not available, so a simpler rescaling will have to suffice, diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index ab7a6fa5fc..2eef171bf5 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -13,7 +13,8 @@ module MOM_hor_visc use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_Leith_viscosity +use MOM_interface_heights, only : thickness_to_dz +use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_slopes, calc_QG_Leith_viscosity use MOM_barotropic, only : barotropic_CS, barotropic_get_tav use MOM_thickness_diffuse, only : thickness_diffuse_CS, thickness_diffuse_get_KH use MOM_io, only : MOM_read_data, slasher @@ -22,9 +23,9 @@ module MOM_hor_visc use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type -use MOM_variables, only : accel_diag_ptrs -use MOM_Zanna_Bolton, only : ZB2020_lateral_stress, ZB2020_init, ZB2020_end, & - ZB2020_CS, ZB2020_copy_gradient_and_thickness +use MOM_variables, only : accel_diag_ptrs, thermo_var_ptrs +use MOM_Zanna_Bolton, only : ZB2020_lateral_stress, ZB2020_init, ZB2020_end +use MOM_Zanna_Bolton, only : ZB2020_CS, ZB2020_copy_gradient_and_thickness implicit none ; private @@ -237,11 +238,11 @@ module MOM_hor_visc !! !! To work, the following fields must be set outside of the usual !! is:ie range before this subroutine is called: -!! u[is-2:ie+2,js-2:je+2] -!! v[is-2:ie+2,js-2:je+2] -!! h[is-1:ie+1,js-1:je+1] +!! u(is-2:ie+2,js-2:je+2) +!! v(is-2:ie+2,js-2:je+2) +!! h(is-1:ie+1,js-1:je+1) or up to h(is-2:ie+2,js-2:je+2) with some Leith options. subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & - CS, OBC, BT, TD, ADp, hu_cont, hv_cont) + CS, tv, dt, OBC, BT, TD, ADp, hu_cont, hv_cont) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & @@ -259,12 +260,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, type(MEKE_type), intent(inout) :: MEKE !< MEKE fields !! related to Mesoscale Eddy Kinetic Energy. type(VarMix_CS), intent(inout) :: VarMix !< Variable mixing control structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure - type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type - type(barotropic_CS), intent(in), optional :: BT !< Barotropic control structure - type(thickness_diffuse_CS), intent(in), optional :: TD !< Thickness diffusion control structure - type(accel_diag_ptrs), intent(in), optional :: ADp !< Acceleration diagnostics + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables + real, intent(in) :: dt !< Time increment [T ~> s] + type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type + type(barotropic_CS), optional, intent(in) :: BT !< Barotropic control structure + type(thickness_diffuse_CS), optional, intent(in) :: TD !< Thickness diffusion control structure + type(accel_diag_ptrs), optional, intent(in) :: ADp !< Acceleration diagnostics real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -341,12 +345,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, GME_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1] ShSt ! A diagnostic array of shear stress [T-1 ~> s-1]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: & - KH_u_GME !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1] + KH_u_GME, & !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1] + slope_x !< Isopycnal slope in i-direction [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: & - KH_v_GME !< Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1] + KH_v_GME, & !< Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1] + slope_y !< Isopycnal slope in j-direction [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & Ah_h, & ! biharmonic viscosity at thickness points [L4 T-1 ~> m4 s-1] Kh_h, & ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1] + dz, & ! Height change across layers [Z ~> m] FrictWork, & ! work done by MKE dissipation mechanisms [R L2 T-3 ~> W m-2] FrictWork_GME, & ! work done by GME [R L2 T-3 ~> W m-2] div_xx_h, & ! horizontal divergence [T-1 ~> s-1] @@ -590,6 +597,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call pass_vector(u_smooth, v_smooth, G%Domain) endif + if (CS%use_QG_Leith_visc .and. ((CS%Leith_Kh) .or. (CS%Leith_Ah))) then + call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=2) + ! Calculate isopycnal slopes that will be used for some forms of viscosity. + call calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, VarMix, OBC) + ! If the following halo update is added, the calculations in calc_QG_slopes could work on just + ! the computational domains, and some halo updates outside of this routine could be smaller. + ! call pass_vector(slope_x, slope_y, G%Domain, halo=2) + endif + !$OMP parallel do default(none) & !$OMP shared( & !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, & @@ -597,7 +613,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP is_vort, ie_vort, js_vort, je_vort, & !$OMP is_Kh, ie_Kh, js_Kh, je_Kh, & !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, & - !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, & + !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, slope_x, slope_y, dz, & !$OMP backscat_subround, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, & !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_GME, & @@ -1008,9 +1024,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, (0.5*(vort_xy_dy(I,j) + vort_xy_dy(I,j+1)))**2 ) enddo ; enddo - ! This accumulates terms, some of which are in VarMix, so rescaling can not be done here. - call calc_QG_Leith_viscosity(VarMix, G, GV, US, h, k, div_xx_dx, div_xx_dy, & - vort_xy_dx, vort_xy_dy) + ! This accumulates terms, some of which are in VarMix. + call calc_QG_Leith_viscosity(VarMix, G, GV, US, h, dz, k, div_xx_dx, div_xx_dy, & + slope_x, slope_y, vort_xy_dx, vort_xy_dy) endif @@ -2207,12 +2223,12 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, & "If true, use QG Leith nonlinear eddy viscosity.", & default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) ) - if (CS%use_QG_Leith_visc) then - call MOM_error(FATAL, "USE_QG_LEITH_VISC=True activates code that is a work-in-progress and "//& - "should not be used until a number of bugs are fixed. Specifically it does not "//& - "reproduce across PE count or layout, and may use arrays that have not been properly "//& - "set or allocated. See github.com/mom-ocean/MOM6/issues/1590 for a discussion.") - endif +! if (CS%use_QG_Leith_visc) then +! call MOM_error(FATAL, "USE_QG_LEITH_VISC=True activates code that is a work-in-progress and "//& +! "should not be used until a number of bugs are fixed. Specifically it does not "//& +! "reproduce across PE count or layout, and may use arrays that have not been properly "//& +! "set or allocated. See github.com/mom-ocean/MOM6/issues/1590 for a discussion.") +! endif if (CS%use_QG_Leith_visc .and. .not. (CS%Leith_Kh .or. CS%Leith_Ah) ) then call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//& "LEITH_KH or LEITH_AH must be True when USE_QG_LEITH_VISC=True.") diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 4f1dbb89ac..defbd78aa7 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -170,7 +170,7 @@ module MOM_lateral_mixing_coeffs end type VarMix_CS public VarMix_init, VarMix_end, calc_slope_functions, calc_resoln_function -public calc_QG_Leith_viscosity, calc_depth_function +public calc_QG_slopes, calc_QG_Leith_viscosity, calc_depth_function contains @@ -474,14 +474,13 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) type(VarMix_CS), intent(inout) :: CS !< Variable mixing control structure type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure ! Local variables - real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & - e ! The interface heights relative to mean sea level [Z ~> m]. - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2] - real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2] - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m] - real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m] - real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] - real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: N2_u ! Square of Brunt-Vaisala freq at u-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: N2_v ! Square of Brunt-Vaisala freq at v-points [L2 Z-2 T-2 ~> s-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzu ! Z-thickness at u-points [Z ~> m] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzv ! Z-thickness at v-points [Z ~> m] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: dzSxN ! |Sx| N times dz at u-points [Z T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: dzSyN ! |Sy| N times dz at v-points [Z T-1 ~> m s-1] if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//& "Module must be initialized before it is used.") @@ -996,18 +995,47 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop end subroutine calc_slope_functions_using_just_e + +!> Calculates and returns isopycnal slopes with wider halos for use in finding QG viscosity. +subroutine calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, CS, OBC) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + real, intent(in) :: dt !< Time increment [T ~> s] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] + type(VarMix_CS), intent(in) :: CS !< Variable mixing control structure + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! The interface heights relative to mean sea level [Z ~> m] + + if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_QG_slopes: "//& + "Module must be initialized before it is used.") + + call find_eta(h, tv, G, GV, US, e, halo_size=3) + call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, CS%use_stanley_iso, & + slope_x, slope_y, halo=2, OBC=OBC) + +end subroutine calc_QG_slopes + !> Calculates the Leith Laplacian and bi-harmonic viscosity coefficients -subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vort_xy_dx, vort_xy_dy) +subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, dz, k, div_xx_dx, div_xx_dy, slope_x, slope_y, & + vort_xy_dx, vort_xy_dy) type(VarMix_CS), intent(inout) :: CS !< Variable mixing coefficients - type(ocean_grid_type), intent(in) :: G !< Ocean 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_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] - integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude - real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence + type(ocean_grid_type), intent(in) :: G !< Ocean 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_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer vertical extents [Z ~> m] + integer, intent(in) :: k !< Layer for which to calculate vorticity magnitude + real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: div_xx_dx !< x-derivative of horizontal divergence !! (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] - real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence + real, dimension(SZI_(G),SZJB_(G)), intent(in) :: div_xx_dy !< y-derivative of horizontal divergence !! (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: vort_xy_dx !< x-derivative of vertical vorticity !! (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1] real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: vort_xy_dy !< y-derivative of vertical vorticity @@ -1030,6 +1058,8 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo real :: h_at_slope_below ! The thickness below [H ~> m or kg m-2] real :: Ih ! The inverse of a combination of thicknesses [H-1 ~> m-1 or m2 kg-1] real :: f ! A copy of the Coriolis parameter [T-1 ~> s-1] + real :: Z_to_H ! A local copy of depth to thickness conversion factors or the inverse of the + ! mass-weighted average specific volumes around an interface [H Z-1 ~> nondim or kg m-3] real :: inv_PI3 ! The inverse of pi cubed [nondim] integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -1038,41 +1068,41 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo nz = GV%ke inv_PI3 = 1.0 / ((4.0*atan(1.0))**3) + Z_to_H = GV%Z_to_H ! This will be replaced with a varying value in non-Boussinesq mode. if ((k > 1) .and. (k < nz)) then - ! With USE_QG_LEITH_VISC=True, this might need to change to - ! do j=js-2,je+2 ; do I=is-2,ie+1 - ! but other arrays used here (e.g., h and CS%slope_x) would also need to have wider valid halos. - do j=js-1,je+1 ; do I=is-2,Ieq+1 + do j=js-2,je+2 ; do I=is-2,ie+1 h_at_slope_above = 2. * ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) * h(i+1,j,k) ) / & ( ( h(i,j,k-1) * h(i+1,j,k-1) ) * ( h(i,j,k) + h(i+1,j,k) ) & - + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff**2 ) + + ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k-1) + h(i+1,j,k-1) ) + GV%H_subroundoff**3 ) h_at_slope_below = 2. * ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) * h(i+1,j,k+1) ) / & ( ( h(i,j,k) * h(i+1,j,k) ) * ( h(i,j,k+1) + h(i+1,j,k+1) ) & - + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff**2 ) - Ih = 1. / ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) - dslopex_dz(I,j) = 2. * ( CS%slope_x(i,j,k) - CS%slope_x(i,j,k+1) ) * (GV%Z_to_H * Ih) + + ( h(i,j,k+1) * h(i+1,j,k+1) ) * ( h(i,j,k) + h(i+1,j,k) ) + GV%H_subroundoff**3 ) + Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) + if (.not.GV%Boussinesq) & + Z_to_H = ( (h(i,j,k-1) + h(i+1,j,k-1)) + (h(i,j,k) + h(i+1,j,k)) ) / & + ( (dz(i,j,k-1) + dz(i+1,j,k-1)) + (dz(i,j,k) + dz(i+1,j,k)) + GV%dZ_subroundoff) + dslopex_dz(I,j) = 2. * ( slope_x(I,j,k) - slope_x(I,j,k+1) ) * (Z_to_H * Ih) h_at_u(I,j) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - ! With USE_QG_LEITH_VISC=True, this might need to change to - ! do J=js-2,je+1 ; do i=is-2,ie+2 - do J=js-2,Jeq+1 ; do i=is-1,ie+1 + do J=js-2,je+1 ; do i=is-2,ie+2 h_at_slope_above = 2. * ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) * h(i,j+1,k) ) / & ( ( h(i,j,k-1) * h(i,j+1,k-1) ) * ( h(i,j,k) + h(i,j+1,k) ) & - + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff**2 ) + + ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k-1) + h(i,j+1,k-1) ) + GV%H_subroundoff**3 ) h_at_slope_below = 2. * ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) * h(i,j+1,k+1) ) / & ( ( h(i,j,k) * h(i,j+1,k) ) * ( h(i,j,k+1) + h(i,j+1,k+1) ) & - + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff**2 ) - Ih = 1. / ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) - dslopey_dz(i,J) = 2. * ( CS%slope_y(i,j,k) - CS%slope_y(i,j,k+1) ) * (GV%Z_to_H * Ih) + + ( h(i,j,k+1) * h(i,j+1,k+1) ) * ( h(i,j,k) + h(i,j+1,k) ) + GV%H_subroundoff**3 ) + Ih = 1./ ( h_at_slope_above + h_at_slope_below + GV%H_subroundoff ) + if (.not.GV%Boussinesq) & + Z_to_H = ( (h(i,j,k-1) + h(i,j+1,k-1)) + (h(i,j,k) + h(i,j+1,k)) ) / & + ( (dz(i,j,k-1) + dz(i,j+1,k-1)) + (dz(i,j,k) + dz(i,j+1,k)) + GV%dZ_subroundoff) + dslopey_dz(i,J) = 2. * ( slope_y(i,J,k) - slope_y(i,J,k+1) ) * (Z_to_H * Ih) h_at_v(i,J) = 2. * ( h_at_slope_above * h_at_slope_below ) * Ih enddo ; enddo - ! With USE_QG_LEITH_VISC=True, this might need to be - ! do J=js-2,je+1 ; do i=is-1,ie+1 - do J=js-1,je ; do i=is-1,Ieq+1 + do J=js-2,je+1 ; do i=is-1,ie+1 f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J) ) vort_xy_dx(i,J) = vort_xy_dx(i,J) - f * & ( ( h_at_u(I,j) * dslopex_dz(I,j) + h_at_u(I-1,j+1) * dslopex_dz(I-1,j+1) ) & @@ -1080,9 +1110,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo ( ( h_at_u(I,j) + h_at_u(I-1,j+1) ) + ( h_at_u(I-1,j) + h_at_u(I,j+1) ) + GV%H_subroundoff) enddo ; enddo - ! With USE_QG_LEITH_VISC=True, this might need to be - ! do j=js-1,je+1 ; do I=is-2,ie+1 - do j=js-1,Jeq+1 ; do I=is-1,ie + do j=js-1,je+1 ; do I=is-2,ie+1 f = 0.5 * ( G%CoriolisBu(I,J) + G%CoriolisBu(I,J-1) ) vort_xy_dy(I,j) = vort_xy_dy(I,j) - f * & ( ( h_at_v(i,J) * dslopey_dz(i,J) + h_at_v(i+1,J-1) * dslopey_dz(i+1,J-1) ) & @@ -1100,7 +1128,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo + (div_xx_dy(i+1,J) + div_xx_dy(i,J-1))))**2) if (CS%use_beta_in_QG_Leith) then beta_u(I,j) = sqrt((0.5*(G%dF_dx(i,j)+G%dF_dx(i+1,j))**2) + & - (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2)) + (0.5*(G%dF_dy(i,j)+G%dF_dy(i+1,j))**2)) CS%KH_u_QG(I,j,k) = MIN(grad_vort_mag_u(I,j) + grad_div_mag_u(I,j), 3.0*beta_u(I,j)) * & CS%Laplac3_const_u(I,j) * inv_PI3 else @@ -1116,7 +1144,7 @@ subroutine calc_QG_Leith_viscosity(CS, G, GV, US, h, k, div_xx_dx, div_xx_dy, vo + (div_xx_dx(I,j+1) + div_xx_dx(I-1,j))))**2) if (CS%use_beta_in_QG_Leith) then beta_v(i,J) = sqrt((0.5*(G%dF_dx(i,j)+G%dF_dx(i,j+1))**2) + & - (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2)) + (0.5*(G%dF_dy(i,j)+G%dF_dy(i,j+1))**2)) CS%KH_v_QG(i,J,k) = MIN(grad_vort_mag_v(i,J) + grad_div_mag_v(i,J), 3.0*beta_v(i,J)) * & CS%Laplac3_const_v(i,J) * inv_PI3 else diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90 index 7f7215c9d8..e7f8a73ab2 100644 --- a/src/parameterizations/lateral/MOM_self_attr_load.F90 +++ b/src/parameterizations/lateral/MOM_self_attr_load.F90 @@ -20,14 +20,17 @@ module MOM_self_attr_load !> The control structure for the MOM_self_attr_load module type, public :: SAL_CS ; private - logical :: use_sal_scalar !< If true, use the scalar approximation to calculate SAL. - logical :: use_sal_sht !< If true, use online spherical harmonics to calculate SAL - logical :: use_tidal_sal_prev !< If true, read the tidal SAL from the previous iteration of - !! the tides to facilitate convergence. + logical :: use_sal_scalar = .false. + !< If true, use the scalar approximation to calculate SAL. + logical :: use_sal_sht = .false. + !< If true, use online spherical harmonics to calculate SAL + logical :: use_tidal_sal_prev = .false. + !< If true, read the tidal SAL from the previous iteration of the tides to + !! facilitate convergence. real :: sal_scalar_value !< The constant of proportionality between sea surface height !! (really it should be bottom pressure) anomalies and bottom !! geopotential anomalies [nondim]. - type(sht_CS) :: sht !< Spherical harmonic transforms (SHT) control structure + type(sht_CS), allocatable :: sht !< Spherical harmonic transforms (SHT) control structure integer :: sal_sht_Nd !< Maximum degree for SHT [nodim] real, allocatable :: Love_Scaling(:) !< Love number for each SHT mode [nodim] real, allocatable :: Snm_Re(:), & !< Real SHT coefficient for SHT SAL [Z ~> m] @@ -218,6 +221,8 @@ subroutine SAL_init(G, US, param_file, CS) allocate(CS%Love_Scaling(lmax)); CS%Love_Scaling(:) = 0.0 call calc_love_scaling(CS%sal_sht_Nd, rhoW, rhoE, CS%Love_Scaling) + + allocate(CS%sht) call spherical_harmonics_init(G, param_file, CS%sht) endif @@ -234,6 +239,7 @@ subroutine SAL_end(CS) if (allocated(CS%Snm_Re)) deallocate(CS%Snm_Re) if (allocated(CS%Snm_Im)) deallocate(CS%Snm_Im) call spherical_harmonics_end(CS%sht) + deallocate(CS%sht) endif end subroutine SAL_end diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index ef2c3125cd..b77949342d 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -380,7 +380,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & real :: dA ! Difference between the reconstruction tracer edge values [conc] real :: mA ! Average of the reconstruction tracer edge values [conc] real :: a6 ! Curvature of the reconstruction tracer values [conc] - logical :: do_i(SZIB_(G),SZJ_(G)) ! If true, work on given points. + logical :: do_i(SZI_(G),SZJ_(G)) ! If true, work on given points. logical :: usePLMslope integer :: i, j, m, n, i_up, stencil, ntr_id type(OBC_segment_type), pointer :: segment=>NULL() @@ -659,7 +659,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & enddo ! diagnostics - if (associated(Tr(m)%ad_x)) then ; do I=is-1,ie ; if (do_i(i,j)) then + if (associated(Tr(m)%ad_x)) then ; do I=is-1,ie ; if (do_i(i,j) .or. do_i(i+1,j)) then Tr(m)%ad_x(I,j,k) = Tr(m)%ad_x(I,j,k) + flux_x(I,j,m)*Idt endif ; enddo ; endif @@ -688,7 +688,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & !$OMP ordered do m=1,ntr ; if (associated(Tr(m)%ad2d_x)) then do j=js,je ; if (domore_u_initial(j,k)) then - do I=is-1,ie ; if (do_i(i,j)) then + do I=is-1,ie ; if (do_i(i,j) .or. do_i(i+1,j)) then Tr(m)%ad2d_x(I,j) = Tr(m)%ad2d_x(I,j) + flux_x(I,j,m)*Idt endif ; enddo endif ; enddo @@ -756,7 +756,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & real :: mA ! Average of the reconstruction tracer edge values [conc] real :: a6 ! Curvature of the reconstruction tracer values [conc] logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles. - logical :: do_i(SZIB_(G), SZJ_(G)) ! If true, work on given points. + logical :: do_i(SZI_(G), SZJ_(G)) ! If true, work on given points. logical :: usePLMslope integer :: i, j, j2, m, n, j_up, stencil, ntr_id type(OBC_segment_type), pointer :: segment=>NULL() @@ -1066,8 +1066,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & !$OMP ordered do m=1,ntr ; if (associated(Tr(m)%ad_y)) then do J=js-1,je ; if (domore_v_initial(J)) then - ! (The logical test could be "do_i(i,j) .or. do_i(i+1,j)" to be clearer, but not needed) - do i=is,ie ; if (do_i(i,j)) then + do i=is,ie ; if (do_i(i,j) .or. do_i(i,j+1)) then Tr(m)%ad_y(i,J,k) = Tr(m)%ad_y(i,J,k) + flux_y(i,m,J)*Idt endif ; enddo endif ; enddo @@ -1075,7 +1074,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & do m=1,ntr ; if (associated(Tr(m)%ad2d_y)) then do J=js-1,je ; if (domore_v_initial(J)) then - do i=is,ie ; if (do_i(i,j)) then + do i=is,ie ; if (do_i(i,j) .or. do_i(i,j+1)) then Tr(m)%ad2d_y(i,J) = Tr(m)%ad2d_y(i,J) + flux_y(i,m,J)*Idt endif ; enddo endif ; enddo diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index 06a781ec94..722a41b7e5 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -26,8 +26,9 @@ module soliton_initialization contains -!> Initialization of thicknesses in Equatorial Rossby soliton test -subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US) +!> Initialization of thicknesses in equatorial Rossby soliton test, as described in section +!! 6.1 of Haidvogel and Beckman (1990) and in Boyd (1980, JPO) and Boyd (1985, JPO). +subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US, param_file, just_read) 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 @@ -35,45 +36,96 @@ subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US) intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, intent(in) :: just_read !< If true, this call will only read + !! parameters without changing h. + ! Local variables + real :: max_depth ! Maximum depth of the model bathymetry [Z ~> m] + real :: cg_max ! The external wave speed based on max_depth [L T-1 ~> m s-1] + real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1] + real :: L_eq ! The equatorial deformation radius used in nondimensionalizing this problem [L ~> m] + real :: scale_pos ! A conversion factor to nondimensionalize the axis units, usually [m-1] + real :: x0 ! Initial x-position of the soliton in the same units as geoLonT, often [m]. + real :: y0 ! Initial y-position of the soliton in the same units as geoLatT, often [m]. + real :: x, y ! Nondimensionalized positions [nondim] + real :: I_nz ! The inverse of the number of layers [nondim] + real :: val1 ! A nondimensionlized zonal decay scale [nondim] + real :: val2 ! An overall surface height anomaly amplitude [L T-1 ~> m s-1] + real :: val3 ! A decay factor [nondim] + real :: val4 ! The local velocity amplitude [L T-1 ~> m s-1] + ! This include declares and sets the variable "version". +# include "version_variable.h" integer :: i, j, k, is, ie, js, je, nz - real :: x, y, x0, y0 - real :: val1, val2, val3, val4 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness") + if (.not.just_read) & + call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness") + + if (.not.just_read) call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, & + units="m", default=-1.e9, scale=US%m_to_Z, do_not_log=.true.) + call get_param(param_file, mdl, "BETA", beta, & + "The northward gradient of the Coriolis parameter with the betaplane option.", & + units="m-1 s-1", default=0.0, scale=US%T_to_s*US%L_to_m, do_not_log=.true.) + + if (just_read) return ! All run-time parameters have been read, so return. + + if (max_depth <= 0.0) call MOM_error(FATAL, & + "soliton_initialization, soliton_initialize_thickness: "//& + "This module requires a positive value of MAXIMUM_DEPTH.") + if (abs(beta) <= 0.0) call MOM_error(FATAL, & + "soliton_initialization, soliton_initialize_thickness: "//& + "This module requires a non-zero value of BETA.") + + cg_max = sqrt(GV%g_Earth * max_depth) + L_eq = sqrt(cg_max / abs(beta)) + scale_pos = US%m_to_L / L_eq + I_nz = 1.0 / real(nz) x0 = 2.0*G%len_lon/3.0 y0 = 0.0 val1 = 0.395 - val2 = US%m_to_Z * 0.771*(val1*val1) + val2 = max_depth * 0.771*(val1*val1) do j = G%jsc,G%jec ; do i = G%isc,G%iec do k = 1, nz - x = G%geoLonT(i,j)-x0 - y = G%geoLatT(i,j)-y0 + x = (G%geoLonT(i,j)-x0) * scale_pos + y = (G%geoLatT(i,j)-y0) * scale_pos val3 = exp(-val1*x) val4 = val2 * ( 2.0*val3 / (1.0 + (val3*val3)) )**2 - h(i,j,k) = (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j)) + h(i,j,k) = (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j)) * I_nz enddo enddo ; enddo end subroutine soliton_initialize_thickness -!> Initialization of u and v in the equatorial Rossby soliton test -subroutine soliton_initialize_velocity(u, v, G, GV, US) +!> Initialization of u and v in the equatorial Rossby soliton test, as described in section +!! 6.1 of Haidvogel and Beckman (1990) and in Boyd (1980, JPO) and Boyd (1985, JPO). +subroutine soliton_initialize_velocity(u, v, G, GV, US, param_file, just_read) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: u !< i-component of velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: v !< j-component of velocity [L T-1 ~> m s-1] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, intent(in) :: just_read !< If true, this call will only read + !! parameters without changing h. ! Local variables - real :: x, x0 ! Positions in the same units as geoLonT. - real :: y, y0 ! Positions in the same units as geoLatT. - real :: val1 ! A zonal decay scale in the inverse of the units of geoLonT. + real :: max_depth ! Maximum depth of the model bathymetry [Z ~> m] + real :: cg_max ! The external wave speed based on max_depth [L T-1 ~> m s-1] + real :: beta ! The meridional gradient of the Coriolis parameter [T-1 L-1 ~> s-1 m-1] + real :: L_eq ! The equatorial deformation radius used in nondimensionalizing this problem [L ~> m] + real :: scale_pos ! A conversion factor to nondimensionalize the axis units, usually [m-1] + real :: x0 ! Initial x-position of the soliton in the same units as geoLonT, often [m]. + real :: y0 ! Initial y-position of the soliton in the same units as geoLatT, often [m]. + real :: x, y ! Nondimensionalized positions [nondim] + real :: val1 ! A nondimensionlized zonal decay scale [nondim] real :: val2 ! An overall velocity amplitude [L T-1 ~> m s-1] real :: val3 ! A decay factor [nondim] real :: val4 ! The local velocity amplitude [L T-1 ~> m s-1] @@ -81,18 +133,40 @@ subroutine soliton_initialize_velocity(u, v, G, GV, US) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + if (.not.just_read) & + call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness") + + call get_param(param_file, mdl, "MAXIMUM_DEPTH", max_depth, & + units="m", default=-1.e9, scale=US%m_to_Z, do_not_log=.true.) + call get_param(param_file, mdl, "BETA", beta, & + "The northward gradient of the Coriolis parameter with the betaplane option.", & + units="m-1 s-1", default=0.0, scale=US%T_to_s*US%L_to_m, do_not_log=.true.) + + if (just_read) return ! All run-time parameters have been read, so return. + + if (max_depth <= 0.0) call MOM_error(FATAL, & + "soliton_initialization, soliton_initialize_velocity: "//& + "This module requires a positive value of MAXIMUM_DEPTH.") + if (abs(beta) <= 0.0) call MOM_error(FATAL, & + "soliton_initialization, soliton_initialize_velocity: "//& + "This module requires a non-zero value of BETA.") + + cg_max = sqrt(GV%g_Earth * max_depth) + L_eq = sqrt(cg_max / abs(beta)) + scale_pos = US%m_to_L / L_eq + x0 = 2.0*G%len_lon/3.0 y0 = 0.0 val1 = 0.395 - val2 = US%m_s_to_L_T * 0.771*(val1*val1) + val2 = cg_max * 0.771*(val1*val1) v(:,:,:) = 0.0 u(:,:,:) = 0.0 do j = G%jsc,G%jec ; do I = G%isc-1,G%iec+1 do k = 1, nz - x = 0.5*(G%geoLonT(i+1,j)+G%geoLonT(i,j))-x0 - y = 0.5*(G%geoLatT(i+1,j)+G%geoLatT(i,j))-y0 + x = (0.5*(G%geoLonT(i+1,j)+G%geoLonT(i,j))-x0) * scale_pos + y = (0.5*(G%geoLatT(i+1,j)+G%geoLatT(i,j))-y0) * scale_pos val3 = exp(-val1*x) val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2) u(I,j,k) = 0.25*val4*(6.0*y*y-9.0) * exp(-0.5*y*y)