From 5fc90ebcf25ce99693e4805f64978c18c69d36d3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 17 Aug 2024 06:58:47 -0400 Subject: [PATCH] Rotate ice shelf forcing and initialization Added the ability to rotate the ice shelf forcing and initialization. The specific changes include correcting the grid passed to a call to initialize_ice_shelf in initialize_MOM, using a temporary mech_forcing type in add_shelf_forces when the grid is rotated, uncommenting and correcting the draft rotate_index code in initialize_ice_shelf, and correcting some error-checking and the grid passed to add_shelf_fluxes in initialize_ice_shelf_forces. In addition, unused hor_index_type elements were removed from the ice_shelf_CS, and some callTree calls were added to save_MOM_restart and initialize_ice_shelf. Without rotation, all answers are bitwise identical, but some cases that would have failed previously will now work. --- src/core/MOM.F90 | 11 ++- src/ice_shelf/MOM_ice_shelf.F90 | 158 +++++++++++++++++--------------- 2 files changed, 92 insertions(+), 77 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2f5ac6e489..47971029cc 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2955,20 +2955,20 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & ! These arrays are not initialized in most solo cases, but are needed ! when using an ice shelf. Passing the ice shelf diagnostics CS from MOM ! for legacy reasons. The actual ice shelf diag CS is internal to the ice shelf - call initialize_ice_shelf(param_file, G_in, Time, ice_shelf_CSp, diag_ptr, & + call initialize_ice_shelf(param_file, G, Time, ice_shelf_CSp, diag_ptr, & Time_init, dirs%output_directory, calve_ice_shelf_bergs=point_calving) allocate(frac_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0) allocate(mass_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0) allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0) allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0) - call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf) + call ice_shelf_query(ice_shelf_CSp, G, CS%frac_shelf_h, CS%mass_shelf) ! MOM_initialize_state is using the unrotated metric call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in) call rotate_array(CS%mass_shelf, -turns, mass_shelf_in) call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in, & - frac_shelf_h=frac_shelf_in, mass_shelf = mass_shelf_in) + frac_shelf_h=frac_shelf_in, mass_shelf=mass_shelf_in) else call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & @@ -4173,9 +4173,14 @@ subroutine save_MOM_restart(CS, directory, time, G, time_stamped, filename, & logical, optional, intent(in) :: write_IC !< If present and true, initial conditions are being written + logical :: showCallTree + showCallTree = callTree_showQuery() + + if (showCallTree) call callTree_waypoint("About to call save_restart (step_MOM)") call save_restart(directory, time, G, CS%restart_CS, & time_stamped=time_stamped, filename=filename, GV=GV, & num_rest_files=num_rest_files, write_IC=write_IC) + if (showCallTree) call callTree_waypoint("Done with call to save_restart (step_MOM)") if (CS%use_particles) call particles_save_restart(CS%particles, CS%h, directory, time, time_stamped) end subroutine save_MOM_restart diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index eb021d8ffb..4f811cac87 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -24,6 +24,8 @@ module MOM_ice_shelf use MOM_domains, only : TO_ALL, CGRID_NE, BGRID_NE, CORNER use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe +use MOM_error_handler, only : callTree_showQuery +use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : read_param, get_param, log_param, log_version, param_file_type use MOM_grid, only : MOM_grid_init, ocean_grid_type use MOM_grid_initialize, only : set_grid_metrics @@ -39,12 +41,12 @@ module MOM_ice_shelf use MOM_restart, only : restart_init, restore_state, MOM_restart_CS, register_restart_pair use MOM_time_manager, only : time_type, time_type_to_real, real_to_time, operator(>), operator(-) use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid -use MOM_transcribe_grid, only : rotate_dyngrid +use MOM_transcribe_grid, only : rotate_dyngrid use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init, fix_restart_unit_scaling use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state use MOM_variables, only : rotate_surface_state use MOM_forcing_type, only : forcing, allocate_forcing_type, deallocate_forcing_type, MOM_forcing_chksum -use MOM_forcing_type, only : mech_forcing, allocate_mech_forcing, MOM_mech_forcing_chksum +use MOM_forcing_type, only : mech_forcing, allocate_mech_forcing, deallocate_mech_forcing, MOM_mech_forcing_chksum use MOM_forcing_type, only : copy_common_forcing_fields, rotate_forcing, rotate_mech_forcing use MOM_get_input, only : directories, Get_MOM_input use MOM_EOS, only : calculate_density, calculate_density_derivs, calculate_TFreeze, EOS_domain @@ -59,7 +61,6 @@ module MOM_ice_shelf use MOM_ice_shelf_state, only : ice_shelf_state, ice_shelf_state_end, ice_shelf_state_init use user_shelf_init, only : USER_initialize_shelf_mass, USER_update_shelf_mass use user_shelf_init, only : user_ice_shelf_CS -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 MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init @@ -90,10 +91,6 @@ module MOM_ice_shelf type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart control !! structure for the ice shelves type(ocean_grid_type), pointer :: Grid_in => NULL() !< un-rotated input grid metric - type(hor_index_type), pointer :: HI_in => NULL() !< Pointer to a horizontal indexing structure for - !! incoming data which has not been rotated. - type(hor_index_type), pointer :: HI => NULL() !< Pointer to a horizontal indexing structure for - !! incoming data which has not been rotated. logical :: rotate_index = .false. !< True if index map is rotated integer :: turns !< The number of quarter turns for rotation testing. type(ocean_grid_type), pointer :: Grid => NULL() !< Grid for the ice-shelf model @@ -1020,18 +1017,18 @@ end subroutine change_thickness_using_melt !> This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on !! the ice state in ice_shelf_CS. -subroutine add_shelf_forces(Ocn_grid, US, CS, forces, do_shelf_area, external_call) +subroutine add_shelf_forces(Ocn_grid, US, CS, forces_in, do_shelf_area, external_call) type(ocean_grid_type), intent(in) :: Ocn_grid !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ice_shelf_CS), pointer :: CS !< This module's control structure. - type(mech_forcing), intent(inout) :: forces !< A structure with the + type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the !! driving mechanical forces logical, optional, intent(in) :: do_shelf_area !< If true find the shelf-covered areas. logical, optional, intent(in) :: external_call !< If true the incoming forcing type - !! is using the input grid metric and needs + !! is using the unrotated input grid and may need !! to be rotated. type(ocean_grid_type), pointer :: G => NULL() !< A pointer to the ocean grid metric. -! type(mech_forcing), target :: forces !< A structure with the driving mechanical forces + type(mech_forcing), pointer :: forces !< A structure with the driving mechanical forces real :: kv_rho_ice ! The viscosity of ice divided by its density [L4 T-1 R-1 Z-2 ~> m5 kg-1 s-1]. real :: press_ice ! The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa]. logical :: find_area ! If true find the shelf areas at u & v points. @@ -1041,29 +1038,25 @@ subroutine add_shelf_forces(Ocn_grid, US, CS, forces, do_shelf_area, external_ca integer :: i, j, is, ie, js, je, isd, ied, jsd, jed - if (present(external_call)) rotate=external_call - - if ((Ocn_grid%isc /= CS%Grid_in%isc) .or. (Ocn_grid%iec /= CS%Grid_in%iec) .or. & - (Ocn_grid%jsc /= CS%Grid_in%jsc) .or. (Ocn_grid%jec /= CS%Grid_in%jec)) & - call MOM_error(FATAL,"add_shelf_forces: Incompatible Ocean and Ice shelf grids.") + rotate = .false. ; if (present(external_call)) rotate = external_call if (CS%rotate_index .and. rotate) then - call MOM_error(FATAL,"add_shelf_forces: Rotation not implemented for ice shelves.") - ! allocate(forces) - ! call allocate_mech_forcing(forces_in, CS%Grid, forces) - ! call rotate_mech_forcing(forces_in, CS%turns, forces) - ! else - ! if ((Ocn_grid%isc /= CS%Grid%isc) .or. (Ocn_grid%iec /= CS%Grid%iec) .or. & - ! (Ocn_grid%jsc /= CS%Grid%jsc) .or. (Ocn_grid%jec /= CS%Grid%jec)) & - ! call MOM_error(FATAL,"add_shelf_forces: Incompatible Ocean and Ice shelf grids.") - - ! forces=>forces_in + if ((Ocn_grid%isc /= CS%Grid_in%isc) .or. (Ocn_grid%iec /= CS%Grid_in%iec) .or. & + (Ocn_grid%jsc /= CS%Grid_in%jsc) .or. (Ocn_grid%jec /= CS%Grid_in%jec)) & + call MOM_error(FATAL,"add_shelf_forces: Incompatible Ocean and external Ice shelf grids.") + allocate(forces) + call allocate_mech_forcing(forces_in, CS%Grid, forces) + call rotate_mech_forcing(forces_in, CS%turns, forces) + else + if ((Ocn_grid%isc /= CS%Grid%isc) .or. (Ocn_grid%iec /= CS%Grid%iec) .or. & + (Ocn_grid%jsc /= CS%Grid%jsc) .or. (Ocn_grid%jec /= CS%Grid%jec)) & + call MOM_error(FATAL,"add_shelf_forces: Incompatible Ocean and internal Ice shelf grids.") + + forces=>forces_in endif G=>CS%Grid - - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed @@ -1125,10 +1118,10 @@ subroutine add_shelf_forces(Ocn_grid, US, CS, forces, do_shelf_area, external_ca scalar_pair=.true.) endif - ! if (CS%rotate_index .and. rotate) then - ! call rotate_mech_forcing(forces, -CS%turns, forces_in) - ! ! TODO: deallocate mech forcing? - ! endif + if (CS%rotate_index .and. rotate) then + call rotate_mech_forcing(forces, -CS%turns, forces_in) + call deallocate_mech_forcing(forces) + endif end subroutine add_shelf_forces @@ -1411,6 +1404,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, !! the ice-shelf state type(directories) :: dirs type(dyn_horgrid_type), pointer :: dG => NULL() + type(dyn_horgrid_type), pointer :: dG_in => NULL() real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic. real :: dz_ocean_min_float ! The minimum ocean thickness above which the ice shelf is considered ! to be floating when CONST_SEA_LEVEL = True [Z ~> m]. @@ -1422,6 +1416,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, character(len=40) :: mdl = "MOM_ice_shelf" ! This module's name. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq integer :: wd_halos(2) + logical :: showCallTree logical :: read_TideAmp, debug logical :: global_indexing character(len=240) :: Tideamp_file ! Input file names @@ -1463,55 +1458,57 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, ! Set up the ice-shelf domain and grid wd_halos(:)=0 - allocate(CS%Grid) - call MOM_domains_init(CS%Grid%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_,& + allocate(CS%Grid_in) + call MOM_domains_init(CS%Grid_in%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_,& domain_name='MOM_Ice_Shelf_in', US=CS%US) !allocate(CS%Grid_in%HI) !call hor_index_init(CS%Grid%Domain, CS%Grid%HI, param_file, & ! local_indexing=.not.global_indexing) - call MOM_grid_init(CS%Grid, param_file, CS%US) + call MOM_grid_init(CS%Grid_in, param_file, CS%US) - ! if (CS%rotate_index) then + if (CS%rotate_index) then ! ! TODO: Index rotation currently only works when index rotation does not ! ! change the MPI rank of each domain. Resolving this will require a ! ! modification to FMS PE assignment. ! ! For now, we only permit single-core runs. - ! if (num_PEs() /= 1) & - ! call MOM_error(FATAL, "Index rotation is only supported on one PE.") - - ! call get_param(param_file, mdl, "INDEX_TURNS", CS%turns, & - ! "Number of counterclockwise quarter-turn index rotations.", & - ! default=1, debuggingParam=.true.) - ! ! NOTE: If indices are rotated, then CS%Grid and CS%Grid_in must both be initialized. - ! ! If not rotated, then CS%Grid_in and CS%Ggrid are the same grid. - ! allocate(CS%Grid) - ! !allocate(CS%HI) - ! call clone_MOM_domain(CS%Grid_in%Domain, CS%Grid%Domain,turns=CS%turns) - ! call rotate_hor_index(CS%Grid_in%HI, CS%turns, CS%Grid%HI) - ! call MOM_grid_init(CS%Grid, param_file, CS%US, CS%HI) - ! call create_dyn_horgrid(dG, CS%Grid%HI) - ! call create_dyn_horgrid(dG_in, CS%Grid_in%HI) - ! call clone_MOM_domain(CS%Grid_in%Domain, dG_in%Domain) - ! ! Set up the bottom depth, G%D either analytically or from file - ! call set_grid_metrics(dG_in,param_file,CS%US) - ! call MOM_initialize_topography(dG_in%bathyT, CS%Grid_in%max_depth, dG_in, param_file) - ! call rescale_dyn_horgrid_bathymetry(dG_in, CS%US%Z_to_m) - ! call rotate_dyngrid(dG_in, dG, CS%US, CS%turns) - ! call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) - ! else - !CS%Grid=>CS%Grid_in - dG => NULL() - !CS%Grid%HI=>CS%Grid_in%HI - call create_dyn_horgrid(dG, CS%Grid%HI) - call clone_MOM_domain(CS%Grid%Domain,dG%Domain) - call set_grid_metrics(dG,param_file,CS%US) - ! Set up the bottom depth, dG%bathyT, either analytically or from file - call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file, CS%US) - call copy_dyngrid_to_MOM_grid(dG, CS%Grid, CS%US) - call destroy_dyn_horgrid(dG) -! endif - G => CS%Grid ; CS%Grid_in => CS%Grid + if (num_PEs() /= 1) call MOM_error(FATAL, "Index rotation is only supported on one PE.") + + call get_param(param_file, mdl, "INDEX_TURNS", CS%turns, & + "Number of counterclockwise quarter-turn index rotations.", & + default=1, debuggingParam=.true.) + ! NOTE: If indices are rotated, then CS%Grid and CS%Grid_in must both be initialized. + ! If not rotated, then CS%Grid_in and CS%Ggrid are the same grid. + call create_dyn_horgrid(dG_in, CS%Grid_in%HI) + call clone_MOM_domain(CS%Grid_in%Domain, dG_in%Domain) + call set_grid_metrics(dG_in, param_file, CS%US) + ! Set up the bottom depth, dG_in%bathyT, either analytically or from file + call MOM_initialize_topography(dG_in%bathyT, CS%Grid_in%max_depth, dG_in, param_file, CS%US) + call copy_dyngrid_to_MOM_grid(dG_in, CS%Grid_in, CS%US) + + ! Now set up the rotated ice-shelf grid. + allocate(CS%Grid) + call clone_MOM_domain(CS%Grid_in%Domain, CS%Grid%Domain, turns=CS%turns) + call rotate_hor_index(CS%Grid_in%HI, CS%turns, CS%Grid%HI) + call MOM_grid_init(CS%Grid, param_file, CS%US, CS%Grid%HI) + call create_dyn_horgrid(dG, CS%Grid%HI) + call rotate_dyngrid(dG_in, dG, CS%US, CS%turns) + call copy_dyngrid_to_MOM_grid(dG, CS%Grid, CS%US) + + call destroy_dyn_horgrid(dG_in) + call destroy_dyn_horgrid(dG) + else + CS%Grid => CS%Grid_in + dG => NULL() + call create_dyn_horgrid(dG, CS%Grid%HI) + call clone_MOM_domain(CS%Grid%Domain, dG%Domain) + call set_grid_metrics(dG, param_file, CS%US) + ! Set up the bottom depth, dG%bathyT, either analytically or from file + call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file, CS%US) + call copy_dyngrid_to_MOM_grid(dG, CS%Grid, CS%US) + call destroy_dyn_horgrid(dG) + endif + G => CS%Grid allocate(CS%diag) call MOM_IS_diag_mediator_init(G, CS%US, param_file, CS%diag, component='MOM_IceShelf') @@ -1979,8 +1976,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, if (save_IC .and. .not.((dirs%input_filename(1:1) == 'r') .and. & (LEN_TRIM(dirs%input_filename) == 1))) then + showCallTree = callTree_showQuery() + if (showCallTree) call callTree_waypoint("About to call save_restart (MOM_ice_shelf)") call save_restart(dirs%output_directory, CS%Time, CS%Grid_in, CS%restart_CSp, & filename=IC_file, write_ic=.true.) + if (showCallTree) call callTree_waypoint("Done with call to save_restart (MOM_ice_shelf)") endif CS%id_area_shelf_h = register_diag_field('ice_shelf_model', 'area_shelf_h', CS%diag%axesT1, CS%Time, & @@ -2130,16 +2130,23 @@ subroutine initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in) end subroutine initialize_ice_shelf_fluxes +!> Allocate and initialize the ice-shelf forcing elements of a mechanical forcing type. +!! This forcing type is on the unrotated grid that is used outside of the MOM6 ice shelf code. subroutine initialize_ice_shelf_forces(CS, ocn_grid, US, forces_in) type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the driving mechanical forces ! Local variables type(mech_forcing), pointer :: forces => NULL() call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating forces.") + + if ((Ocn_grid%isc /= CS%Grid_in%isc) .or. (Ocn_grid%iec /= CS%Grid_in%iec) .or. & + (Ocn_grid%jsc /= CS%Grid_in%jsc) .or. (Ocn_grid%jec /= CS%Grid_in%jec)) & + call MOM_error(FATAL,"initialize_ice_shelf_forces: Incompatible ocean and external ice shelf grids.") + call allocate_mech_forcing(CS%Grid_in, forces_in, ustar=.true., shelf=.true., press=.true., tau_mag=.true.) if (CS%rotate_index) then allocate(forces) @@ -2149,10 +2156,13 @@ subroutine initialize_ice_shelf_forces(CS, ocn_grid, US, forces_in) forces=>forces_in endif - call add_shelf_forces(ocn_grid, US, CS, forces, do_shelf_area=.not.CS%solo_ice_sheet) + call add_shelf_forces(CS%grid, US, CS, forces, do_shelf_area=.not.CS%solo_ice_sheet, & + external_call=.false.) - if (CS%rotate_index) & + if (CS%rotate_index) then call rotate_mech_forcing(forces, -CS%turns, forces_in) + call deallocate_mech_forcing(forces) + endif end subroutine initialize_ice_shelf_forces