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 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/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/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/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/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. 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/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 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 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)