Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into f2023
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward committed Sep 12, 2024
2 parents 691ceab + 5fc90eb commit 6b632ac
Show file tree
Hide file tree
Showing 2 changed files with 92 additions and 77 deletions.
11 changes: 8 additions & 3 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down
158 changes: 84 additions & 74 deletions src/ice_shelf/MOM_ice_shelf.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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].
Expand All @@ -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
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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)
Expand All @@ -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

Expand Down

0 comments on commit 6b632ac

Please sign in to comment.