From df5cb9b5d199ba6b517fcada9113563f93c976ee Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 31 May 2023 15:18:27 -0800 Subject: [PATCH 01/14] Salt data structures --- config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 | 1 + docs/zotero.bib | 9 +++++++++ src/core/MOM_forcing_type.F90 | 8 +++++++- 3 files changed, 17 insertions(+), 1 deletion(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index 26ab6269ef..c04cf1750b 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -172,6 +172,7 @@ module MOM_surface_forcing_gfdl real, pointer, dimension(:,:) :: t_flux =>NULL() !< sensible heat flux [W m-2] real, pointer, dimension(:,:) :: q_flux =>NULL() !< specific humidity flux [kg m-2 s-1] real, pointer, dimension(:,:) :: salt_flux =>NULL() !< salt flux [kg m-2 s-1] + real, pointer, dimension(:,:) :: excess_salt =>NULL() !< salt left behind by brine rejection [kg m-2 s-1] real, pointer, dimension(:,:) :: lw_flux =>NULL() !< long wave radiation [W m-2] real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() !< direct visible sw radiation [W m-2] real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() !< diffuse visible sw radiation [W m-2] diff --git a/docs/zotero.bib b/docs/zotero.bib index c0c7ee3bd9..5acaee968a 100644 --- a/docs/zotero.bib +++ b/docs/zotero.bib @@ -2738,3 +2738,12 @@ @article{kraus1967 journal = {Tellus} } +@article{Nguyen2009, + doi = {10.1029/2008JC005121}, + year = {2009}, + journal = {JGR Oceans}, + volume = {114}, + author = {A. T. Nguyen and D. Menemenlis and R. Kwok}, + title = {Improved modeling of the Arctic halocline with a subgrid-scale brine rejection parameterization}, + pages = {C11014} +} diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 9a5e1f48f5..e6d4efd424 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -131,8 +131,10 @@ module MOM_forcing_type real, pointer, dimension(:,:) :: & salt_flux => NULL(), & !< net salt flux into the ocean [R Z T-1 ~> kgSalt m-2 s-1] salt_flux_in => NULL(), & !< salt flux provided to the ocean from coupler [R Z T-1 ~> kgSalt m-2 s-1] - salt_flux_added => NULL() !< additional salt flux from restoring or flux adjustment before adjustment + salt_flux_added => NULL(), & !< additional salt flux from restoring or flux adjustment before adjustment !! to net zero [R Z T-1 ~> kgSalt m-2 s-1] + salt_left_behind => NULL() !< salt left in ocean at the surface from brine rejection + !! [R Z T-1 ~> kgSalt m-2 s-1] ! applied surface pressure from other component models (e.g., atmos, sea ice, land ice) real, pointer, dimension(:,:) :: p_surf_full => NULL() @@ -1929,6 +1931,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, diag%axesT1,Time,'Salt flux into ocean at surface due to restoring or flux adjustment', & units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) + handles%id_saltFluxAdded = register_diag_field('ocean_model', 'salt_left_behind', & + diag%axesT1,Time,'Salt left in ocean at surface due to ice formation', & + units='kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) + handles%id_saltFluxGlobalAdj = register_scalar_field('ocean_model', & 'salt_flux_global_restoring_adjustment', Time, diag, & 'Adjustment needed to balance net global salt flux into ocean at surface', & From 1a5e9ad0d0567aab67f9fb37540b77db9c971697 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 5 Jun 2023 09:25:40 -0800 Subject: [PATCH 02/14] First steps at brine plume: pass info from SIS2 --- .../drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index c04cf1750b..a9724c5964 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -61,7 +61,7 @@ module MOM_surface_forcing_gfdl !! from MOM_domains) to indicate the staggering of !! the winds that are being provided in calls to !! update_ocean_model. - logical :: use_temperature !< If true, temp and saln used as state variables + logical :: use_temperature !< If true, temp and saln used as state variables. real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress [nondim]. real :: Rho0 !< Boussinesq reference density [R ~> kg m-3] @@ -302,6 +302,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) + if (associated(IOB%excess_salt)) call safe_alloc_ptr(fluxes%salt_left_behind,isd,ied,jsd,jed) + do j=js-2,je+2 ; do i=is-2,ie+2 fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) @@ -574,6 +576,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, call check_mask_val_consistency(IOB%salt_flux(i-i0,j-j0), G%mask2dT(i,j), i, j, 'salt_flux', G) enddo ; enddo endif + if (associated(IOB%excess_salt)) then + do j=js,je ; do i=is,ie + fluxes%salt_left_behind(i,j) = G%mask2dT(i,j)*( - kg_m2_s_conversion*IOB%excess_salt(i-i0,j-j0)) + enddo ; enddo + endif !#CTRL# if (associated(CS%ctrl_forcing_CSp)) then !#CTRL# do j=js,je ; do i=is,ie @@ -1713,6 +1720,9 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) if (associated(iobt%mass_berg)) then chks = field_chksum( iobt%mass_berg ) ; if (root) write(outunit,100) 'iobt%mass_berg ', chks endif + if (associated(iobt%excess_salt)) then + chks = field_chksum( iobt%excess_salt ) ; if (root) write(outunit,100) 'iobt%excess_salt ', chks + endif 100 FORMAT(" CHECKSUM::",A20," = ",Z20) call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') From b1cf0d20f835bf0e12093dbc3c1422a24fee75d9 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 31 May 2023 14:05:39 -0400 Subject: [PATCH 03/14] FMS2 interpolation ID replaced with derived type All instances of an FMS ID to the internal interpolation content is replaced with a derived type containing additional metadata recording the field's origin filename and fieldname. This additional information is required in order to replicate the axis data from the field, which is no longer provided by FMS2. The abstraction of this type also allows us to either extend it or redefine it in other frameworks as needed in the future. This primarily affects the usage of the following functions: - init_external_field - time_interp_external - horiz_interp_and_extrap_tracer The following solvers are updated: - MOM_open_boundary - MOM_ice_shelf - MOM_oda_driver - MOM_MEKE - MOM_ALE_sponge - MOM_diabatic_aux Of these, OBC was the most significant. The integer handle (fid) was previously used to determine if each segment field was constant or (if negative) read from a file. After being replaced by the derived type, a new flag was added to make this determination. All of the coupled drivers have been modified, since they support time interpolation of T and S fields. - FMS - MCT - NUOPC The NUOPC driver also includes modifications to its CFC11 and CFC12 fields. Changes to the MOM CFC modules replaces an `id == -1`-like test, which is not used by the derived type. This check has been removed, and we now solely rely on the `present(cfc_handle)` test. While this could change behavior, there does not seem to be any scenario where init_external_field would return -1 but would be passed to the function. (But I may eat these words.) --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 15 +++-- .../mct_cap/mom_surface_forcing_mct.F90 | 15 +++-- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 28 +++++---- config_src/infra/FMS1/MOM_interp_infra.F90 | 54 +++++++++++------- config_src/infra/FMS2/MOM_interp_infra.F90 | 57 +++++++++++-------- src/core/MOM_open_boundary.F90 | 31 +++++----- src/framework/MOM_horizontal_regridding.F90 | 11 ++-- src/framework/MOM_interpolate.F90 | 27 +++++---- src/ice_shelf/MOM_ice_shelf.F90 | 17 +++--- src/ocean_data_assim/MOM_oda_driver.F90 | 15 ++--- src/parameterizations/lateral/MOM_MEKE.F90 | 7 ++- .../vertical/MOM_ALE_sponge.F90 | 32 +++++------ .../vertical/MOM_diabatic_aux.F90 | 3 +- src/tracer/MOM_CFC_cap.F90 | 20 ++++--- 14 files changed, 187 insertions(+), 145 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index 26ab6269ef..f70cd34012 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -27,6 +27,7 @@ module MOM_surface_forcing_gfdl use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : read_netCDF_data use MOM_io, only : stdout_if_root @@ -153,8 +154,10 @@ module MOM_surface_forcing_gfdl !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' real, pointer, dimension(:,:) :: trestore_mask => NULL() !< Mask for SST restoring [nondim] - integer :: id_srestore = -1 !< An id number for time_interp_external. - integer :: id_trestore = -1 !< An id number for time_interp_external. + type(external_field) :: srestore_handle + !< Handle for time-interpolated salt restoration field + type(external_field) :: trestore_handle + !< Handle for time-interpolated temperature restoration field type(forcing_diags), public :: handles !< Diagnostics handles @@ -345,7 +348,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, scale=US%ppt_to_S) + call time_interp_external(CS%srestore_handle, Time, data_restore, scale=US%ppt_to_S) ! 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 @@ -403,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, scale=US%degC_to_C) + call time_interp_external(CS%trestore_handle, Time, data_restore, scale=US%degC_to_C) if ( CS%trestore_SPEAR_ECDA ) then do j=js,je ; do i=is,ie if (abs(data_restore(i,j)+1.8*US%degC_to_C) < 0.0001*US%degC_to_C) then @@ -1610,7 +1613,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) if (CS%restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) - CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, MOM_domain=G%Domain) + CS%srestore_handle = init_external_field(salt_file, CS%salt_restore_var_name, MOM_domain=G%Domain) call safe_alloc_ptr(CS%srestore_mask,isd,ied,jsd,jed); CS%srestore_mask(:,:) = 1.0 if (CS%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'salt_restore_mask.nc' @@ -1620,7 +1623,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) if (CS%restore_temp) then temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) - CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, MOM_domain=G%Domain) + CS%trestore_handle = init_external_field(temp_file, CS%temp_restore_var_name, MOM_domain=G%Domain) call safe_alloc_ptr(CS%trestore_mask,isd,ied,jsd,jed); CS%trestore_mask(:,:) = 1.0 if (CS%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'temp_restore_mask.nc' diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 index 0364d46ddc..9b858af94e 100644 --- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 @@ -25,6 +25,7 @@ module MOM_surface_forcing_mct use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : stdout use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS @@ -134,8 +135,10 @@ module MOM_surface_forcing_mct !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. + type(external_field) :: srestore_handle + !< Handle for time-interpolated salt restoration field + type(external_field) :: trestore_handle + !< Handle for time-interpolated temperature restoration field type(forcing_diags), public :: handles !< diagnostics handles type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< restart pointer @@ -348,7 +351,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (restore_salinity) then - call time_interp_external(CS%id_srestore, Time, data_restore, scale=US%ppt_to_S) + call time_interp_external(CS%srestore_handle, Time, data_restore, scale=US%ppt_to_S) ! 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 @@ -405,7 +408,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (restore_sst) then - call time_interp_external(CS%id_trestore, Time, data_restore, scale=US%degC_to_C) + call time_interp_external(CS%trestore_handle, Time, data_restore, scale=US%degC_to_C) 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) @@ -1292,7 +1295,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_salt)) then ; if (restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) - CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) + CS%srestore_handle = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%srestore_mask,isd,ied,jsd,jed); CS%srestore_mask(:,:) = 1.0 if (CS%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'salt_restore_mask.nc' @@ -1302,7 +1305,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_temp)) then ; if (restore_temp) then temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) - CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) + CS%trestore_handle = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%trestore_mask,isd,ied,jsd,jed); CS%trestore_mask(:,:) = 1.0 if (CS%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'temp_restore_mask.nc' diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 2c8e3db8bd..b8162b4b59 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -26,6 +26,7 @@ module MOM_surface_forcing_nuopc use MOM_grid, only : ocean_grid_type use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_CFC_cap, only : CFC_cap_fluxes use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_io, only : stdout @@ -146,10 +147,14 @@ module MOM_surface_forcing_nuopc character(len=30) :: cfc11_var_name !< name of cfc11 in CFC_BC_file character(len=30) :: cfc12_var_name !< name of cfc11 in CFC_BC_file real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. - integer :: id_cfc11_atm = -1 !< id number for time_interp_external. - integer :: id_cfc12_atm = -1 !< id number for time_interp_external. + type(external_field) :: srestore_handle + !< Handle for time-interpolated salt restoration field + type(external_field) :: trestore_handle + !< Handle for time-interpolated temperature restoration field + type(external_field) :: cfc11_atm_handle + !< Handle for time-interpolated CFC11 restoration field + type(external_field) :: cfc12_atm_handle + !< Handle for time-interpolated CFC12 restoration field ! Diagnostics handles type(forcing_diags), public :: handles @@ -377,7 +382,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (restore_salinity) then - call time_interp_external(CS%id_srestore, Time, data_restore, scale=US%ppt_to_S) + call time_interp_external(CS%srestore_handle, Time, data_restore, scale=US%ppt_to_S) ! 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 @@ -434,7 +439,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (restore_sst) then - call time_interp_external(CS%id_trestore, Time, data_restore, scale=US%degC_to_C) + call time_interp_external(CS%trestore_handle, Time, data_restore, scale=US%degC_to_C) 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) @@ -596,7 +601,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! CFCs if (CS%use_CFC) then - call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) + call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, & + CS%cfc11_atm_handle, CS%cfc11_atm_handle) endif if (associated(IOB%salt_flux)) then @@ -1394,7 +1400,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_salt)) then ; if (restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) - CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) + CS%srestore_handle = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%srestore_mask,isd,ied,jsd,jed); CS%srestore_mask(:,:) = 1.0 if (CS%mask_srestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'salt_restore_mask.nc' @@ -1404,7 +1410,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, if (present(restore_temp)) then ; if (restore_temp) then temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) - CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) + CS%trestore_handle = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%trestore_mask,isd,ied,jsd,jed); CS%trestore_mask(:,:) = 1.0 if (CS%mask_trestore) then ! read a 2-d file containing a mask for restoring fluxes flnam = trim(CS%inputdir) // 'temp_restore_mask.nc' @@ -1430,8 +1436,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "The name of the variable representing CFC-12 in "//& "CFC_BC_FILE.", default="CFC_12", do_not_log=.true.) - CS%id_cfc11_atm = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) - CS%id_cfc12_atm = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) + CS%cfc11_atm_handle = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) + CS%cfc12_atm_handle = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) endif endif diff --git a/config_src/infra/FMS1/MOM_interp_infra.F90 b/config_src/infra/FMS1/MOM_interp_infra.F90 index 224e26a051..e14233a64b 100644 --- a/config_src/infra/FMS1/MOM_interp_infra.F90 +++ b/config_src/infra/FMS1/MOM_interp_infra.F90 @@ -18,6 +18,18 @@ module MOM_interp_infra public :: time_interp_extern, init_extern_field, time_interp_extern_init public :: get_external_field_info, axistype, get_axis_data public :: run_horiz_interp, build_horiz_interp_weights +public :: external_field + +!< Handle of an external field for interpolation +type :: external_field + private + integer :: id + !< FMS ID for the interpolated field + character(len=:), allocatable :: filename + !< Filename containing the field values + character(len=:), allocatable :: label + !< Field name in the file +end type external_field !> Read a field based on model time, and rotate to the model domain. interface time_interp_extern @@ -167,8 +179,8 @@ end function get_extern_field_missing !> 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 +subroutine get_external_field_info(field, size, axes, missing) + type(external_field), intent(in) :: field !< Handle for time interpolated external !! field returned from a previous !! call to init_external_field() integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data @@ -176,37 +188,35 @@ subroutine get_external_field_info(field_id, size, axes, missing) real, optional, intent(inout) :: missing !< Missing value for the input data if (present(size)) then - size(1:4) = get_extern_field_size(field_id) + size(1:4) = get_extern_field_size(field%id) endif if (present(axes)) then - axes(1:4) = get_extern_field_axes(field_id) + axes(1:4) = get_extern_field_axes(field%id) endif if (present(missing)) then - missing = get_extern_field_missing(field_id) + missing = get_extern_field_missing(field%id) endif end subroutine get_external_field_info !> Read a scalar field based on model time. -subroutine time_interp_extern_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() +subroutine time_interp_extern_0d(field, time, data_in, verbose) + type(external_field), intent(in) :: field !< Handle for time interpolated 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(field_id, time, data_in, verbose=verbose) + call time_interp_external(field%id, time, data_in, verbose=verbose) end subroutine time_interp_extern_0d !> Read a 2d field from an external based on model time, potentially including horizontal !! interpolation and rotation of the data -subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_2d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated 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 @@ -216,15 +226,14 @@ subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_2d !> Read a 3d field based on model time, and rotate to the model grid -subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_3d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated 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 @@ -234,14 +243,15 @@ subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_3d !> initialize an external field -integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & - threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency ) +function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & + threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency) & + result(field) character(len=*), intent(in) :: file !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -261,17 +271,17 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, !! is in use, and (2) the modulo time period of the !! data is an integer number of years, then map !! a model date of Feb 29. onto a common year on Feb. 28. + type(external_field) :: field !< Handle to external field if (present(MOM_Domain)) then - init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & + field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) else - init_extern_field = init_external_field(file, fieldname, domain=domain, & + field%id = init_external_field(file, fieldname, domain=domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) endif - end function init_extern_field end module MOM_interp_infra diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index c29459aad1..7964b3537f 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -18,6 +18,18 @@ module MOM_interp_infra public :: time_interp_extern, init_extern_field, time_interp_extern_init public :: get_external_field_info, axistype, get_axis_data public :: run_horiz_interp, build_horiz_interp_weights +public :: external_field + +!< Handle of an external field for interpolation +type :: external_field + private + integer :: id + !< FMS ID for the interpolated field + character(len=:), allocatable :: filename + !< Filename containing the field values + character(len=:), allocatable :: label + !< Field name in the file +end type external_field !> Read a field based on model time, and rotate to the model domain. interface time_interp_extern @@ -166,8 +178,8 @@ end function get_extern_field_missing !> 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 +subroutine get_external_field_info(field, size, axes, missing) + type(external_field), intent(in) :: field !< Handle for time interpolated external !! field returned from a previous !! call to init_external_field() integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data @@ -175,37 +187,35 @@ subroutine get_external_field_info(field_id, size, axes, missing) real, optional, intent(inout) :: missing !< Missing value for the input data if (present(size)) then - size(1:4) = get_extern_field_size(field_id) + size(1:4) = get_extern_field_size(field%id) endif if (present(axes)) then - axes(1:4) = get_extern_field_axes(field_id) + axes(1:4) = get_extern_field_axes(field%id) endif if (present(missing)) then - missing = get_extern_field_missing(field_id) + missing = get_extern_field_missing(field%id) endif end subroutine get_external_field_info !> Read a scalar field based on model time. -subroutine time_interp_extern_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() +subroutine time_interp_extern_0d(field, time, data_in, verbose) + type(external_field), intent(in) :: field !< Handle for time interpolated 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(field_id, time, data_in, verbose=verbose) + call time_interp_external(field%id, time, data_in, verbose=verbose) end subroutine time_interp_extern_0d !> Read a 2d field from an external based on model time, potentially including horizontal !! interpolation and rotation of the data -subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_2d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated 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 @@ -215,15 +225,14 @@ subroutine time_interp_extern_2d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_2d !> Read a 3d field based on model time, and rotate to the model grid -subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_interp, mask_out) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_extern_3d(field, time, data_in, interp, verbose, horz_interp, mask_out) + type(external_field), intent(in) :: field !< Handle for time interpolated 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 @@ -233,14 +242,15 @@ subroutine time_interp_extern_3d(field_id, time, data_in, interp, verbose, horz_ logical, dimension(:,:,:), & optional, intent(out) :: mask_out !< An array that is true where there is valid data - call time_interp_external(field_id, time, data_in, interp=interp, verbose=verbose, & + call time_interp_external(field%id, time, data_in, interp=interp, verbose=verbose, & horz_interp=horz_interp, mask_out=mask_out) end subroutine time_interp_extern_3d !> initialize an external field -integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & - threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency ) +function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & + threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency) & + result(field) character(len=*), intent(in) :: file !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -260,19 +270,20 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, !! is in use, and (2) the modulo time period of the !! data is an integer number of years, then map !! a model date of Feb 29. onto a common year on Feb. 28. + type(external_field) :: field !< Handle to external field - + field%filename = file + field%label = fieldname if (present(MOM_Domain)) then - init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & + field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) else - init_extern_field = init_external_field(file, fieldname, domain=domain, & + field%id = init_external_field(file, fieldname, domain=domain, & verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency) endif - end function init_extern_field end module MOM_interp_infra diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9bd292e796..ba8b8ce818 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -24,6 +24,7 @@ module MOM_open_boundary 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 MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init +use MOM_interpolate, only : external_field 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 @@ -81,8 +82,9 @@ module MOM_open_boundary !> Open boundary segment data from files (mostly). type, public :: OBC_segment_data_type - integer :: fid !< handle from FMS associated with segment data on disk - integer :: fid_dz !< handle from FMS associated with segment thicknesses on disk + type(external_field) :: handle !< handle from FMS associated with segment data on disk + type(external_field) :: dz_handle !< handle from FMS associated with segment thicknesses on disk + logical :: use_IO = .false. !< True if segment data is based on file input character(len=32) :: name !< a name identifier for the segment data character(len=8) :: genre !< an identifier for the segment data real :: scale !< A scaling factor for converting input data to @@ -96,7 +98,7 @@ module MOM_open_boundary real, allocatable :: buffer_dst(:,:,:) !< buffer src data remapped to the target vertical grid. !! The values for tracers should have the same units as the field !! they are being applied to? - real :: value !< constant value if fid is equal to -1 + real :: value !< constant value if not read from file real :: resrv_lfac_in = 1. !< reservoir inverse length scale factor for IN direction per field !< the general 1/Lscale_IN is multiplied by this factor for each tracer real :: resrv_lfac_out= 1. !< reservoir inverse length scale factor for OUT direction per field @@ -842,6 +844,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) ! The scale factor for tracers may also be set in register_segment_tracer, and a constant input ! value is rescaled there. segment%field(m)%scale = scale_factor_from_name(fields(m), GV, US, segment%tr_Reg) + segment%field(m)%use_IO = .true. if (segment%field(m)%name == 'TEMP') then segment%temp_segment_data_exists = .true. segment%t_values_needed = .false. @@ -957,7 +960,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) endif endif segment%field(m)%buffer_src(:,:,:) = 0.0 - segment%field(m)%fid = init_external_field(trim(filename), trim(fieldname), & + segment%field(m)%handle = init_external_field(trim(filename), trim(fieldname), & ignore_axis_atts=.true., threading=SINGLE_FILE) if (siz(3) > 1) then if ((index(segment%field(m)%name, 'phase') > 0) .or. (index(segment%field(m)%name, 'amp') > 0)) then @@ -988,7 +991,7 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) endif segment%field(m)%dz_src(:,:,:)=0.0 segment%field(m)%nk_src=siz(3) - segment%field(m)%fid_dz = init_external_field(trim(filename), trim(fieldname), & + segment%field(m)%dz_handle = init_external_field(trim(filename), trim(fieldname), & ignore_axis_atts=.true., threading=SINGLE_FILE) endif else @@ -996,12 +999,12 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF) endif endif else - segment%field(m)%fid = -1 segment%field(m)%name = trim(fields(m)) ! The scale factor for tracers may also be set in register_segment_tracer, and a constant input ! value is rescaled there. segment%field(m)%scale = scale_factor_from_name(fields(m), GV, US, segment%tr_Reg) segment%field(m)%value = segment%field(m)%scale * value + segment%field(m)%use_IO = .false. ! Check if this is a tidal field. If so, the number ! of expected constituents must be 1. @@ -3892,7 +3895,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) !a less frequent update as set by the parameter update_OBC_period_max in MOM.F90. !Cycle if it is not the time to update OBC segment data for this field. if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle - if (segment%field(m)%fid > 0) then + if (segment%field(m)%use_IO) then siz(1)=size(segment%field(m)%buffer_src,1) siz(2)=size(segment%field(m)%buffer_src,2) siz(3)=size(segment%field(m)%buffer_src,3) @@ -3972,7 +3975,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif ! This is where the data values are actually read in. - call time_interp_external(segment%field(m)%fid, Time, tmp_buffer_in, scale=segment%field(m)%scale) + call time_interp_external(segment%field(m)%handle, Time, tmp_buffer_in, scale=segment%field(m)%scale) ! NOTE: Rotation of face-points require that we skip the final value if (turns /= 0) then @@ -4045,7 +4048,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (segment%field(m)%nk_src > 1 .and.& (index(segment%field(m)%name, 'phase') <= 0 .and. index(segment%field(m)%name, 'amp') <= 0)) then ! This is where the 2-d tidal data values are actually read in. - call time_interp_external(segment%field(m)%fid_dz, Time, tmp_buffer_in, scale=US%m_to_Z) + call time_interp_external(segment%field(m)%dz_handle, Time, tmp_buffer_in, scale=US%m_to_Z) if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. if (segment%is_E_or_W & @@ -4211,7 +4214,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) deallocate(tmp_buffer) if (turns /= 0) & deallocate(tmp_buffer_in) - else ! fid <= 0 (Uniform value) + else ! use_IO = .false. (Uniform value) if (.not. allocated(segment%field(m)%buffer_dst)) then if (segment%is_E_or_W) then if (segment%field(m)%name == 'V') then @@ -4257,7 +4260,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) do m = 1,segment%num_fields !cycle if it is not the time to update OBGC tracers from source if (trim(segment%field(m)%genre) == 'obgc' .and. (.not. OBC%update_OBC_seg_data)) cycle - ! if (segment%field(m)%fid>0) then + ! if (segment%field(m)%use_IO) then ! calculate external BT velocity and transport if needed if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then if (trim(segment%field(m)%name) == 'U' .and. segment%is_E_or_W) then @@ -4684,7 +4687,7 @@ subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, & ! rescale the previously stored input values. Note that calls to register_segment_tracer ! can come before or after calls to initialize_segment_data. if (uppercase(segment%field(m)%name) == uppercase(segment%tr_Reg%Tr(ntseg)%name)) then - if (segment%field(m)%fid == -1) then + if (.not. segment%field(m)%use_IO) then rescale = scale if ((segment%field(m)%scale /= 0.0) .and. (segment%field(m)%scale /= 1.0)) & rescale = scale / segment%field(m)%scale @@ -5948,8 +5951,8 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%num_fields = segment_in%num_fields do n = 1, num_fields - segment%field(n)%fid = segment_in%field(n)%fid - segment%field(n)%fid_dz = segment_in%field(n)%fid_dz + segment%field(n)%handle = segment_in%field(n)%handle + segment%field(n)%dz_handle = segment_in%field(n)%dz_handle if (modulo(turns, 2) /= 0) then select case (segment_in%field(n)%name) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 83e7718311..bedf710582 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -17,6 +17,7 @@ module MOM_horizontal_regridding use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights use MOM_interp_infra, only : horiz_interp_type, horizontal_interp_init use MOM_interp_infra, only : axistype, get_external_field_info, get_axis_data +use MOM_interp_infra, only : external_field use MOM_time_manager, only : time_type use MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data use MOM_io, only : read_attribute, read_variable @@ -598,12 +599,12 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr end subroutine horiz_interp_and_extrap_tracer_record !> Extrapolate and interpolate using a FMS time interpolation handle -subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, & +subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & z_in, z_edges_in, missing_value, scale, & homogenize, spongeOngrid, m_to_Z, & answers_2018, tr_iter_tol, answer_date) - integer, intent(in) :: fms_id !< A unique id used by the FMS time interpolator + type(external_field), intent(in) :: field !< Handle for the time interpolated field type(time_type), intent(in) :: Time !< A FMS time type type(ocean_grid_type), intent(inout) :: G !< Grid object real, allocatable, dimension(:,:,:), intent(out) :: tr_z @@ -716,7 +717,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, call cpu_clock_begin(id_clock_read) - call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_val_in) + call get_external_field_info(field, size=fld_sz, axes=axes_data, missing=missing_val_in) missing_value = scale*missing_val_in verbosity = MOM_get_verbosity() @@ -790,7 +791,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, if (.not.is_ongrid) then if (is_root_pe()) & - call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>5), turns=turns) + call time_interp_external(field, Time, data_in, verbose=(verbosity>5), turns=turns) ! Loop through each data level and interpolate to model grid. ! After interpolating, fill in points which will be needed to define the layers. @@ -897,7 +898,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, G, tr_z, mask_z, enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=(verbosity>5), turns=turns) + call time_interp_external(field, Time, data_in, verbose=(verbosity>5), turns=turns) do k=1,kd do j=js,je do i=is,ie diff --git a/src/framework/MOM_interpolate.F90 b/src/framework/MOM_interpolate.F90 index 38a786e593..e131e8db9d 100644 --- a/src/framework/MOM_interpolate.F90 +++ b/src/framework/MOM_interpolate.F90 @@ -9,12 +9,14 @@ module MOM_interpolate use MOM_interp_infra, only : time_interp_external_init=>time_interp_extern_init use MOM_interp_infra, only : horiz_interp_type, get_external_field_info use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights +use MOM_interp_infra, only : external_field use MOM_time_manager, only : time_type implicit none ; private public :: time_interp_external, init_external_field, time_interp_external_init, get_external_field_info public :: horiz_interp_type, run_horiz_interp, build_horiz_interp_weights +public :: external_field !> Read a field based on model time, and rotate to the model domain. interface time_interp_external @@ -26,9 +28,8 @@ module MOM_interpolate contains !> Read a scalar field based on model time. -subroutine time_interp_external_0d(field_id, time, data_in, verbose, scale) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() +subroutine time_interp_external_0d(field, time, data_in, verbose, scale) + type(external_field), intent(in) :: field !< Handle for time interpolated 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 @@ -48,7 +49,7 @@ subroutine time_interp_external_0d(field_id, time, data_in, verbose, scale) data_in = data_in * I_scale endif ; endif - call time_interp_extern(field_id, time, data_in, verbose=verbose) + call time_interp_extern(field, time, data_in, verbose=verbose) if (present(scale)) then ; if (scale /= 1.0) then ! Rescale data that has been newly set and restore the scaling of unset data. @@ -63,10 +64,9 @@ 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, & +subroutine time_interp_external_2d(field, time, data_in, interp, & verbose, horz_interp, mask_out, turns, scale) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() + type(external_field), intent(in) :: field !< Handle for time interpolated 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 @@ -105,11 +105,11 @@ subroutine time_interp_external_2d(field_id, time, data_in, interp, & qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) if (qturns == 0) then - call time_interp_extern(field_id, time, data_in, interp=interp, & + call time_interp_extern(field, 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_extern(field_id, time, data_pre_rot, interp=interp, & + call time_interp_extern(field, 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) @@ -136,10 +136,9 @@ 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, & +subroutine time_interp_external_3d(field, time, data_in, interp, & verbose, horz_interp, mask_out, turns, scale) - integer, intent(in) :: field_id !< The integer index of the external field returned - !! from a previous call to init_external_field() + type(external_field), intent(in) :: field !< Handle for time interpolated 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 @@ -178,11 +177,11 @@ subroutine time_interp_external_3d(field_id, time, data_in, interp, & qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) if (qturns == 0) then - call time_interp_extern(field_id, time, data_in, interp=interp, & + call time_interp_extern(field, 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_extern(field_id, time, data_pre_rot, interp=interp, & + call time_interp_extern(field, 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) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 113b6c045b..8e0e58c1b6 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -61,6 +61,7 @@ module MOM_ice_shelf use MOM_spatial_means, only : global_area_integral use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init +use MOM_interpolate, only : external_field implicit none ; private @@ -196,10 +197,10 @@ module MOM_ice_shelf id_shelf_sfc_mass_flux = -1 !>@} - integer :: id_read_mass !< An integer handle used in time interpolation of - !! the ice shelf mass read from a file - integer :: id_read_area !< An integer handle used in time interpolation of - !! the ice shelf mass read from a file + type(external_field) :: mass_handle + !< Handle for reading the time interpolated ice shelf mass from a file + type(external_field) :: area_handle + !< Handle for reading the time interpolated ice shelf area from a file type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to control diagnostic output. type(user_ice_shelf_CS), pointer :: user_CS => NULL() !< A pointer to the control structure for @@ -1118,7 +1119,7 @@ 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_external(CS%mass_handle, Time0, last_mass_shelf) do j=js,je ; do i=is,ie ! 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 @@ -1937,7 +1938,7 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) filename = trim(slasher(inputdir))//trim(shelf_file) call log_param(param_file, mdl, "INPUTDIR/SHELF_FILE", filename) - CS%id_read_mass = init_external_field(filename, shelf_mass_var, & + CS%mass_handle = init_external_field(filename, shelf_mass_var, & MOM_domain=CS%Grid_in%Domain, verbose=CS%debug) if (read_shelf_area) then @@ -1945,7 +1946,7 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) "The variable in SHELF_FILE with the shelf area.", & default="shelf_area") - CS%id_read_area = init_external_field(filename, shelf_area_var, & + CS%area_handle = init_external_field(filename, shelf_area_var, & MOM_domain=CS%Grid_in%Domain) endif @@ -2040,7 +2041,7 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) allocate(tmp2d(is:ie,js:je), source=0.0) endif - call time_interp_external(CS%id_read_mass, Time, tmp2d) + call time_interp_external(CS%mass_handle, Time, tmp2d) call rotate_array(tmp2d, CS%turns, ISS%mass_shelf) deallocate(tmp2d) diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 8a1aab3328..53615b0063 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -17,6 +17,7 @@ module MOM_oda_driver_mod use MOM_io, only : SINGLE_FILE use MOM_interp_infra, only : init_extern_field, get_external_field_info use MOM_interp_infra, only : time_interp_extern +use MOM_interpolate, only : external_field use MOM_remapping, only : remappingSchemesDoc use MOM_time_manager, only : time_type, real_to_time, get_date use MOM_time_manager, only : operator(+), operator(>=), operator(/=) @@ -80,8 +81,8 @@ module MOM_oda_driver_mod !> A structure containing integer handles for bias adjustment of tracers type :: INC_CS integer :: fldno = 0 !< The number of tracers - integer :: T_id !< The integer handle for the temperature file - integer :: S_id !< The integer handle for the salinity file + type(external_field) :: T !< The handle for the temperature file + type(external_field) :: S !< The handle for the salinity file end type INC_CS !> Control structure that contains a transpose of the ocean state across ensemble members. @@ -391,11 +392,11 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) "tendency adjustments", default='temp_salt_adjustment.nc') inc_file = trim(inputdir) // trim(bias_correction_file) - CS%INC_CS%T_id = init_extern_field(inc_file, "temp_increment", & + CS%INC_CS%T = init_extern_field(inc_file, "temp_increment", & correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) - CS%INC_CS%S_id = init_extern_field(inc_file, "salt_increment", & + CS%INC_CS%S = init_extern_field(inc_file, "salt_increment", & correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) - call get_external_field_info(CS%INC_CS%T_id,size=fld_sz) + call get_external_field_info(CS%INC_CS%T, size=fld_sz) CS%INC_CS%fldno = 2 if (CS%nk /= fld_sz(3)) call MOM_error(FATAL,'Increment levels /= ODA levels') @@ -578,9 +579,9 @@ subroutine get_bias_correction_tracer(Time, US, CS) call cpu_clock_begin(id_clock_bias_adjustment) - call horiz_interp_and_extrap_tracer(CS%INC_CS%T_id, Time, CS%G, T_bias, & + call horiz_interp_and_extrap_tracer(CS%INC_CS%T, Time, CS%G, T_bias, & valid_flag, z_in, z_edges_in, missing_value, scale=US%degC_to_C*US%s_to_T, spongeOngrid=.true.) - call horiz_interp_and_extrap_tracer(CS%INC_CS%S_id, Time, CS%G, S_bias, & + call horiz_interp_and_extrap_tracer(CS%INC_CS%S, Time, CS%G, S_bias, & valid_flag, z_in, z_edges_in, missing_value, scale=US%ppt_to_S*US%s_to_T, spongeOngrid=.true.) ! This should be replaced to use mask_z instead of the following lines diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 2a5cef5974..6a439dfd22 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -21,6 +21,7 @@ module MOM_MEKE use MOM_interface_heights, only : find_eta use MOM_interpolate, only : init_external_field, time_interp_external use MOM_interpolate, only : time_interp_external_init +use MOM_interpolate, only : external_field use MOM_io, only : vardesc, var_desc, slasher use MOM_isopycnal_slopes, only : calc_isoneutral_slopes use MOM_restart, only : MOM_restart_CS, register_restart_field, query_initialized @@ -129,7 +130,7 @@ module MOM_MEKE integer :: id_Lrhines = -1, id_Leady = -1 integer :: id_MEKE_equilibrium = -1 !>@} - integer :: id_eke = -1 !< Handle for reading in EKE from a file + type(external_field) :: eke_handle !< Handle for reading in EKE from a file ! Infrastructure integer :: id_clock_pass !< Clock for group pass calls type(group_pass_type) :: pass_MEKE !< Group halo pass handle for MEKE%MEKE and maybe MEKE%Kh_diff @@ -627,7 +628,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif case(EKE_FILE) - call time_interp_external(CS%id_eke, Time, data_eke, scale=US%m_s_to_L_T**2) + call time_interp_external(CS%eke_handle, Time, data_eke, scale=US%m_s_to_L_T**2) do j=js,je ; do i=is,ie MEKE%MEKE(i,j) = data_eke(i,j) * G%mask2dT(i,j) enddo; enddo @@ -1153,7 +1154,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, dbcomms_CS, CS, MEKE, inputdir = slasher(inputdir) eke_filename = trim(inputdir) // trim(eke_filename) - CS%id_eke = init_external_field(eke_filename, eke_varname, domain=G%Domain%mpp_domain) + CS%eke_handle = init_external_field(eke_filename, eke_varname, domain=G%Domain%mpp_domain) case("prog") CS%eke_src = EKE_PROG ! Read all relevant parameters and write them to the model log. diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 584ccccc93..2a30f68b42 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -22,6 +22,7 @@ module MOM_ALE_sponge use MOM_grid, only : ocean_grid_type use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer use MOM_interpolate, only : init_external_field, get_external_field_info, time_interp_external_init +use MOM_interpolate, only : external_field use MOM_remapping, only : remapping_cs, remapping_core_h, initialize_remapping use MOM_spatial_means, only : global_i_mean use MOM_time_manager, only : time_type @@ -66,7 +67,7 @@ module MOM_ALE_sponge !> A structure for creating arrays of pointers to 3D arrays with extra gridding information type :: p3d - integer :: id !< id for FMS external time interpolator + !integer :: id !< id for FMS external time interpolator integer :: nz_data !< The number of vertical levels in the input field. integer :: num_tlevs !< The number of time records contained in the file real, dimension(:,:,:), pointer :: p => NULL() !< pointer to the data [various] @@ -75,7 +76,7 @@ module MOM_ALE_sponge !> A structure for creating arrays of pointers to 2D arrays with extra gridding information type :: p2d - integer :: id !< id for FMS external time interpolator + type(external_field) :: field !< Time interpolator field handle integer :: nz_data !< The number of vertical levels in the input field integer :: num_tlevs !< The number of time records contained in the file real :: scale = 1.0 !< A multiplicative factor by which to rescale input data [various] @@ -771,7 +772,6 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, !! if not given, use 'none' real, optional, intent(in) :: scale !< A factor by which to rescale the input data, including any !! contributions due to dimensional rescaling [various ~> 1]. - !! The default is 1. ! Local variables integer :: isd, ied, jsd, jed @@ -798,15 +798,15 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, ! get a unique time interp id for this field. If sponge data is on-grid, then setup ! to only read on the computational domain if (CS%spongeDataOngrid) then - CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname, MOM_domain=G%Domain) + CS%Ref_val(CS%fldno)%field = init_external_field(filename, fieldname, MOM_domain=G%Domain) else - CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname) + CS%Ref_val(CS%fldno)%field = init_external_field(filename, fieldname) endif CS%Ref_val(CS%fldno)%name = sp_name CS%Ref_val(CS%fldno)%long_name = long_name CS%Ref_val(CS%fldno)%unit = unit fld_sz(1:4) = -1 - call get_external_field_info(CS%Ref_val(CS%fldno)%id, size=fld_sz) + call get_external_field_info(CS%Ref_val(CS%fldno)%field, size=fld_sz) nz_data = fld_sz(3) CS%Ref_val(CS%fldno)%nz_data = nz_data !< individual sponge fields may reside on a different vertical grid CS%Ref_val(CS%fldno)%num_tlevs = fld_sz(4) @@ -899,23 +899,23 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename ! containing time-interpolated values from an external file corresponding ! to the current model date. if (CS%spongeDataOngrid) then - CS%Ref_val_u%id = init_external_field(filename_u, fieldname_u, domain=G%Domain%mpp_domain) + CS%Ref_val_u%field = init_external_field(filename_u, fieldname_u, domain=G%Domain%mpp_domain) else - CS%Ref_val_u%id = init_external_field(filename_u, fieldname_u) + CS%Ref_val_u%field = init_external_field(filename_u, fieldname_u) endif fld_sz(1:4) = -1 - call get_external_field_info(CS%Ref_val_u%id, size=fld_sz) + call get_external_field_info(CS%Ref_val_u%field, size=fld_sz) CS%Ref_val_u%nz_data = fld_sz(3) CS%Ref_val_u%num_tlevs = fld_sz(4) CS%Ref_val_u%scale = US%m_s_to_L_T ; if (present(scale)) CS%Ref_val_u%scale = scale if (CS%spongeDataOngrid) then - CS%Ref_val_v%id = init_external_field(filename_v, fieldname_v, domain=G%Domain%mpp_domain) + CS%Ref_val_v%field = init_external_field(filename_v, fieldname_v, domain=G%Domain%mpp_domain) else - CS%Ref_val_v%id = init_external_field(filename_v, fieldname_v) + CS%Ref_val_v%field = init_external_field(filename_v, fieldname_v) endif fld_sz(1:4) = -1 - call get_external_field_info(CS%Ref_val_v%id, size=fld_sz) + call get_external_field_info(CS%Ref_val_v%field, size=fld_sz) CS%Ref_val_v%nz_data = fld_sz(3) CS%Ref_val_v%num_tlevs = fld_sz(4) CS%Ref_val_v%scale = US%m_s_to_L_T ; if (present(scale)) CS%Ref_val_v%scale = scale @@ -989,7 +989,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%time_varying_sponges) then do m=1,CS%fldno nz_data = CS%Ref_val(m)%nz_data - call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, G, sp_val, & + call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%field, Time, G, sp_val, & mask_z, z_in, z_edges_in, missing_value, & scale=CS%Ref_val(m)%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & answer_date=CS%hor_regrid_answer_date) @@ -1073,7 +1073,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) if (CS%time_varying_sponges) then nz_data = CS%Ref_val_u%nz_data ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, G, sp_val, & + call horiz_interp_and_extrap_tracer(CS%Ref_val_u%field, Time, G, sp_val, & mask_z, z_in, z_edges_in, missing_value, & scale=CS%Ref_val_u%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & answer_date=CS%hor_regrid_answer_date) @@ -1121,7 +1121,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) deallocate(sp_val, mask_u, mask_z, hsrc) nz_data = CS%Ref_val_v%nz_data ! Interpolate from the external horizontal grid and in time - call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, G, sp_val, & + call horiz_interp_and_extrap_tracer(CS%Ref_val_v%field, Time, G, sp_val, & mask_z, z_in, z_edges_in, missing_value, & scale=CS%Ref_val_v%scale, spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,& answer_date=CS%hor_regrid_answer_date) @@ -1341,7 +1341,7 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, GV, turns, param_file) ! We don't want to repeat FMS init in set_up_ALE_sponge_field_varying() ! (time_interp_external_init, init_external_field, etc), so we manually ! do a portion of this function below. - sponge%Ref_val(n)%id = sponge_in%Ref_val(n)%id + sponge%Ref_val(n)%field = sponge_in%Ref_val(n)%field sponge%Ref_val(n)%num_tlevs = sponge_in%Ref_val(n)%num_tlevs nz_data = sponge_in%Ref_val(n)%nz_data diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 5f7acd982b..3096fe72cd 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -16,6 +16,7 @@ module MOM_diabatic_aux 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_external, time_interp_external_init +use MOM_interpolate, only : external_field 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 @@ -64,7 +65,7 @@ module MOM_diabatic_aux !! is added with a temperature of the local SST. logical :: var_pen_sw !< If true, use one of the CHL_A schemes to determine the !! e-folding depth of incoming shortwave radiation. - integer :: sbc_chl !< An integer handle used in time interpolation of + type(external_field) :: sbc_chl !< A handle used in time interpolation of !! chlorophyll read from a file. logical :: chl_from_file !< If true, chl_a is read from a file. diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 2a5e3f8854..ef8e712b7a 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -1,4 +1,4 @@ -!> Simulates CFCs using atmospheric pressure, wind speed and sea ice cover + !> Simulates CFCs using atmospheric pressure, wind speed and sea ice cover !! provided via cap (only NUOPC cap is implemented so far). module MOM_CFC_cap @@ -19,7 +19,8 @@ module MOM_CFC_cap use MOM_restart, only : query_initialized, set_initialized, MOM_restart_CS use MOM_spatial_means, only : global_mass_int_EFP use MOM_time_manager, only : time_type -use time_interp_external_mod, only : init_external_field, time_interp_external +use MOM_interpolate, only : time_interp_external +use MOM_interpolate, only : external_field use MOM_tracer_registry, only : register_tracer use MOM_tracer_types, only : tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut @@ -428,7 +429,8 @@ end subroutine CFC_cap_surface_state !> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM !! concentration, and calculating the solubility, Schmidt number, and gas exchange. -subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id_cfc12_atm) +subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, & + cfc11_atm_handle, cfc12_atm_handle) type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type type(surface), intent(in ) :: sfc_state !< A structure containing fields @@ -439,8 +441,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id real, intent(in ) :: Rho0 !< The mean ocean density [R ~> kg m-3] type(time_type), intent(in ) :: Time !< The time of the fluxes, used for interpolating the !! CFC's concentration in the atmosphere. - integer, optional, intent(inout):: id_cfc11_atm !< id number for time_interp_external. - integer, optional, intent(inout):: id_cfc12_atm !< id number for time_interp_external. + type(external_field), optional, intent(inout) :: cfc11_atm_handle !< Handle for time-interpolated CFC11 + type(external_field), optional, intent(inout) :: cfc12_atm_handle !< Handle for time-interpolated CFC12 ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & @@ -463,8 +465,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ! CFC11 ATM concentration - if (present(id_cfc11_atm) .and. (id_cfc11_atm /= -1)) then - call time_interp_external(id_cfc11_atm, Time, cfc11_atm) + if (present(cfc11_atm_handle)) then + call time_interp_external(cfc11_atm_handle, Time, cfc11_atm) ! convert from ppt (pico mol/mol) to mol/mol cfc11_atm = cfc11_atm * 1.0e-12 else @@ -474,8 +476,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id endif ! CFC12 ATM concentration - if (present(id_cfc12_atm) .and. (id_cfc12_atm /= -1)) then - call time_interp_external(id_cfc12_atm, Time, cfc12_atm) + if (present(cfc12_atm_handle)) then + call time_interp_external(cfc12_atm_handle, Time, cfc12_atm) ! convert from ppt (pico mol/mol) to mol/mol cfc12_atm = cfc12_atm * 1.0e-12 else From 1268c9758e008e26ce2650a19d895b99c8bdf654 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 6 Jun 2023 20:35:24 -0400 Subject: [PATCH 04/14] FMS2: Remove MPP-based axis data access With removal of axis-based operations in FMS2 I/O, this patch removes references to these calls and replaces them with MOM `axes_info` types. References to FMS1 read into an `axistype`, but the contents are transferred to an `axis_info`. FMS2 directly populates the `axis_info` content. The `get_external_field_info` calls are modified to return `axis_info` rather than `axistype`. The redundant `get_axis_data` function is also removed from `MOM_interp_infra`, since `get_axis_info` provides an equivalent operation. Generally speaking, this is not an improvement of the codebase. The FMS1 layer does a redundant copy of data from `axistype` to `axis_info`. The FMS2 layer is significantly worse, and re-opens the file to read the axis data for each field! But if the intention is to leverage the existing API, then I don't think we have any choice at the moment. Assuming this is a relatively infrequent operation, this should not cause any measureable issues, but it needs to be watched carefully. --- config_src/infra/FMS1/MOM_interp_infra.F90 | 44 +++++++++++++++------ config_src/infra/FMS2/MOM_interp_infra.F90 | 32 ++++++--------- src/framework/MOM_horizontal_regridding.F90 | 10 ++--- 3 files changed, 49 insertions(+), 37 deletions(-) diff --git a/config_src/infra/FMS1/MOM_interp_infra.F90 b/config_src/infra/FMS1/MOM_interp_infra.F90 index e14233a64b..70bc99827e 100644 --- a/config_src/infra/FMS1/MOM_interp_infra.F90 +++ b/config_src/infra/FMS1/MOM_interp_infra.F90 @@ -4,9 +4,11 @@ module MOM_interp_infra ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domain_infra, only : MOM_domain_type, domain2d +use MOM_io, only : axis_info +use MOM_io, only : set_axis_info use MOM_time_manager, only : time_type use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type -use mpp_io_mod, only : axistype, mpp_get_axis_data +use mpp_io_mod, only : axistype, mpp_get_axis_data, mpp_get_atts use time_interp_external_mod, only : 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 @@ -157,13 +159,33 @@ end function get_extern_field_size !> get axes of an external field from field index -function get_extern_field_axes(index) +function get_extern_field_axes(index) result(axes) - integer, intent(in) :: index !< field index - type(axistype), dimension(4) :: get_extern_field_axes !< field axes + integer, intent(in) :: index !< FMS interpolation field index + type(axis_info) :: axes(4) !< MOM IO field axes handle - get_extern_field_axes = get_external_field_axes(index) + type(axistype), dimension(4) :: fms_axes(4) + ! FMS axis handles + character(len=32) :: name + ! Axis name + real, allocatable :: points(:) + ! Axis line points + integer :: length + ! Axis line point length + integer :: i + ! Loop index + fms_axes = get_external_field_axes(index) + + do i = 1, 4 + call mpp_get_atts(fms_axes(i), name=name, len=length) + + allocate(points(length)) + call mpp_get_axis_data(fms_axes(i), points) + call set_axis_info(axes(i), name=name, ax_data=points) + + deallocate(points) + enddo end function get_extern_field_axes @@ -180,12 +202,12 @@ end function get_extern_field_missing !> Get information about the external fields. subroutine get_external_field_info(field, size, axes, missing) - type(external_field), intent(in) :: field !< Handle for time interpolated 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 + type(external_field), intent(in) :: field !< Handle for time interpolated external + !! field returned from a previous + !! call to init_external_field() + integer, optional, intent(inout) :: size(4) !< Dimension sizes for the input data + type(axis_info), optional, intent(inout) :: axes(4) !< 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_extern_field_size(field%id) diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index 7964b3537f..db471d0fc1 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -4,9 +4,10 @@ module MOM_interp_infra ! This file is part of MOM6. See LICENSE.md for the license. use MOM_domain_infra, only : MOM_domain_type, domain2d +use MOM_io, only : axis_info +use MOM_io, only : get_var_axes_info use MOM_time_manager, only : time_type use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type -use mpp_io_mod, only : axistype, mpp_get_axis_data use time_interp_external_mod, only : 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 @@ -16,7 +17,7 @@ module MOM_interp_infra public :: horiz_interp_type, horizontal_interp_init public :: time_interp_extern, init_extern_field, time_interp_extern_init -public :: get_external_field_info, axistype, get_axis_data +public :: get_external_field_info public :: run_horiz_interp, build_horiz_interp_weights public :: external_field @@ -135,15 +136,6 @@ subroutine build_horiz_interp_weights_2d_to_2d(Interp, lon_in, lat_in, lon_out, end subroutine build_horiz_interp_weights_2d_to_2d -!> Extracts and returns the axis data stored in an axistype. -subroutine get_axis_data( axis, dat ) - type(axistype), intent(in) :: axis !< An axis type - real, dimension(:), intent(out) :: dat !< The data in the axis variable - - call mpp_get_axis_data( axis, dat ) -end subroutine get_axis_data - - !> get size of an external field from field index function get_extern_field_size(index) @@ -156,13 +148,11 @@ end function get_extern_field_size !> get axes of an external field from field index -function get_extern_field_axes(index) - - integer, intent(in) :: index !< field index - type(axistype), dimension(4) :: get_extern_field_axes !< field axes - - get_extern_field_axes = get_external_field_axes(index) +function get_extern_field_axes(field) result(axes) + type(external_field), intent(in) :: field !< Field handle + type(axis_info), dimension(4) :: axes !< Field axes + call get_var_axes_info(field%filename, field%label, axes) end function get_extern_field_axes @@ -182,16 +172,16 @@ subroutine get_external_field_info(field, size, axes, missing) type(external_field), intent(in) :: field !< Handle for time interpolated 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 + integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data + type(axis_info), 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_extern_field_size(field%id) endif if (present(axes)) then - axes(1:4) = get_extern_field_axes(field%id) + axes(1:4) = get_extern_field_axes(field) endif if (present(missing)) then diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index bedf710582..883653d715 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -16,7 +16,7 @@ module MOM_horizontal_regridding use MOM_interpolate, only : time_interp_external use MOM_interp_infra, only : run_horiz_interp, build_horiz_interp_weights use MOM_interp_infra, only : horiz_interp_type, horizontal_interp_init -use MOM_interp_infra, only : axistype, get_external_field_info, get_axis_data +use MOM_interp_infra, only : get_external_field_info use MOM_interp_infra, only : external_field use MOM_time_manager, only : time_type use MOM_io, only : axis_info, get_axis_info, get_var_axes_info, MOM_read_data @@ -668,7 +668,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim] logical :: add_np type(horiz_interp_type) :: Interp - type(axistype), dimension(4) :: axes_data + type(axis_info), dimension(4) :: axes_data integer :: is, ie, js, je ! compute domain indices integer :: isg, ieg, jsg, jeg ! global extent integer :: isd, ied, jsd, jed ! data domain indices @@ -728,8 +728,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & if (PRESENT(spongeOngrid)) is_ongrid = spongeOngrid if (.not. is_ongrid) then 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) + call get_axis_info(axes_data(1), ax_data=lon_in) + call get_axis_info(axes_data(2), ax_data=lat_in) endif allocate(z_in(kd), z_edges_in(kd+1)) @@ -737,7 +737,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, & allocate(tr_z(isd:ied,jsd:jed,kd), source=0.0) allocate(mask_z(isd:ied,jsd:jed,kd), source=0.0) - call get_axis_data(axes_data(3), z_in) + call get_axis_info(axes_data(3), ax_data=z_in) if (present(m_to_Z)) then ; do k=1,kd ; z_in(k) = m_to_Z * z_in(k) ; enddo ; endif From 35e3642dc222572672c78310b35ff9456816d76e Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 7 Jun 2023 00:24:30 -0400 Subject: [PATCH 05/14] FMS2: Update time_interp_external functions This patch shifts all remaining time_interp_external functions from time_interp_external to equivalent ones in time_interp_external2. Internally, time-interpolated fields are initialized with `ongrid` set to `.true.`, and such fields are assumed to be on-grid. This seems to hold for all existing instances of `time_interp_external`, but needs to be monitored in the future somehow. --- config_src/infra/FMS2/MOM_interp_infra.F90 | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index db471d0fc1..09a9eb0421 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -8,10 +8,10 @@ module MOM_interp_infra use MOM_io, only : get_var_axes_info use MOM_time_manager, only : time_type 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 -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_interp_external2_mod, only : time_interp_external +use time_interp_external2_mod, only : init_external_field, time_interp_external_init +use time_interp_external2_mod, only : get_external_field_size +use time_interp_external2_mod, only : get_external_field_missing implicit none ; private @@ -267,12 +267,14 @@ function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & if (present(MOM_Domain)) then field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & - correct_leap_year_inconsistency=correct_leap_year_inconsistency) + verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency, & + ongrid=.true.) else field%id = init_external_field(file, fieldname, domain=domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & - correct_leap_year_inconsistency=correct_leap_year_inconsistency) + verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency, & + ongrid=.true.) endif end function init_extern_field From d5ce3364437c151326bc74a4fea632455ebc5162 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 16 Jun 2023 10:40:57 -0400 Subject: [PATCH 06/14] FMS2: Case-insensitive init_external_field The FMS1 implementation of init_external_field is case-insensitive, but the FMS2 implementation is case-sensitive, which can cause errors in older established input files. This patch sweeps through the fields of the input files and checks for a case-insensitive match (using lowercase()). This requires an additional open/close of the file. --- config_src/infra/FMS2/MOM_interp_infra.F90 | 56 ++++++++++++++++++++-- 1 file changed, 52 insertions(+), 4 deletions(-) diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index 09a9eb0421..4420108b8d 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -7,7 +7,11 @@ module MOM_interp_infra use MOM_io, only : axis_info use MOM_io, only : get_var_axes_info use MOM_time_manager, only : time_type -use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use MOM_error_handler, only : MOM_error, FATAL +use MOM_string_functions, only : lowercase +use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use netcdf_io_mod, only : FmsNetcdfFile_t, netcdf_file_open, netcdf_file_close +use netcdf_io_mod, only : get_num_variables, get_variable_names use time_interp_external2_mod, only : time_interp_external use time_interp_external2_mod, only : init_external_field, time_interp_external_init use time_interp_external2_mod, only : get_external_field_size @@ -262,16 +266,60 @@ function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & !! a model date of Feb 29. onto a common year on Feb. 28. type(external_field) :: field !< Handle to external field + type(FmsNetcdfFile_t) :: extern_file + ! Local instance of netCDF file used to locate case-insensitive field name + integer :: num_fields + ! Number of fields in external file + character(len=256), allocatable :: extern_fieldnames(:) + ! List of field names in file + ! NOTE: length should NF90_MAX_NAME, but I don't know how to read it + character(len=:), allocatable :: label + ! Case-insensitive match to fieldname in file + logical :: rc + ! Return status + integer :: i + ! Loop index + field%filename = file - field%label = fieldname + + ! FMS2's init_external_field is case sensitive, so we must replicate the + ! case-insensitivity of FMS1. This requires opening the file twice. + + rc = netcdf_file_open(extern_file, file, 'read') + if (.not. rc) then + call MOM_error(FATAL, 'init_extern_file: file ' // trim(file) & + // ' could not be opened.') + endif + + ! TODO: broadcast = .false.? + num_fields = get_num_variables(extern_file) + + allocate(extern_fieldnames(num_fields)) + call get_variable_names(extern_file, extern_fieldnames) + + do i = 1, num_fields + if (lowercase(extern_fieldnames(i)) == lowercase(fieldname)) then + field%label = extern_fieldnames(i) + exit + endif + enddo + + call netcdf_file_close(extern_file) + + if (.not. allocated(field%label)) then + call MOM_error(FATAL, 'init_extern_field: field ' // trim(fieldname) & + // ' not found in ' // trim(file) // '.') + endif + + ! Pass to FMS2 implementation of init_external_field if (present(MOM_Domain)) then - field%id = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & + field%id = init_external_field(file, field%label, domain=MOM_domain%mpp_domain, & verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency, & ongrid=.true.) else - field%id = init_external_field(file, fieldname, domain=domain, & + field%id = init_external_field(file, field%label, domain=domain, & verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & correct_leap_year_inconsistency=correct_leap_year_inconsistency, & ongrid=.true.) From cd5faea0e71f85608831994efc4ec9876dbf5a23 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Mon, 19 Jun 2023 19:11:48 -0800 Subject: [PATCH 07/14] The brine plume parameterization, - including now passing the dimensional scaling tests. --- .../FMS_cap/MOM_surface_forcing_gfdl.F90 | 2 +- config_src/infra/FMS2/MOM_interp_infra.F90 | 4 ++ src/core/MOM_forcing_type.F90 | 8 +-- .../vertical/MOM_diabatic_aux.F90 | 72 +++++++++++++++++-- .../vertical/MOM_diabatic_driver.F90 | 4 +- 5 files changed, 77 insertions(+), 13 deletions(-) diff --git a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 index 022054c54e..6fb237eddf 100644 --- a/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 +++ b/config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90 @@ -581,7 +581,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, endif if (associated(IOB%excess_salt)) then do j=js,je ; do i=is,ie - fluxes%salt_left_behind(i,j) = G%mask2dT(i,j)*( - kg_m2_s_conversion*IOB%excess_salt(i-i0,j-j0)) + fluxes%salt_left_behind(i,j) = G%mask2dT(i,j)*(kg_m2_s_conversion*IOB%excess_salt(i-i0,j-j0)) enddo ; enddo endif diff --git a/config_src/infra/FMS2/MOM_interp_infra.F90 b/config_src/infra/FMS2/MOM_interp_infra.F90 index 4420108b8d..0b45b752ae 100644 --- a/config_src/infra/FMS2/MOM_interp_infra.F90 +++ b/config_src/infra/FMS2/MOM_interp_infra.F90 @@ -313,6 +313,10 @@ function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & ! Pass to FMS2 implementation of init_external_field + ! NOTE: external fields are currently assumed to be on-grid, which holds + ! across the current codebase. In the future, we may need to either enforce + ! this or somehow relax this requirement. + if (present(MOM_Domain)) then field%id = init_external_field(file, field%label, domain=MOM_domain%mpp_domain, & verbose=verbose, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index e6d4efd424..c78c8532e7 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -741,15 +741,15 @@ subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & endif ! Salt fluxes - Net_salt(i) = 0.0 - if (do_NSR) Net_salt_rate(i) = 0.0 + net_salt(i) = 0.0 + if (do_NSR) net_salt_rate(i) = 0.0 ! Convert salt_flux from kg (salt)/(m^2 * s) to ! Boussinesq: (ppt * m) ! non-Bouss: (g/m^2) if (associated(fluxes%salt_flux)) then - Net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H + net_salt(i) = (scale * dt * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H !Repeat above code for 'rate' term - if (do_NSR) Net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H + if (do_NSR) net_salt_rate(i) = (scale * 1. * (1000.0*US%ppt_to_S * fluxes%salt_flux(i,j))) * GV%RZ_to_H endif ! Diagnostics follow... diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 3096fe72cd..e82ca247be 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -67,7 +67,10 @@ module MOM_diabatic_aux !! e-folding depth of incoming shortwave radiation. type(external_field) :: sbc_chl !< A handle used in time interpolation of !! chlorophyll read from a file. - logical :: chl_from_file !< If true, chl_a is read from a file. + logical :: chl_from_file !< If true, chl_a is read from a file. + logical :: do_brine_plume !< If true, insert salt flux below the surface according to + !! a parameterization by \cite Nguyen2009. + integer :: brine_plume_n !< The exponent in the brine plume parameterization. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output @@ -1032,7 +1035,7 @@ end subroutine diagnoseMLDbyEnergy subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, & aggregate_FW_forcing, evap_CFL_limit, & minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, & - SkinBuoyFlux ) + SkinBuoyFlux, MLD) type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -1062,6 +1065,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !! salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3]. + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: MLD !< Mixed layer depth for brine plumes [Z ~> m] ! Local variables integer, parameter :: maxGroundings = 5 @@ -1100,7 +1105,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t netheat_rate, & ! netheat but for dt=1 [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate) ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - netMassInOut_rate! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1] + netMassInOut_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1] + mixing_depth ! Mixed layer depth [Z -> m] real, dimension(SZI_(G), SZK_(GV)) :: & h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2] T2d, & ! A 2-d copy of the layer temperatures [C ~> degC] @@ -1130,6 +1136,11 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! and rejected brine are initially applied in vanishingly thin layers at the ! top of the layer before being mixed throughout the layer. logical :: calculate_buoyancy ! If true, calculate the surface buoyancy flux. + real, dimension(SZI_(G)) :: dK ! Depth [Z ~> m]. + real, dimension(SZI_(G)) :: A_brine ! Constant [Z-(n+1) ~> m-(n+1)]. + real :: fraction_left_brine ! Sum for keeping track of the fraction of brine so far (in depth) + real :: plume_fraction ! Sum for keeping track of the fraction of brine so far (in depth) + real :: plume_flux ! Brine flux to move downwards [S H ~> ppt m or ppt kg m-2] integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, k, nz, nb character(len=45) :: mesg @@ -1156,6 +1167,17 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 endif + if (CS%do_brine_plume .and. .not. present(MLD)) then + call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& + "Brine plume parameterization requires a mixed-layer depth,\n"//& + "currently coming from the energetic PBL scheme.") + endif + if (CS%do_brine_plume .and. .not. associated(fluxes%salt_left_behind)) then + call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& + "Brine plume parameterization requires DO_BRINE_PLUME\n"//& + "to be turned on in SIS2 as well as MOM6.") + endif + ! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total ! depth of the ocean is vanishing. It does not (yet) handle a value of zero. ! To accommodate vanishing upper layers, we need to allow for an instantaneous @@ -1173,7 +1195,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, & !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho, & !$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, & - !$OMP EnthalpyConst) & + !$OMP EnthalpyConst,MLD) & !$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, & !$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, & !$OMP IforcingDepthScale, & @@ -1298,6 +1320,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! ocean (and corresponding outward heat content), and ignoring penetrative SW. ! B/ update mass, salt, temp from mass leaving ocean. ! C/ update temp due to penetrative SW + if (CS%do_brine_plume) then + do i=is,ie + mixing_depth(i) = max(MLD(i,j) - minimum_forcing_depth * GV%H_to_Z, minimum_forcing_depth * GV%H_to_Z) + mixing_depth(i) = min(mixing_depth(i), sum(h(i,j,:)) * GV%H_to_Z) + A_brine(i) = (CS%brine_plume_n + 1) / (mixing_depth(i) ** (CS%brine_plume_n + 1)) + enddo + endif + do i=is,ie if (G%mask2dT(i,j) > 0.) then @@ -1370,8 +1400,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t enddo ! k=1,1 ! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt. + fraction_left_brine = 1.0 do k=1,nz - ! Place forcing into this layer if this layer has nontrivial thickness. ! For layers thin relative to 1/IforcingDepthScale, then distribute ! forcing into deeper layers. @@ -1386,11 +1416,34 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t fractionOfForcing = -evap_CFL_limit*h2d(i,k)/netMassOut(i) endif + if (CS%do_brine_plume .and. associated(fluxes%salt_left_behind)) then + if (fluxes%salt_left_behind(i,j) > 0 .and. fraction_left_brine > 0.0) then + ! Place forcing into this layer by depth for brine plume parameterization. + if (k == 1) then + dK(i) = 0.5 * h(i,j,k) * GV%H_to_Z ! Depth of center of layer K + plume_flux = - (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H + plume_fraction = 1.0 + else + dK(i) = dK(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * GV%H_to_Z ! Depth of center of layer K + plume_flux = 0.0 + endif + if (dK(i) <= mixing_depth(i) .and. fraction_left_brine > 0.0) then + plume_fraction = min(fraction_left_brine, A_brine(i) * dK(i) ** CS%brine_plume_n * h(i,j,k) * GV%H_to_Z) + else + IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) ) + ! plume_fraction = fraction_left_brine, unless h2d is less than IforcingDepthScale. + plume_fraction = min(fraction_left_brine, h2d(i,k)*IforcingDepthScale) + endif + fraction_left_brine = fraction_left_brine - plume_fraction + plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H + endif + endif + ! Change in state due to forcing dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) ) dTemp = fractionOfForcing*netHeat(i) - dSalt = max( fractionOfForcing*netSalt(i), -CS%dSalt_frac_max * h2d(i,k) * tv%S(i,j,k)) + dSalt = max( fractionOfForcing*netSalt(i) + plume_flux, -CS%dSalt_frac_max * h2d(i,k) * tv%S(i,j,k)) ! Update the forcing by the part to be consumed within the present k-layer. ! If fractionOfForcing = 1, then new netMassOut vanishes. @@ -1698,6 +1751,13 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori CS%use_calving_heat_content = .false. endif + call get_param(param_file, mdl, "DO_BRINE_PLUME", CS%do_brine_plume, & + "If true, use a brine plume parameterization from "//& + "Nguyen et al., 2009.", default=.false.) + call get_param(param_file, mdl, "BRINE_PLUME_EXPONENT", CS%brine_plume_n, & + "If using the brine plume parameterization, set the integer exponent.", & + default=5, do_not_log=.not.CS%do_brine_plume) + if (useALEalgorithm) then CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, & Time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index b0d04e434c..2136c65a61 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -835,7 +835,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=Hml) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_m) @@ -1380,7 +1380,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=Hml) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_m) From 874684ef7b84169461a07e87c4e6610a07c09536 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 28 Jun 2023 12:04:33 -0800 Subject: [PATCH 08/14] Fix problem when running Tidal_bay case with gnu. --- .../vertical/MOM_set_viscosity.F90 | 34 ++++++++++--------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 47d4dffef6..481aa5e9fc 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1027,22 +1027,24 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv) endif endif - if (CS%body_force_drag .and. (h_bbl_drag(i) > 0.0)) then - ! Increment the Rayleigh drag as a way introduce the bottom drag as a body force. - h_sum = 0.0 - I_hwtot = 1.0 / h_bbl_drag(i) - do k=nz,1,-1 - h_bbl_fr = min(h_bbl_drag(i) - h_sum, h_at_vel(i,k)) * I_hwtot - if (m==1) then - visc%Ray_u(I,j,k) = visc%Ray_u(I,j,k) + (CS%cdrag*US%L_to_Z*umag_avg(I)) * h_bbl_fr - else - visc%Ray_v(i,J,k) = visc%Ray_v(i,J,k) + (CS%cdrag*US%L_to_Z*umag_avg(i)) * h_bbl_fr - endif - h_sum = h_sum + h_at_vel(i,k) - if (h_sum >= h_bbl_drag(i)) exit ! The top of this layer is above the drag zone. - enddo - ! Do not enhance the near-bottom viscosity in this case. - Kv_bbl = CS%Kv_BBL_min + if (CS%body_force_drag) then + if (h_bbl_drag(i) > 0.0) then + ! Increment the Rayleigh drag as a way introduce the bottom drag as a body force. + h_sum = 0.0 + I_hwtot = 1.0 / h_bbl_drag(i) + do k=nz,1,-1 + h_bbl_fr = min(h_bbl_drag(i) - h_sum, h_at_vel(i,k)) * I_hwtot + if (m==1) then + visc%Ray_u(I,j,k) = visc%Ray_u(I,j,k) + (CS%cdrag*US%L_to_Z*umag_avg(I)) * h_bbl_fr + else + visc%Ray_v(i,J,k) = visc%Ray_v(i,J,k) + (CS%cdrag*US%L_to_Z*umag_avg(i)) * h_bbl_fr + endif + h_sum = h_sum + h_at_vel(i,k) + if (h_sum >= h_bbl_drag(i)) exit ! The top of this layer is above the drag zone. + enddo + ! Do not enhance the near-bottom viscosity in this case. + Kv_bbl = CS%Kv_BBL_min + endif endif kv_bbl = max(CS%Kv_BBL_min, kv_bbl) From 236f71b92ea94f46cf72f46d6d74d5324fca55ed Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 30 Jun 2023 11:17:58 -0800 Subject: [PATCH 09/14] Avoiding visc_rem issues inside land mask. Tweaking the brine plume code. --- src/core/MOM_continuity_PPM.F90 | 8 ++++---- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 8 +++++--- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 090d1ee0fb..73c6503242 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -378,9 +378,9 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are dx_E = ratio_max(G%areaT(i+1,j), G%dy_Cu(I,j), 1000.0*G%dxT(i+1,j)) else ; dx_W = G%dxT(i,j) ; dx_E = G%dxT(i+1,j) ; endif - if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)) & + if (du_max_CFL(I) * visc_rem(I,k) > dx_W*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) & du_max_CFL(I) = (dx_W*CFL_dt - u(I,j,k)) / visc_rem(I,k) - if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)) & + if (du_min_CFL(I) * visc_rem(I,k) < -dx_E*CFL_dt - u(I,j,k)*G%mask2dCu(I,j)) & du_min_CFL(I) = -(dx_E*CFL_dt + u(I,j,k)) / visc_rem(I,k) enddo ; enddo endif @@ -1201,9 +1201,9 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac dy_S = ratio_max(G%areaT(i,j), G%dx_Cv(I,j), 1000.0*G%dyT(i,j)) dy_N = ratio_max(G%areaT(i,j+1), G%dx_Cv(I,j), 1000.0*G%dyT(i,j+1)) else ; dy_S = G%dyT(i,j) ; dy_N = G%dyT(i,j+1) ; endif - if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)) & + if (dv_max_CFL(i) * visc_rem(i,k) > dy_S*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) & dv_max_CFL(i) = (dy_S*CFL_dt - v(i,J,k)) / visc_rem(i,k) - if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)) & + if (dv_min_CFL(i) * visc_rem(i,k) < -dy_N*CFL_dt - v(i,J,k)*G%mask2dCv(i,J)) & dv_min_CFL(i) = -(dy_N*CFL_dt + v(i,J,k)) / visc_rem(i,k) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index e82ca247be..1f0a42dd03 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1323,7 +1323,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t if (CS%do_brine_plume) then do i=is,ie mixing_depth(i) = max(MLD(i,j) - minimum_forcing_depth * GV%H_to_Z, minimum_forcing_depth * GV%H_to_Z) - mixing_depth(i) = min(mixing_depth(i), sum(h(i,j,:)) * GV%H_to_Z) + mixing_depth(i) = min(mixing_depth(i), max(sum(h(i,j,:)), GV%angstrom_h) * GV%H_to_Z) A_brine(i) = (CS%brine_plume_n + 1) / (mixing_depth(i) ** (CS%brine_plume_n + 1)) enddo endif @@ -1437,13 +1437,15 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t fraction_left_brine = fraction_left_brine - plume_fraction plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H endif + else + plume_flux = 0.0 endif ! Change in state due to forcing dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) ) dTemp = fractionOfForcing*netHeat(i) - dSalt = max( fractionOfForcing*netSalt(i) + plume_flux, -CS%dSalt_frac_max * h2d(i,k) * tv%S(i,j,k)) + dSalt = max( fractionOfForcing*netSalt(i), -CS%dSalt_frac_max * h2d(i,k) * tv%S(i,j,k)) ! Update the forcing by the part to be consumed within the present k-layer. ! If fractionOfForcing = 1, then new netMassOut vanishes. @@ -1483,7 +1485,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t endif Ithickness = 1.0/h2d(i,k) ! Inverse of new thickness T2d(i,k) = (hOld*T2d(i,k) + dTemp)*Ithickness - tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt)*Ithickness + tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt + plume_flux)*Ithickness elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (h<0)') write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j) From 82e6f18a5fc8382bde0b24521c5f3dcc9eeb83eb Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 6 Jul 2023 17:41:01 -0800 Subject: [PATCH 10/14] Using the proper MLD in the brine plumes - it now works better on restart --- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 9 ++++----- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 4 ++-- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 1f0a42dd03..b572f8248b 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1065,8 +1065,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !! salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3]. - real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: MLD !< Mixed layer depth for brine plumes [Z ~> m] + real, pointer, dimension(:,:), optional :: MLD!< Mixed layer depth for brine plumes [Z ~> m] ! Local variables integer, parameter :: maxGroundings = 5 @@ -1167,7 +1166,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 endif - if (CS%do_brine_plume .and. .not. present(MLD)) then + if (CS%do_brine_plume .and. .not. associated(MLD)) then call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//& "Brine plume parameterization requires a mixed-layer depth,\n"//& "currently coming from the energetic PBL scheme.") @@ -1436,9 +1435,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t endif fraction_left_brine = fraction_left_brine - plume_fraction plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H + else + plume_flux = 0.0 endif - else - plume_flux = 0.0 endif ! Change in state due to forcing diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 42c7379a37..3bfa93481a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -815,7 +815,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=Hml) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=visc%MLD) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_m) @@ -1360,7 +1360,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, skinbuoyflux(:,:) = 0.0 call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=Hml) + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD=visc%MLD) if (CS%debug) then call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, scale=GV%H_to_m) From b44d65807827952304a245511986a73ee595d7da Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 7 Jul 2023 16:35:11 -0800 Subject: [PATCH 11/14] Always including MLD in call to applyBoundary... - I could move it up and make it not optional. --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 3bfa93481a..f1ac39fb0d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -875,7 +875,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim else call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth) + CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD=visc%MLD) endif ! endif for CS%use_energetic_PBL @@ -1414,7 +1414,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, else call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, & optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth) + CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD=visc%MLD) endif ! endif for CS%use_energetic_PBL From 85f3e0d0b4b03e61c0385a3968326b54c33e70fd Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 11 Jul 2023 22:00:05 -0800 Subject: [PATCH 12/14] Adding some OpenMP directives to brine plumes --- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index b572f8248b..8f60e7a979 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1192,9 +1192,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP H_limit_fluxes,numberOfGroundings,iGround,jGround,& !$OMP nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing, & !$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, & - !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho, & + !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,& !$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, & - !$OMP EnthalpyConst,MLD) & + !$OMP EnthalpyConst,MLD) & !$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, & !$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, & !$OMP IforcingDepthScale, & @@ -1202,7 +1202,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t !$OMP netMassIn,pres,d_pres,p_lay,dSV_dT_2d, & !$OMP netmassinout_rate,netheat_rate,netsalt_rate, & !$OMP drhodt,drhods,pen_sw_bnd_rate, & - !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst) & + !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, & + !$OMP mixing_depth,A_brine,fraction_left_brine, & + !$OMP plume_flux,plume_fraction,dK) & !$OMP firstprivate(SurfPressure) do j=js,je ! Work in vertical slices for efficiency From 977ac37b0d5a023013979e2653ce39404bbedc20 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 13 Jul 2023 11:46:33 -0800 Subject: [PATCH 13/14] Need to initialize plume_flux. --- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 8f60e7a979..52de3d1663 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1147,6 +1147,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Idt = 1.0 / dt + plume_flux = 0.0 calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS)) calculate_buoyancy = present(SkinBuoyFlux) From 9129d1f19d2b26a3236d312b8a6e9fa34f4affe8 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 13 Jul 2023 11:55:00 -0800 Subject: [PATCH 14/14] Fixing the line length on one line --- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 52de3d1663..689df230d1 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1430,14 +1430,16 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t plume_flux = 0.0 endif if (dK(i) <= mixing_depth(i) .and. fraction_left_brine > 0.0) then - plume_fraction = min(fraction_left_brine, A_brine(i) * dK(i) ** CS%brine_plume_n * h(i,j,k) * GV%H_to_Z) + plume_fraction = min(fraction_left_brine, A_brine(i) * dK(i) ** CS%brine_plume_n & + * h(i,j,k) * GV%H_to_Z) else IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) ) ! plume_fraction = fraction_left_brine, unless h2d is less than IforcingDepthScale. plume_fraction = min(fraction_left_brine, h2d(i,k)*IforcingDepthScale) endif fraction_left_brine = fraction_left_brine - plume_fraction - plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H + plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) & + * GV%RZ_to_H else plume_flux = 0.0 endif