From b1cf0d20f835bf0e12093dbc3c1422a24fee75d9 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 31 May 2023 14:05:39 -0400 Subject: [PATCH 1/4] 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 2/4] 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 3/4] 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 c010f0e3007d8cd6122c0f96246180ba8b6f274d Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 16 Jun 2023 10:40:57 -0400 Subject: [PATCH 4/4] 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 | 60 ++++++++++++++++++++-- 1 file changed, 56 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..0b45b752ae 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,64 @@ 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 + + ! 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, 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.)