From 7a2483e8466687c837a4aa8a159e47738e1bcd77 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Wed, 23 Aug 2023 18:19:39 -0400 Subject: [PATCH 1/2] Ice-shelf solo driver and MISMIP+ updates -Several edits to the ice shelf solo driver so that it works with the rest of the current MOM6 -Added capability to initialize a surface mass balance (SMB) that is held contstant over time when running from the ice-shelf solo driver (see new subroutine initialize_ice_SMB). This is required for MISMIP+. A constant SMB can also be used from the MOM driver for coupled ice-shelf/ocean experiments (e.g. MISOMIP). -The new, constant SMB is passed into solo_step_ice_shelf, where change_thickness_using_precip is called -Added capability to save both non-time-stamped and time-stamped restart files when using the ice shelf solo driver. This is useful for debugging. -slight reorganization to when ice shelf post_data calls are made --- .../ice_solo_driver/ice_shelf_driver.F90 | 61 +++++++++--- config_src/drivers/solo_driver/MOM_driver.F90 | 6 +- src/framework/MOM_diag_mediator.F90 | 95 +++++++++++++++++++ src/ice_shelf/MOM_ice_shelf.F90 | 20 ++-- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 2 +- src/ice_shelf/MOM_ice_shelf_initialize.F90 | 49 +++++++++- 6 files changed, 209 insertions(+), 24 deletions(-) diff --git a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 index 8ea0867d03..b864903ddc 100644 --- a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 @@ -24,8 +24,8 @@ program Shelf_main use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT use MOM_debugging, only : MOM_debugging_init - use MOM_diag_mediator, only : diag_mediator_init, diag_mediator_infrastructure_init - use MOM_diag_mediator, only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration + use MOM_diag_mediator, only : diag_mediator_init, diag_mediator_infrastructure_init, set_axes_info + use MOM_diag_mediator, only : solo_ice_shelf_diag_mediator_end, diag_ctrl, diag_mediator_close_registration use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_domains, only : MOM_domains_init, clone_MOM_domain, pass_var use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid @@ -54,6 +54,8 @@ program Shelf_main use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd use MOM_write_cputime, only : write_cputime, MOM_write_cputime_init use MOM_write_cputime, only : write_cputime_start_clock, write_cputime_CS + use MOM_forcing_type, only : forcing + use MOM_ice_shelf_initialize, only : initialize_ice_SMB use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS use MOM_ice_shelf, only : ice_shelf_save_restart, solo_step_ice_shelf @@ -75,7 +77,9 @@ program Shelf_main ! CPU time limit. nmax is determined by evaluating the CPU time used between successive calls to ! write_cputime. Initially it is set to be very large. integer :: nmax=2000000000 - + ! A structure containing pointers to the thermodynamic forcing fields + ! at the ocean surface. + type(forcing) :: fluxes ! A structure containing several relevant directory paths. type(directories) :: dirs @@ -104,7 +108,7 @@ program Shelf_main real :: time_step ! The time step [T ~> s] ! A pointer to a structure containing metrics and related information. - type(ocean_grid_type), pointer :: ocn_grid + type(ocean_grid_type), pointer :: ocn_grid => NULL() type(dyn_horgrid_type), pointer :: dG => NULL() ! A dynamic version of the horizontal grid type(hor_index_type), pointer :: HI => NULL() ! A hor_index_type for array extents @@ -114,7 +118,7 @@ program Shelf_main type(ocean_OBC_type), pointer :: OBC => NULL() ! A pointer to a structure containing dimensional unit scaling factors. - type(unit_scale_type), pointer :: US + type(unit_scale_type), pointer :: US => NULL() type(diag_ctrl), pointer :: & diag => NULL() ! A pointer to the diagnostic regulatory structure @@ -125,6 +129,7 @@ program Shelf_main ! files and +2 (bit 1) for time-stamped files. A ! restart file is saved at the end of a run segment ! unless Restart_control is negative. + logical :: save_both_restarts ! To save both time-stamped and non-time-stamped restart files real :: Time_unit ! The time unit for the following input fields [s]. type(time_type) :: restint ! The time between saves of the restart file. @@ -138,8 +143,9 @@ program Shelf_main integer :: yr, mon, day, hr, mins, sec ! Temp variables for writing the date. type(param_file_type) :: param_file ! The structure indicating the file(s) ! containing all run-time parameters. + real :: smb !A constant surface mass balance that can be specified in the param_file character(len=9) :: month - character(len=16) :: calendar = 'julian' + character(len=16) :: calendar = 'noleap' integer :: calendar_type=-1 integer :: unit, io_status, ierr @@ -184,6 +190,8 @@ program Shelf_main endif endif + ! Get the names of the I/O directories and initialization file. + ! Also calls the subroutine that opens run-time parameter files. call Get_MOM_Input(param_file, dirs) ! Read ocean_solo restart, which can override settings from the namelist. @@ -252,8 +260,11 @@ program Shelf_main ! Set up the ocean model domain and grid; the ice model grid is set in initialize_ice_shelf, ! but the grids have strong commonalities in this configuration, and the ocean grid is required ! to set up the diag mediator control structure. - call MOM_domains_init(ocn_grid%domain, param_file) + allocate(ocn_grid) + call MOM_domains_init(ocn_grid%domain, param_file) !, domain_name='MOM') + allocate(HI) call hor_index_init(ocn_grid%Domain, HI, param_file) + allocate(dG) call create_dyn_horgrid(dG, HI) call clone_MOM_domain(ocn_grid%Domain, dG%Domain) @@ -266,11 +277,16 @@ program Shelf_main ! Initialize the diag mediator. The ocean's vertical grid is not really used here, but at ! present the interface to diag_mediator_init assumes the presence of ocean-specific information. call verticalGridInit(param_file, GV, US) + allocate(diag) call diag_mediator_init(ocn_grid, GV, US, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory) call callTree_waypoint("returned from diag_mediator_init()") - call initialize_ice_shelf(param_file, ocn_grid, Time, ice_shelf_CSp, diag) + call set_axes_info(ocn_grid, GV, US, param_file, diag) + + call initialize_ice_shelf(param_file, ocn_grid, Time, ice_shelf_CSp, diag, fluxes_in=fluxes, solo_ice_sheet_in=.true.) + + call initialize_ice_SMB(fluxes%shelf_sfc_mass_flux, ocn_grid, US, param_file) ! This is the end of the code that is the counterpart of MOM_initialization. call callTree_waypoint("End of ice shelf initialization.") @@ -318,6 +334,9 @@ program Shelf_main "(bit 0) for a non-time-stamped file. A non-time-stamped "//& "restart file is saved at the end of the run segment "//& "for any non-negative value.", default=1) + call get_param(param_file, mod_name, "SAVE_BOTH_RESTARTS", save_both_restarts, & + "If true, save both a time-stamped and non-time-stamped "//& + "restart file after each restart interval.", default=.false.) call get_param(param_file, mod_name, "RESTINT", restint, & "The interval between saves of the restart file in units "//& "of TIMEUNIT. Use 0 (the default) to not save "//& @@ -378,7 +397,7 @@ program Shelf_main ! This call steps the model over a time time_step. Time1 = Master_Time ; Time = Master_Time - call solo_step_ice_shelf(ice_shelf_CSp, Time_step_shelf, ns_ice, Time) + call solo_step_ice_shelf(ice_shelf_CSp, Time_step_shelf, ns_ice, Time, fluxes_in=fluxes) ! Time = Time + Time_step_shelf ! This is here to enable fractional-second time steps. @@ -406,11 +425,26 @@ program Shelf_main ! See if it is time to write out a restart file - timestamped or not. if ((permit_incr_restart) .and. (Time + (Time_step_shelf/2) > restart_time)) then - if (BTEST(Restart_control,1)) then + if (BTEST(Restart_control,1) .or. save_both_restarts) then call ice_shelf_save_restart(ice_shelf_CSp, Time, dirs%restart_output_dir, .true.) endif - if (BTEST(Restart_control,0)) then + if (BTEST(Restart_control,0) .or. save_both_restarts) then call ice_shelf_save_restart(ice_shelf_CSp, Time, dirs%restart_output_dir) + + ! Write ice shelf solo restart file. + if (is_root_pe())then + call open_ASCII_file(unit, trim(dirs%restart_output_dir)//'shelf.res') + write(unit, '(i6,8x,a)') calendar_type, & + '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' + + call get_date(Start_time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Model start time: year, month, day, hour, minute, second' + call get_date(Time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Current model time: year, month, day, hour, minute, second' + call close_file(unit) + endif endif restart_time = restart_time + restint endif @@ -456,12 +490,11 @@ program Shelf_main endif call callTree_waypoint("End Shelf_main") - call diag_mediator_end(Time, diag, end_diag_manager=.true.) + call ice_shelf_end(ice_shelf_CSp) + call solo_ice_shelf_diag_mediator_end(Time, diag, end_diag_manager=.true.) if (cpu_steps > 0) call write_cputime(Time, ns-1, write_CPU_CSp, call_end=.true.) call cpu_clock_end(termClock) call io_infra_end ; call MOM_infra_end - call ice_shelf_end(ice_shelf_CSp) - end program Shelf_main diff --git a/config_src/drivers/solo_driver/MOM_driver.F90 b/config_src/drivers/solo_driver/MOM_driver.F90 index 84c2eec5b5..0e355f8638 100644 --- a/config_src/drivers/solo_driver/MOM_driver.F90 +++ b/config_src/drivers/solo_driver/MOM_driver.F90 @@ -49,6 +49,7 @@ program MOM6 use MOM_ice_shelf, only : shelf_calc_flux, add_shelf_forces, ice_shelf_save_restart use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces use MOM_ice_shelf, only : ice_shelf_query + use MOM_ice_shelf_initialize, only : initialize_ice_SMB use MOM_interpolate, only : time_interp_external_init use MOM_io, only : file_exists, open_ASCII_file, close_file use MOM_io, only : check_nml_error, io_infra_init, io_infra_end @@ -134,7 +135,7 @@ program MOM6 real :: dtdia ! The diabatic timestep [T ~> s] real :: t_elapsed_seg ! The elapsed time in this run segment [T ~> s] integer :: n, ns, n_max, nts, n_last_thermo - logical :: diabatic_first, single_step_call + logical :: diabatic_first, single_step_call, initialize_smb type(time_type) :: Time2, time_chg ! Temporary time variables integer :: Restart_control ! An integer that is bit-tested to determine whether @@ -302,6 +303,9 @@ program MOM6 call initialize_ice_shelf_forces(ice_shelf_CSp, grid, US, forces) call ice_shelf_query(ice_shelf_CSp, grid, data_override_shelf_fluxes=override_shelf_fluxes) if (override_shelf_fluxes) call data_override_init(Ocean_Domain_in=grid%domain%mpp_domain) + call get_param(param_file, mod_name, "INITIALIZE_ICE_SHEET_SMB", & + initialize_smb, "Read in a constant SMB for the ice sheet", default=.false.) + if (initialize_smb) call initialize_ice_SMB(fluxes%shelf_sfc_mass_flux, grid, US, param_file) endif diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 61290cb579..9d756adeff 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -47,6 +47,7 @@ module MOM_diag_mediator public safe_alloc_ptr, safe_alloc_alloc public enable_averaging, enable_averages, disable_averaging, query_averaging_enabled public diag_mediator_init, diag_mediator_end, set_diag_mediator_grid +public solo_ice_shelf_diag_mediator_end public diag_mediator_infrastructure_init public diag_mediator_close_registration, get_diag_time_end public diag_axis_init, ocean_register_diag, register_static_field @@ -3630,6 +3631,100 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) end subroutine diag_mediator_end +!> A copy of diag_mediator_end, but modified to be called from the ice shelf solo-driver by +!! removing deallocation of fields that were never allocated +subroutine solo_ice_shelf_diag_mediator_end(time, diag_CS, end_diag_manager) + type(time_type), intent(in) :: time !< The current model time + type(diag_ctrl), intent(inout) :: diag_CS !< Structure used to regulate diagnostic output + logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end() + + ! Local variables + type(diag_type), pointer :: diag, next_diag + integer :: i, dl + + if (diag_CS%available_diag_doc_unit > -1) then + close(diag_CS%available_diag_doc_unit) ; diag_CS%available_diag_doc_unit = -3 + endif + if (diag_CS%chksum_iounit > -1) then + close(diag_CS%chksum_iounit) ; diag_CS%chksum_iounit = -3 + endif + + do i=1, diag_cs%next_free_diag_id - 1 + if (associated(diag_cs%diags(i)%next)) then + next_diag => diag_cs%diags(i)%next + do while (associated(next_diag)) + diag => next_diag + next_diag => diag%next + deallocate(diag) + enddo + endif + enddo + + deallocate(diag_cs%diags) + + do i=1, diag_cs%num_diag_coords + call diag_remap_end(diag_cs%diag_remap_cs(i)) + enddo + + call diag_grid_storage_end(diag_cs%diag_grid_temp) + + ! axes_grp masks may point to diag_cs masks, so do these after mask dealloc + do i=1, diag_cs%num_diag_coords + call axes_grp_end(diag_cs%remap_axesZL(i)) + call axes_grp_end(diag_cs%remap_axesZi(i)) + call axes_grp_end(diag_cs%remap_axesTL(i)) + call axes_grp_end(diag_cs%remap_axesTi(i)) + call axes_grp_end(diag_cs%remap_axesBL(i)) + call axes_grp_end(diag_cs%remap_axesBi(i)) + call axes_grp_end(diag_cs%remap_axesCuL(i)) + call axes_grp_end(diag_cs%remap_axesCui(i)) + call axes_grp_end(diag_cs%remap_axesCvL(i)) + call axes_grp_end(diag_cs%remap_axesCvi(i)) + enddo + + if (diag_cs%num_diag_coords > 0) then + deallocate(diag_cs%remap_axesZL) + deallocate(diag_cs%remap_axesZi) + deallocate(diag_cs%remap_axesTL) + deallocate(diag_cs%remap_axesTi) + deallocate(diag_cs%remap_axesBL) + deallocate(diag_cs%remap_axesBi) + deallocate(diag_cs%remap_axesCuL) + deallocate(diag_cs%remap_axesCui) + deallocate(diag_cs%remap_axesCvL) + deallocate(diag_cs%remap_axesCvi) + endif + + do dl=2,MAX_DSAMP_LEV + if (allocated(diag_cs%dsamp(dl)%remap_axesTL)) & + deallocate(diag_cs%dsamp(dl)%remap_axesTL) + if (allocated(diag_cs%dsamp(dl)%remap_axesTi)) & + deallocate(diag_cs%dsamp(dl)%remap_axesTi) + if (allocated(diag_cs%dsamp(dl)%remap_axesBL)) & + deallocate(diag_cs%dsamp(dl)%remap_axesBL) + if (allocated(diag_cs%dsamp(dl)%remap_axesBi)) & + deallocate(diag_cs%dsamp(dl)%remap_axesBi) + if (allocated(diag_cs%dsamp(dl)%remap_axesCuL)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCuL) + if (allocated(diag_cs%dsamp(dl)%remap_axesCui)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCui) + if (allocated(diag_cs%dsamp(dl)%remap_axesCvL)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCvL) + if (allocated(diag_cs%dsamp(dl)%remap_axesCvi)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCvi) + enddo + + +#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) + deallocate(diag_cs%h_old) +#endif + + if (present(end_diag_manager)) then + if (end_diag_manager) call MOM_diag_manager_end(time) + endif + +end subroutine solo_ice_shelf_diag_mediator_end + !> Convert the first n elements (up to 3) of an integer array to an underscore delimited string. function i2s(a,n_in) ! "Convert the first n elements of an integer array to a string." diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index f5a85da95a..7176e3ccdf 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1873,7 +1873,8 @@ subroutine initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in) tau_mag=.true.) else call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.") - call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., press=.true., tau_mag=.true.) + call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., & + press=.true., shelf_sfc_accumulation = CS%active_shelf_dynamics, tau_mag=.true.) endif if (CS%rotate_index) then allocate(fluxes) @@ -2178,13 +2179,14 @@ subroutine ice_shelf_end(CS) end subroutine ice_shelf_end !> This routine is for stepping a stand-alone ice shelf model without an ocean. -subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in) +subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in, fluxes_in) type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure type(time_type), intent(in) :: time_interval !< The time interval for this update [s]. integer, intent(inout) :: nsteps !< The running number of ice shelf steps. type(time_type), intent(inout) :: Time !< The current model time real, optional, intent(in) :: min_time_step_in !< The minimum permitted time step [T ~> s]. - + type(forcing), optional, target, intent(inout) :: fluxes_in !< A structure containing pointers to any + !! possible thermodynamic or mass-flux forcing fields. type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the ocean's grid structure type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing ! various unit conversion factors @@ -2192,6 +2194,7 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in !! the ice-shelf state real :: remaining_time ! The remaining time in this call [T ~> s] real :: time_step ! The internal time step during this call [T ~> s] + real :: full_time_step ! The external time step (sum of internal time steps) during this call [T ~> s] real :: min_time_step ! The minimal required timestep that would indicate a fatal problem [T ~> s] character(len=240) :: mesg logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities. @@ -2205,6 +2208,7 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in is = G%isc ; iec = G%iec ; js = G%jsc ; jec = G%jec remaining_time = US%s_to_T*time_type_to_real(time_interval) + full_time_step = remaining_time if (present (min_time_step_in)) then min_time_step = min_time_step_in @@ -2228,6 +2232,8 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in call MOM_mesg("solo_step_ice_shelf: "//mesg, 5) endif + call change_thickness_using_precip(CS, ISS, G, US, fluxes_in, time_step, Time) + remaining_time = remaining_time - time_step ! If the last mini-timestep is a day or less, we cannot expect velocities to change by much. @@ -2237,13 +2243,13 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in call update_ice_shelf(CS%dCS, ISS, G, US, time_step, Time, must_update_vel=update_ice_vel) - call enable_averages(time_step, Time, CS%diag) + enddo + + call enable_averages(full_time_step, Time, CS%diag) if (CS%id_area_shelf_h > 0) call post_data(CS%id_area_shelf_h, ISS%area_shelf_h, CS%diag) if (CS%id_h_shelf > 0) call post_data(CS%id_h_shelf, ISS%h_shelf, CS%diag) if (CS%id_h_mask > 0) call post_data(CS%id_h_mask, ISS%hmask, CS%diag) - call disable_averaging(CS%diag) - - enddo + call disable_averaging(CS%diag) end subroutine solo_step_ice_shelf diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 81a4c7e21b..25f6b9f73f 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -768,7 +768,7 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled ! call ice_shelf_temp(CS, ISS, G, US, time_step, ISS%water_flux, Time) - if (update_ice_vel) then + if (CS%elapsed_velocity_time >= CS%velocity_update_time_step) then call enable_averages(CS%elapsed_velocity_time, Time, CS%diag) if (CS%id_col_thick > 0) call post_data(CS%id_col_thick, CS%OD_av, CS%diag) if (CS%id_u_shelf > 0) call post_data(CS%id_u_shelf, CS%u_shelf, CS%diag) diff --git a/src/ice_shelf/MOM_ice_shelf_initialize.F90 b/src/ice_shelf/MOM_ice_shelf_initialize.F90 index 20a48730f3..1e2076f889 100644 --- a/src/ice_shelf/MOM_ice_shelf_initialize.F90 +++ b/src/ice_shelf/MOM_ice_shelf_initialize.F90 @@ -22,6 +22,7 @@ module MOM_ice_shelf_initialize public initialize_ice_shelf_boundary_from_file public initialize_ice_C_basal_friction public initialize_ice_AGlen +public initialize_ice_SMB ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units @@ -657,5 +658,51 @@ subroutine initialize_ice_AGlen(AGlen, G, US, PF) call MOM_read_data(filename,trim(varname), AGlen, G%Domain) endif -end subroutine +end subroutine initialize_ice_AGlen + +!> Initialize ice surface mass balance field that is held constant over time +subroutine initialize_ice_SMB(SMB, G, US, PF) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZDI_(G),SZDJ_(G)), & + intent(inout) :: SMB !< Ice surface mass balance parameter, often in [kg m-2 s-1] + type(unit_scale_type), intent(in) :: US !< A structure containing unit conversion factors + type(param_file_type), intent(in) :: PF !< A structure to parse for run-time parameters + + real :: SMB_val ! Constant ice surface mass balance parameter, often in [kg m-2 s-1] + character(len=40) :: mdl = "initialize_ice_SMB" ! This subroutine's name. + character(len=200) :: config + character(len=200) :: varname + character(len=200) :: inputdir, filename, SMB_file + + call get_param(PF, mdl, "ICE_SMB_CONFIG", config, & + "This specifies how the initial ice surface mass balance parameter is specified. "//& + "Valid values are: CONSTANT and FILE.", & + default="CONSTANT") + + if (trim(config)=="CONSTANT") then + call get_param(PF, mdl, "SMB", SMB_val, & + "Surface mass balance.", units="kg m-2 s-1", default=0.0) + + SMB(:,:) = SMB_val + + elseif (trim(config)=="FILE") then + call MOM_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading SMB parameter") + call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + + call get_param(PF, mdl, "ICE_SMB_FILE", SMB_file, & + "The file from which the ice surface mass balance is read.", & + default="ice_SMB.nc") + filename = trim(inputdir)//trim(SMB_file) + call log_param(PF, mdl, "INPUTDIR/ICE_SMB_FILE", filename) + call get_param(PF, mdl, "ICE_SMB_VARNAME", varname, & + "The variable to use as surface mass balance.", & + default="SMB") + + if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, & + " initialize_ice_SMV_from_file: Unable to open "//trim(filename)) + call MOM_read_data(filename,trim(varname), SMB, G%Domain) + + endif +end subroutine initialize_ice_SMB end module MOM_ice_shelf_initialize From 11a2b257f9a27b0b71c52daf8507176972f49942 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Fri, 13 Oct 2023 14:28:20 -0400 Subject: [PATCH 2/2] Added safety checks to diag_mediator_end() so that it works with the ice shelf solo-driver, which now calls it instead of (now removed) solo_ice_shelf_diag_mediator_end() routine. Removed the runtime parameter SAVE_BOTH_RESTARTS from the ice shelf solo-driver, which is no longer needed. --- .../ice_solo_driver/ice_shelf_driver.F90 | 41 ++--- src/framework/MOM_diag_mediator.F90 | 159 ++++-------------- 2 files changed, 54 insertions(+), 146 deletions(-) diff --git a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 index b864903ddc..f91595bd51 100644 --- a/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 +++ b/config_src/drivers/ice_solo_driver/ice_shelf_driver.F90 @@ -25,7 +25,7 @@ program Shelf_main use MOM_cpu_clock, only : CLOCK_COMPONENT use MOM_debugging, only : MOM_debugging_init use MOM_diag_mediator, only : diag_mediator_init, diag_mediator_infrastructure_init, set_axes_info - use MOM_diag_mediator, only : solo_ice_shelf_diag_mediator_end, diag_ctrl, diag_mediator_close_registration + use MOM_diag_mediator, only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration use MOM_domains, only : MOM_infra_init, MOM_infra_end use MOM_domains, only : MOM_domains_init, clone_MOM_domain, pass_var use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid @@ -129,7 +129,6 @@ program Shelf_main ! files and +2 (bit 1) for time-stamped files. A ! restart file is saved at the end of a run segment ! unless Restart_control is negative. - logical :: save_both_restarts ! To save both time-stamped and non-time-stamped restart files real :: Time_unit ! The time unit for the following input fields [s]. type(time_type) :: restint ! The time between saves of the restart file. @@ -334,9 +333,6 @@ program Shelf_main "(bit 0) for a non-time-stamped file. A non-time-stamped "//& "restart file is saved at the end of the run segment "//& "for any non-negative value.", default=1) - call get_param(param_file, mod_name, "SAVE_BOTH_RESTARTS", save_both_restarts, & - "If true, save both a time-stamped and non-time-stamped "//& - "restart file after each restart interval.", default=.false.) call get_param(param_file, mod_name, "RESTINT", restint, & "The interval between saves of the restart file in units "//& "of TIMEUNIT. Use 0 (the default) to not save "//& @@ -425,26 +421,25 @@ program Shelf_main ! See if it is time to write out a restart file - timestamped or not. if ((permit_incr_restart) .and. (Time + (Time_step_shelf/2) > restart_time)) then - if (BTEST(Restart_control,1) .or. save_both_restarts) then + if (BTEST(Restart_control,1)) then call ice_shelf_save_restart(ice_shelf_CSp, Time, dirs%restart_output_dir, .true.) endif - if (BTEST(Restart_control,0) .or. save_both_restarts) then + if (BTEST(Restart_control,0)) then call ice_shelf_save_restart(ice_shelf_CSp, Time, dirs%restart_output_dir) - - ! Write ice shelf solo restart file. - if (is_root_pe())then - call open_ASCII_file(unit, trim(dirs%restart_output_dir)//'shelf.res') - write(unit, '(i6,8x,a)') calendar_type, & - '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' - - call get_date(Start_time, yr, mon, day, hr, mins, sec) - write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & - 'Model start time: year, month, day, hour, minute, second' - call get_date(Time, yr, mon, day, hr, mins, sec) - write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & - 'Current model time: year, month, day, hour, minute, second' - call close_file(unit) - endif + endif + ! Write ice shelf solo restart file. + if (is_root_pe())then + call open_ASCII_file(unit, trim(dirs%restart_output_dir)//'shelf.res') + write(unit, '(i6,8x,a)') calendar_type, & + '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)' + + call get_date(Start_time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Model start time: year, month, day, hour, minute, second' + call get_date(Time, yr, mon, day, hr, mins, sec) + write(unit, '(6i6,8x,a)') yr, mon, day, hr, mins, sec, & + 'Current model time: year, month, day, hour, minute, second' + call close_file(unit) endif restart_time = restart_time + restint endif @@ -491,7 +486,7 @@ program Shelf_main call callTree_waypoint("End Shelf_main") call ice_shelf_end(ice_shelf_CSp) - call solo_ice_shelf_diag_mediator_end(Time, diag, end_diag_manager=.true.) + call diag_mediator_end(Time, diag, end_diag_manager=.true.) if (cpu_steps > 0) call write_cputime(Time, ns-1, write_CPU_CSp, call_end=.true.) call cpu_clock_end(termClock) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 9d756adeff..2c71a93e42 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -47,7 +47,6 @@ module MOM_diag_mediator public safe_alloc_ptr, safe_alloc_alloc public enable_averaging, enable_averages, disable_averaging, query_averaging_enabled public diag_mediator_init, diag_mediator_end, set_diag_mediator_grid -public solo_ice_shelf_diag_mediator_end public diag_mediator_infrastructure_init public diag_mediator_close_registration, get_diag_time_end public diag_axis_init, ocean_register_diag, register_static_field @@ -3540,37 +3539,45 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) enddo call diag_grid_storage_end(diag_cs%diag_grid_temp) - deallocate(diag_cs%mask3dTL) - deallocate(diag_cs%mask3dBL) - deallocate(diag_cs%mask3dCuL) - deallocate(diag_cs%mask3dCvL) - deallocate(diag_cs%mask3dTi) - deallocate(diag_cs%mask3dBi) - deallocate(diag_cs%mask3dCui) - deallocate(diag_cs%mask3dCvi) + if (associated(diag_cs%mask3dTL)) deallocate(diag_cs%mask3dTL) + if (associated(diag_cs%mask3dBL)) deallocate(diag_cs%mask3dBL) + if (associated(diag_cs%mask3dCuL)) deallocate(diag_cs%mask3dCuL) + if (associated(diag_cs%mask3dCvL)) deallocate(diag_cs%mask3dCvL) + if (associated(diag_cs%mask3dTi)) deallocate(diag_cs%mask3dTi) + if (associated(diag_cs%mask3dBi)) deallocate(diag_cs%mask3dBi) + if (associated(diag_cs%mask3dCui)) deallocate(diag_cs%mask3dCui) + if (associated(diag_cs%mask3dCvi)) deallocate(diag_cs%mask3dCvi) do dl=2,MAX_DSAMP_LEV - deallocate(diag_cs%dsamp(dl)%mask2dT) - deallocate(diag_cs%dsamp(dl)%mask2dBu) - deallocate(diag_cs%dsamp(dl)%mask2dCu) - deallocate(diag_cs%dsamp(dl)%mask2dCv) - deallocate(diag_cs%dsamp(dl)%mask3dTL) - deallocate(diag_cs%dsamp(dl)%mask3dBL) - deallocate(diag_cs%dsamp(dl)%mask3dCuL) - deallocate(diag_cs%dsamp(dl)%mask3dCvL) - deallocate(diag_cs%dsamp(dl)%mask3dTi) - deallocate(diag_cs%dsamp(dl)%mask3dBi) - deallocate(diag_cs%dsamp(dl)%mask3dCui) - deallocate(diag_cs%dsamp(dl)%mask3dCvi) + if (associated(diag_cs%dsamp(dl)%mask2dT)) deallocate(diag_cs%dsamp(dl)%mask2dT) + if (associated(diag_cs%dsamp(dl)%mask2dBu)) deallocate(diag_cs%dsamp(dl)%mask2dBu) + if (associated(diag_cs%dsamp(dl)%mask2dCu)) deallocate(diag_cs%dsamp(dl)%mask2dCu) + if (associated(diag_cs%dsamp(dl)%mask2dCv)) deallocate(diag_cs%dsamp(dl)%mask2dCv) + if (associated(diag_cs%dsamp(dl)%mask3dTL)) deallocate(diag_cs%dsamp(dl)%mask3dTL) + if (associated(diag_cs%dsamp(dl)%mask3dBL)) deallocate(diag_cs%dsamp(dl)%mask3dBL) + if (associated(diag_cs%dsamp(dl)%mask3dCuL)) deallocate(diag_cs%dsamp(dl)%mask3dCuL) + if (associated(diag_cs%dsamp(dl)%mask3dCvL)) deallocate(diag_cs%dsamp(dl)%mask3dCvL) + if (associated(diag_cs%dsamp(dl)%mask3dTi)) deallocate(diag_cs%dsamp(dl)%mask3dTi) + if (associated(diag_cs%dsamp(dl)%mask3dBi)) deallocate(diag_cs%dsamp(dl)%mask3dBi) + if (associated(diag_cs%dsamp(dl)%mask3dCui)) deallocate(diag_cs%dsamp(dl)%mask3dCui) + if (associated(diag_cs%dsamp(dl)%mask3dCvi)) deallocate(diag_cs%dsamp(dl)%mask3dCvi) do i=1,diag_cs%num_diag_coords - deallocate(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d) - deallocate(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d) - deallocate(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d) - deallocate(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d) - deallocate(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d) - deallocate(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d) - deallocate(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d) - deallocate(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesTL(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCuL(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCvL(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesBL(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesTi(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCui(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesCvi(i)%dsamp(dl)%mask3d) + if (associated(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d)) & + deallocate(diag_cs%dsamp(dl)%remap_axesBi(i)%dsamp(dl)%mask3d) enddo enddo @@ -3631,100 +3638,6 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) end subroutine diag_mediator_end -!> A copy of diag_mediator_end, but modified to be called from the ice shelf solo-driver by -!! removing deallocation of fields that were never allocated -subroutine solo_ice_shelf_diag_mediator_end(time, diag_CS, end_diag_manager) - type(time_type), intent(in) :: time !< The current model time - type(diag_ctrl), intent(inout) :: diag_CS !< Structure used to regulate diagnostic output - logical, optional, intent(in) :: end_diag_manager !< If true, call diag_manager_end() - - ! Local variables - type(diag_type), pointer :: diag, next_diag - integer :: i, dl - - if (diag_CS%available_diag_doc_unit > -1) then - close(diag_CS%available_diag_doc_unit) ; diag_CS%available_diag_doc_unit = -3 - endif - if (diag_CS%chksum_iounit > -1) then - close(diag_CS%chksum_iounit) ; diag_CS%chksum_iounit = -3 - endif - - do i=1, diag_cs%next_free_diag_id - 1 - if (associated(diag_cs%diags(i)%next)) then - next_diag => diag_cs%diags(i)%next - do while (associated(next_diag)) - diag => next_diag - next_diag => diag%next - deallocate(diag) - enddo - endif - enddo - - deallocate(diag_cs%diags) - - do i=1, diag_cs%num_diag_coords - call diag_remap_end(diag_cs%diag_remap_cs(i)) - enddo - - call diag_grid_storage_end(diag_cs%diag_grid_temp) - - ! axes_grp masks may point to diag_cs masks, so do these after mask dealloc - do i=1, diag_cs%num_diag_coords - call axes_grp_end(diag_cs%remap_axesZL(i)) - call axes_grp_end(diag_cs%remap_axesZi(i)) - call axes_grp_end(diag_cs%remap_axesTL(i)) - call axes_grp_end(diag_cs%remap_axesTi(i)) - call axes_grp_end(diag_cs%remap_axesBL(i)) - call axes_grp_end(diag_cs%remap_axesBi(i)) - call axes_grp_end(diag_cs%remap_axesCuL(i)) - call axes_grp_end(diag_cs%remap_axesCui(i)) - call axes_grp_end(diag_cs%remap_axesCvL(i)) - call axes_grp_end(diag_cs%remap_axesCvi(i)) - enddo - - if (diag_cs%num_diag_coords > 0) then - deallocate(diag_cs%remap_axesZL) - deallocate(diag_cs%remap_axesZi) - deallocate(diag_cs%remap_axesTL) - deallocate(diag_cs%remap_axesTi) - deallocate(diag_cs%remap_axesBL) - deallocate(diag_cs%remap_axesBi) - deallocate(diag_cs%remap_axesCuL) - deallocate(diag_cs%remap_axesCui) - deallocate(diag_cs%remap_axesCvL) - deallocate(diag_cs%remap_axesCvi) - endif - - do dl=2,MAX_DSAMP_LEV - if (allocated(diag_cs%dsamp(dl)%remap_axesTL)) & - deallocate(diag_cs%dsamp(dl)%remap_axesTL) - if (allocated(diag_cs%dsamp(dl)%remap_axesTi)) & - deallocate(diag_cs%dsamp(dl)%remap_axesTi) - if (allocated(diag_cs%dsamp(dl)%remap_axesBL)) & - deallocate(diag_cs%dsamp(dl)%remap_axesBL) - if (allocated(diag_cs%dsamp(dl)%remap_axesBi)) & - deallocate(diag_cs%dsamp(dl)%remap_axesBi) - if (allocated(diag_cs%dsamp(dl)%remap_axesCuL)) & - deallocate(diag_cs%dsamp(dl)%remap_axesCuL) - if (allocated(diag_cs%dsamp(dl)%remap_axesCui)) & - deallocate(diag_cs%dsamp(dl)%remap_axesCui) - if (allocated(diag_cs%dsamp(dl)%remap_axesCvL)) & - deallocate(diag_cs%dsamp(dl)%remap_axesCvL) - if (allocated(diag_cs%dsamp(dl)%remap_axesCvi)) & - deallocate(diag_cs%dsamp(dl)%remap_axesCvi) - enddo - - -#if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__) - deallocate(diag_cs%h_old) -#endif - - if (present(end_diag_manager)) then - if (end_diag_manager) call MOM_diag_manager_end(time) - endif - -end subroutine solo_ice_shelf_diag_mediator_end - !> Convert the first n elements (up to 3) of an integer array to an underscore delimited string. function i2s(a,n_in) ! "Convert the first n elements of an integer array to a string."