diff --git a/.github/actions/testing-setup/action.yml b/.github/actions/testing-setup/action.yml index c6fae4ad58..c47270af0d 100644 --- a/.github/actions/testing-setup/action.yml +++ b/.github/actions/testing-setup/action.yml @@ -1,5 +1,14 @@ name: 'Build-.testing-prerequisites' description: 'Build pre-requisites for .testing including FMS and a symmetric MOM6 executable' +inputs: + build_symmetric: + description: 'If true, will build the symmetric MOM6 executable' + required: false + default: 'true' + install_python: + description: 'If true, will install the local python env needed for .testing' + required: false + default: 'true' runs: using: 'composite' steps: @@ -51,7 +60,7 @@ runs: run: | echo "::group::Compile MOM6 in symmetric memory mode" cd .testing - make build/symmetric/MOM6 -j + test ${{ inputs.build_symmetric }} == true && make build/symmetric/MOM6 -j echo "::endgroup::" - name: Install local python venv for generating input data @@ -59,7 +68,7 @@ runs: run: | echo "::group::Create local python env for input data generation" cd .testing - make work/local-env + test ${{ inputs.install_python }} == true && make work/local-env echo "::endgroup::" - name: Set flags diff --git a/.github/workflows/coupled-api.yml b/.github/workflows/coupled-api.yml new file mode 100644 index 0000000000..86d7262548 --- /dev/null +++ b/.github/workflows/coupled-api.yml @@ -0,0 +1,33 @@ +name: API for coupled drivers + +on: [push, pull_request] + +jobs: + test-top-api: + + runs-on: ubuntu-latest + defaults: + run: + working-directory: .testing + + steps: + - uses: actions/checkout@v2 + with: + submodules: recursive + + - uses: ./.github/actions/testing-setup + with: + build_symmetric: 'false' + install_python: 'false' + + - name: Compile MOM6 for the GFDL coupled driver + shell: bash + run: make check_mom6_api_coupled -j + + - name: Compile MOM6 for the NUOPC driver + shell: bash + run: make check_mom6_api_nuopc -j + + - name: Compile MOM6 for the MCT driver + shell: bash + run: make check_mom6_api_mct -j diff --git a/.testing/Makefile b/.testing/Makefile index 33d3ea4c4e..2806d54130 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -226,7 +226,9 @@ build/asymmetric/Makefile: MOM_ENV=$(PATH_FMS) $(ASYMMETRIC_FCFLAGS) $(MOM_LDFLA build/repro/Makefile: MOM_ENV=$(PATH_FMS) $(REPRO_FCFLAGS) $(MOM_LDFLAGS) build/openmp/Makefile: MOM_ENV=$(PATH_FMS) $(OPENMP_FCFLAGS) $(MOM_LDFLAGS) build/target/Makefile: MOM_ENV=$(PATH_FMS) $(TARGET_FCFLAGS) $(MOM_LDFLAGS) - +build/coupled/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) +build/nuopc/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) +build/mct/Makefile: MOM_ENV=$(PATH_FMS) $(SYMMETRIC_FCFLAGS) $(SYMMETRIC_LDFLAGS) # Configure script flags build/symmetric/Makefile: MOM_ACFLAGS= @@ -234,7 +236,9 @@ build/asymmetric/Makefile: MOM_ACFLAGS=--enable-asymmetric build/repro/Makefile: MOM_ACFLAGS= build/openmp/Makefile: MOM_ACFLAGS=--enable-openmp build/target/Makefile: MOM_ACFLAGS= - +build/coupled/Makefile: MOM_ACFLAGS=--with-driver=coupled_driver +build/nuopc/Makefile: MOM_ACFLAGS=--with-driver=nuopc_driver +build/mct/Makefile: MOM_ACFLAGS=--with-driver=mct_driver # Fetch regression target source code build/target/Makefile: | $(TARGET_CODEBASE) @@ -328,6 +332,25 @@ $(DEPS)/Makefile: ../ac/deps/Makefile mkdir -p $(@D) cp $< $@ +#--- +# The following block does a non-library build of a coupled driver interface to MOM, along with everything below it. +# This simply checks that we have not broken the ability to compile. This is not a means to build a complete coupled executable. +# Todo: +# - avoid re-building FMS and MOM6 src by re-using existing object/mod files +# - use autoconf rather than mkmf templates +MK_TEMPLATE ?= ../../$(DEPS)/mkmf/templates/ncrc-gnu.mk +# NUOPC driver +build/nuopc/mom_ocean_model_nuopc.o: build/nuopc/Makefile + cd $(@D) && make $(@F) +check_mom6_api_nuopc: build/nuopc/mom_ocean_model_nuopc.o +# GFDL coupled driver +build/coupled/ocean_model_MOM.o: build/coupled/Makefile + cd $(@D) && make $(@F) +check_mom6_api_coupled: build/coupled/ocean_model_MOM.o +# MCT driver +build/mct/mom_ocean_model_mct.o: build/mct/Makefile + cd $(@D) && make $(@F) +check_mom6_api_mct: build/mct/mom_ocean_model_mct.o #--- # Python preprocessing diff --git a/ac/configure.ac b/ac/configure.ac index ad8ed83603..487230beb8 100644 --- a/ac/configure.ac +++ b/ac/configure.ac @@ -48,6 +48,12 @@ AC_ARG_ENABLE([asymmetric], AS_IF([test "$enable_asymmetric" = yes], [MEM_LAYOUT=${srcdir}/config_src/dynamic]) +# Default to solo_driver +DRIVER_DIR=${srcdir}/config_src/solo_driver +AC_ARG_WITH([driver], + AS_HELP_STRING([--with-driver=coupled_driver|solo_driver], [Select directory for driver source code])) +AS_IF([test "x$with_driver" != "x"], + [DRIVER_DIR=${srcdir}/config_src/${with_driver}]) # TODO: Rather than point to a pre-configured header file, autoconf could be # used to configure a header based on a template. @@ -210,10 +216,10 @@ AS_IF([test -z "$MKMF"], [ AC_CONFIG_COMMANDS([path_names], [list_paths -l \ ${srcdir}/src \ - ${srcdir}/config_src/solo_driver \ ${srcdir}/config_src/ext* \ + ${DRIVER_DIR} \ ${MEM_LAYOUT} -], [MEM_LAYOUT=$MEM_LAYOUT]) +], [MEM_LAYOUT=$MEM_LAYOUT DRIVER_DIR=$DRIVER_DIR]) AC_CONFIG_COMMANDS([Makefile.mkmf], diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index be960accd6..67f6643a42 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -5,7 +5,7 @@ module MOM_surface_forcing_gfdl !#CTRL# use MOM_controlled_forcing, only : apply_ctrl_forcing, register_ctrl_forcing_restarts !#CTRL# use MOM_controlled_forcing, only : controlled_forcing_init, controlled_forcing_end !#CTRL# use MOM_controlled_forcing, only : ctrl_forcing_CS -use MOM_coms, only : reproducing_sum +use MOM_coms, only : reproducing_sum, field_chksum use MOM_constants, only : hlv, hlf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT @@ -21,6 +21,8 @@ module MOM_surface_forcing_gfdl use MOM_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : init_external_field, time_interp_extern +use MOM_interpolate, only : time_interp_external_init use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS use MOM_restart, only : restart_init_end, save_restart, restore_state @@ -36,9 +38,6 @@ module MOM_surface_forcing_gfdl use coupler_types_mod, only : coupler_type_copy_data use data_override_mod, only : data_override_init, data_override use fms_mod, only : stdout -use mpp_mod, only : mpp_chksum -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init implicit none ; private @@ -350,7 +349,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (CS%restore_salt) then - call time_interp_external(CS%id_srestore,Time,data_restore) + call time_interp_extern(CS%id_srestore, Time, data_restore) ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) open_ocn_mask(:,:) = 1.0 if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice @@ -407,7 +406,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (CS%restore_temp) then - call time_interp_external(CS%id_trestore,Time,data_restore) + call time_interp_extern(CS%id_trestore, Time, data_restore) do j=js,je ; do i=is,ie delta_sst = data_restore(i,j)- sfc_state%SST(i,j) delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) @@ -1486,7 +1485,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) enddo ; enddo endif - call time_interp_external_init + call time_interp_external_init() ! Optionally read a x-y gustiness field in place of a global constant. call get_param(param_file, mdl, "READ_GUST_2D", CS%read_gust_2d, & @@ -1632,27 +1631,27 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) outunit = stdout() write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep - write(outunit,100) 'iobt%u_flux ', mpp_chksum( iobt%u_flux ) - write(outunit,100) 'iobt%v_flux ', mpp_chksum( iobt%v_flux ) - write(outunit,100) 'iobt%t_flux ', mpp_chksum( iobt%t_flux ) - write(outunit,100) 'iobt%q_flux ', mpp_chksum( iobt%q_flux ) - write(outunit,100) 'iobt%salt_flux ', mpp_chksum( iobt%salt_flux ) - write(outunit,100) 'iobt%lw_flux ', mpp_chksum( iobt%lw_flux ) - write(outunit,100) 'iobt%sw_flux_vis_dir', mpp_chksum( iobt%sw_flux_vis_dir) - write(outunit,100) 'iobt%sw_flux_vis_dif', mpp_chksum( iobt%sw_flux_vis_dif) - write(outunit,100) 'iobt%sw_flux_nir_dir', mpp_chksum( iobt%sw_flux_nir_dir) - write(outunit,100) 'iobt%sw_flux_nir_dif', mpp_chksum( iobt%sw_flux_nir_dif) - write(outunit,100) 'iobt%lprec ', mpp_chksum( iobt%lprec ) - write(outunit,100) 'iobt%fprec ', mpp_chksum( iobt%fprec ) - write(outunit,100) 'iobt%runoff ', mpp_chksum( iobt%runoff ) - write(outunit,100) 'iobt%calving ', mpp_chksum( iobt%calving ) - write(outunit,100) 'iobt%p ', mpp_chksum( iobt%p ) + write(outunit,100) 'iobt%u_flux ', field_chksum( iobt%u_flux ) + write(outunit,100) 'iobt%v_flux ', field_chksum( iobt%v_flux ) + write(outunit,100) 'iobt%t_flux ', field_chksum( iobt%t_flux ) + write(outunit,100) 'iobt%q_flux ', field_chksum( iobt%q_flux ) + write(outunit,100) 'iobt%salt_flux ', field_chksum( iobt%salt_flux ) + write(outunit,100) 'iobt%lw_flux ', field_chksum( iobt%lw_flux ) + write(outunit,100) 'iobt%sw_flux_vis_dir', field_chksum( iobt%sw_flux_vis_dir) + write(outunit,100) 'iobt%sw_flux_vis_dif', field_chksum( iobt%sw_flux_vis_dif) + write(outunit,100) 'iobt%sw_flux_nir_dir', field_chksum( iobt%sw_flux_nir_dir) + write(outunit,100) 'iobt%sw_flux_nir_dif', field_chksum( iobt%sw_flux_nir_dif) + write(outunit,100) 'iobt%lprec ', field_chksum( iobt%lprec ) + write(outunit,100) 'iobt%fprec ', field_chksum( iobt%fprec ) + write(outunit,100) 'iobt%runoff ', field_chksum( iobt%runoff ) + write(outunit,100) 'iobt%calving ', field_chksum( iobt%calving ) + write(outunit,100) 'iobt%p ', field_chksum( iobt%p ) if (associated(iobt%ustar_berg)) & - write(outunit,100) 'iobt%ustar_berg ', mpp_chksum( iobt%ustar_berg ) + write(outunit,100) 'iobt%ustar_berg ', field_chksum( iobt%ustar_berg ) if (associated(iobt%area_berg)) & - write(outunit,100) 'iobt%area_berg ', mpp_chksum( iobt%area_berg ) + write(outunit,100) 'iobt%area_berg ', field_chksum( iobt%area_berg ) if (associated(iobt%mass_berg)) & - write(outunit,100) 'iobt%mass_berg ', mpp_chksum( iobt%mass_berg ) + write(outunit,100) 'iobt%mass_berg ', field_chksum( iobt%mass_berg ) 100 FORMAT(" CHECKSUM::",A20," = ",Z20) call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index b429da649b..12f803a970 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -15,6 +15,7 @@ module ocean_model_mod use MOM, only : extract_surface_state, allocate_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized use MOM, only : get_ocean_stocks, step_offline +use MOM_coms, only : field_chksum use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end @@ -22,6 +23,7 @@ module ocean_model_mod use MOM_domains, only : TO_ALL, Omit_Corners use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave +use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type use MOM_forcing_type, only : forcing, mech_forcing, allocate_forcing_type use MOM_forcing_type, only : fluxes_accumulate, get_net_mass_forcing @@ -48,6 +50,8 @@ module ocean_model_mod use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart +use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init +use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data @@ -56,10 +60,6 @@ module ocean_model_mod use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux use fms_mod, only : stdout -use mpp_mod, only : mpp_chksum -use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct -use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init -use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves #include @@ -1130,13 +1130,13 @@ subroutine ocean_public_type_chksum(id, timestep, ocn) outunit = stdout() write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep - write(outunit,100) 'ocean%t_surf ',mpp_chksum(ocn%t_surf ) - write(outunit,100) 'ocean%s_surf ',mpp_chksum(ocn%s_surf ) - write(outunit,100) 'ocean%u_surf ',mpp_chksum(ocn%u_surf ) - write(outunit,100) 'ocean%v_surf ',mpp_chksum(ocn%v_surf ) - write(outunit,100) 'ocean%sea_lev ',mpp_chksum(ocn%sea_lev) - write(outunit,100) 'ocean%frazil ',mpp_chksum(ocn%frazil ) - write(outunit,100) 'ocean%melt_potential ',mpp_chksum(ocn%melt_potential) + write(outunit,100) 'ocean%t_surf ', field_chksum(ocn%t_surf ) + write(outunit,100) 'ocean%s_surf ', field_chksum(ocn%s_surf ) + write(outunit,100) 'ocean%u_surf ', field_chksum(ocn%u_surf ) + write(outunit,100) 'ocean%v_surf ', field_chksum(ocn%v_surf ) + write(outunit,100) 'ocean%sea_lev ', field_chksum(ocn%sea_lev) + write(outunit,100) 'ocean%frazil ', field_chksum(ocn%frazil ) + write(outunit,100) 'ocean%melt_potential ', field_chksum(ocn%melt_potential) call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%') 100 FORMAT(" CHECKSUM::",A20," = ",Z20) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 9726aa1281..c9383a4287 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -32,6 +32,7 @@ program MOM_main use MOM, only : extract_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized use MOM, only : step_offline + use MOM_coms, only : Set_PElist use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_error_handler, only : MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint @@ -41,6 +42,7 @@ program MOM_main use MOM_forcing_type, only : mech_forcing_diags, MOM_forcing_chksum, MOM_mech_forcing_chksum use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type + use MOM_interpolate, only : time_interp_external_init use MOM_io, only : file_exists, open_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE @@ -64,8 +66,6 @@ program MOM_main use MOM_get_input, only : get_MOM_input use ensemble_manager_mod, only : ensemble_manager_init, get_ensemble_size use ensemble_manager_mod, only : ensemble_pelist_setup - use mpp_mod, only : set_current_pelist => mpp_set_current_pelist - use time_interp_external_mod, only : time_interp_external_init use fms_affinity_mod, only : fms_affinity_init, fms_affinity_set,fms_affinity_get use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS @@ -229,7 +229,7 @@ program MOM_main allocate(ocean_pelist(nPEs_per)) call ensemble_pelist_setup(.true., 0, nPEs_per, 0, 0, atm_pelist, ocean_pelist, & land_pelist, ice_pelist) - call set_current_pelist(ocean_pelist) + call Set_PElist(ocean_pelist) deallocate(ocean_pelist) endif @@ -286,17 +286,17 @@ program MOM_main if (sum(date_init) > 0) then - Start_time = set_date(date_init(1),date_init(2), date_init(3), & - date_init(4),date_init(5),date_init(6)) + Start_time = set_date(date_init(1), date_init(2), date_init(3), & + date_init(4), date_init(5), date_init(6)) else Start_time = real_to_time(0.0) endif - call time_interp_external_init + call time_interp_external_init() if (sum(date) >= 0) then ! In this case, the segment starts at a time fixed by ocean_solo.res - segment_start_time = set_date(date(1),date(2),date(3),date(4),date(5),date(6)) + segment_start_time = set_date(date(1), date(2), date(3), date(4), date(5), date(6)) Time = segment_start_time else ! In this case, the segment starts at a time read from the MOM restart file diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 0673d7ca5b..c06cbfeb11 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -24,8 +24,7 @@ module MOM_open_boundary use MOM_tidal_forcing, only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) use MOM_tracer_registry, only : tracer_type, tracer_registry_type, tracer_name_lookup -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init +use MOM_interpolate, only : init_external_field, time_interp_extern, time_interp_external_init use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping use MOM_regridding, only : regridding_CS @@ -3897,7 +3896,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! TODO: Since we conditionally rotate a subset of tmp_buffer_in after ! reading the value, it is currently not possible to use the rotated - ! implementation of time_interp_external. + ! implementation of time_interp_extern. ! For now, we must explicitly allocate and rotate this array. if (turns /= 0) then if (modulo(turns, 2) /= 0) then @@ -3909,7 +3908,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) tmp_buffer_in => tmp_buffer endif - call time_interp_external(segment%field(m)%fid,Time, tmp_buffer_in) + call time_interp_extern(segment%field(m)%fid,Time, tmp_buffer_in) ! NOTE: Rotation of face-points require that we skip the final value if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. @@ -3976,7 +3975,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! no dz for tidal variables if (segment%field(m)%nk_src > 1 .and.& (index(segment%field(m)%name, 'phase') .le. 0 .and. index(segment%field(m)%name, 'amp') .le. 0)) then - call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer_in) + call time_interp_extern(segment%field(m)%fid_dz,Time, tmp_buffer_in) if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. if (segment%is_E_or_W & diff --git a/src/diagnostics/MOM_obsolete_diagnostics.F90 b/src/diagnostics/MOM_obsolete_diagnostics.F90 index e30749984d..bba8379bbb 100644 --- a/src/diagnostics/MOM_obsolete_diagnostics.F90 +++ b/src/diagnostics/MOM_obsolete_diagnostics.F90 @@ -4,10 +4,10 @@ module MOM_obsolete_diagnostics ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_diag_manager, only : register_static_field_fms +use MOM_diag_mediator, only : diag_ctrl use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : param_file_type, log_version, get_param -use MOM_diag_mediator, only : diag_ctrl -use diag_manager_mod, only : register_static_field_fms=>register_static_field implicit none ; private diff --git a/src/framework/MOM_coms.F90 b/src/framework/MOM_coms.F90 index a84540a33f..dec59a717d 100644 --- a/src/framework/MOM_coms.F90 +++ b/src/framework/MOM_coms.F90 @@ -5,25 +5,19 @@ module MOM_coms ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING -use fms_mod, only : fms_end, MOM_infra_init => fms_init -use memutils_mod, only : print_memuse_stats -use mpp_mod, only : PE_here => mpp_pe, root_PE => mpp_root_pe, num_PEs => mpp_npes -use mpp_mod, only : Set_PElist => mpp_set_current_pelist, Get_PElist => mpp_get_current_pelist -use mpp_mod, only : broadcast => mpp_broadcast -use mpp_domains_mod, only : broadcast_domain => mpp_broadcast_domain -use mpp_mod, only : set_rootPE => mpp_set_root_pe -use mpp_mod, only : field_chksum => mpp_chksum -use mpp_mod, only : sum_across_PEs => mpp_sum, max_across_PEs => mpp_max, min_across_PEs => mpp_min +use MOM_coms_wrapper, only : PE_here, root_PE, num_PEs, set_rootPE, Set_PElist, Get_PElist +use MOM_coms_wrapper, only : broadcast, field_chksum, MOM_infra_init, MOM_infra_end +use MOM_coms_wrapper, only : sum_across_PEs, max_across_PEs, min_across_PEs implicit none ; private public :: PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end -public :: broadcast, broadcast_domain, sum_across_PEs, min_across_PEs, max_across_PEs, field_chksum +public :: broadcast, sum_across_PEs, min_across_PEs, max_across_PEs, field_chksum +public :: set_PElist, Get_PElist, Set_rootPE public :: reproducing_sum, reproducing_sum_EFP, EFP_sum_across_PEs, EFP_list_sum_across_PEs public :: EFP_plus, EFP_minus, EFP_to_real, real_to_EFP, EFP_real_diff public :: operator(+), operator(-), assignment(=) public :: query_EFP_overflow_error, reset_EFP_overflow_error -public :: Set_PElist, Get_PElist, Set_rootPE ! This module provides interfaces to the non-domain-oriented communication subroutines. @@ -883,12 +877,4 @@ subroutine EFP_val_sum_across_PEs(EFP, error) end subroutine EFP_val_sum_across_PEs - -!> This subroutine carries out all of the calls required to close out the infrastructure cleanly. -!! This should only be called in ocean-only runs, as the coupler takes care of this in coupled runs. -subroutine MOM_infra_end - call print_memuse_stats( 'Memory HiWaterMark', always=.TRUE. ) - call fms_end -end subroutine MOM_infra_end - end module MOM_coms diff --git a/src/framework/MOM_coms_wrapper.F90 b/src/framework/MOM_coms_wrapper.F90 new file mode 100644 index 0000000000..034381c9f3 --- /dev/null +++ b/src/framework/MOM_coms_wrapper.F90 @@ -0,0 +1,161 @@ +!> Thin interfaces to non-domain-oriented mpp communication subroutines +module MOM_coms_wrapper + +! This file is part of MOM6. See LICENSE.md for the license. + +use fms_mod, only : fms_end, MOM_infra_init => fms_init +use memutils_mod, only : print_memuse_stats +use mpp_mod, only : PE_here => mpp_pe, root_PE => mpp_root_pe, num_PEs => mpp_npes +use mpp_mod, only : set_rootPE => mpp_set_root_pe +use mpp_mod, only : Set_PElist => mpp_set_current_pelist, Get_PElist => mpp_get_current_pelist +use mpp_mod, only : mpp_broadcast, mpp_sync, mpp_sync_self, field_chksum => mpp_chksum +use mpp_mod, only : sum_across_PEs => mpp_sum, max_across_PEs => mpp_max, min_across_PEs => mpp_min + +implicit none ; private + +public :: PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end, Set_PElist, Get_PElist +public :: set_rootPE, broadcast, sum_across_PEs, min_across_PEs, max_across_PEs, field_chksum + +! This module provides interfaces to the non-domain-oriented communication subroutines. + +!> Communicate an array, string or scalar from one PE to others +interface broadcast + module procedure broadcast_char, broadcast_int0D, broadcast_int1D + module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D +end interface broadcast + +contains + +!> Communicate a 1-D array of character strings from one PE to others +subroutine broadcast_char(dat, length, from_PE, PElist, blocking) + character(len=*), intent(inout) :: dat(:) !< The data to communicate and destination + integer, intent(in) :: length !< The length of each string + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_char + +!> Communicate an integer from one PE to others +subroutine broadcast_int0D(dat, from_PE, PElist, blocking) + integer, intent(inout) :: dat !< The data to communicate and destination + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_int0D + +!> Communicate a 1-D array of integers from one PE to others +subroutine broadcast_int1D(dat, length, from_PE, PElist, blocking) + integer, dimension(:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_int1D + +!> Communicate a real number from one PE to others +subroutine broadcast_real0D(dat, from_PE, PElist, blocking) + real, intent(inout) :: dat !< The data to communicate and destination + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real0D + +!> Communicate a 1-D array of reals from one PE to others +subroutine broadcast_real1D(dat, length, from_PE, PElist, blocking) + real, dimension(:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real1D + +!> Communicate a 2-D array of reals from one PE to others +subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking) + real, dimension(:,:), intent(inout) :: dat !< The data to communicate and destination + integer, intent(in) :: length !< The total number of data elements + integer, optional, intent(in) :: from_PE !< The source PE, by default the root PE + integer, optional, intent(in) :: PElist(:) !< The list of participating PEs, by default the + !! active PE set as previously set via Set_PElist. + logical, optional, intent(in) :: blocking !< If true, barriers are added around the call + + integer :: src_PE ! The processor that is sending the data + logical :: do_block ! If true add synchronizing barriers + + do_block = .false. ; if (present(blocking)) do_block = blocking + if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif + + if (do_block) call mpp_sync(PElist) + call mpp_broadcast(dat, length, src_PE, PElist) + if (do_block) call mpp_sync_self(PElist) + +end subroutine broadcast_real2D + + +!> This subroutine carries out all of the calls required to close out the infrastructure cleanly. +!! This should only be called in ocean-only runs, as the coupler takes care of this in coupled runs. +subroutine MOM_infra_end + call print_memuse_stats( 'Memory HiWaterMark', always=.TRUE. ) + call fms_end() +end subroutine MOM_infra_end + +end module MOM_coms_wrapper diff --git a/src/framework/MOM_diag_manager_wrapper.F90 b/src/framework/MOM_diag_manager.F90 similarity index 88% rename from src/framework/MOM_diag_manager_wrapper.F90 rename to src/framework/MOM_diag_manager.F90 index 47dc701798..0c9f875bcd 100644 --- a/src/framework/MOM_diag_manager_wrapper.F90 +++ b/src/framework/MOM_diag_manager.F90 @@ -1,14 +1,23 @@ -!> A simple (very thin) wrapper for register_diag_field to avoid a compiler bug with PGI -module MOM_diag_manager_wrapper +!> A simple (very thin) wrapper for the FMS diag_manager routines, with some name changes +module MOM_diag_manager ! This file is part of MOM6. See LICENSE.md for the license. use MOM_time_manager, only : time_type +use diag_axis_mod, only : diag_axis_init, get_diag_axis_name, EAST, NORTH +use diag_data_mod, only : null_axis_id +use diag_manager_mod, only : diag_manager_init, diag_manager_end +use diag_manager_mod, only : send_data, diag_field_add_attribute, DIAG_FIELD_NOT_FOUND use diag_manager_mod, only : 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 implicit none ; private -public register_diag_field_fms +public diag_manager_init, diag_manager_end +public diag_axis_init, get_diag_axis_name, EAST, NORTH, null_axis_id +public send_data, diag_field_add_attribute, DIAG_FIELD_NOT_FOUND +public register_diag_field_fms, register_static_field_fms, get_diag_field_id_fms !> A wrapper for register_diag_field_array() interface register_diag_field_fms @@ -85,11 +94,11 @@ integer function register_diag_field_scalar_fms(module_name, field_name, init_ti end function register_diag_field_scalar_fms -!> \namespace mom_diag_manager_wrapper +!> \namespace mom_diag_manager !! !! This module simply wraps register_diag_field() from FMS's diag_manager_mod. !! We used to be able to import register_diag_field and rename it to register_diag_field_fms !! with a simple "use, only : register_diag_field_fms => register_diag_field" but PGI 16.5 !! has a bug that refuses to compile this - earlier versions did work. -end module MOM_diag_manager_wrapper +end module MOM_diag_manager diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index fa4a4a2701..071585a951 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -9,37 +9,28 @@ module MOM_diag_mediator use MOM_coms, only : PE_here use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_diag_manager, only : diag_manager_init, diag_manager_end +use MOM_diag_manager, only : diag_axis_init, get_diag_axis_name, null_axis_id +use MOM_diag_manager, only : send_data, diag_field_add_attribute, EAST, NORTH +use MOM_diag_manager, only : register_diag_field_fms, register_static_field_fms +use MOM_diag_manager, only : get_diag_field_id_fms, DIAG_FIELD_NOT_FOUND +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 +use MOM_diag_remap, only : horizontally_average_diag_field, diag_remap_get_axes_info +use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured +use MOM_diag_remap, only : diag_remap_diag_registration_closed, diag_remap_set_active +use MOM_EOS, only : EOS_type use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, assert use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : slasher, vardesc, query_vardesc, mom_read_data +use MOM_io, only : slasher, vardesc, query_vardesc, MOM_read_data use MOM_io, only : get_filename_appendix use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : EOS_type -use MOM_diag_remap, only : diag_remap_ctrl -use MOM_diag_remap, only : diag_remap_update -use MOM_diag_remap, only : 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 -use MOM_diag_remap, only : diag_remap_configure_axes, diag_remap_axes_configured -use MOM_diag_remap, only : diag_remap_get_axes_info, diag_remap_set_active -use MOM_diag_remap, only : diag_remap_diag_registration_closed -use MOM_diag_remap, only : horizontally_average_diag_field - -use diag_axis_mod, only : get_diag_axis_name -use diag_data_mod, only : null_axis_id -use diag_manager_mod, only : diag_manager_init, diag_manager_end -use diag_manager_mod, only : send_data, diag_axis_init, EAST, NORTH, diag_field_add_attribute -! The following module is needed for PGI since the following line does not compile with PGI 6.5.0 -! was: use diag_manager_mod, only : register_diag_field_fms=>register_diag_field -use MOM_diag_manager_wrapper, only : register_diag_field_fms -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 diag_manager_mod, only : DIAG_FIELD_NOT_FOUND implicit none ; private @@ -482,10 +473,9 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) call define_axes_group(diag_cs, (/ id_xh, id_yq /), diag_cs%axesCv1, & x_cell_method='mean', y_cell_method='point', is_v_point=.true.) - ! Axis group for special null axis from diag manager + ! Axis group for special null axis from diag manager. (Could null_axis_id be made MOM specific?) call define_axes_group(diag_cs, (/ null_axis_id /), diag_cs%axesNull) - !Non-native Non-downsampled if (diag_cs%num_diag_coords>0) then allocate(diag_cs%remap_axesZL(diag_cs%num_diag_coords)) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index f27a153a2b..462ba7bf5e 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -60,6 +60,8 @@ module MOM_diag_remap use MOM_coms, only : reproducing_sum_EFP, EFP_to_real use MOM_coms, only : EFP_type, assignment(=), EFP_sum_across_PEs use MOM_error_handler, only : MOM_error, FATAL, assert, WARNING +use MOM_debugging, only : check_column_integrals +use MOM_diag_manager, only : diag_axis_init use MOM_diag_vkernels, only : interpolate_column, reintegrate_column use MOM_file_parser, only : get_param, log_param, param_file_type use MOM_io, only : slasher, mom_read_data @@ -80,9 +82,7 @@ module MOM_diag_remap use coord_sigma, only : build_sigma_column use coord_rho, only : build_rho_column -use diag_manager_mod, only : diag_axis_init -use MOM_debugging, only : check_column_integrals implicit none ; private public diag_remap_ctrl diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index fc29cd154f..69fa617a6c 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -33,6 +33,7 @@ module MOM_domains use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE => NORTH, EAST_FACE => EAST use mpp_domains_mod, only : global_field => mpp_global_field use mpp_domains_mod, only : mpp_redistribute +use mpp_domains_mod, only : broadcast_domain => mpp_broadcast_domain use fms_io_mod, only : file_exist, parse_mask_table use fms_affinity_mod, only : fms_affinity_init, fms_affinity_set, fms_affinity_get @@ -54,7 +55,7 @@ module MOM_domains public :: compute_block_extent, get_global_shape, get_layout_extents public :: MOM_thread_affinity_set, set_MOM_thread_affinity public :: get_simple_array_i_ind, get_simple_array_j_ind -public :: domain2D, domain1D, global_field, redistribute_array +public :: domain2D, domain1D, global_field, redistribute_array, broadcast_domain !> Do a halo update on an array interface pass_var diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 44df470928..b91509cc1d 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -3,34 +3,21 @@ module MOM_horizontal_regridding ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_debugging, only : hchksum -use MOM_coms, only : max_across_PEs, min_across_PEs -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP -use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast -use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID +use MOM_debugging, only : hchksum +use MOM_coms, only : max_across_PEs, min_across_PEs, sum_across_PEs, broadcast +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_LOOP +use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint -use MOM_file_parser, only : get_param, read_param, log_param, param_file_type -use MOM_file_parser, only : log_version -use MOM_get_input, only : directories -use MOM_grid, only : ocean_grid_type, isPointInCell -use MOM_io, only : close_file, fieldtype, file_exists -use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE -use MOM_io, only : slasher, vardesc, write_field -use MOM_string_functions, only : uppercase -use MOM_time_manager, only : time_type, get_external_field_size -use MOM_time_manager, only : init_external_field -use MOM_time_manager, only : get_external_field_axes, get_external_field_missing -use MOM_transform_FMS, only : time_interp_external => rotated_time_interp_external -use MOM_variables, only : thermo_var_ptrs - -use mpp_io_mod, only : axistype, mpp_get_axis_data -use mpp_mod, only : mpp_broadcast, mpp_sync, mpp_sync_self, mpp_max -use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type -use horiz_interp_mod, only : horiz_interp_init, horiz_interp_del - -use netcdf +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : time_interp_extern, get_external_field_info, horiz_interp_init +use MOM_interpolate, only : horiz_interp_new, horiz_interp, horiz_interp_type +use MOM_io_wrapper, only : axistype, get_axis_data +use MOM_time_manager, only : time_type + +use netcdf, only : NF90_OPEN, NF90_NOWRITE, NF90_GET_ATT, NF90_GET_VAR +use netcdf, only : NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, NF90_INQUIRE_DIMENSION implicit none ; private @@ -38,14 +25,6 @@ module MOM_horizontal_regridding public :: horiz_interp_and_extrap_tracer, myStats -! character(len=40) :: mdl = "MOM_horizontal_regridding" ! This module's name. - -!> Fill grid edges -interface fill_boundaries - module procedure fill_boundaries_real - module procedure fill_boundaries_int -end interface - !> Extrapolate and interpolate data interface horiz_interp_and_extrap_tracer module procedure horiz_interp_and_extrap_tracer_record @@ -94,10 +73,9 @@ subroutine myStats(array, missing, is, ie, js, je, k, mesg) end subroutine myStats -!> Use ICE-9 algorithm to populate points (fill=1) with -!! valid data (good=1). If no information is available, -!! Then use a previous guess (prev). Optionally (smooth) -!! blend the filled points to achieve a more desirable result. +!> Use ICE-9 algorithm to populate points (fill=1) with valid data (good=1). If no information +!! is available, use a previous guess (prev). Optionally (smooth) blend the filled points to +!! achieve a more desirable result. subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, debug, answers_2018) use MOM_coms, only : sum_across_PEs @@ -122,15 +100,20 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, !! answers as the code did in late 2018. Otherwise !! add parentheses for rotational symmetry. - real, dimension(SZI_(G),SZJ_(G)) :: b,r - real, dimension(SZI_(G),SZJ_(G)) :: fill_pts, good_, good_new - + real, dimension(SZI_(G),SZJ_(G)) :: a_filled ! The aout with missing values filled in + real, dimension(SZI_(G),SZJ_(G)) :: a_chg ! The change in aout due to an iteration of smoothing + real, dimension(SZI_(G),SZJ_(G)) :: fill_pts ! 1 for points that still need to be filled + real, dimension(SZI_(G),SZJ_(G)) :: good_ ! The values that are valid for the current iteration + real, dimension(SZI_(G),SZJ_(G)) :: good_new ! The values of good_ to use for the next iteration + + real :: east, west, north, south ! Valid neighboring values or 0 for invalid values + real :: ge, gw, gn, gs ! Flags indicating which neighbors have valid values + real :: ngood ! The number of valid values in neighboring points + logical :: do_smooth ! Indicates whether to do smoothing of the array + real :: nfill ! The remaining number of points to fill + real :: nfill_prev ! The previous value of nfill character(len=256) :: mesg ! The text of an error message - integer :: i,j,k - real :: east,west,north,south,sor - real :: ge,gw,gn,gs,ngood - logical :: do_smooth,siena_bug - real :: nfill, nfill_prev + integer :: i, j, k integer, parameter :: num_pass_default = 10000 real, parameter :: relc_default = 0.25, crit_default = 1.e-3 @@ -165,15 +148,15 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, nfill_prev = nfill good_(:,:) = good(:,:) - r(:,:) = 0.0 + a_chg(:,:) = 0.0 do while (nfill > 0.0) call pass_var(good_,G%Domain) call pass_var(aout,G%Domain) - b(:,:)=aout(:,:) - good_new(:,:)=good_(:,:) + a_filled(:,:) = aout(:,:) + good_new(:,:) = good_(:,:) do j=js,je ; do i=is,ie @@ -194,16 +177,16 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, endif if (ngood > 0.) then if (ans_2018) then - b(i,j)=(east+west+north+south)/ngood + a_filled(i,j) = (east+west+north+south)/ngood else - b(i,j) = ((east+west) + (north+south))/ngood + a_filled(i,j) = ((east+west) + (north+south))/ngood endif fill_pts(i,j) = 0.0 good_new(i,j) = 1.0 endif enddo ; enddo - aout(is:ie,js:je) = b(is:ie,js:je) + aout(is:ie,js:je) = a_filled(is:ie,js:je) good_(is:ie,js:je) = good_new(is:ie,js:je) nfill_prev = nfill nfill = sum(fill_pts(is:ie,js:je)) @@ -223,11 +206,13 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, call MOM_error(WARNING, mesg, .true.) endif + ! Determine the number of remaining points to fill globally. nfill = sum(fill_pts(is:ie,js:je)) call sum_across_PEs(nfill) - enddo + enddo ! while block for remaining points to fill. + ! Do Laplacian smoothing for the points that have been filled in. if (do_smooth) then ; do k=1,npass call pass_var(aout,G%Domain) do j=js,je ; do i=is,ie @@ -235,22 +220,22 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, smooth, num_pass, relc, crit, east = max(good(i+1,j),fill(i+1,j)) ; west = max(good(i-1,j),fill(i-1,j)) north = max(good(i,j+1),fill(i,j+1)) ; south = max(good(i,j-1),fill(i,j-1)) if (ans_2018) then - r(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + & - west*aout(i-1,j)+east*aout(i+1,j) - & - (south+north+west+east)*aout(i,j)) + a_chg(i,j) = relax_coeff*(south*aout(i,j-1)+north*aout(i,j+1) + & + west*aout(i-1,j)+east*aout(i+1,j) - & + (south+north+west+east)*aout(i,j)) else - r(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + & + a_chg(i,j) = relax_coeff*( ((south*aout(i,j-1) + north*aout(i,j+1)) + & (west*aout(i-1,j)+east*aout(i+1,j))) - & ((south+north)+(west+east))*aout(i,j) ) endif else - r(i,j) = 0. + a_chg(i,j) = 0. endif enddo ; enddo ares = 0.0 do j=js,je ; do i=is,ie - aout(i,j) = r(i,j) + aout(i,j) - ares = max(ares, abs(r(i,j))) + aout(i,j) = a_chg(i,j) + aout(i,j) + ares = max(ares, abs(a_chg(i,j))) enddo ; enddo call max_across_PEs(ares) if (ares <= acrit) exit @@ -299,18 +284,18 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, !! extrapolation is performed by this routine ! Local variables - real, dimension(:,:), allocatable :: tr_in, tr_inp ! A 2-d array for holding input data on - ! native horizontal grid and extended grid - ! with poles. + real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its + !! native horizontal grid. + real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles. real, dimension(:,:), allocatable :: mask_in ! A 2-d mask for extended input grid. real :: PI_180 integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp - integer :: i,j,k + integer :: i, j, k integer, dimension(4) :: start, count, dims, dim_id real, dimension(:,:), allocatable :: x_in, y_in - real, dimension(:), allocatable :: lon_in, lat_in - real, dimension(:), allocatable :: lat_inp, last_row + real, dimension(:), allocatable :: lon_in, lat_in ! The longitude and latitude in the input file + real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole real :: max_lat, min_lat, pole, max_depth, npole real :: roundoff ! The magnitude of roundoff, usually ~2e-16. real :: add_offset, scale_factor @@ -319,18 +304,21 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, character(len=8) :: laynum type(horiz_interp_type) :: Interp integer :: is, ie, js, je ! compute domain indices - integer :: isc,iec,jsc,jec ! global compute domain indices + integer :: isc, iec, jsc, jec ! global compute domain indices integer :: isg, ieg, jsg, jeg ! global extent integer :: isd, ied, jsd, jed ! data domain indices integer :: id_clock_read character(len=12) :: dim_name(4) logical :: debug=.false. - real :: npoints,varAvg - real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out - real, dimension(SZI_(G),SZJ_(G)) :: good, fill - real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev - real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2 - real, dimension(SZI_(G),SZJ_(G)) :: nlevs + real :: npoints, varAvg + real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out ! The longitude and latitude of points on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: tr_out, mask_out ! The tracer and mask on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1. + real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in + real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above + real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -338,14 +326,14 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP) - is_ongrid=.false. - if (present(ongrid)) is_ongrid=ongrid + is_ongrid = .false. + if (present(ongrid)) is_ongrid = ongrid if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) if (allocated(z_edges_in)) deallocate(z_edges_in) - PI_180=atan(1.0)/45. + PI_180 = atan(1.0)/45. ! Open NetCDF file and if present, extract data and spatial coordinate information ! The convention adopted here requires that the data be written in (i,j,k) ordering. @@ -401,7 +389,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) - allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1)) + allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1)) allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd)) start = 1 ; count = 1 ; count(1) = id @@ -421,7 +409,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (present(m_to_Z)) then ; do k=1,kd ; z_in(k) = m_to_Z * z_in(k) ; enddo ; endif - ! extrapolate the input data to the north pole using the northerm-most latitude + ! extrapolate the input data to the north pole using the northern-most latitude add_np = .false. jdp = jd if (.not. is_ongrid) then @@ -459,11 +447,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, allocate(tr_in(id,jd)) ; tr_in(:,:) = 0.0 allocate(tr_inp(id,jdp)) ; tr_inp(:,:) = 0.0 allocate(mask_in(id,jdp)) ; mask_in(:,:) = 0.0 - allocate(last_row(id)) ; last_row(:) = 0.0 endif max_depth = maxval(G%bathyT) - call mpp_max(max_depth) + call max_across_PEs(max_depth) if (z_edges_in(kd+1) abs(roundoff*missing_value)) then - pole = pole+last_row(i) - npole = npole+1.0 + pole = pole + tr_in(i,jd) + npole = npole + 1.0 endif enddo if (npole > 0) then @@ -521,9 +508,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, endif endif - call mpp_sync() - call mpp_broadcast(tr_inp, id*jdp, root_PE()) - call mpp_sync_self() + call broadcast(tr_inp, id*jdp, blocking=.true.) do j=1,jdp do i=1,id @@ -538,17 +523,15 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, endif - - ! call fms routine horiz_interp to interpolate input level data to model horizontal grid if (.not. is_ongrid) then if (k == 1) then - call horiz_interp_new(Interp,x_in,y_in,lon_out(is:ie,js:je),lat_out(is:ie,js:je), & - interp_method='bilinear',src_modulo=.true.) + call horiz_interp_new(Interp, x_in, y_in, lon_out(is:ie,js:je), lat_out(is:ie,js:je), & + interp_method='bilinear', src_modulo=.true.) endif if (debug) then - call myStats(tr_inp,missing_value, is,ie,js,je,k,'Tracer from file') + call myStats(tr_inp,missing_value, is, ie, js, je, k,'Tracer from file') endif endif @@ -556,7 +539,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (is_ongrid) then tr_out(is:ie,js:je)=tr_in(is:ie,js:je) else - call horiz_interp(Interp,tr_inp,tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.) + call horiz_interp(Interp, tr_inp, tr_out(is:ie,js:je), missing_value=missing_value, new_missing_handle=.true.) endif mask_out=1.0 @@ -609,7 +592,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, fill2(:,:) = fill(:,:) call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, smooth=.true., answers_2018=answers_2018) - call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()') + if (debug) then + call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()') + endif tr_z(:,:,k) = tr_outf(:,:) * G%mask2dT(:,:) mask_z(:,:,k) = good2(:,:) + fill2(:,:) @@ -652,9 +637,9 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t !! add parentheses for rotational symmetry. ! Local variables - real, dimension(:,:), allocatable :: tr_in,tr_inp !< A 2-d array for holding input data on - !! native horizontal grid and extended grid - !! with poles. + real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its + !! native horizontal grid. + real, dimension(:,:), allocatable :: tr_inp !< Native horizontal grid data extended to the poles. real, dimension(:,:,:), allocatable :: data_in !< A buffer for storing the full 3-d time-interpolated array !! on the original grid real, dimension(:,:), allocatable :: mask_in !< A 2-d mask for extended input grid. @@ -664,8 +649,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t integer :: i,j,k integer, dimension(4) :: start, count, dims, dim_id real, dimension(:,:), allocatable :: x_in, y_in - real, dimension(:), allocatable :: lon_in, lat_in - real, dimension(:), allocatable :: lat_inp, last_row + real, dimension(:), allocatable :: lon_in, lat_in ! The longitude and latitude in the input file + real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole real :: max_lat, min_lat, pole, max_depth, npole real :: roundoff ! The magnitude of roundoff, usually ~2e-16. logical :: add_np @@ -681,12 +666,15 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t character(len=12) :: dim_name(4) logical :: debug=.false. logical :: spongeDataOngrid - real :: npoints,varAvg - real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out - real, dimension(SZI_(G),SZJ_(G)) :: good, fill - real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev - real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2 - real, dimension(SZI_(G),SZJ_(G)) :: nlevs + real :: npoints, varAvg + real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out ! The longitude and latitude of points on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: tr_out, mask_out ! The tracer and mask on the model grid + real, dimension(SZI_(G),SZJ_(G)) :: good ! Where the data is valid, this is 1. + real, dimension(SZI_(G),SZJ_(G)) :: fill ! 1 where the data needs to be filled in + real, dimension(SZI_(G),SZJ_(G)) :: tr_outf ! The tracer concentrations after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: tr_prev ! The tracer concentrations in the layer above + real, dimension(SZI_(G),SZJ_(G)) :: good2 ! 1 where the data is valid after Ice-9 + real, dimension(SZI_(G),SZJ_(G)) :: fill2 ! 1 for points that still need to be filled after Ice-9 integer :: turns turns = G%HI%turns @@ -697,14 +685,14 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP) - PI_180=atan(1.0)/45. + PI_180 = atan(1.0)/45. ! Open NetCDF file and if present, extract data and spatial coordinate information ! The convention adopted here requires that the data be written in (i,j,k) ordering. call cpu_clock_begin(id_clock_read) - fld_sz = get_external_field_size(fms_id) + call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_value) if (allocated(lon_in)) deallocate(lon_in) if (allocated(lat_in)) deallocate(lat_in) if (allocated(z_in)) deallocate(z_in) @@ -712,34 +700,30 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) - axes_data = get_external_field_axes(fms_id) - id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3) - spongeDataOngrid=.false. - if (PRESENT(spongeOngrid)) spongeDataOngrid=spongeOngrid + spongeDataOngrid = .false. + if (PRESENT(spongeOngrid)) spongeDataOngrid = spongeOngrid if (.not. spongeDataOngrid) then - allocate(lon_in(id),lat_in(jd)) - call mpp_get_axis_data(axes_data(1), lon_in) - call mpp_get_axis_data(axes_data(2), lat_in) + allocate(lon_in(id), lat_in(jd)) + call get_axis_data(axes_data(1), lon_in) + call get_axis_data(axes_data(2), lat_in) endif - allocate(z_in(kd),z_edges_in(kd+1)) + allocate(z_in(kd), z_edges_in(kd+1)) allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd)) - call mpp_get_axis_data(axes_data(3), z_in) + call get_axis_data(axes_data(3), z_in) if (present(m_to_Z)) then ; do k=1,kd ; z_in(k) = m_to_Z * z_in(k) ; enddo ; endif call cpu_clock_end(id_clock_read) - missing_value = get_external_field_missing(fms_id) - if (.not. spongeDataOngrid) then - ! extrapolate the input data to the north pole using the northerm-most latitude + ! Extrapolate the input data to the north pole using the northerm-most latitude. max_lat = maxval(lat_in) - add_np=.false. + add_np = .false. if (max_lat < 90.0) then add_np = .true. jdp = jd+1 @@ -750,7 +734,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t allocate(lat_in(1:jdp)) lat_in(:) = lat_inp(:) else - jdp=jd + jdp = jd endif call horiz_interp_init() lon_in = lon_in*PI_180 @@ -763,7 +747,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t allocate(tr_in(id,jd)) ; tr_in(:,:)=0.0 allocate(tr_inp(id,jdp)) ; tr_inp(:,:)=0.0 allocate(mask_in(id,jdp)) ; mask_in(:,:)=0.0 - allocate(last_row(id)) ; last_row(:)=0.0 else allocate(data_in(isd:ied,jsd:jed,kd)) endif @@ -776,7 +759,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t max_depth = maxval(G%bathyT) - call mpp_max(max_depth) + call max_across_PEs(max_depth) if (z_edges_in(kd+1) abs(roundoff*missing_value)) then - pole = pole+last_row(i) - npole = npole+1.0 - endif - enddo - if (npole > 0) then - pole=pole/npole - else - pole=missing_value - endif - tr_inp(:,1:jd) = tr_in(:,:) - tr_inp(:,jdp) = pole + pole = 0.0 ; npole = 0.0 + do i=1,id + if (abs(tr_in(i,jd)-missing_value) > abs(roundoff*missing_value)) then + pole = pole + tr_in(i,jd) + npole = npole+1.0 + endif + enddo + if (npole > 0) then + pole = pole / npole + else + pole = missing_value + endif + tr_inp(:,1:jd) = tr_in(:,:) + tr_inp(:,jdp) = pole else - tr_inp(:,:) = tr_in(:,:) + tr_inp(:,:) = tr_in(:,:) endif endif - call mpp_sync() - call mpp_broadcast(tr_inp, id*jdp, root_PE()) - call mpp_sync_self() + call broadcast(tr_inp, id*jdp, blocking=.true.) - mask_in=0.0 + mask_in(:,:) = 0.0 do j=1,jdp ; do i=1,id if (abs(tr_inp(i,j)-missing_value) > abs(roundoff*missing_value)) then - mask_in(i,j)=1.0 + mask_in(i,j) = 1.0 tr_inp(i,j) = tr_inp(i,j) * conversion else tr_inp(i,j) = missing_value @@ -861,7 +841,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t endif if ((G%mask2dT(i,j) == 1.0) .and. (z_edges_in(k) <= G%bathyT(i,j)) .and. & (mask_out(i,j) < 1.0)) & - fill(i,j)=1.0 + fill(i,j) = 1.0 enddo ; enddo call pass_var(fill, G%Domain) call pass_var(good, G%Domain) @@ -905,15 +885,15 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) - do k=1,kd - do j=js,je - do i=is,ie - tr_z(i,j,k)=data_in(i,j,k) - if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. - enddo + call time_interp_extern(fms_id, Time, data_in, verbose=.true., turns=turns) + do k=1,kd + do j=js,je + do i=is,ie + tr_z(i,j,k) = data_in(i,j,k) + if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. enddo enddo + enddo endif end subroutine horiz_interp_and_extrap_tracer_fms_id @@ -925,9 +905,9 @@ subroutine meshgrid(x, y, x_T, y_T) real, dimension(size(x,1),size(y,1)), intent(inout) :: x_T !< output 2-dimensional array real, dimension(size(x,1),size(y,1)), intent(inout) :: y_T !< output 2-dimensional array - integer :: ni,nj,i,j + integer :: ni, nj, i, j - ni=size(x,1) ; nj=size(y,1) + ni = size(x,1) ; nj = size(y,1) do j=1,nj ; do i=1,ni x_T(i,j) = x(i) @@ -939,122 +919,4 @@ subroutine meshgrid(x, y, x_T, y_T) end subroutine meshgrid -! None of the subsequent code appears to be used at all. - -!> Fill grid edges for integer data -function fill_boundaries_int(m,cyclic_x,tripolar_n) result(mp) - integer, dimension(:,:), intent(in) :: m !< input array (ND) - logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant - logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold - integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - - real, dimension(size(m,1),size(m,2)) :: m_real - real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real - - m_real = real(m) - - mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n) - - mp = int(mp_real) - -end function fill_boundaries_int - -!> Fill grid edges for real data -function fill_boundaries_real(m,cyclic_x,tripolar_n) result(mp) - real, dimension(:,:), intent(in) :: m !< input array (ND) - logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant - logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold - real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - - integer :: ni,nj,i,j - - ni=size(m,1); nj=size(m,2) - - mp(1:ni,1:nj)=m(:,:) - - if (cyclic_x) then - mp(0,1:nj)=m(ni,1:nj) - mp(ni+1,1:nj)=m(1,1:nj) - else - mp(0,1:nj)=m(1,1:nj) - mp(ni+1,1:nj)=m(ni,1:nj) - endif - - mp(1:ni,0)=m(1:ni,1) - if (tripolar_n) then - do i=1,ni - mp(i,nj+1)=m(ni-i+1,nj) - enddo - else - mp(1:ni,nj+1)=m(1:ni,nj) - endif - -end function fill_boundaries_real - -!> Solve del2 (zi) = 0 using successive iterations -!! with a 5 point stencil. Only points fill==1 are -!! modified. Except where bad==1, information propagates -!! isotropically in index space. The resulting solution -!! in each region is an approximation to del2(zi)=0 subject to -!! boundary conditions along the valid points curve bounding this region. -subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) - real, dimension(:,:), intent(inout) :: zi !< input and output array (ND) - integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< same shape as zi, 1=fill - integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< same shape as zi, 1=bad data - real, intent(in) :: sor !< relaxation coefficient (ND) - integer, intent(in) :: niter !< maximum number of iterations - logical, intent(in) :: cyclic_x !< true if domain is zonally reentrant - logical, intent(in) :: tripolar_n !< true if domain has an Arctic fold - - ! Local variables - real, dimension(size(zi,1),size(zi,2)) :: res, m - integer, dimension(size(zi,1),size(zi,2),4) :: B - real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp - integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm - integer :: i,j,k,n - integer :: ni,nj - real :: Isum, bsum - - ni=size(zi,1) ; nj=size(zi,2) - - - mp(:,:) = fill_boundaries(zi,cyclic_x,tripolar_n) - - B(:,:,:) = 0.0 - nm(:,:) = fill_boundaries(bad,cyclic_x,tripolar_n) - - do j=1,nj - do i=1,ni - if (fill(i,j) == 1) then - B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) - B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) - endif - enddo - enddo - - do n=1,niter - do j=1,nj - do i=1,ni - if (fill(i,j) == 1) then - bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) - Isum = 1.0/bsum - res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& - B(i,j,3)*mp(i,j+1)+B(i,j,4)*mp(i,j-1)) - mp(i,j) - endif - enddo - enddo - res(:,:)=res(:,:)*sor - - do j=1,nj - do i=1,ni - mp(i,j)=mp(i,j)+res(i,j) - enddo - enddo - - zi(:,:)=mp(1:ni,1:nj) - mp = fill_boundaries(zi,cyclic_x,tripolar_n) - enddo - -end subroutine smooth_heights - end module MOM_horizontal_regridding diff --git a/src/framework/MOM_interpolate.F90 b/src/framework/MOM_interpolate.F90 new file mode 100644 index 0000000000..c63c847e55 --- /dev/null +++ b/src/framework/MOM_interpolate.F90 @@ -0,0 +1,143 @@ +!> This module wraps the FMS temporal and spatial interpolation routines +module MOM_interpolate + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_array_transform, only : allocate_rotated_array, rotate_array +use MOM_error_handler, only : MOM_error, FATAL +use MOM_io_wrapper, only : axistype +use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use time_interp_external_mod, only : time_interp_external_fms=>time_interp_external +use time_interp_external_mod, only : init_external_field, time_interp_external_init +use time_interp_external_mod, only : get_external_field_size +use time_interp_external_mod, only : get_external_field_axes, get_external_field_missing +use time_manager_mod, only : time_type + +implicit none ; private + +public :: time_interp_extern, init_external_field, time_interp_external_init +public :: get_external_field_info +public :: horiz_interp_type, horiz_interp_init, horiz_interp, horiz_interp_new + +!> Read a field based on model time, and rotate to the model domain. +! This inerface does not share the name time_interp_external with the module it primarily +! wraps because of errors (perhaps a bug) that arise with the PGI 19.10.0 compiler. +interface time_interp_extern + module procedure time_interp_external_0d + module procedure time_interp_external_2d + module procedure time_interp_external_3d +end interface time_interp_extern + +contains + +!> Get information about the external fields. +subroutine get_external_field_info(field_id, size, axes, missing) + integer, intent(in) :: field_id !< The integer index of the external + !! field returned from a previous + !! call to init_external_field() + integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data + type(axistype), dimension(4), optional, intent(inout) :: axes !< Axis types for the input data + real, optional, intent(inout) :: missing !< Missing value for the input data + + if (present(size)) then + size(1:4) = get_external_field_size(field_id) + endif + + if (present(axes)) then + axes(1:4) = get_external_field_axes(field_id) + endif + + if (present(missing)) then + missing = get_external_field_missing(field_id) + endif + +end subroutine get_external_field_info + + +!> Read a scalar field based on model time. +subroutine time_interp_external_0d(field_id, time, data_in, verbose) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, intent(inout) :: data_in !< The interpolated value + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + + call time_interp_external_fms(field_id, time, data_in, verbose=verbose) +end subroutine time_interp_external_0d + +!> Read a 2d field from an external based on model time, potentially including horizontal +!! interpolation and rotation of the data +subroutine time_interp_external_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out, turns) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, dimension(:,:), intent(inout) :: data_in !< The array in which to store the interpolated values + integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + type(horiz_interp_type), & + optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation + logical, dimension(:,:), & + optional, intent(out) :: mask_out !< An array that is true where there is valid data + integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data + + real, allocatable :: data_pre_rot(:,:) ! The input data before rotation + integer :: qturns ! The number of quarter turns to rotate the data + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call time_interp_external_fms(field_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + else + call allocate_rotated_array(data_in, [1,1], -qturns, data_pre_rot) + call time_interp_external_fms(field_id, time, data_pre_rot, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + call rotate_array(data_pre_rot, turns, data_in) + deallocate(data_pre_rot) + endif +end subroutine time_interp_external_2d + + +!> Read a 3d field based on model time, and rotate to the model grid +subroutine time_interp_external_3d(field_id, time, data_in, interp, & + verbose, horz_interp, mask_out, turns) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, dimension(:,:,:), intent(inout) :: data_in !< The array in which to store the interpolated values + integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + type(horiz_interp_type), & + optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation + logical, dimension(:,:,:), & + optional, intent(out) :: mask_out !< An array that is true where there is valid data + integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data + + real, allocatable :: data_pre_rot(:,:,:) ! The input data before rotation + integer :: qturns ! The number of quarter turns to rotate the data + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call time_interp_external_fms(field_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + else + call allocate_rotated_array(data_in, [1,1,1], -qturns, data_pre_rot) + call time_interp_external_fms(field_id, time, data_pre_rot, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + call rotate_array(data_pre_rot, turns, data_in) + deallocate(data_pre_rot) + endif +end subroutine time_interp_external_3d + +end module MOM_interpolate diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index b9b42d0cbf..aadbf1ace0 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -57,8 +57,8 @@ module MOM_ice_shelf use MOM_coms, only : reproducing_sum use MOM_spatial_means, only : global_area_integral use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init +use MOM_interpolate, only : init_external_field, time_interp_extern, time_interp_external_init + implicit none ; private #include @@ -1084,9 +1084,9 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) do j=js,je ; do i=is,ie last_hmask(i,j) = ISS%hmask(i,j) ; last_area_shelf_h(i,j) = ISS%area_shelf_h(i,j) enddo ; enddo - call time_interp_external(CS%id_read_mass, Time0, last_mass_shelf) + call time_interp_extern(CS%id_read_mass, Time0, last_mass_shelf) do j=js,je ; do i=is,ie - ! This should only be done if time_interp_external did an update. + ! This should only be done if time_interp_extern did an update. last_mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * last_mass_shelf(i,j) ! Rescale after time_interp last_h_shelf(i,j) = last_mass_shelf(i,j) / CS%density_ice enddo ; enddo @@ -1510,9 +1510,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, inputdir = slasher(inputdir) TideAmp_file = trim(inputdir) // trim(TideAmp_file) if (CS%rotate_index) then - allocate(tmp2d(CS%Grid_in%isd:CS%Grid_in%ied,CS%Grid_in%jsd:CS%Grid_in%jed)) ; tmp2d(:,:)=0.0 + allocate(tmp2d(CS%Grid_in%isd:CS%Grid_in%ied,CS%Grid_in%jsd:CS%Grid_in%jed)) ; tmp2d(:,:) = 0.0 call MOM_read_data(TideAmp_file, 'tideamp', tmp2d, CS%Grid_in%domain, timelevel=1, scale=US%m_s_to_L_T) - call rotate_array(tmp2d,CS%turns, CS%utide) + call rotate_array(tmp2d, CS%turns, CS%utide) deallocate(tmp2d) else call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) @@ -1984,8 +1984,8 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) allocate(tmp2d(is:ie,js:je)) ; tmp2d(:,:) = 0.0 endif - call time_interp_external(CS%id_read_mass, Time, tmp2d) - call rotate_array(tmp2d,CS%turns, ISS%mass_shelf) + call time_interp_extern(CS%id_read_mass, Time, tmp2d) + call rotate_array(tmp2d, CS%turns, ISS%mass_shelf) deallocate(tmp2d) ! This should only be done if time_interp_external did an update. diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index 90ae47450d..397696c0ba 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -3,19 +3,15 @@ module MOM_IS_diag_mediator ! This file is a part of SIS2. See LICENSE.md for the license. -use MOM_grid, only : ocean_grid_type - -use MOM_coms, only : PE_here +use MOM_coms, only : PE_here +use MOM_diag_manager, only : diag_manager_init, send_data, diag_axis_init, EAST, NORTH +use MOM_diag_manager, only : register_diag_field_fms, register_static_field_fms use MOM_error_handler, only : MOM_error, FATAL, is_root_pe, assert -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase, uppercase, slasher -use MOM_time_manager, only : time_type - -use diag_manager_mod, only : diag_manager_init -use diag_manager_mod, only : send_data, diag_axis_init,EAST,NORTH -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 MOM_time_manager, only : time_type implicit none ; private diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index a85f9c8484..42c4894dd3 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -6,16 +6,16 @@ module MOM_oda_driver_mod ! MOM infrastructure use MOM_error_handler, only : stdout, stdlog, MOM_error use MOM_coms, only : PE_here, num_PEs -use MOM_coms, only : set_PElist, set_rootPE, Get_PElist, broadcast, broadcast_domain +use MOM_coms, only : set_PElist, set_rootPE, Get_PElist, broadcast use MOM_io, only : SINGLE_FILE use MOM_domains, only : domain2d, global_field, get_domain_extent -use MOM_domains, only : pass_var, redistribute_array +use MOM_domains, only : pass_var, redistribute_array, broadcast_domain use MOM_diag_mediator, only : register_diag_field, diag_axis_init, post_data use MOM_ensemble_manager, only : get_ensemble_id, get_ensemble_size use MOM_ensemble_manager, only : get_ensemble_pelist, get_ensemble_filter_pelist use MOM_time_manager, only : time_type, real_to_time, get_date -use MOM_time_manager, only : operator(+), operator(>=), operator(/=) -use MOM_time_manager, only : operator(==),operator(<) +use MOM_time_manager, only : operator(+), operator(>=), operator(/=) +use MOM_time_manager, only : operator(==), operator(<) ! ODA Modules use ocean_da_types_mod, only : grid_type, ocean_profile_type, ocean_control_struct use ocean_da_core_mod, only : ocean_da_core_init, get_profiles diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index ed3ef7173e..003b134b2a 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -196,7 +196,6 @@ module MOM_hor_visc integer :: id_normstress = -1, id_shearstress = -1 !>@} - end type hor_visc_CS contains @@ -265,8 +264,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME [H L2 T-2 ~> m3 s-2 or kg s-2] bhstr_xx, & ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R L2 T-3 ~> W m-2] - ! Leith_Kh_h, & ! Leith Laplacian viscosity at h-points [L2 T-1 ~> m2 s-1] - ! Leith_Ah_h, & ! Leith bi-harmonic viscosity at h-points [L4 T-1 ~> m4 s-1] grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1] Del2vort_h, & ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1] @@ -289,7 +286,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xy, & ! str_xy is the cross term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2] str_xy_GME, & ! smoothed cross term in the stress tensor from GME [H L2 T-2 ~> m3 s-2 or kg s-2] bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2] - vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1] + vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1] Leith_Kh_q, & ! Leith Laplacian viscosity at q-points [L2 T-1 ~> m2 s-1] Leith_Ah_q, & ! Leith bi-harmonic viscosity at q-points [L4 T-1 ~> m4 s-1] grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1] @@ -297,8 +294,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, Del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1] grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1] grad_vel_mag_q, & ! Magnitude of the velocity gradient tensor squared at q-points [T-2 ~> s-2] - hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] - ! This form guarantees that hq/hu < 4. + hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2] + ! This form guarantees that hq/hu < 4. grad_vel_mag_bt_q, & ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2] boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim] @@ -323,30 +320,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, div_xx_h, & ! horizontal divergence [T-1 ~> s-1] sh_xx_h, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1] NoSt ! A diagnostic array of normal stress [T-1 ~> s-1]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & - grid_Re_Kh, & !< Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] - grid_Re_Ah, & !< Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] - GME_coeff_h !< GME coeff. at h-points [L2 T-1 ~> m2 s-1] - real :: Ah ! biharmonic viscosity [L4 T-1 ~> m4 s-1] - real :: Kh ! Laplacian viscosity [L2 T-1 ~> m2 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & + grid_Re_Kh, & ! Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim] + grid_Re_Ah, & ! Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim] + GME_coeff_h ! GME coeff. at h-points [L2 T-1 ~> m2 s-1] real :: AhSm ! Smagorinsky biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: AhLth ! 2D Leith biharmonic viscosity [L4 T-1 ~> m4 s-1] real :: mod_Leith ! nondimensional coefficient for divergence part of modified Leith ! viscosity. Here set equal to nondimensional Laplacian Leith constant. ! This is set equal to zero if modified Leith is not used. - real :: Shear_mag ! magnitude of the shear [T-1 ~> s-1] - real :: vert_vort_mag ! magnitude of the vertical vorticity gradient [L-1 T-1 ~> m-1 s-1] + real :: Shear_mag_bc ! Shear_mag value in backscatter [T-1 ~> s-1] + real :: sh_xx_sq ! Square of tension (sh_xx) [T-2 ~> s-2] + real :: sh_xy_sq ! Square of shearing strain (sh_xy) [T-2 ~> s-2] real :: h2uq, h2vq ! temporary variables [H2 ~> m2 or kg2 m-4]. real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner ! points; these are first interpolated to u or v velocity ! points where masks are applied [H ~> m or kg m-2]. real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2] real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6] - real :: hrat_min ! minimum thicknesses at the 4 neighboring - ! velocity points divided by the thickness at the stress - ! point (h or q point) [nondim] - real :: visc_bound_rem ! fraction of overall viscous bounds that - ! remain to be applied [nondim] + real :: h_min ! Minimum h at the 4 neighboring velocity points [H ~> m] real :: Kh_scale ! A factor between 0 and 1 by which the horizontal ! Laplacian viscosity is rescaled [nondim] real :: RoScl ! The scaling function for MEKE source term [nondim] @@ -365,6 +357,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! calculation gives the same value as if f were 0 [nondim]. real :: H0_GME ! Depth used to scale down GME coefficient in shallow areas [Z ~> m] real :: KE ! Local kinetic energy [L2 T-2 ~> m2 s-2] + real :: d_del2u ! dy-weighted Laplacian(u) diff in x [L-2 T-1 ~> m-2 s-1] + real :: d_del2v ! dx-weighted Laplacian(v) diff in y [L-2 T-1 ~> m-2 s-1] + real :: d_str ! Stress tensor update [H L2 T-2 ~> m3 s-2 or kg s-2] + real :: grad_vort ! Vorticity gradient magnitude [L-1 T-1 ~> m-1 s-1] + real :: grad_vort_qg ! QG-based vorticity gradient magnitude [L-1 T-1 ~> m-1 s-1] + real :: grid_Kh ! Laplacian viscosity bound by grid [L2 T-1 ~> m2 s-1] + real :: grid_Ah ! Biharmonic viscosity bound by grid [L4 T-1 ~> m4 s-1] logical :: rescale_Kh, legacy_bound logical :: find_FrictWork @@ -374,6 +373,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: i, j, k, n real :: inv_PI3, inv_PI2, inv_PI6 + + ! Fields evaluated on active layers, used for constructing 3D stress fields + ! NOTE: The position of these declarations can impact performance, due to the + ! very large number of stack arrays in this function. Move with caution! + real, dimension(SZIB_(G),SZJB_(G)) :: & + Ah, & ! biharmonic viscosity (h or q) [L4 T-1 ~> m4 s-1] + Kh, & ! Laplacian viscosity [L2 T-1 ~> m2 s-1] + Shear_mag, & ! magnitude of the shear [T-1 ~> s-1] + vert_vort_mag, & ! magnitude of the vertical vorticity gradient [L-1 T-1 ~> m-1 s-1] + hrat_min, & ! h_min divided by the thickness at the stress point (h or q) [nondim] + visc_bound_rem ! fraction of overall viscous bounds that remain to be applied [nondim] + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -383,9 +394,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, inv_PI2 = 1.0/((4.0*atan(1.0))**2) inv_PI6 = inv_PI3 * inv_PI3 - Ah_h(:,:,:) = 0.0 - Kh_h(:,:,:) = 0.0 - if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally apply_OBC = .true. @@ -505,10 +513,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !$OMP str_xx, str_xy, bhstr_xx, bhstr_xy, str_xx_GME, str_xy_GME, & !$OMP vort_xy, vort_xy_dx, vort_xy_dy, div_xx, div_xx_dx, div_xx_dy, & !$OMP grad_div_mag_h, grad_div_mag_q, grad_vort_mag_h, grad_vort_mag_q, & - !$OMP grad_vort_mag_h_2d, grad_vort_mag_q_2d, & + !$OMP grad_vort, grad_vort_qg, grad_vort_mag_h_2d, grad_vort_mag_q_2d, & !$OMP grad_vel_mag_h, grad_vel_mag_q, & !$OMP grad_vel_mag_bt_h, grad_vel_mag_bt_q, grad_d2vel_mag_h, & - !$OMP meke_res_fn, Shear_mag, vert_vort_mag, hrat_min, visc_bound_rem, & + !$OMP meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, h_min, hrat_min, visc_bound_rem, & + !$OMP sh_xx_sq, sh_xy_sq, grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, & !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, & !$OMP dDel2vdx, dDel2udy, DY_dxCv, DX_dyCu, Del2vort_q, Del2vort_h, KE, & !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff & @@ -519,6 +528,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, ! shearing strain advocated by Smagorinsky (1993) and discussed in ! Griffies and Hallberg (2000). + ! NOTE: There is a ~1% speedup when the tension and shearing loops below + ! are fused (presumably due to shared access of Id[xy]C[uv]). However, + ! this breaks the center/vertex index case convention, and also evaluates + ! the dudx and dvdy terms beyond their valid bounds. + ! TODO: Explore methods for retaining both the syntax and speedup. + ! Calculate horizontal tension do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 dudx(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u(I,j,k) - & @@ -526,7 +541,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, dvdy(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v(i,J,k) - & G%IdxCv(i,J-1) * v(i,J-1,k)) sh_xx(i,j) = dudx(i,j) - dvdy(i,j) - if (CS%id_normstress > 0) NoSt(i,j,k) = sh_xx(i,j) enddo ; enddo ! Components for the shearing strain @@ -535,6 +549,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j)) enddo ; enddo + if (CS%id_normstress > 0) then + do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2 + NoSt(i,j,k) = sh_xx(i,j) + enddo ; enddo + endif + ! Interpolate the thicknesses to velocity points. ! The extra wide halos are to accommodate the cross-corner-point projections ! in OBCs, which are not ordinarily be necessary, and might not be necessary @@ -608,7 +628,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif - if (OBC%segment(n)%direction == OBC_DIRECTION_N) then ! There are extra wide halos here to accommodate the cross-corner-point ! OBC projections, but they might not be necessary if the accelerations @@ -765,7 +784,6 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_div_mag_h(i,j) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2) enddo ; enddo - !do J=js-1,Jeq ; do I=is-1,Ieq do j=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1 grad_div_mag_q(I,J) =sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + & (0.5 * (div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2) @@ -828,133 +846,249 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, meke_res_fn = 1. - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then - Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + & - 0.25*((sh_xy(I-1,J-1)*sh_xy(I-1,J-1) + sh_xy(I,J)*sh_xy(I,J)) + & - (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1)))) + if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + sh_xx_sq = sh_xx(i,j) * sh_xx(i,j) + sh_xy_sq = 0.25 * ( & + (sh_xy(I-1,J-1) * sh_xy(I-1,J-1) + sh_xy(I,J) * sh_xy(I,J)) & + + (sh_xy(I-1,J) * sh_xy(I-1,J) + sh_xy(I,J-1) * sh_xy(I,J-1)) & + ) + Shear_mag(i,j) = sqrt(sh_xx_sq + sh_xy_sq) + enddo ; enddo + endif + + if (CS%better_bound_Ah .or. CS%better_bound_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + h_min = min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) + hrat_min(i,j) = min(1.0, h_min / (h(i,j,k) + h_neglect)) + enddo ; enddo + + if (CS%better_bound_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + visc_bound_rem(i,j) = 1.0 + enddo ; enddo endif + endif + + if (CS%Laplacian) then if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then - vert_vort_mag = MIN(grad_vort_mag_h(i,j) + grad_div_mag_h(i,j),3.*grad_vort_mag_h_2d(i,j)) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + grad_vort = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j) + grad_vort_qg = 3. * grad_vort_mag_h_2d(i,j) + vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg) + enddo ; enddo else - vert_vort_mag = (grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)) + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + vert_vort_mag(i,j) = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j) + enddo ; enddo endif endif - if (CS%better_bound_Ah .or. CS%better_bound_Kh) then - hrat_min = min(1.0, min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) / & - (h(i,j,k) + h_neglect) ) - visc_bound_rem = 1.0 - endif - if (CS%Laplacian) then - ! Determine the Laplacian viscosity at h points, using the - ! largest value from several parameterizations. - Kh = CS%Kh_bg_xx(i,j) ! Static (pre-computed) background viscosity + ! Determine the Laplacian viscosity at h points, using the + ! largest value from several parameterizations. + + ! Static (pre-computed) background viscosity + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = CS%Kh_bg_xx(i,j) + enddo ; enddo + + ! NOTE: The following do-block can be decomposed and vectorized after the + ! stack size has been reduced. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if (CS%add_LES_viscosity) then - if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag - if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + if (CS%Smagorinsky_Kh) & + Kh(i,j) = Kh(i,j) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j) + if (CS%Leith_Kh) & + Kh(i,j) = Kh(i,j) + CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3 else - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xx(i,j) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3) + if (CS%Smagorinsky_Kh) & + Kh(i,j) = max(Kh(i,j), CS%Laplac2_const_xx(i,j) * Shear_mag(i,j)) + if (CS%Leith_Kh) & + Kh(i,j) = max(Kh(i,j), CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3) endif - ! All viscosity contributions above are subject to resolution scaling - if (rescale_Kh) Kh = VarMix%Res_fn_h(i,j) * Kh - if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j) + enddo ; enddo + + ! All viscosity contributions above are subject to resolution scaling + + if (rescale_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = VarMix%Res_fn_h(i,j) * Kh(i,j) + enddo ; enddo + endif + + if (legacy_bound) then ! Older method of bounding for stability - if (legacy_bound) Kh = min(Kh, CS%Kh_Max_xx(i,j)) - Kh = max( Kh, CS%Kh_bg_min ) ! Place a floor on the viscosity, if desired. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = min(Kh(i,j), CS%Kh_Max_xx(i,j)) + enddo ; enddo + endif + + ! Place a floor on the viscosity, if desired. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) + enddo ; enddo + + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_h(i,j) + if (use_MEKE_Ku) & - Kh = Kh + MEKE%Ku(i,j) * meke_res_fn ! *Add* the MEKE contribution (might be negative) - if (CS%anisotropic) Kh = Kh + CS%Kh_aniso * ( 1. - CS%n1n2_h(i,j)**2 ) ! *Add* the tension component - ! of anisotropic viscosity + ! *Add* the MEKE contribution (might be negative) + Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * meke_res_fn + enddo ; enddo - ! Newer method of bounding for stability + if (CS%anisotropic) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + ! *Add* the tension component of anisotropic viscosity + Kh(i,j) = Kh(i,j) + CS%Kh_aniso * (1. - CS%n1n2_h(i,j)**2) + enddo ; enddo + endif + + ! Newer method of bounding for stability + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if (CS%better_bound_Kh) then - if (Kh >= hrat_min*CS%Kh_Max_xx(i,j)) then - visc_bound_rem = 0.0 - Kh = hrat_min*CS%Kh_Max_xx(i,j) + if (Kh(i,j) >= hrat_min(i,j) * CS%Kh_Max_xx(i,j)) then + visc_bound_rem(i,j) = 0.0 + Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xx(i,j) else - visc_bound_rem = 1.0 - Kh / (hrat_min*CS%Kh_Max_xx(i,j)) + visc_bound_rem(i,j) = 1.0 - Kh(i,j) / (hrat_min(i,j) * CS%Kh_Max_xx(i,j)) endif endif + enddo ; enddo - if ((CS%id_Kh_h>0) .or. find_FrictWork .or. CS%debug) Kh_h(i,j,k) = Kh + if (CS%id_Kh_h>0 .or. CS%debug) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Kh_h(i,j,k) = Kh(i,j) + enddo ; enddo + endif - if (CS%id_grid_Re_Kh>0) then + if (CS%id_grid_Re_Kh>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j))) & - / max(Kh, CS%min_grid_Kh) - endif + grid_Kh = max(Kh(i,j), CS%min_grid_Kh) + grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j))) / grid_Kh + enddo ; enddo + endif - if (CS%id_div_xx_h>0) div_xx_h(i,j,k) = div_xx(i,j) - if (CS%id_sh_xx_h>0) sh_xx_h(i,j,k) = sh_xx(i,j) + if (CS%id_div_xx_h>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + div_xx_h(i,j,k) = div_xx(i,j) + enddo ; enddo + endif + + if (CS%id_sh_xx_h>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + sh_xx_h(i,j,k) = sh_xx(i,j) + enddo ; enddo + endif - str_xx(i,j) = -Kh * sh_xx(i,j) - else ! not Laplacian + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + str_xx(i,j) = -Kh(i,j) * sh_xx(i,j) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 str_xx(i,j) = 0.0 - endif ! Laplacian + enddo ; enddo + endif - if (CS%anisotropic) then + if (CS%anisotropic) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ! Shearing-strain averaged to h-points local_strain = 0.25 * ( (sh_xy(I,J) + sh_xy(I-1,J-1)) + (sh_xy(I-1,J) + sh_xy(I,J-1)) ) ! *Add* the shear-strain contribution to the xx-component of stress str_xx(i,j) = str_xx(i,j) - CS%Kh_aniso * CS%n1n2_h(i,j) * CS%n1n1_m_n2n2_h(i,j) * local_strain - endif + enddo ; enddo + endif - if (CS%biharmonic) then - ! Determine the biharmonic viscosity at h points, using the - ! largest value from several parameterizations. - AhSm = 0.0; AhLth = 0.0 - if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then - if (CS%Smagorinsky_Ah) then - if (CS%bound_Coriolis) then - AhSm = Shear_mag * (CS%Biharm_const_xx(i,j) + & - CS%Biharm_const2_xx(i,j)*Shear_mag) - else - AhSm = CS%Biharm_const_xx(i,j) * Shear_mag - endif - endif - if (CS%Leith_Ah) AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h(i,j)) * inv_PI6 - Ah = MAX(MAX(CS%Ah_bg_xx(i,j), AhSm), AhLth) - if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & - Ah = MIN(Ah, CS%Ah_Max_xx(i,j)) - else - Ah = CS%Ah_bg_xx(i,j) - endif ! Smagorinsky_Ah or Leith_Ah + if (CS%biharmonic) then + ! Determine the biharmonic viscosity at h points, using the + ! largest value from several parameterizations. + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = CS%Ah_bg_xx(i,j) + enddo ; enddo - if (use_MEKE_Au) Ah = Ah + MEKE%Au(i,j) ! *Add* the MEKE contribution + if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah)) then + if (CS%Smagorinsky_Ah) then + if (CS%bound_Coriolis) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + AhSm = Shear_mag(i,j) * (CS%Biharm_const_xx(i,j) & + + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j) & + ) + Ah(i,j) = max(Ah(i,j), AhSm) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + AhSm = CS%Biharm_const_xx(i,j) * Shear_mag(i,j) + Ah(i,j) = max(Ah(i,j), AhSm) + enddo ; enddo + endif + endif - if (CS%Re_Ah > 0.0) then - KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - Ah = sqrt(KE) * CS%Re_Ah_const_xx(i,j) + if (CS%Leith_Ah) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h(i,j)) * inv_PI6 + Ah(i,j) = max(Ah(i,j), AhLth) + enddo ; enddo endif - if (CS%better_bound_Ah) then - Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xx(i,j)) + if (CS%bound_Ah .and. .not. CS%better_bound_Ah) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xx(i,j)) + enddo ; enddo endif + endif ! Smagorinsky_Ah or Leith_Ah - if ((CS%id_Ah_h>0) .or. find_FrictWork .or. CS%debug) Ah_h(i,j,k) = Ah + if (use_MEKE_Au) then + ! *Add* the MEKE contribution + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = Ah(i,j) + MEKE%Au(i,j) + enddo ; enddo + endif - if (CS%id_grid_Re_Ah>0) then + if (CS%Re_Ah > 0.0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2) - grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j)) & - / max(Ah, CS%min_grid_Ah) + Ah(i,j) = sqrt(KE) * CS%Re_Ah_const_xx(i,j) + enddo ; enddo + endif + + if (CS%better_bound_Ah) then + if (CS%better_bound_Kh) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xx(i,j)) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah(i,j) = min(Ah(i,j), hrat_min(i,j) * CS%Ah_Max_xx(i,j)) + enddo ; enddo endif + endif - str_xx(i,j) = str_xx(i,j) + Ah * & - (CS%DY_dxT(i,j) * (G%IdyCu(I,j)*Del2u(I,j) - G%IdyCu(I-1,j)*Del2u(I-1,j)) - & - CS%DX_dyT(i,j) * (G%IdxCv(i,J)*Del2v(i,J) - G%IdxCv(i,J-1)*Del2v(i,J-1))) + if ((CS%id_Ah_h>0) .or. CS%debug) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + Ah_h(i,j,k) = Ah(i,j) + enddo ; enddo + endif - ! Keep a copy of the biharmonic contribution for backscatter parameterization - bhstr_xx(i,j) = Ah * & - (CS%DY_dxT(i,j) * (G%IdyCu(I,j)*Del2u(I,j) - G%IdyCu(I-1,j)*Del2u(I-1,j)) - & - CS%DX_dyT(i,j) * (G%IdxCv(i,J)*Del2v(i,J) - G%IdxCv(i,J-1)*Del2v(i,J-1))) - bhstr_xx(i,j) = bhstr_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j)) + if (CS%id_grid_Re_Ah>0) then + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + KE = 0.125 * ((u(I,j,k) + u(I-1,j,k))**2 + (v(i,J,k) + v(i,J-1,k))**2) + grid_Ah = max(Ah(i,j), CS%min_grid_Ah) + grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j)) / grid_Ah + enddo ; enddo + endif - endif ! biharmonic + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + d_del2u = G%IdyCu(I,j) * Del2u(I,j) - G%IdyCu(I-1,j) * Del2u(I-1,j) + d_del2v = G%IdxCv(i,J) * Del2v(i,J) - G%IdxCv(i,J-1) * Del2v(i,J-1) + d_str = Ah(i,j) * (CS%DY_dxT(i,j) * d_del2u - CS%DX_dyT(i,j) * d_del2v) - enddo ; enddo + str_xx(i,j) = str_xx(i,j) + d_str + + ! Keep a copy of the biharmonic contribution for backscatter parameterization + bhstr_xx(i,j) = d_str * (h(i,j,k) * CS%reduction_xx(i,j)) + enddo ; enddo + endif if (CS%biharmonic) then ! Gradient of Laplacian, for use in bi-harmonic term @@ -989,147 +1123,256 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, meke_res_fn = 1. + if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then + do J=js-1,Jeq ; do I=is-1,Ieq + sh_xy_sq = sh_xy(I,J) * sh_xy(I,J) + sh_xx_sq = 0.25 * ( & + (sh_xx(i,j) * sh_xx(i,j) + sh_xx(i+1,j+1) * sh_xx(i+1,j+1)) & + + (sh_xx(i,j+1) * sh_xx(i,j+1) + sh_xx(i+1,j) * sh_xx(i+1,j)) & + ) + Shear_mag(i,j) = sqrt(sh_xy_sq + sh_xx_sq) + enddo ; enddo + endif + do J=js-1,Jeq ; do I=is-1,Ieq - if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then - Shear_mag = sqrt(sh_xy(I,J)*sh_xy(I,J) + & - 0.25*((sh_xx(i,j)*sh_xx(i,j) + sh_xx(i+1,j+1)*sh_xx(i+1,j+1)) + & - (sh_xx(i,j+1)*sh_xx(i,j+1) + sh_xx(i+1,j)*sh_xx(i+1,j)))) + h2uq = 4.0 * (h_u(I,j) * h_u(I,j+1)) + h2vq = 4.0 * (h_v(i,J) * h_v(i+1,J)) + hq(I,J) = (2.0 * (h2uq * h2vq)) & + / (h_neglect3 + (h2uq + h2vq) * ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J)))) + enddo ; enddo + + if (CS%better_bound_Ah .or. CS%better_bound_Kh) then + do J=js-1,Jeq ; do I=is-1,Ieq + h_min = min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J)) + hrat_min(i,j) = min(1.0, h_min / (hq(I,J) + h_neglect)) + enddo ; enddo + + if (CS%better_bound_Kh) then + do J=js-1,Jeq ; do I=is-1,Ieq + visc_bound_rem(i,j) = 1.0 + enddo ; enddo endif + endif + + if (CS%no_slip) then + do J=js-1,Jeq ; do I=is-1,Ieq + if (CS%no_slip .and. (G%mask2dBu(I,J) < 0.5)) then + if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + & + (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 0.0) then + ! This is a coastal vorticity point, so modify hq and hrat_min. + + hu = G%mask2dCu(I,j) * h_u(I,j) + G%mask2dCu(I,j+1) * h_u(I,j+1) + hv = G%mask2dCv(i,J) * h_v(i,J) + G%mask2dCv(i+1,J) * h_v(i+1,J) + if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) * & + (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) == 0.0) then + ! Only one of hu and hv is nonzero, so just add them. + hq(I,J) = hu + hv + hrat_min(i,j) = 1.0 + else + ! Both hu and hv are nonzero, so take the harmonic mean. + hq(I,J) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect) + hrat_min(i,j) = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) ) + endif + endif + endif + enddo ; enddo + endif + + if (CS%Laplacian) then if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then if (CS%use_QG_Leith_visc) then - vert_vort_mag = MIN(grad_vort_mag_q(I,J) + grad_div_mag_q(I,J), 3.*grad_vort_mag_q_2d(I,J)) + do J=js-1,Jeq ; do I=is-1,Ieq + grad_vort = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J) + grad_vort_qg = 3. * grad_vort_mag_q_2d(I,J) + vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg) + enddo ; enddo else - vert_vort_mag = (grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)) + do J=js-1,Jeq ; do I=is-1,Ieq + vert_vort_mag(i,j) = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J) + enddo ; enddo endif endif - h2uq = 4.0 * h_u(I,j) * h_u(I,j+1) - h2vq = 4.0 * h_v(i,J) * h_v(i+1,J) - hq(I,J) = 2.0 * h2uq * h2vq / (h_neglect3 + (h2uq + h2vq) * & - ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J)))) - - if (CS%better_bound_Ah .or. CS%better_bound_Kh) then - hrat_min = min(1.0, min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J)) / & - (hq(I,J) + h_neglect) ) - visc_bound_rem = 1.0 - endif - if (CS%no_slip .and. (G%mask2dBu(I,J) < 0.5)) then - if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + & - (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 0.0) then - ! This is a coastal vorticity point, so modify hq and hrat_min. - - hu = G%mask2dCu(I,j) * h_u(I,j) + G%mask2dCu(I,j+1) * h_u(I,j+1) - hv = G%mask2dCv(i,J) * h_v(i,J) + G%mask2dCv(i+1,J) * h_v(i+1,J) - if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) * & - (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) == 0.0) then - ! Only one of hu and hv is nonzero, so just add them. - hq(I,J) = hu + hv - hrat_min = 1.0 - else - ! Both hu and hv are nonzero, so take the harmonic mean. - hq(I,J) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect) - hrat_min = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) ) - endif + ! Determine the Laplacian viscosity at q points, using the + ! largest value from several parameterizations. + + ! Static (pre-computed) background viscosity + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = CS%Kh_bg_xy(i,j) + enddo ; enddo + + if (CS%Smagorinsky_Kh) then + if (CS%add_LES_viscosity) then + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = Kh(i,j) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = max(Kh(i,j), CS%Laplac2_const_xy(I,J) * Shear_mag(i,j) ) + enddo ; enddo endif endif - if (CS%Laplacian) then - ! Determine the Laplacian viscosity at q points, using the - ! largest value from several parameterizations. - Kh = CS%Kh_bg_xy(i,j) ! Static (pre-computed) background viscosity + if (CS%Leith_Kh) then if (CS%add_LES_viscosity) then - if (CS%Smagorinsky_Kh) Kh = Kh + CS%Laplac2_const_xx(i,j) * Shear_mag - if (CS%Leith_Kh) Kh = Kh + CS%Laplac3_const_xx(i,j) * vert_vort_mag*inv_PI3 + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = Kh(i,j) + CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3 + enddo ; enddo else - if (CS%Smagorinsky_Kh) Kh = max( Kh, CS%Laplac2_const_xy(I,J) * Shear_mag ) - if (CS%Leith_Kh) Kh = max( Kh, CS%Laplac3_const_xy(I,J) * vert_vort_mag*inv_PI3) + do J=js-1,Jeq ; do I=is-1,Ieq + Kh(i,j) = max(Kh(i,j), CS%Laplac3_const_xy(I,J) * vert_vort_mag(i,j) * inv_PI3) + enddo ; enddo endif - ! All viscosity contributions above are subject to resolution scaling - if (rescale_Kh) Kh = VarMix%Res_fn_q(i,j) * Kh - if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(i,j) + endif + + ! All viscosity contributions above are subject to resolution scaling + + ! NOTE: The following do-block can be decomposed and vectorized after the + ! stack size has been reduced. + do J=js-1,Jeq ; do I=is-1,Ieq + if (rescale_Kh) & + Kh(i,j) = VarMix%Res_fn_q(i,j) * Kh(i,j) + + if (CS%res_scale_MEKE) & + meke_res_fn = VarMix%Res_fn_q(i,j) + ! Older method of bounding for stability - if (legacy_bound) Kh = min(Kh, CS%Kh_Max_xy(i,j)) - Kh = max( Kh, CS%Kh_bg_min ) ! Place a floor on the viscosity, if desired. - if (use_MEKE_Ku) then ! *Add* the MEKE contribution (might be negative) - Kh = Kh + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + & + if (legacy_bound) & + Kh(i,j) = min(Kh(i,j), CS%Kh_Max_xy(i,j)) + + Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired. + + if (use_MEKE_Ku) then + ! *Add* the MEKE contribution (might be negative) + Kh(i,j) = Kh(i,j) + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + & (MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) ) * meke_res_fn endif + ! Older method of bounding for stability - if (CS%anisotropic) Kh = Kh + CS%Kh_aniso * CS%n1n2_q(I,J)**2 ! *Add* the shear component - ! of anisotropic viscosity + if (CS%anisotropic) & + ! *Add* the shear component of anisotropic viscosity + Kh(i,j) = Kh(i,j) + CS%Kh_aniso * CS%n1n2_q(I,J)**2 ! Newer method of bounding for stability if (CS%better_bound_Kh) then - if (Kh >= hrat_min*CS%Kh_Max_xy(I,J)) then - visc_bound_rem = 0.0 - Kh = hrat_min*CS%Kh_Max_xy(I,J) + if (Kh(i,j) >= hrat_min(i,j) * CS%Kh_Max_xy(I,J)) then + visc_bound_rem(i,j) = 0.0 + Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xy(I,J) elseif (CS%Kh_Max_xy(I,J)>0.) then - visc_bound_rem = 1.0 - Kh / (hrat_min*CS%Kh_Max_xy(I,J)) + visc_bound_rem(i,j) = 1.0 - Kh(i,j) / (hrat_min(i,j) * CS%Kh_Max_xy(I,J)) endif endif - if (CS%id_Kh_q>0 .or. CS%debug) Kh_q(I,J,k) = Kh - if (CS%id_vort_xy_q>0) vort_xy_q(I,J,k) = vort_xy(I,J) - if (CS%id_sh_xy_q>0) sh_xy_q(I,J,k) = sh_xy(I,J) + if (CS%id_Kh_q>0 .or. CS%debug) & + Kh_q(I,J,k) = Kh(i,j) - str_xy(I,J) = -Kh * sh_xy(I,J) - else ! not Laplacian - str_xy(I,J) = 0.0 - endif ! Laplacian + if (CS%id_vort_xy_q>0) & + vort_xy_q(I,J,k) = vort_xy(I,J) - if (CS%anisotropic) then + if (CS%id_sh_xy_q>0) & + sh_xy_q(I,J,k) = sh_xy(I,J) + enddo ; enddo + + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = -Kh(i,j) * sh_xy(I,J) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + str_xy(I,J) = 0. + enddo ; enddo + endif + + if (CS%anisotropic) then + do J=js-1,Jeq ; do I=is-1,Ieq ! Horizontal-tension averaged to q-points local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) ) ! *Add* the tension contribution to the xy-component of stress str_xy(I,J) = str_xy(I,J) - CS%Kh_aniso * CS%n1n2_q(i,j) * CS%n1n1_m_n2n2_q(i,j) * local_strain - endif + enddo ; enddo + endif - if (CS%biharmonic) then + if (CS%biharmonic) then ! Determine the biharmonic viscosity at q points, using the ! largest value from several parameterizations. - AhSm = 0.0 ; AhLth = 0.0 - if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then - if (CS%Smagorinsky_Ah) then - if (CS%bound_Coriolis) then - AhSm = Shear_mag * (CS%Biharm_const_xy(I,J) + & - CS%Biharm_const2_xy(I,J)*Shear_mag) - else - AhSm = CS%Biharm_const_xy(I,J) * Shear_mag - endif - endif - if (CS%Leith_Ah) AhLth = CS%Biharm6_const_xy(I,J) * abs(Del2vort_q(I,J)) * inv_PI6 - Ah = MAX(MAX(CS%Ah_bg_xy(I,J), AhSm), AhLth) - if (CS%bound_Ah .and. .not.CS%better_bound_Ah) & - Ah = MIN(Ah, CS%Ah_Max_xy(I,J)) - else - Ah = CS%Ah_bg_xy(I,J) - endif ! Smagorinsky_Ah or Leith_Ah + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = CS%Ah_bg_xy(I,J) + enddo ; enddo - if (use_MEKE_Au) then ! *Add* the MEKE contribution - Ah = Ah + 0.25*( (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + & - (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) ) + if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then + if (CS%Smagorinsky_Ah) then + if (CS%bound_Coriolis) then + do J=js-1,Jeq ; do I=is-1,Ieq + AhSm = Shear_mag(i,j) * (CS%Biharm_const_xy(I,J) & + + CS%Biharm_const2_xy(I,J) * Shear_mag(i,j) & + ) + Ah(i,j) = max(Ah(I,J), AhSm) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + AhSm = CS%Biharm_const_xy(I,J) * Shear_mag(i,j) + Ah(i,j) = max(Ah(I,J), AhSm) + enddo ; enddo + endif endif - if (CS%Re_Ah > 0.0) then - KE = 0.125*((u(I,j,k)+u(I,j+1,k))**2 + (v(i,J,k)+v(i+1,J,k))**2) - Ah = sqrt(KE) * CS%Re_Ah_const_xy(i,j) + if (CS%Leith_Ah) then + do J=js-1,Jeq ; do I=is-1,Ieq + AhLth = CS%Biharm6_const_xy(I,J) * abs(Del2vort_q(I,J)) * inv_PI6 + Ah(i,j) = max(Ah(I,J), AhLth) + enddo ; enddo endif - if (CS%better_bound_Ah) then - Ah = MIN(Ah, visc_bound_rem*hrat_min*CS%Ah_Max_xy(I,J)) + if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xy(I,J)) + enddo ; enddo endif + endif ! Smagorinsky_Ah or Leith_Ah - if (CS%id_Ah_q>0 .or. CS%debug) Ah_q(I,J,k) = Ah + if (use_MEKE_Au) then + ! *Add* the MEKE contribution + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = Ah(i,j) + 0.25 * ( & + (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) & + ) + enddo ; enddo + endif - str_xy(I,J) = str_xy(I,J) + Ah * ( dDel2vdx(I,J) + dDel2udy(I,J) ) + if (CS%Re_Ah > 0.0) then + do J=js-1,Jeq ; do I=is-1,Ieq + KE = 0.125 * ((u(I,j,k) + u(I,j+1,k))**2 + (v(i,J,k) + v(i+1,J,k))**2) + Ah(i,j) = sqrt(KE) * CS%Re_Ah_const_xy(i,j) + enddo ; enddo + endif - ! Keep a copy of the biharmonic contribution for backscatter parameterization - bhstr_xy(I,J) = Ah * ( dDel2vdx(I,J) + dDel2udy(I,J) ) * & - (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) + if (CS%better_bound_Ah) then + if (CS%better_bound_Kh) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xy(I,J)) + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq + Ah(i,j) = min(Ah(i,j), hrat_min(i,j) * CS%Ah_Max_xy(I,J)) + enddo ; enddo + endif + endif - endif ! biharmonic + if (CS%id_Ah_q>0 .or. CS%debug) then + do J=js-1,Jeq ; do I=is-1,Ieq + Ah_q(I,J,k) = Ah(i,j) + enddo ; enddo + endif - enddo ; enddo + ! Again, need to initialize str_xy as if its biharmonic + do J=js-1,Jeq ; do I=is-1,Ieq + d_str = Ah(i,j) * (dDel2vdx(I,J) + dDel2udy(I,J)) + + str_xy(I,J) = str_xy(I,J) + d_str + + ! Keep a copy of the biharmonic contribution for backscatter parameterization + bhstr_xy(I,J) = d_str * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) + enddo ; enddo + endif if (CS%use_GME) then call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G, GV) @@ -1197,14 +1440,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j)) enddo ; enddo - do J=js-1,Jeq ; do I=is-1,Ieq - if (CS%no_slip) then + if (CS%no_slip) then + do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = str_xy(I,J) * (hq(I,J) * CS%reduction_xy(I,J)) - else + enddo ; enddo + else + do J=js-1,Jeq ; do I=is-1,Ieq str_xy(I,J) = str_xy(I,J) * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J)) - endif - enddo ; enddo - + enddo ; enddo + endif endif ! use_GME ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent. @@ -1214,8 +1458,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, G%IdxCu(I,j)*(CS%dx2q(I,J-1)*str_xy(I,J-1) - & CS%dx2q(I,J) *str_xy(I,J))) * & G%IareaCu(I,j)) / (h_u(I,j) + h_neglect) - enddo ; enddo + if (apply_OBC) then ! This is not the right boundary condition. If all the masking of tendencies are done ! correctly later then eliminating this block should not change answers. @@ -1237,6 +1481,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, CS%dx2h(i,j+1)*str_xx(i,j+1))) * & G%IareaCv(i,J)) / (h_v(i,J) + h_neglect) enddo ; enddo + if (apply_OBC) then ! This is not the right boundary condition. If all the masking of tendencies are done ! correctly later then eliminating this block should not change answers. @@ -1288,28 +1533,27 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=js,je ; do i=is,ie FatH = 0.25*( (abs(G%CoriolisBu(I-1,J-1)) + abs(G%CoriolisBu(I,J))) + & (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1))) ) - Shear_mag = sqrt(sh_xx(i,j)*sh_xx(i,j) + & + Shear_mag_bc = sqrt(sh_xx(i,j) * sh_xx(i,j) + & 0.25*((sh_xy(I-1,J-1)*sh_xy(I-1,J-1) + sh_xy(I,J)*sh_xy(I,J)) + & (sh_xy(I-1,J)*sh_xy(I-1,J) + sh_xy(I,J-1)*sh_xy(I,J-1)))) if (CS%answers_2018) then FatH = (US%s_to_T*FatH)**MEKE%backscatter_Ro_pow ! f^n ! Note the hard-coded dimensional constant in the following line that can not ! be rescaled for dimensional consistency. - Shear_mag = ( ( (US%s_to_T*Shear_mag)**MEKE%backscatter_Ro_pow ) + 1.e-30 ) & + Shear_mag_bc = (((US%s_to_T * Shear_mag_bc)**MEKE%backscatter_Ro_pow) + 1.e-30) & * MEKE%backscatter_Ro_c ! c * D^n ! The Rossby number function is g(Ro) = 1/(1+c.Ro^n) ! RoScl = 1 - g(Ro) - RoScl = Shear_mag / ( FatH + Shear_mag ) ! = 1 - f^n/(f^n+c*D^n) + RoScl = Shear_mag_bc / (FatH + Shear_mag_bc) ! = 1 - f^n/(f^n+c*D^n) else - if (FatH <= backscat_subround*Shear_mag) then + if (FatH <= backscat_subround*Shear_mag_bc) then RoScl = 1.0 else - Sh_F_pow = MEKE%backscatter_Ro_c * (Shear_mag / FatH)**MEKE%backscatter_Ro_pow + Sh_F_pow = MEKE%backscatter_Ro_c * (Shear_mag_bc / FatH)**MEKE%backscatter_Ro_pow RoScl = Sh_F_pow / (1.0 + Sh_F_pow) ! = 1 - f^n/(f^n+c*D^n) endif endif - MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + GV%H_to_RZ * ( & ((str_xx(i,j)-RoScl*bhstr_xx(i,j))*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j) & -(str_xx(i,j)-RoScl*bhstr_xx(i,j))*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j)) & diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 470098a08a..8ba9dd959a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -15,15 +15,15 @@ module MOM_diabatic_aux use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : init_external_field, time_interp_extern +use MOM_interpolate, only : time_interp_external_init use MOM_io, only : slasher use MOM_opacity, only : set_opacity, opacity_CS, extract_optics_slice, extract_optics_fields use MOM_opacity, only : optics_type, optics_nbands, absorbRemainingSW, sumSWoverBands use MOM_tracer_flow_control, only : get_chl_from_model, tracer_flow_control_CS use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs ! , vertvisc_type, accel_diag_ptrs +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init implicit none ; private @@ -621,7 +621,7 @@ subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity_CSp, tracer_ if (CS%chl_from_file) then ! Only the 2-d surface chlorophyll can be read in from a file. The ! same value is assumed for all layers. - call time_interp_external(CS%sbc_chl, CS%Time, chl_2d) + call time_interp_extern(CS%sbc_chl, CS%Time, chl_2d) do j=js,je ; do i=is,ie if ((G%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&