From 1417dceb4e73b6de7c9d18808b14b5e128c09393 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 13 Jan 2021 07:10:45 -0500 Subject: [PATCH] +Added MOM_interpolate.F90 Added the new module MOM_interpolate to wrap the time_interp_external and horiz_interp modules. The overloaded interface time_interp_external had to be renamed to time_interp_extern in the version offered by MOM_interpolate because of a weird compile time problem with the PGI 19.10.0 compiler. Some large blocks of unused code in MOM_horizontal_regridding were commented out. Also modified 6 files outside of the framework directory to use these new interfaces, but they will all work without these changes. All answers are bitwise identical. --- .../MOM_surface_forcing_gfdl.F90 | 49 ++-- config_src/coupled_driver/ocean_model_MOM.F90 | 22 +- config_src/solo_driver/MOM_driver.F90 | 14 +- src/core/MOM_open_boundary.F90 | 9 +- src/framework/MOM_horizontal_regridding.F90 | 234 +++++++++--------- src/framework/MOM_interpolate.F90 | 143 +++++++++++ src/ice_shelf/MOM_ice_shelf.F90 | 14 +- .../vertical/MOM_diabatic_aux.F90 | 8 +- 8 files changed, 313 insertions(+), 180 deletions(-) create mode 100644 src/framework/MOM_interpolate.F90 diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index be960accd6..67f6643a42 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -5,7 +5,7 @@ module MOM_surface_forcing_gfdl !#CTRL# use MOM_controlled_forcing, only : apply_ctrl_forcing, register_ctrl_forcing_restarts !#CTRL# use MOM_controlled_forcing, only : controlled_forcing_init, controlled_forcing_end !#CTRL# use MOM_controlled_forcing, only : ctrl_forcing_CS -use MOM_coms, only : reproducing_sum +use MOM_coms, only : reproducing_sum, field_chksum use MOM_constants, only : hlv, hlf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT @@ -21,6 +21,8 @@ module MOM_surface_forcing_gfdl use MOM_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : init_external_field, time_interp_extern +use MOM_interpolate, only : time_interp_external_init use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS use MOM_restart, only : restart_init_end, save_restart, restore_state @@ -36,9 +38,6 @@ module MOM_surface_forcing_gfdl use coupler_types_mod, only : coupler_type_copy_data use data_override_mod, only : data_override_init, data_override use fms_mod, only : stdout -use mpp_mod, only : mpp_chksum -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init implicit none ; private @@ -350,7 +349,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! Salinity restoring logic if (CS%restore_salt) then - call time_interp_external(CS%id_srestore,Time,data_restore) + call time_interp_extern(CS%id_srestore, Time, data_restore) ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) open_ocn_mask(:,:) = 1.0 if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice @@ -407,7 +406,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! SST restoring logic if (CS%restore_temp) then - call time_interp_external(CS%id_trestore,Time,data_restore) + call time_interp_extern(CS%id_trestore, Time, data_restore) do j=js,je ; do i=is,ie delta_sst = data_restore(i,j)- sfc_state%SST(i,j) delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) @@ -1486,7 +1485,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger) enddo ; enddo endif - call time_interp_external_init + call time_interp_external_init() ! Optionally read a x-y gustiness field in place of a global constant. call get_param(param_file, mdl, "READ_GUST_2D", CS%read_gust_2d, & @@ -1632,27 +1631,27 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) outunit = stdout() write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep - write(outunit,100) 'iobt%u_flux ', mpp_chksum( iobt%u_flux ) - write(outunit,100) 'iobt%v_flux ', mpp_chksum( iobt%v_flux ) - write(outunit,100) 'iobt%t_flux ', mpp_chksum( iobt%t_flux ) - write(outunit,100) 'iobt%q_flux ', mpp_chksum( iobt%q_flux ) - write(outunit,100) 'iobt%salt_flux ', mpp_chksum( iobt%salt_flux ) - write(outunit,100) 'iobt%lw_flux ', mpp_chksum( iobt%lw_flux ) - write(outunit,100) 'iobt%sw_flux_vis_dir', mpp_chksum( iobt%sw_flux_vis_dir) - write(outunit,100) 'iobt%sw_flux_vis_dif', mpp_chksum( iobt%sw_flux_vis_dif) - write(outunit,100) 'iobt%sw_flux_nir_dir', mpp_chksum( iobt%sw_flux_nir_dir) - write(outunit,100) 'iobt%sw_flux_nir_dif', mpp_chksum( iobt%sw_flux_nir_dif) - write(outunit,100) 'iobt%lprec ', mpp_chksum( iobt%lprec ) - write(outunit,100) 'iobt%fprec ', mpp_chksum( iobt%fprec ) - write(outunit,100) 'iobt%runoff ', mpp_chksum( iobt%runoff ) - write(outunit,100) 'iobt%calving ', mpp_chksum( iobt%calving ) - write(outunit,100) 'iobt%p ', mpp_chksum( iobt%p ) + write(outunit,100) 'iobt%u_flux ', field_chksum( iobt%u_flux ) + write(outunit,100) 'iobt%v_flux ', field_chksum( iobt%v_flux ) + write(outunit,100) 'iobt%t_flux ', field_chksum( iobt%t_flux ) + write(outunit,100) 'iobt%q_flux ', field_chksum( iobt%q_flux ) + write(outunit,100) 'iobt%salt_flux ', field_chksum( iobt%salt_flux ) + write(outunit,100) 'iobt%lw_flux ', field_chksum( iobt%lw_flux ) + write(outunit,100) 'iobt%sw_flux_vis_dir', field_chksum( iobt%sw_flux_vis_dir) + write(outunit,100) 'iobt%sw_flux_vis_dif', field_chksum( iobt%sw_flux_vis_dif) + write(outunit,100) 'iobt%sw_flux_nir_dir', field_chksum( iobt%sw_flux_nir_dir) + write(outunit,100) 'iobt%sw_flux_nir_dif', field_chksum( iobt%sw_flux_nir_dif) + write(outunit,100) 'iobt%lprec ', field_chksum( iobt%lprec ) + write(outunit,100) 'iobt%fprec ', field_chksum( iobt%fprec ) + write(outunit,100) 'iobt%runoff ', field_chksum( iobt%runoff ) + write(outunit,100) 'iobt%calving ', field_chksum( iobt%calving ) + write(outunit,100) 'iobt%p ', field_chksum( iobt%p ) if (associated(iobt%ustar_berg)) & - write(outunit,100) 'iobt%ustar_berg ', mpp_chksum( iobt%ustar_berg ) + write(outunit,100) 'iobt%ustar_berg ', field_chksum( iobt%ustar_berg ) if (associated(iobt%area_berg)) & - write(outunit,100) 'iobt%area_berg ', mpp_chksum( iobt%area_berg ) + write(outunit,100) 'iobt%area_berg ', field_chksum( iobt%area_berg ) if (associated(iobt%mass_berg)) & - write(outunit,100) 'iobt%mass_berg ', mpp_chksum( iobt%mass_berg ) + write(outunit,100) 'iobt%mass_berg ', field_chksum( iobt%mass_berg ) 100 FORMAT(" CHECKSUM::",A20," = ",Z20) call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%') diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index b429da649b..12f803a970 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -15,6 +15,7 @@ module ocean_model_mod use MOM, only : extract_surface_state, allocate_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized use MOM, only : get_ocean_stocks, step_offline +use MOM_coms, only : field_chksum use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end @@ -22,6 +23,7 @@ module ocean_model_mod use MOM_domains, only : TO_ALL, Omit_Corners use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave +use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type use MOM_forcing_type, only : forcing, mech_forcing, allocate_forcing_type use MOM_forcing_type, only : fluxes_accumulate, get_net_mass_forcing @@ -48,6 +50,8 @@ module ocean_model_mod use MOM_verticalGrid, only : verticalGrid_type use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart +use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init +use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type use coupler_types_mod, only : coupler_type_spawn, coupler_type_write_chksums use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data @@ -56,10 +60,6 @@ module ocean_model_mod use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux use fms_mod, only : stdout -use mpp_mod, only : mpp_chksum -use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct -use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init -use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves #include @@ -1130,13 +1130,13 @@ subroutine ocean_public_type_chksum(id, timestep, ocn) outunit = stdout() write(outunit,*) "BEGIN CHECKSUM(ocean_type):: ", id, timestep - write(outunit,100) 'ocean%t_surf ',mpp_chksum(ocn%t_surf ) - write(outunit,100) 'ocean%s_surf ',mpp_chksum(ocn%s_surf ) - write(outunit,100) 'ocean%u_surf ',mpp_chksum(ocn%u_surf ) - write(outunit,100) 'ocean%v_surf ',mpp_chksum(ocn%v_surf ) - write(outunit,100) 'ocean%sea_lev ',mpp_chksum(ocn%sea_lev) - write(outunit,100) 'ocean%frazil ',mpp_chksum(ocn%frazil ) - write(outunit,100) 'ocean%melt_potential ',mpp_chksum(ocn%melt_potential) + write(outunit,100) 'ocean%t_surf ', field_chksum(ocn%t_surf ) + write(outunit,100) 'ocean%s_surf ', field_chksum(ocn%s_surf ) + write(outunit,100) 'ocean%u_surf ', field_chksum(ocn%u_surf ) + write(outunit,100) 'ocean%v_surf ', field_chksum(ocn%v_surf ) + write(outunit,100) 'ocean%sea_lev ', field_chksum(ocn%sea_lev) + write(outunit,100) 'ocean%frazil ', field_chksum(ocn%frazil ) + write(outunit,100) 'ocean%melt_potential ', field_chksum(ocn%melt_potential) call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%') 100 FORMAT(" CHECKSUM::",A20," = ",Z20) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 9726aa1281..c9383a4287 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -32,6 +32,7 @@ program MOM_main use MOM, only : extract_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized use MOM, only : step_offline + use MOM_coms, only : Set_PElist use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_error_handler, only : MOM_error, MOM_mesg, WARNING, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint @@ -41,6 +42,7 @@ program MOM_main use MOM_forcing_type, only : mech_forcing_diags, MOM_forcing_chksum, MOM_mech_forcing_chksum use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type + use MOM_interpolate, only : time_interp_external_init use MOM_io, only : file_exists, open_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end use MOM_io, only : APPEND_FILE, ASCII_FILE, READONLY_FILE, SINGLE_FILE @@ -64,8 +66,6 @@ program MOM_main use MOM_get_input, only : get_MOM_input use ensemble_manager_mod, only : ensemble_manager_init, get_ensemble_size use ensemble_manager_mod, only : ensemble_pelist_setup - use mpp_mod, only : set_current_pelist => mpp_set_current_pelist - use time_interp_external_mod, only : time_interp_external_init use fms_affinity_mod, only : fms_affinity_init, fms_affinity_set,fms_affinity_get use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS @@ -229,7 +229,7 @@ program MOM_main allocate(ocean_pelist(nPEs_per)) call ensemble_pelist_setup(.true., 0, nPEs_per, 0, 0, atm_pelist, ocean_pelist, & land_pelist, ice_pelist) - call set_current_pelist(ocean_pelist) + call Set_PElist(ocean_pelist) deallocate(ocean_pelist) endif @@ -286,17 +286,17 @@ program MOM_main if (sum(date_init) > 0) then - Start_time = set_date(date_init(1),date_init(2), date_init(3), & - date_init(4),date_init(5),date_init(6)) + Start_time = set_date(date_init(1), date_init(2), date_init(3), & + date_init(4), date_init(5), date_init(6)) else Start_time = real_to_time(0.0) endif - call time_interp_external_init + call time_interp_external_init() if (sum(date) >= 0) then ! In this case, the segment starts at a time fixed by ocean_solo.res - segment_start_time = set_date(date(1),date(2),date(3),date(4),date(5),date(6)) + segment_start_time = set_date(date(1), date(2), date(3), date(4), date(5), date(6)) Time = segment_start_time else ! In this case, the segment starts at a time read from the MOM restart file diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 3fc0c9bcba..4faab21177 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -24,8 +24,7 @@ module MOM_open_boundary use MOM_tidal_forcing, only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency use MOM_time_manager, only : set_date, time_type, time_type_to_real, operator(-) use MOM_tracer_registry, only : tracer_type, tracer_registry_type, tracer_name_lookup -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init +use MOM_interpolate, only : init_external_field, time_interp_extern, time_interp_external_init use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS use MOM_remapping, only : initialize_remapping, remapping_core_h, end_remapping use MOM_regridding, only : regridding_CS @@ -3897,7 +3896,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! TODO: Since we conditionally rotate a subset of tmp_buffer_in after ! reading the value, it is currently not possible to use the rotated - ! implementation of time_interp_external. + ! implementation of time_interp_extern. ! For now, we must explicitly allocate and rotate this array. if (turns /= 0) then if (modulo(turns, 2) /= 0) then @@ -3909,7 +3908,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) tmp_buffer_in => tmp_buffer endif - call time_interp_external(segment%field(m)%fid,Time, tmp_buffer_in) + call time_interp_extern(segment%field(m)%fid,Time, tmp_buffer_in) ! NOTE: Rotation of face-points require that we skip the final value if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. @@ -3976,7 +3975,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ! no dz for tidal variables if (segment%field(m)%nk_src > 1 .and.& (index(segment%field(m)%name, 'phase') .le. 0 .and. index(segment%field(m)%name, 'amp') .le. 0)) then - call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer_in) + call time_interp_extern(segment%field(m)%fid_dz,Time, tmp_buffer_in) if (turns /= 0) then ! TODO: This is hardcoded for 90 degrees, and needs to be generalized. if (segment%is_E_or_W & diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 9b340f3aa7..6e09d391d9 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -13,12 +13,8 @@ module MOM_horizontal_regridding use MOM_grid, only : ocean_grid_type use MOM_io_wrapper, only : axistype, get_axis_data use MOM_time_manager, only : time_type -use MOM_time_manager, only : init_external_field, get_external_field_size -use MOM_time_manager, only : get_external_field_axes, get_external_field_missing -use MOM_transform_FMS, only : time_interp_external => rotated_time_interp_external - -use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type -use horiz_interp_mod, only : horiz_interp_init, horiz_interp_del +use MOM_interpolate, only : time_interp_extern, get_external_field_info, horiz_interp_init +use MOM_interpolate, only : horiz_interp_new, horiz_interp, horiz_interp_type use netcdf @@ -31,10 +27,10 @@ module MOM_horizontal_regridding ! character(len=40) :: mdl = "MOM_horizontal_regridding" ! This module's name. !> Fill grid edges -interface fill_boundaries - module procedure fill_boundaries_real - module procedure fill_boundaries_int -end interface +! interface fill_boundaries +! module procedure fill_boundaries_real +! module procedure fill_boundaries_int +! end interface !> Extrapolate and interpolate data interface horiz_interp_and_extrap_tracer @@ -296,7 +292,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real :: PI_180 integer :: rcode, ncid, varid, ndims, id, jd, kd, jdp - integer :: i,j,k + integer :: i, j, k integer, dimension(4) :: start, count, dims, dim_id real, dimension(:,:), allocatable :: x_in, y_in real, dimension(:), allocatable :: lon_in, lat_in @@ -309,7 +305,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, character(len=8) :: laynum type(horiz_interp_type) :: Interp integer :: is, ie, js, je ! compute domain indices - integer :: isc,iec,jsc,jec ! global compute domain indices + integer :: isc, iec, jsc, jec ! global compute domain indices integer :: isg, ieg, jsg, jeg ! global extent integer :: isd, ied, jsd, jed ! data domain indices integer :: id_clock_read @@ -318,9 +314,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real :: npoints,varAvg real, dimension(SZI_(G),SZJ_(G)) :: lon_out, lat_out, tr_out, mask_out real, dimension(SZI_(G),SZJ_(G)) :: good, fill - real, dimension(SZI_(G),SZJ_(G)) :: tr_outf,tr_prev - real, dimension(SZI_(G),SZJ_(G)) :: good2,fill2 - real, dimension(SZI_(G),SZJ_(G)) :: nlevs + real, dimension(SZI_(G),SZJ_(G)) :: tr_outf, tr_prev + real, dimension(SZI_(G),SZJ_(G)) :: good2, fill2 + real, dimension(SZI_(G),SZJ_(G)) :: nlevs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -328,14 +324,14 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, id_clock_read = cpu_clock_id('(Initialize tracer from Z) read', grain=CLOCK_LOOP) - is_ongrid=.false. - if (present(ongrid)) is_ongrid=ongrid + is_ongrid = .false. + if (present(ongrid)) is_ongrid = ongrid if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) if (allocated(z_edges_in)) deallocate(z_edges_in) - PI_180=atan(1.0)/45. + PI_180 = atan(1.0)/45. ! Open NetCDF file and if present, extract data and spatial coordinate information ! The convention adopted here requires that the data be written in (i,j,k) ordering. @@ -391,7 +387,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) - allocate(lon_in(id),lat_in(jd),z_in(kd),z_edges_in(kd+1)) + allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1)) allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd)) start = 1 ; count = 1 ; count(1) = id @@ -692,7 +688,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t call cpu_clock_begin(id_clock_read) - fld_sz = get_external_field_size(fms_id) + call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_value) if (allocated(lon_in)) deallocate(lon_in) if (allocated(lat_in)) deallocate(lat_in) if (allocated(z_in)) deallocate(z_in) @@ -700,8 +696,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (allocated(tr_z)) deallocate(tr_z) if (allocated(mask_z)) deallocate(mask_z) - axes_data = get_external_field_axes(fms_id) - id = fld_sz(1) ; jd = fld_sz(2) ; kd = fld_sz(3) spongeDataOngrid=.false. @@ -722,8 +716,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t call cpu_clock_end(id_clock_read) - missing_value = get_external_field_missing(fms_id) - if (.not. spongeDataOngrid) then ! extrapolate the input data to the north pole using the northerm-most latitude max_lat = maxval(lat_in) @@ -773,7 +765,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t if (.not.spongeDataOngrid) then if (is_root_pe()) & - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) + call time_interp_extern(fms_id, Time, data_in, verbose=.true., 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 @@ -891,7 +883,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t enddo ! kd else - call time_interp_external(fms_id, Time, data_in, verbose=.true., turns=turns) + call time_interp_extern(fms_id, Time, data_in, verbose=.true., turns=turns) do k=1,kd do j=js,je do i=is,ie @@ -927,55 +919,55 @@ end subroutine meshgrid ! None of the subsequent code appears to be used at all. -!> Fill grid edges for integer data -function fill_boundaries_int(m,cyclic_x,tripolar_n) result(mp) - integer, dimension(:,:), intent(in) :: m !< input array (ND) - logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant - logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold - integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp +! !> Fill grid edges for integer data +! function fill_boundaries_int(m,cyclic_x,tripolar_n) result(mp) +! integer, dimension(:,:), intent(in) :: m !< input array (ND) +! logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant +! logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold +! integer, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - real, dimension(size(m,1),size(m,2)) :: m_real - real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real +! real, dimension(size(m,1),size(m,2)) :: m_real +! real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp_real - m_real = real(m) +! m_real = real(m) - mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n) +! mp_real = fill_boundaries_real(m_real,cyclic_x,tripolar_n) - mp = int(mp_real) +! mp = int(mp_real) -end function fill_boundaries_int +! end function fill_boundaries_int !> Fill grid edges for real data -function fill_boundaries_real(m,cyclic_x,tripolar_n) result(mp) - real, dimension(:,:), intent(in) :: m !< input array (ND) - logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant - logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold - real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp +! function fill_boundaries_real(m,cyclic_x,tripolar_n) result(mp) +! real, dimension(:,:), intent(in) :: m !< input array (ND) +! logical, intent(in) :: cyclic_x !< True if domain is zonally re-entrant +! logical, intent(in) :: tripolar_n !< True if domain has an Arctic fold +! real, dimension(0:size(m,1)+1,0:size(m,2)+1) :: mp - integer :: ni,nj,i,j +! integer :: ni,nj,i,j - ni=size(m,1); nj=size(m,2) +! ni=size(m,1); nj=size(m,2) - mp(1:ni,1:nj)=m(:,:) +! mp(1:ni,1:nj)=m(:,:) - if (cyclic_x) then - mp(0,1:nj)=m(ni,1:nj) - mp(ni+1,1:nj)=m(1,1:nj) - else - mp(0,1:nj)=m(1,1:nj) - mp(ni+1,1:nj)=m(ni,1:nj) - endif +! if (cyclic_x) then +! mp(0,1:nj)=m(ni,1:nj) +! mp(ni+1,1:nj)=m(1,1:nj) +! else +! mp(0,1:nj)=m(1,1:nj) +! mp(ni+1,1:nj)=m(ni,1:nj) +! endif - mp(1:ni,0)=m(1:ni,1) - if (tripolar_n) then - do i=1,ni - mp(i,nj+1)=m(ni-i+1,nj) - enddo - else - mp(1:ni,nj+1)=m(1:ni,nj) - endif +! mp(1:ni,0)=m(1:ni,1) +! if (tripolar_n) then +! do i=1,ni +! mp(i,nj+1)=m(ni-i+1,nj) +! enddo +! else +! mp(1:ni,nj+1)=m(1:ni,nj) +! endif -end function fill_boundaries_real +! end function fill_boundaries_real !> Solve del2 (zi) = 0 using successive iterations !! with a 5 point stencil. Only points fill==1 are @@ -983,64 +975,64 @@ end function fill_boundaries_real !! isotropically in index space. The resulting solution !! in each region is an approximation to del2(zi)=0 subject to !! boundary conditions along the valid points curve bounding this region. -subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) - real, dimension(:,:), intent(inout) :: zi !< input and output array (ND) - integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< same shape as zi, 1=fill - integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< same shape as zi, 1=bad data - real, intent(in) :: sor !< relaxation coefficient (ND) - integer, intent(in) :: niter !< maximum number of iterations - logical, intent(in) :: cyclic_x !< true if domain is zonally reentrant - logical, intent(in) :: tripolar_n !< true if domain has an Arctic fold - - ! Local variables - real, dimension(size(zi,1),size(zi,2)) :: res, m - integer, dimension(size(zi,1),size(zi,2),4) :: B - real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp - integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm - integer :: i,j,k,n - integer :: ni,nj - real :: Isum, bsum - - ni=size(zi,1) ; nj=size(zi,2) - - - mp(:,:) = fill_boundaries(zi,cyclic_x,tripolar_n) - - B(:,:,:) = 0.0 - nm(:,:) = fill_boundaries(bad,cyclic_x,tripolar_n) - - do j=1,nj - do i=1,ni - if (fill(i,j) == 1) then - B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) - B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) - endif - enddo - enddo - - do n=1,niter - do j=1,nj - do i=1,ni - if (fill(i,j) == 1) then - bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) - Isum = 1.0/bsum - res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& - B(i,j,3)*mp(i,j+1)+B(i,j,4)*mp(i,j-1)) - mp(i,j) - endif - enddo - enddo - res(:,:)=res(:,:)*sor - - do j=1,nj - do i=1,ni - mp(i,j)=mp(i,j)+res(i,j) - enddo - enddo - - zi(:,:)=mp(1:ni,1:nj) - mp = fill_boundaries(zi,cyclic_x,tripolar_n) - enddo - -end subroutine smooth_heights +! subroutine smooth_heights(zi,fill,bad,sor,niter,cyclic_x, tripolar_n) +! real, dimension(:,:), intent(inout) :: zi !< input and output array (ND) +! integer, dimension(size(zi,1),size(zi,2)), intent(in) :: fill !< same shape as zi, 1=fill +! integer, dimension(size(zi,1),size(zi,2)), intent(in) :: bad !< same shape as zi, 1=bad data +! real, intent(in) :: sor !< relaxation coefficient (ND) +! integer, intent(in) :: niter !< maximum number of iterations +! logical, intent(in) :: cyclic_x !< true if domain is zonally reentrant +! logical, intent(in) :: tripolar_n !< true if domain has an Arctic fold + +! ! Local variables +! real, dimension(size(zi,1),size(zi,2)) :: res, m +! integer, dimension(size(zi,1),size(zi,2),4) :: B +! real, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: mp +! integer, dimension(0:size(zi,1)+1,0:size(zi,2)+1) :: nm +! integer :: i,j,k,n +! integer :: ni,nj +! real :: Isum, bsum + +! ni=size(zi,1) ; nj=size(zi,2) + + +! mp(:,:) = fill_boundaries(zi,cyclic_x,tripolar_n) + +! B(:,:,:) = 0.0 +! nm(:,:) = fill_boundaries(bad,cyclic_x,tripolar_n) + +! do j=1,nj +! do i=1,ni +! if (fill(i,j) == 1) then +! B(i,j,1)=1-nm(i+1,j);B(i,j,2)=1-nm(i-1,j) +! B(i,j,3)=1-nm(i,j+1);B(i,j,4)=1-nm(i,j-1) +! endif +! enddo +! enddo + +! do n=1,niter +! do j=1,nj +! do i=1,ni +! if (fill(i,j) == 1) then +! bsum = real(B(i,j,1)+B(i,j,2)+B(i,j,3)+B(i,j,4)) +! Isum = 1.0/bsum +! res(i,j)=Isum*(B(i,j,1)*mp(i+1,j)+B(i,j,2)*mp(i-1,j)+& +! B(i,j,3)*mp(i,j+1)+B(i,j,4)*mp(i,j-1)) - mp(i,j) +! endif +! enddo +! enddo +! res(:,:)=res(:,:)*sor + +! do j=1,nj +! do i=1,ni +! mp(i,j)=mp(i,j)+res(i,j) +! enddo +! enddo + +! zi(:,:)=mp(1:ni,1:nj) +! mp = fill_boundaries(zi,cyclic_x,tripolar_n) +! enddo + +! end subroutine smooth_heights end module MOM_horizontal_regridding diff --git a/src/framework/MOM_interpolate.F90 b/src/framework/MOM_interpolate.F90 new file mode 100644 index 0000000000..c63c847e55 --- /dev/null +++ b/src/framework/MOM_interpolate.F90 @@ -0,0 +1,143 @@ +!> This module wraps the FMS temporal and spatial interpolation routines +module MOM_interpolate + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_array_transform, only : allocate_rotated_array, rotate_array +use MOM_error_handler, only : MOM_error, FATAL +use MOM_io_wrapper, only : axistype +use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_init, horiz_interp_type +use time_interp_external_mod, only : time_interp_external_fms=>time_interp_external +use time_interp_external_mod, only : init_external_field, time_interp_external_init +use time_interp_external_mod, only : get_external_field_size +use time_interp_external_mod, only : get_external_field_axes, get_external_field_missing +use time_manager_mod, only : time_type + +implicit none ; private + +public :: time_interp_extern, init_external_field, time_interp_external_init +public :: get_external_field_info +public :: horiz_interp_type, horiz_interp_init, horiz_interp, horiz_interp_new + +!> Read a field based on model time, and rotate to the model domain. +! This inerface does not share the name time_interp_external with the module it primarily +! wraps because of errors (perhaps a bug) that arise with the PGI 19.10.0 compiler. +interface time_interp_extern + module procedure time_interp_external_0d + module procedure time_interp_external_2d + module procedure time_interp_external_3d +end interface time_interp_extern + +contains + +!> Get information about the external fields. +subroutine get_external_field_info(field_id, size, axes, missing) + integer, intent(in) :: field_id !< The integer index of the external + !! field returned from a previous + !! call to init_external_field() + integer, dimension(4), optional, intent(inout) :: size !< Dimension sizes for the input data + type(axistype), dimension(4), optional, intent(inout) :: axes !< Axis types for the input data + real, optional, intent(inout) :: missing !< Missing value for the input data + + if (present(size)) then + size(1:4) = get_external_field_size(field_id) + endif + + if (present(axes)) then + axes(1:4) = get_external_field_axes(field_id) + endif + + if (present(missing)) then + missing = get_external_field_missing(field_id) + endif + +end subroutine get_external_field_info + + +!> Read a scalar field based on model time. +subroutine time_interp_external_0d(field_id, time, data_in, verbose) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, intent(inout) :: data_in !< The interpolated value + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + + call time_interp_external_fms(field_id, time, data_in, verbose=verbose) +end subroutine time_interp_external_0d + +!> Read a 2d field from an external based on model time, potentially including horizontal +!! interpolation and rotation of the data +subroutine time_interp_external_2d(field_id, time, data_in, interp, verbose, horz_interp, mask_out, turns) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, dimension(:,:), intent(inout) :: data_in !< The array in which to store the interpolated values + integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + type(horiz_interp_type), & + optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation + logical, dimension(:,:), & + optional, intent(out) :: mask_out !< An array that is true where there is valid data + integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data + + real, allocatable :: data_pre_rot(:,:) ! The input data before rotation + integer :: qturns ! The number of quarter turns to rotate the data + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call time_interp_external_fms(field_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + else + call allocate_rotated_array(data_in, [1,1], -qturns, data_pre_rot) + call time_interp_external_fms(field_id, time, data_pre_rot, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + call rotate_array(data_pre_rot, turns, data_in) + deallocate(data_pre_rot) + endif +end subroutine time_interp_external_2d + + +!> Read a 3d field based on model time, and rotate to the model grid +subroutine time_interp_external_3d(field_id, time, data_in, interp, & + verbose, horz_interp, mask_out, turns) + integer, intent(in) :: field_id !< The integer index of the external field returned + !! from a previous call to init_external_field() + type(time_type), intent(in) :: time !< The target time for the data + real, dimension(:,:,:), intent(inout) :: data_in !< The array in which to store the interpolated values + integer, optional, intent(in) :: interp !< A flag indicating the temporal interpolation method + logical, optional, intent(in) :: verbose !< If true, write verbose output for debugging + type(horiz_interp_type), & + optional, intent(in) :: horz_interp !< A structure to control horizontal interpolation + logical, dimension(:,:,:), & + optional, intent(out) :: mask_out !< An array that is true where there is valid data + integer, optional, intent(in) :: turns !< Number of quarter turns to rotate the data + + real, allocatable :: data_pre_rot(:,:,:) ! The input data before rotation + integer :: qturns ! The number of quarter turns to rotate the data + + ! TODO: Mask rotation requires logical array rotation support + if (present(mask_out)) & + call MOM_error(FATAL, "Rotation of masked output not yet support") + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call time_interp_external_fms(field_id, time, data_in, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + else + call allocate_rotated_array(data_in, [1,1,1], -qturns, data_pre_rot) + call time_interp_external_fms(field_id, time, data_pre_rot, interp=interp, & + verbose=verbose, horz_interp=horz_interp) + call rotate_array(data_pre_rot, turns, data_in) + deallocate(data_pre_rot) + endif +end subroutine time_interp_external_3d + +end module MOM_interpolate diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index a5dd92b640..634ff80325 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -57,8 +57,8 @@ module MOM_ice_shelf use MOM_coms, only : reproducing_sum use MOM_spatial_means, only : global_area_integral use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init +use MOM_interpolate, only : init_external_field, time_interp_extern, time_interp_external_init + implicit none ; private #include @@ -1084,9 +1084,9 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) do j=js,je ; do i=is,ie last_hmask(i,j) = ISS%hmask(i,j) ; last_area_shelf_h(i,j) = ISS%area_shelf_h(i,j) enddo ; enddo - call time_interp_external(CS%id_read_mass, Time0, last_mass_shelf) + call time_interp_extern(CS%id_read_mass, Time0, last_mass_shelf) do j=js,je ; do i=is,ie - ! This should only be done if time_interp_external did an update. + ! This should only be done if time_interp_extern did an update. last_mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * last_mass_shelf(i,j) ! Rescale after time_interp last_h_shelf(i,j) = last_mass_shelf(i,j) / CS%density_ice enddo ; enddo @@ -1512,7 +1512,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, if (CS%rotate_index) then allocate(tmp2d(CS%Grid_in%isd:CS%Grid_in%ied,CS%Grid_in%jsd:CS%Grid_in%jed));tmp2d(:,:)=0.0 call MOM_read_data(TideAmp_file, 'tideamp', tmp2d, CS%Grid_in%domain, timelevel=1, scale=US%m_s_to_L_T) - call rotate_array(tmp2d,CS%turns, CS%utide) + call rotate_array(tmp2d, CS%turns, CS%utide) deallocate(tmp2d) else call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) @@ -1984,8 +1984,8 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) allocate(tmp2d(is:ie,js:je)) ; tmp2d(:,:) = 0.0 endif - call time_interp_external(CS%id_read_mass, Time, tmp2d) - call rotate_array(tmp2d,CS%turns, ISS%mass_shelf) + call time_interp_extern(CS%id_read_mass, Time, tmp2d) + call rotate_array(tmp2d, CS%turns, ISS%mass_shelf) deallocate(tmp2d) ! This should only be done if time_interp_external did an update. diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 470098a08a..8ba9dd959a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -15,15 +15,15 @@ module MOM_diabatic_aux use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint use MOM_grid, only : ocean_grid_type +use MOM_interpolate, only : init_external_field, time_interp_extern +use MOM_interpolate, only : time_interp_external_init use MOM_io, only : slasher use MOM_opacity, only : set_opacity, opacity_CS, extract_optics_slice, extract_optics_fields use MOM_opacity, only : optics_type, optics_nbands, absorbRemainingSW, sumSWoverBands use MOM_tracer_flow_control, only : get_chl_from_model, tracer_flow_control_CS use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs ! , vertvisc_type, accel_diag_ptrs +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use time_interp_external_mod, only : init_external_field, time_interp_external -use time_interp_external_mod, only : time_interp_external_init implicit none ; private @@ -621,7 +621,7 @@ subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity_CSp, tracer_ if (CS%chl_from_file) then ! Only the 2-d surface chlorophyll can be read in from a file. The ! same value is assumed for all layers. - call time_interp_external(CS%sbc_chl, CS%Time, chl_2d) + call time_interp_extern(CS%sbc_chl, CS%Time, chl_2d) do j=js,je ; do i=is,ie if ((G%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&