From b32aea7bf3f9e2a774afa23d3386c88156cd1182 Mon Sep 17 00:00:00 2001 From: Daniel Sarmiento <42810219+dpsarmie@users.noreply.github.com> Date: Tue, 21 May 2024 13:46:30 -0400 Subject: [PATCH 01/12] Add end of run restart functionality to MOM6 (#133) * Update mom_cap.F90 to add end of run restart file functionality controlled by write_restart_at_endofrun configuration option in CMEPS. Author: Daniel Sarmient --- config_src/drivers/nuopc_cap/mom_cap.F90 | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 3574943918..d2dcf96067 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,12 @@ 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 +1644,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 +1785,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 +1815,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 From bfe3e4f03d2795b42d557a2ebe2713f6cee3bff3 Mon Sep 17 00:00:00 2001 From: jiandewang Date: Wed, 22 May 2024 15:51:27 -0400 Subject: [PATCH 02/12] fix line length too long issue in mom_cap.F90 --- config_src/drivers/nuopc_cap/mom_cap.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index d2dcf96067..3e3abba674 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -383,7 +383,8 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc) 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) + 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. From 7384dba37cb3ba82bb9eec204b7e8bd1c8031322 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 23 May 2024 16:14:22 -0400 Subject: [PATCH 03/12] Autoconf: Detect FMS IO implementation This patch allows autoconf to detect the version of FMS IO that is implemented in its library, and select the appropriate "infra" in the MOM6 source. The allows for removal of the --with-framework flag, and should prevent future FMS library/infra mismatches. --- .testing/Makefile | 1 - ac/configure.ac | 21 ++++++++------------- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index a1461fe7ab..2efead315d 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -286,7 +286,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 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 From c7d3268b65292e6cfcdca595cfbfb9d5f1d2f519 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 24 May 2024 12:53:36 -0400 Subject: [PATCH 04/12] .testing: Target build uses its own Makefile The build/target/MOM6 rule was modified so that build/target is now a symlink into build/target_codebase/.testing/build/symmetric. In other words, build/target now uses the entire build system of target_codebase (i.e. the reference codebase) rather than borrowing from bits and pieces of the current codebase. This should prevent potential false positives in future regression tests related to buildsystem changes. Several apparently unused blocks of rules and other content in .testing/Makefile have also been removed. --- .testing/Makefile | 60 ++++++++++------------------------------------- 1 file changed, 12 insertions(+), 48 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index 2efead315d..8ca3fe21eb 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -133,9 +133,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 +204,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 @@ -297,11 +266,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: @@ -333,27 +312,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) From 8ecb2603119b8c9d4dee9ee6b174e1864415d1cc Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 28 May 2024 14:48:54 -0400 Subject: [PATCH 05/12] CI: Fix FMS_COMMIT, remove FRAMEWORK This patch fixes FMS_COMMIT and FMS_URL so that they are properly exported to ``ac/deps/Makefile`` and used in the .testing FMS build. It also removes any references to the FRAMEWORK flag, which is no longer needed to identify the FMS API. --- .github/workflows/macos-regression.yml | 1 - .github/workflows/macos-stencil.yml | 1 - .testing/Makefile | 8 +++++--- .testing/README.rst | 8 +++++--- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/.github/workflows/macos-regression.yml b/.github/workflows/macos-regression.yml index 5e0380bcd2..d769e15131 100644 --- a/.github/workflows/macos-regression.yml +++ b/.github/workflows/macos-regression.yml @@ -11,7 +11,6 @@ jobs: CC: gcc FC: gfortran FMS_COMMIT: 2019.01.03 - FRAMEWORK: fms1 defaults: run: diff --git a/.github/workflows/macos-stencil.yml b/.github/workflows/macos-stencil.yml index e54912b607..6e77a5c4a6 100644 --- a/.github/workflows/macos-stencil.yml +++ b/.github/workflows/macos-stencil.yml @@ -11,7 +11,6 @@ jobs: CC: gcc FC: gfortran FMS_COMMIT: 2019.01.03 - FRAMEWORK: fms1 defaults: run: diff --git a/.testing/Makefile b/.testing/Makefile index 8ca3fe21eb..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 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 From f21ec031ecb5a78caaae5a6d153feba8abf5d907 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 24 May 2024 04:24:15 -0400 Subject: [PATCH 06/12] *Fix ALE_remap_scalar call with h_in_Z_units Corrected an argument to the call to ALE_remap_scalar() when h_in_Z_units is true in MOM_initialize_tracer_from_Z(), to avoid the problems documented in github.com/NOAA-GFDL/MOM6/issues/589. A comment was also added explaining the logic of what is going on in this fork of the code. This commit will change answers with some generic tracers that are initialized from a Z-space input file, restoring them to previous values that worked previously (before about Feb. 1, 2024 on dev/gfdl) in Boussinesq configurations without dimensional consistency testing, but in a new form that does pass the dimensional consistency testing for depths and thicknesses. All answers are bitwise identical in any cases that do not use generic tracers. --- src/initialization/MOM_tracer_initialization_from_Z.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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, From 705a38493fef60cd4d577315fe1d0ff9609a7d58 Mon Sep 17 00:00:00 2001 From: Uriel Ramirez Date: Mon, 1 Apr 2024 13:27:35 -0400 Subject: [PATCH 07/12] Add diag send complete calls to MOM --- config_src/infra/FMS1/MOM_diag_manager_infra.F90 | 6 ++++++ config_src/infra/FMS2/MOM_diag_manager_infra.F90 | 12 +++++++++++- src/framework/MOM_diag_mediator.F90 | 2 ++ 3 files changed, 19 insertions(+), 1 deletion(-) diff --git a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 index 18ccdaae67..cd9544d6e6 100644 --- a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 @@ -57,6 +57,7 @@ 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 MOM_diag_field_add_attribute public register_diag_field_infra public register_static_field_infra @@ -451,4 +452,9 @@ 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 +!! Needed for backwards compatibility, does nothing +subroutine diag_send_complete_infra () +end subroutine diag_send_complete_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..230ee79566 100644 --- a/config_src/infra/FMS2/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS2/MOM_diag_manager_infra.F90 @@ -14,13 +14,14 @@ 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 : 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 +58,7 @@ 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 MOM_diag_field_add_attribute public register_diag_field_infra public register_static_field_infra @@ -451,4 +453,12 @@ 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 + end module MOM_diag_manager_infra diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 5289d266df..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 @@ -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. From 7305528fd99bc09bda648693cd7130e4627de25c Mon Sep 17 00:00:00 2001 From: Uriel Ramirez Date: Tue, 14 May 2024 16:43:16 -0400 Subject: [PATCH 08/12] Add diag_manager_set_time_end calls to the solo drivers --- config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 | 3 +++ config_src/drivers/solo_driver/MOM_driver.F90 | 3 +++ config_src/infra/FMS1/MOM_diag_manager_infra.F90 | 9 +++++++-- config_src/infra/FMS2/MOM_diag_manager_infra.F90 | 9 +++++++++ 4 files changed, 22 insertions(+), 2 deletions(-) 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/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/infra/FMS1/MOM_diag_manager_infra.F90 b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 index cd9544d6e6..232986f480 100644 --- a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS1/MOM_diag_manager_infra.F90 @@ -58,6 +58,7 @@ module MOM_diag_manager_infra 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 @@ -452,9 +453,13 @@ 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 -!! Needed for backwards compatibility, does nothing +!> 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 230ee79566..f05baa4474 100644 --- a/config_src/infra/FMS2/MOM_diag_manager_infra.F90 +++ b/config_src/infra/FMS2/MOM_diag_manager_infra.F90 @@ -15,6 +15,7 @@ module MOM_diag_manager_infra 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 @@ -59,6 +60,7 @@ module MOM_diag_manager_infra 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 @@ -461,4 +463,11 @@ subroutine diag_send_complete_infra () 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 From d90ff6b09acf421695d73d4c9206c07152b2a5d9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 19 Apr 2024 15:46:33 -0400 Subject: [PATCH 09/12] (*)Refactor and document soliton_initialization Refactored the code in soliton_initialization.F90 to more accurately reflect the nondimensionalization that was applied in developing this test case. This change includes reading in the maximum depth and beta, and using them to calculate the equatorial deformation radius and the external gravity wave speed. The references to the papers describing the test case in the module were added to the doxygen comments describing the routines in this module. This commit also includes adding comments documenting the nature and units of all the internal variables in this module. There are two new arguments each (param_file and just_read) to soliton_initialize_thickness and soliton_initialize_velocity to accommodate these changes, bringing them into line with the interfaces for other similar user initialization routines, and MOM_initialize_state was changed accordingly. This change could change answers in general, but in the specific examples that use this code, both beta and the external wave speed are deliberately set to 1 in MKS units, so this commit does not change answers for that specific case. --- .../MOM_state_initialization.F90 | 5 +- src/user/soliton_initialization.F90 | 108 +++++++++++++++--- 2 files changed, 94 insertions(+), 19 deletions(-) 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/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) From 3c87d981dc4a895fa7528a1b805078085370c27a Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 11 May 2024 03:53:38 -0400 Subject: [PATCH 10/12] (*)Correct rescaling of KPP pseudo-salt flux Calculate the rescaled pseudo-salt flux when KPP is in use inside of pseudo_salt_tracer_column_physics(), rather than using fluxes%KPP_salt_flux. This commit also corrects for the fact that salinity and pseudo-salt have scaling that differs by a factor of US%S_to_ppt, which was not previously being taken into account. All answers are bitwise identical when no dimensional rescaling is being used, but they will change (and be corrected) when salinity is being rescaled. --- src/tracer/pseudo_salt_tracer.F90 | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/tracer/pseudo_salt_tracer.F90 b/src/tracer/pseudo_salt_tracer.F90 index 1185f70d22..ade36bad19 100644 --- a/src/tracer/pseudo_salt_tracer.F90 +++ b/src/tracer/pseudo_salt_tracer.F90 @@ -196,6 +196,8 @@ subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) ! Local variables + real :: net_salt_rate(SZI_(G),SZJ_(G)) ! Net salt flux into the ocean + ! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] real :: net_salt(SZI_(G),SZJ_(G)) ! Net salt flux into the ocean integrated over ! a timestep [ppt H ~> ppt m or ppt kg m-2] real :: htot(SZI_(G)) ! Total ocean depth [H ~> m or kg m-2] @@ -216,11 +218,27 @@ subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G call hchksum(CS%ps,"pseudo_salt pre pseudo-salt vertdiff", G%HI) endif + FluxRescaleDepth = max( GV%Angstrom_H, 1.e-30*GV%m_to_H ) + Ih_limit = 0.0 ; if (FluxRescaleDepth > 0.0) Ih_limit = 1.0 / FluxRescaleDepth + ! Compute KPP nonlocal term if necessary if (present(KPP_CSp)) then - if (associated(KPP_CSp) .and. present(nonLocalTrans)) & - call KPP_NonLocalTransport(KPP_CSp, G, GV, h_old, nonLocalTrans, fluxes%KPP_salt_flux(:,:), & + if (associated(KPP_CSp) .and. present(nonLocalTrans)) then + ! Determine the salt flux, including limiting for small total ocean depths. + net_salt_rate(:,:) = 0.0 + if (associated(fluxes%salt_flux)) then + do j=js,je + do i=is,ie ; htot(i) = h_old(i,j,1) ; enddo + do k=2,nz ; do i=is,ie ; htot(i) = htot(i) + h_old(i,j,k) ; enddo ; enddo + do i=is,ie + scale = 1.0 ; if ((Ih_limit > 0.0) .and. (htot(i)*Ih_limit < 1.0)) scale = htot(i)*Ih_limit + net_salt_rate(i,j) = (scale * (1000.0 * fluxes%salt_flux(i,j))) * GV%RZ_to_H + enddo + enddo + endif + call KPP_NonLocalTransport(KPP_CSp, G, GV, h_old, nonLocalTrans, net_salt_rate, & dt, CS%diag, CS%tr_ptr, CS%ps(:,:,:)) + endif endif ! This uses applyTracerBoundaryFluxesInOut, usually in ALE mode @@ -229,8 +247,6 @@ subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G ! Determine the time-integrated salt flux, including limiting for small total ocean depths. net_Salt(:,:) = 0.0 - FluxRescaleDepth = max( GV%Angstrom_H, 1.e-30*GV%m_to_H ) - Ih_limit = 0.0 ; if (FluxRescaleDepth > 0.0) Ih_limit = 1.0 / FluxRescaleDepth do j=js,je do i=is,ie ; htot(i) = h_old(i,j,1) ; enddo do k=2,nz ; do i=is,ie ; htot(i) = htot(i) + h_old(i,j,k) ; enddo ; enddo From 506d1864ca51452cce32b85a326d0f7fd8949677 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 11 May 2024 03:54:21 -0400 Subject: [PATCH 11/12] +Eliminate fluxes%KPP_salt_flux Eliminated the KPP_salt_flux element of the forcing type, which is no longer being used after the revisions to pseudo_salt_tracer_column_physics. Also moved 5 allocatable arrays (KPP_NLTheat, KPP_NLTscalar, KPP_buoy_flux, KPP_temp_flux and KPP_salt_flux) out of diabatic_CS, replacing them with ordinary arrays in diabatic_ALE(), diabatic_ALE_legacy() and layered_diabatic(). This change could reduce the high-water memory footprint of the model when KPP is in use. All answers are bitwise identical, but one element of a publicly visible type has been eliminated. --- src/core/MOM_forcing_type.F90 | 3 +- .../vertical/MOM_diabatic_driver.F90 | 149 ++++++++---------- 2 files changed, 70 insertions(+), 82 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 452161c6ca..72c67253ed 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -126,9 +126,8 @@ module MOM_forcing_type real, pointer, dimension(:,:) :: & netMassIn => NULL(), & !< Sum of water mass fluxes into the ocean integrated over a !! forcing timestep [H ~> m or kg m-2] - netMassOut => NULL(), & !< Net water mass flux out of the ocean integrated over a forcing timestep, + netMassOut => NULL() !< Net water mass flux out of the ocean integrated over a forcing timestep, !! with negative values for water leaving the ocean [H ~> m or kg m-2] - KPP_salt_flux => NULL() !< KPP effective salt flux [ppt m s-1] ! heat associated with water crossing ocean surface real, pointer, dimension(:,:) :: & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 081f065f3e..010f17b978 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -244,15 +244,6 @@ module MOM_diabatic_driver type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass type(group_pass_type) :: pass_Kv !< For group halo pass type(diag_grid_storage) :: diag_grids_prev!< Stores diagnostic grids at some previous point in the algorithm - ! Data arrays for communicating between components - !### Why are these arrays in this control structure, and not local variables in the various routines? - real, allocatable, dimension(:,:,:) :: KPP_NLTheat !< KPP non-local transport for heat [nondim] - real, allocatable, dimension(:,:,:) :: KPP_NLTscalar !< KPP non-local transport for scalars [nondim] - real, allocatable, dimension(:,:,:) :: KPP_buoy_flux !< KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3] - real, allocatable, dimension(:,:) :: KPP_temp_flux !< KPP effective temperature flux - !! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] - real, allocatable, dimension(:,:) :: KPP_salt_flux !< KPP effective salt flux - !! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] type(time_type), pointer :: Time !< Pointer to model time (needed for sponges) end type diabatic_CS @@ -558,11 +549,16 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_ePBL, & ! test array of diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + KPP_NLTheat, & ! KPP non-local transport for heat [nondim] + KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim] + KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3] Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] real, dimension(SZI_(G),SZJ_(G)) :: & U_star, & ! The friction velocity [Z T-1 ~> m s-1]. + KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] + KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] SkinBuoyFlux ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL logical, dimension(SZI_(G)) :: & @@ -677,11 +673,11 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! KPP needs the surface buoyancy flux but does not update state variables. ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes. - ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux - ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second) + ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux + ! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second) ! unlike other instances where the fluxes are integrated in time over a time-step. call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, & - CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) + KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux) ! Determine the friction velocity, perhaps using the evovling surface density. call find_ustar(fluxes, tv, U_star, G, GV, US) @@ -689,16 +685,16 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! The KPP scheme calculates boundary layer diffusivities and non-local transport. if ( associated(fluxes%lamult) ) then call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - U_star, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) + U_star, KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) else call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - U_star, CS%KPP_buoy_flux, Waves=Waves) + U_star, KPP_buoy_flux, Waves=Waves) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves) endif if (associated(Hml)) then @@ -708,7 +704,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) endif if (associated(visc%sfc_buoy_flx)) then - visc%sfc_buoy_flx(:,:) = CS%KPP_buoy_flux(:,:,1) + visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1) call pass_var(visc%sfc_buoy_flx, G%domain, halo=1) endif @@ -733,18 +729,18 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call MOM_thermovar_chksum("after KPP", tv, G, US) call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) - call hchksum(CS%KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, & + call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, & scale=US%C_to_degC*GV%H_to_m*US%s_to_T) - call hchksum(CS%KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, & + call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, & scale=US%S_to_ppt*GV%H_to_m*US%s_to_T) - call hchksum(CS%KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0) - call hchksum(CS%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0) + call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0) + call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0) endif ! Apply non-local transport of heat and salt ! Changes: tv%T, tv%S - call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, CS%KPP_NLTheat, CS%KPP_temp_flux, & + call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, KPP_NLTheat, KPP_temp_flux, & dt, tv%tr_T, tv%T, tv%C_p) - call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, CS%KPP_NLTscalar, CS%KPP_salt_flux, & + call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, KPP_NLTscalar, KPP_salt_flux, & dt, tv%tr_S, tv%S) call cpu_clock_end(id_clock_kpp) if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)") @@ -755,7 +751,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0) call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US) endif - if (.not.associated(fluxes%KPP_salt_flux)) fluxes%KPP_salt_flux => CS%KPP_salt_flux endif ! endif for KPP ! This is the "old" method for applying differential diffusion. @@ -1084,7 +1079,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call call_tracer_column_fns(h_orig, h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), fluxes, Hml, dt, & G, GV, US, tv, CS%optics, CS%tracer_flow_CSp, CS%debug, & KPP_CSp=CS%KPP_CSp, & - nonLocalTrans=CS%KPP_NLTscalar, & + nonLocalTrans=KPP_NLTscalar, & evap_CFL_limit=CS%evap_CFL_limit, & minimum_forcing_depth=CS%minimum_forcing_depth) @@ -1174,11 +1169,16 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int returned from set_diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_ePBL, & ! boundary layer or convective diapycnal diffusivities at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + KPP_NLTheat, & ! KPP non-local transport for heat [nondim] + KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim] + KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3] Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] real, dimension(SZI_(G),SZJ_(G)) :: & U_star, & ! The friction velocity [Z T-1 ~> m s-1]. + KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] + KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] SkinBuoyFlux ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL logical, dimension(SZI_(G)) :: & @@ -1298,11 +1298,11 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, ! KPP needs the surface buoyancy flux but does not update state variables. ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes. - ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux - ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second) + ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux + ! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second) ! unlike other instances where the fluxes are integrated in time over a time-step. call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, & - CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) + KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux) ! Determine the friction velocity, perhaps using the evovling surface density. call find_ustar(fluxes, tv, U_star, G, GV, US) @@ -1310,16 +1310,16 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, ! The KPP scheme calculates boundary layer diffusivities and non-local transport. if ( associated(fluxes%lamult) ) then call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - U_star, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) + U_star, KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) else call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - U_star, CS%KPP_buoy_flux, Waves=Waves) + U_star, KPP_buoy_flux, Waves=Waves) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves) endif if (associated(Hml)) then @@ -1329,7 +1329,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) endif if (associated(visc%sfc_buoy_flx)) then - visc%sfc_buoy_flx(:,:) = CS%KPP_buoy_flux(:,:,1) + visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1) call pass_var(visc%sfc_buoy_flx, G%domain, halo=1) endif @@ -1340,18 +1340,18 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call MOM_thermovar_chksum("after KPP", tv, G, US) call hchksum(Kd_heat, "after KPP Kd_heat", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) call hchksum(Kd_salt, "after KPP Kd_salt", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) - call hchksum(CS%KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, & + call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, & scale=US%C_to_degC*GV%H_to_m*US%s_to_T) - call hchksum(CS%KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, & + call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, & scale=US%S_to_ppt*GV%H_to_m*US%s_to_T) - call hchksum(CS%KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0) - call hchksum(CS%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0) + call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0) + call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0) endif ! Apply non-local transport of heat and salt ! Changes: tv%T, tv%S - call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, CS%KPP_NLTheat, CS%KPP_temp_flux, & + call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, KPP_NLTheat, KPP_temp_flux, & dt, tv%tr_T, tv%T, tv%C_p) - call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, CS%KPP_NLTscalar, CS%KPP_salt_flux, & + call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, KPP_NLTscalar, KPP_salt_flux, & dt, tv%tr_S, tv%S) call cpu_clock_end(id_clock_kpp) if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)") @@ -1362,7 +1362,6 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, US, haloshift=0) call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G, US) endif - if (.not.associated(fluxes%KPP_salt_flux)) fluxes%KPP_salt_flux => CS%KPP_salt_flux endif ! endif for KPP ! Calculate vertical mixing due to convection (computed via CVMix) @@ -1607,7 +1606,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call call_tracer_column_fns(h_orig, h, ent_s(:,:,1:nz), ent_s(:,:,2:nz+1), fluxes, Hml, dt, & G, GV, US, tv, CS%optics, CS%tracer_flow_CSp, CS%debug, & KPP_CSp=CS%KPP_CSp, & - nonLocalTrans=CS%KPP_NLTscalar, & + nonLocalTrans=KPP_NLTscalar, & evap_CFL_limit=CS%evap_CFL_limit, & minimum_forcing_depth=CS%minimum_forcing_depth) @@ -1693,6 +1692,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e saln_diag ! Diagnostic array of previous salinity [S ~> ppt] real, dimension(SZI_(G),SZJ_(G)) :: & U_star, & ! The friction velocity [Z T-1 ~> m s-1]. + KPP_temp_flux, & ! KPP effective temperature flux [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] + KPP_salt_flux, & ! KPP effective salt flux [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] Rcv_ml ! Coordinate density of mixed layer [R ~> kg m-3], used for applying sponges real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: & @@ -1710,6 +1711,9 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_extra_S , & ! The extra diffusivity of salinity due to double diffusion relative to ! Kd_int [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + KPP_NLTheat, & ! KPP non-local transport for heat [nondim] + KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim] + KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3] Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] Tadv_flx, & ! advective diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] @@ -1902,11 +1906,11 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e call cpu_clock_begin(id_clock_kpp) ! KPP needs the surface buoyancy flux but does not update state variables. ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes. - ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux - ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second) + ! Sets: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux + ! NOTE: KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux are returned as rates (i.e. stuff per second) ! unlike other instances where the fluxes are integrated in time over a time-step. call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, & - CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) + KPP_buoy_flux, KPP_temp_flux, KPP_salt_flux) ! The KPP scheme calculates boundary layer diffusivities and non-local transport. ! Set diffusivities for heat and salt separately @@ -1931,16 +1935,16 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e if ( associated(fluxes%lamult) ) then call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - U_star, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) + U_star, KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) else call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - U_star, CS%KPP_buoy_flux, Waves=Waves) + U_star, KPP_buoy_flux, Waves=Waves) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, KPP_NLTheat, KPP_NLTscalar, Waves=Waves) endif if (associated(Hml)) then @@ -1950,7 +1954,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e if (associated(visc%MLD)) visc%MLD(:,:) = Hml(:,:) endif if (associated(visc%sfc_buoy_flx)) then - visc%sfc_buoy_flx(:,:) = CS%KPP_buoy_flux(:,:,1) + visc%sfc_buoy_flx(:,:) = KPP_buoy_flux(:,:,1) call pass_var(visc%sfc_buoy_flx, G%domain, halo=1) endif @@ -1977,7 +1981,6 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e call hchksum(Kd_lay, "after KPP Kd_lay", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) call hchksum(Kd_Int, "after KPP Kd_Int", G%HI, haloshift=0, scale=GV%HZ_T_to_m2_s) endif - if (.not.associated(fluxes%KPP_salt_flux)) fluxes%KPP_salt_flux => CS%KPP_salt_flux endif ! endif for KPP ! Add vertical diff./visc. due to convection (computed via CVMix) @@ -1988,18 +1991,18 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e if (CS%useKPP) then call cpu_clock_begin(id_clock_kpp) if (CS%debug) then - call hchksum(CS%KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, & + call hchksum(KPP_temp_flux, "before KPP_applyNLT netHeat", G%HI, haloshift=0, & scale=US%C_to_degC*GV%H_to_m*US%s_to_T) - call hchksum(CS%KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, & + call hchksum(KPP_salt_flux, "before KPP_applyNLT netSalt", G%HI, haloshift=0, & scale=US%S_to_ppt*GV%H_to_m*US%s_to_T) - call hchksum(CS%KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0) - call hchksum(CS%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0) + call hchksum(KPP_NLTheat, "before KPP_applyNLT NLTheat", G%HI, haloshift=0) + call hchksum(KPP_NLTscalar, "before KPP_applyNLT NLTscalar", G%HI, haloshift=0) endif ! Apply non-local transport of heat and salt ! Changes: tv%T, tv%S - call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, CS%KPP_NLTheat, CS%KPP_temp_flux, & + call KPP_NonLocalTransport_temp(CS%KPP_CSp, G, GV, h, KPP_NLTheat, KPP_temp_flux, & dt, tv%tr_T, tv%T, tv%C_p) - call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, CS%KPP_NLTscalar, CS%KPP_salt_flux, & + call KPP_NonLocalTransport_saln(CS%KPP_CSp, G, GV, h, KPP_NLTscalar, KPP_salt_flux, & dt, tv%tr_S, tv%S) call cpu_clock_end(id_clock_kpp) if (showCallTree) call callTree_waypoint("done with KPP_applyNonLocalTransport (diabatic)") @@ -2401,7 +2404,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, US, tv, & CS%optics, CS%tracer_flow_CSp, CS%debug, & KPP_CSp=CS%KPP_CSp, & - nonLocalTrans=CS%KPP_NLTscalar) + nonLocalTrans=KPP_NLTscalar) elseif (CS%double_diffuse) then ! extra diffusivity for passive tracers @@ -2423,13 +2426,13 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, US, tv, & CS%optics, CS%tracer_flow_CSp, CS%debug, & KPP_CSp=CS%KPP_CSp, & - nonLocalTrans=CS%KPP_NLTscalar) + nonLocalTrans=KPP_NLTscalar) else call call_tracer_column_fns(hold, h, ea, eb, fluxes, Hml, dt, G, GV, US, tv, & CS%optics, CS%tracer_flow_CSp, CS%debug, & KPP_CSp=CS%KPP_CSp, & - nonLocalTrans=CS%KPP_NLTscalar) + nonLocalTrans=KPP_NLTscalar) endif ! (CS%mix_boundary_tracers) @@ -3315,14 +3318,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise. ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive CS%useKPP = KPP_init(param_file, G, GV, US, diag, Time, CS%KPP_CSp, passive=CS%KPPisPassive) - if (CS%useKPP) then - allocate(CS%KPP_NLTheat(isd:ied,jsd:jed,nz+1), source=0.0) - allocate(CS%KPP_NLTscalar(isd:ied,jsd:jed,nz+1), source=0.0) - allocate(CS%KPP_buoy_flux(isd:ied,jsd:jed,nz+1), source=0.0) - allocate(CS%KPP_temp_flux(isd:ied,jsd:jed), source=0.0) - allocate(CS%KPP_salt_flux(isd:ied,jsd:jed), source=0.0) - endif - ! Diagnostics for tendencies of temperature and salinity due to diabatic processes, ! available only for ALE algorithm. @@ -3616,14 +3611,8 @@ subroutine diabatic_driver_end(CS) if (CS%use_geothermal) & call geothermal_end(CS%geothermal) - if (CS%useKPP) then - deallocate( CS%KPP_buoy_flux ) - deallocate( CS%KPP_temp_flux ) - deallocate( CS%KPP_salt_flux ) - deallocate( CS%KPP_NLTheat ) - deallocate( CS%KPP_NLTscalar ) + if (CS%useKPP) & call KPP_end(CS%KPP_CSp) - endif ! GMM, the following is commented out because arrays in ! CS%diag_grids_prev are neither pointers or allocatables From c121215cbfbfc0944418abb4ea3490649f97c946 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 24 May 2024 05:23:22 -0400 Subject: [PATCH 12/12] Fix TEOS10 argument descriptions Corrected the descriptions of several conservative temperature and absolute salinity arguments in the comments for arguments to 4 TEOS10 equation of state routines, and eliminated the commented out code to mask out negative salinities in 8 TEOS10 routines. Only comments are changed, and all answers are bitwise identical. --- src/equation_of_state/MOM_EOS_TEOS10.F90 | 139 +++++++++-------------- 1 file changed, 53 insertions(+), 86 deletions(-) diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 index 3f138e20bb..6e4aaa762f 100644 --- a/src/equation_of_state/MOM_EOS_TEOS10.F90 +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -50,13 +50,6 @@ real elemental function density_elem_TEOS10(this, T, S, pressure) real, intent(in) :: S !< Absolute salinity [g kg-1]. real, intent(in) :: pressure !< pressure [Pa]. - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA -! if (S < -1.0e-10) then ! Can we assume safely that this is a missing value? -! density_elem_TEOS10 = 1000.0 -! else -! density_elem_TEOS10 = gsw_rho(S, T, pressure * Pa2db) -! endif - density_elem_TEOS10 = gsw_rho(S, T, pressure * Pa2db) end function density_elem_TEOS10 @@ -69,13 +62,6 @@ real elemental function density_anomaly_elem_TEOS10(this, T, S, pressure, rho_re real, intent(in) :: pressure !< pressure [Pa]. real, intent(in) :: rho_ref !< A reference density [kg m-3]. - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA -! if (S < -1.0e-10) then ! Can we assume safely that this is a missing value? -! density_elem_TEOS10 = 1000.0 -! else -! density_elem_TEOS10 = gsw_rho(S, T, pressure * Pa2db) -! endif - density_anomaly_elem_TEOS10 = gsw_rho(S, T, pressure * Pa2db) density_anomaly_elem_TEOS10 = density_anomaly_elem_TEOS10 - rho_ref @@ -88,13 +74,6 @@ real elemental function spec_vol_elem_TEOS10(this, T, S, pressure) real, intent(in) :: S !< Absolute salinity [g kg-1]. real, intent(in) :: pressure !< pressure [Pa]. - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA -! if (S < -1.0e-10) then ! Can we assume safely that this is a missing value? -! spec_vol_elem_TEOS10 = 0.001 -! else -! spec_vol_elem_TEOS10 = gsw_rho(S, T, pressure * Pa2db) -! endif - spec_vol_elem_TEOS10 = gsw_specvol(S, T, pressure * Pa2db) end function spec_vol_elem_TEOS10 @@ -107,13 +86,6 @@ real elemental function spec_vol_anomaly_elem_TEOS10(this, T, S, pressure, spv_r real, intent(in) :: pressure !< pressure [Pa]. real, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1]. - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA -! if (S < -1.0e-10) then ! Can we assume safely that this is a missing value? -! spec_vol_elem_TEOS10 = 0.001 -! else -! spec_vol_elem_TEOS10 = gsw_rho(S, T, pressure * Pa2db) -! endif - spec_vol_anomaly_elem_TEOS10 = gsw_specvol(S, T, pressure * Pa2db) - spv_ref end function spec_vol_anomaly_elem_TEOS10 @@ -122,28 +94,27 @@ end function spec_vol_anomaly_elem_TEOS10 !! temperature and absolute salinity, using the TEOS10 expressions. elemental subroutine calculate_density_derivs_elem_TEOS10(this, T, S, pressure, drho_dT, drho_dS) class(TEOS10_EOS), intent(in) :: this !< This EOS - real, intent(in) :: T !< Potential temperature relative to the surface [degC] - real, intent(in) :: S !< Salinity [PSU] + real, intent(in) :: T !< Conservative temperature [degC] + real, intent(in) :: S !< Absolute salinity [g kg-1] = [ppt] real, intent(in) :: pressure !< Pressure [Pa] - real, intent(out) :: drho_dT !< The partial derivative of density with potential + real, intent(out) :: drho_dT !< The partial derivative of density with conservative !! temperature [kg m-3 degC-1] real, intent(out) :: drho_dS !< The partial derivative of density with salinity, - !! in [kg m-3 PSU-1] + !! in [kg m-3 ppt-1] ! Local variables real :: zs ! Absolute salinity [g kg-1] real :: zt ! Conservative temperature [degC] real :: zp ! Pressure converted to decibars [dbar] - !Conversions - zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity - zt = T !gsw_ct_from_pt(S,T) !Convert potential temp to conservative temp - zp = pressure* Pa2db !Convert pressure from Pascal to decibar - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA - !if (S < -1.0e-10) then !Can we assume safely that this is a missing value? - ! drho_dT = 0.0 ; drho_dS = 0.0 - !else - call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS, drho_dct=drho_dT) - !endif + ! Conversions + zs = S + zt = T + zp = pressure * Pa2db ! Convert pressure from Pascal to decibar + ! The following conversions are unnecessary because the arguments are already the right variables. + ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity + ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp + + call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS, drho_dct=drho_dT) end subroutine calculate_density_derivs_elem_TEOS10 @@ -151,17 +122,17 @@ end subroutine calculate_density_derivs_elem_TEOS10 elemental subroutine calculate_density_second_derivs_elem_TEOS10(this, T, S, pressure, & drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP) class(TEOS10_EOS), intent(in) :: this !< This EOS - real, intent(in) :: T !< Potential temperature referenced to 0 dbar [degC] - real, intent(in) :: S !< Salinity [PSU] + real, intent(in) :: T !< Conservative temperature [degC] + real, intent(in) :: S !< Absolute salinity [g kg-1] = [ppt] real, intent(in) :: pressure !< Pressure [Pa] real, intent(inout) :: drho_ds_ds !< Partial derivative of beta with respect - !! to S [kg m-3 PSU-2] + !! to S [kg m-3 ppt-2] real, intent(inout) :: drho_ds_dt !< Partial derivative of beta with respect - !! to T [kg m-3 PSU-1 degC-1] + !! to T [kg m-3 ppt-1 degC-1] real, intent(inout) :: drho_dt_dt !< Partial derivative of alpha with respect !! to T [kg m-3 degC-2] real, intent(inout) :: drho_ds_dp !< Partial derivative of beta with respect - !! to pressure [kg m-3 PSU-1 Pa-1] = [s2 m-2 PSU-1] + !! to pressure [kg m-3 ppt-1 Pa-1] = [s2 m-2 ppt-1] real, intent(inout) :: drho_dt_dp !< Partial derivative of alpha with respect !! to pressure [kg m-3 degC-1 Pa-1] = [s2 m-2 degC-1] ! Local variables @@ -169,18 +140,16 @@ elemental subroutine calculate_density_second_derivs_elem_TEOS10(this, T, S, pre real :: zt ! Conservative temperature [degC] real :: zp ! Pressure converted to decibars [dbar] - !Conversions - zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity - zt = T !gsw_ct_from_pt(S,T) !Convert potential temp to conservative temp - zp = pressure* Pa2db !Convert pressure from Pascal to decibar - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA - !if (S < -1.0e-10) then !Can we assume safely that this is a missing value? - ! drho_dS_dS = 0.0 ; drho_dS_dT = 0.0 ; drho_dT_dT = 0.0 - ! drho_dS_dP = 0.0 ; drho_dT_dP = 0.0 - !else - call gsw_rho_second_derivatives(zs, zt, zp, rho_sa_sa=drho_dS_dS, rho_sa_ct=drho_dS_dT, & - rho_ct_ct=drho_dT_dT, rho_sa_p=drho_dS_dP, rho_ct_p=drho_dT_dP) - !endif + ! Conversions + zs = S + zt = T + zp = pressure * Pa2db ! Convert pressure from Pascal to decibar + ! The following conversions are unnecessary because the arguments are already the right variables. + ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity + ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp + + call gsw_rho_second_derivatives(zs, zt, zp, rho_sa_sa=drho_dS_dS, rho_sa_ct=drho_dS_dT, & + rho_ct_ct=drho_dT_dT, rho_sa_p=drho_dS_dP, rho_ct_p=drho_dT_dP) end subroutine calculate_density_second_derivs_elem_TEOS10 @@ -188,28 +157,27 @@ end subroutine calculate_density_second_derivs_elem_TEOS10 !! temperature and absolute salinity, using the TEOS10 expressions. elemental subroutine calculate_specvol_derivs_elem_TEOS10(this, T, S, pressure, dSV_dT, dSV_dS) class(TEOS10_EOS), intent(in) :: this !< This EOS - real, intent(in) :: T !< Potential temperature [degC] - real, intent(in) :: S !< Salinity [PSU] + real, intent(in) :: T !< Conservative temperature [degC] + real, intent(in) :: S !< Absolute salinity [g kg-1] = [ppt] real, intent(in) :: pressure !< Pressure [Pa] real, intent(inout) :: dSV_dT !< The partial derivative of specific volume with - !! potential temperature [m3 kg-1 degC-1] + !! conservative temperature [m3 kg-1 degC-1] real, intent(inout) :: dSV_dS !< The partial derivative of specific volume with - !! salinity [m3 kg-1 PSU-1] + !! absolute salinity [m3 kg-1 ppt-1] ! Local variables real :: zs ! Absolute salinity [g kg-1] real :: zt ! Conservative temperature [degC] real :: zp ! Pressure converted to decibars [dbar] - !Conversions - zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity - zt = T !gsw_ct_from_pt(S,T) !Convert potential temp to conservative temp - zp = pressure* Pa2db !Convert pressure from Pascal to decibar - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA - !if (S < -1.0e-10) then !Can we assume safely that this is a missing value? - ! dSV_dT = 0.0 ; dSV_dS = 0.0 - !else - call gsw_specvol_first_derivatives(zs,zt,zp, v_sa=dSV_dS, v_ct=dSV_dT) - !endif + ! Conversions + zs = S + zt = T + zp = pressure * Pa2db ! Convert pressure from Pascal to decibar + ! The following conversions are unnecessary because the arguments are already the right variables. + ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity + ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp + + call gsw_specvol_first_derivatives(zs, zt, zp, v_sa=dSV_dS, v_ct=dSV_dT) end subroutine calculate_specvol_derivs_elem_TEOS10 @@ -220,8 +188,8 @@ end subroutine calculate_specvol_derivs_elem_TEOS10 !! subroutines from TEOS10 website elemental subroutine calculate_compress_elem_TEOS10(this, T, S, pressure, rho, drho_dp) class(TEOS10_EOS), intent(in) :: this !< This EOS - real, intent(in) :: T !< Potential temperature relative to the surface [degC] - real, intent(in) :: S !< Salinity [PSU] + real, intent(in) :: T !< Conservative temperature [degC] + real, intent(in) :: S !< Absolute salinity [g kg-1] real, intent(in) :: pressure !< Pressure [Pa] real, intent(out) :: rho !< In situ density [kg m-3] real, intent(out) :: drho_dp !< The partial derivative of density with pressure @@ -233,17 +201,16 @@ elemental subroutine calculate_compress_elem_TEOS10(this, T, S, pressure, rho, d real :: zt ! Conservative temperature [degC] real :: zp ! Pressure converted to decibars [dbar] - !Conversions - zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity - zt = T !gsw_ct_from_pt(S,T) !Convert potential temp to conservative temp - zp = pressure* Pa2db !Convert pressure from Pascal to decibar - !!! #### This code originally had this "masking" line. The answer to the question below is "no" -AJA - !if (S < -1.0e-10) then !Can we assume safely that this is a missing value? - ! rho = 1000.0 ; drho_dp = 0.0 - !else - rho = gsw_rho(zs,zt,zp) - call gsw_rho_first_derivatives(zs,zt,zp, drho_dp=drho_dp) - !endif + ! Conversions + zs = S + zt = T + zp = pressure * Pa2db ! Convert pressure from Pascal to decibar + ! The following conversions are unnecessary because the arguments are already the right variables. + ! zs = gsw_sr_from_sp(S) ! Uncomment to convert practical salinity to absolute salinity + ! zt = gsw_ct_from_pt(S,T) ! Uncomment to convert potential temp to conservative temp + + rho = gsw_rho(zs, zt, zp) + call gsw_rho_first_derivatives(zs, zt, zp, drho_dp=drho_dp) end subroutine calculate_compress_elem_TEOS10