From 2c93933a971e9054058ff3512c79622ff3a9cd92 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 21 Dec 2020 12:43:08 -0500 Subject: [PATCH 01/11] Revert changes in the MOM drivers for ice shelf initialization - Call ice shelf initialization from within initialize_MOM. --- config_src/solo_driver/MOM_driver.F90 | 28 +++++++++++++-------------- src/core/MOM.F90 | 11 ++++++++--- src/ice_shelf/MOM_ice_shelf.F90 | 15 ++++++++------ 3 files changed, 31 insertions(+), 23 deletions(-) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 9981f291b1..3e11299cea 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -307,20 +307,8 @@ program MOM_main Time = Start_time endif - ! Read paths and filenames from namelist and store in "dirs". - ! Also open the parsed input parameter file(s) and setup param_file. - call get_MOM_input(param_file, dirs) - - call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, & - "If true, enables the ice shelf model.", default=.false.) - if (use_ice_shelf) then - ! These arrays are not initialized in most solo cases, but are needed - ! when using an ice shelf - call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, & - diag_IS, forces, fluxes, sfc_state) - endif - call close_param_file(param_file) - + ! Call initialize MOM with an optional Ice Shelf CS which, if present triggers + ! initialization of ice shelf parameters and arrays. if (sum(date) >= 0) then call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, & segment_start_time, offline_tracer_mode=offline_tracer_mode, & @@ -331,6 +319,18 @@ program MOM_main tracer_flow_CSp=tracer_flow_CSp,ice_shelf_CSp=ice_shelf_CSp) endif + call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, & + "If true, enables the ice shelf model.", default=.false.) + + if (use_ice_shelf) then + ! These arrays are not initialized in most solo cases, but are needed + ! when using an ice shelf + ice_shelf_CSp => NULL() ! Reset the pointer and make another call to reinitialize + ! the ice shelf and associated forcing types + call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, & + diag, forces, fluxes, sfc_state) + endif + call get_MOM_state_elements(MOM_CSp, G=grid, GV=GV, US=US, C_p_scaled=fluxes%C_p) Master_Time = Time diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 3f4507f546..fa7c3ca565 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -141,7 +141,7 @@ module MOM use MOM_offline_main, only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline -use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query +use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf implicit none ; private @@ -2037,7 +2037,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & use_ice_shelf=.false. if (present(ice_shelf_CSp)) then - if (associated(ice_shelf_CSp)) use_ice_shelf=.true. + call get_param(param_file, "MOM", "ICE_SHELF", use_ice_shelf, & + "If true, enables the ice shelf model.", default=.false.) endif CS%ensemble_ocean=.false. @@ -2381,6 +2382,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif if (use_ice_shelf) then + ! 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) allocate(frac_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed)) frac_shelf_in(:,:) = 0.0 allocate(CS%frac_shelf_h(isd:ied, jsd:jed)) @@ -2431,10 +2436,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & deallocate(frac_shelf_in) else if (use_ice_shelf) then + call initialize_ice_shelf(param_file, G, Time, ice_shelf_CSp, diag_ptr) allocate(CS%frac_shelf_h(isd:ied, jsd:jed)) CS%frac_shelf_h(:,:) = 0.0 call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h) - call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in, & diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 92545d6dc4..88ad6c7123 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -9,6 +9,7 @@ module MOM_ice_shelf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_ROUTINE use MOM_coms, only : num_PEs +use MOM_diag_mediator, only : MOM_diag_ctrl=>diag_ctrl use MOM_IS_diag_mediator, only : post_data, register_diag_field=>register_MOM_IS_diag_field, safe_alloc_ptr use MOM_IS_diag_mediator, only : set_axes_info use MOM_IS_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid, diag_ctrl, time_type @@ -1156,7 +1157,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure type(time_type), intent(inout) :: Time !< The clock that that will indicate the model time type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure - type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the diagnostic output. + type(MOM_diag_ctrl), pointer :: diag !< This is a pointer to the MOM diag CS + !! which will be discarded + type(mech_forcing), optional, target, intent(inout) :: forces_in !< A structure with the driving mechanical forces type(forcing), optional, target, intent(inout) :: fluxes_in !< A structure containing pointers to any !! possible thermodynamic or mass-flux forcing fields. @@ -1284,11 +1287,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif G=>CS%Grid - allocate(diag) - call diag_mediator_init(G, param_file,diag,component='MOM_IceShelf') + allocate(CS%diag) + call diag_mediator_init(G, param_file,CS%diag,component='MOM_IceShelf') ! This call sets up the diagnostic axes. These are needed, ! e.g. to generate the target grids below. - call set_axes_info(G, param_file, diag) + call set_axes_info(G, param_file, CS%diag) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -1302,7 +1305,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, ! Convenience pointers OG => CS%ocn_grid US => CS%US - CS%diag=>diag + !CS%diag=>diag ! Are we being called from the solo ice-sheet driver? When called by the ocean ! model solo_ice_sheet_in is not preset. @@ -1777,7 +1780,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif if (shelf_mass_is_dynamic) & - call initialize_ice_shelf_dyn(param_file, Time, ISS, CS%dCS, G, US, diag, new_sim, solo_ice_sheet_in) + call initialize_ice_shelf_dyn(param_file, Time, ISS, CS%dCS, G, US, CS%diag, new_sim, solo_ice_sheet_in) call fix_restart_unit_scaling(US) From f3dc73f7b78770c2e857e31a1441ebb5ca88ba2e Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 21 Dec 2020 13:10:19 -0500 Subject: [PATCH 02/11] Revert forcing,fluxes and sfc_state from pointers. - This removes suppport for rotation. --- config_src/solo_driver/MOM_driver.F90 | 9 +++-- src/ice_shelf/MOM_ice_shelf.F90 | 51 ++++++++++++++------------- 2 files changed, 31 insertions(+), 29 deletions(-) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 3e11299cea..d61d6c986a 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -80,13 +80,12 @@ program MOM_main #include ! A structure with the driving mechanical surface forces - type(mech_forcing), pointer :: forces => NULL() + type(mech_forcing) :: forces ! A structure containing pointers to the thermodynamic forcing fields ! at the ocean surface. - type(forcing), pointer :: fluxes => NULL() - + type(forcing) :: fluxes ! A structure containing pointers to the ocean surface state fields. - type(surface), pointer :: sfc_state => NULL() + type(surface) :: sfc_state ! A pointer to a structure containing metrics and related information. type(ocean_grid_type), pointer :: grid => NULL() @@ -221,7 +220,7 @@ program MOM_main call MOM_infra_init() ; call io_infra_init() - allocate(forces,fluxes,sfc_state) + !allocate(forces,fluxes,sfc_state) ! Initialize the ensemble manager. If there are no settings for ensemble_size ! in input.nml(ensemble.nml), these should not do anything. In coupled diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 88ad6c7123..3638652bfc 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -807,18 +807,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_in, do_shelf_area, external_call) +subroutine add_shelf_forces(Ocn_grid, US, CS, forces, 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),target, intent(inout) :: forces_in !< A structure with the + type(mech_forcing), intent(inout) :: forces !< 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 !! to be rotated. type(ocean_grid_type), pointer :: G => NULL() !< A pointer to the ocean grid metric. - type(mech_forcing), pointer :: forces => NULL() !< A structure with the driving mechanical forces +! type(mech_forcing), target :: 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. @@ -830,20 +830,25 @@ subroutine add_shelf_forces(Ocn_grid, US, CS, forces_in, do_shelf_area, external if (present(external_call)) rotate=external_call - if (CS%rotate_index .and. rotate) then - 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.") - - 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.") + 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.") - forces=>forces_in + if (CS%rotate_index .and. rotate) then + call MOM_error(FATAL,"add_shelf_forces: Rotation not implemented for ice shelves.") + ! 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.") + + ! 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 endif G=>CS%Grid @@ -911,10 +916,10 @@ subroutine add_shelf_forces(Ocn_grid, US, CS, forces_in, do_shelf_area, external 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) + ! ! TODO: deallocate mech forcing? + ! endif end subroutine add_shelf_forces @@ -923,8 +928,7 @@ subroutine add_shelf_pressure(Ocn_grid, US, CS, fluxes) 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), intent(in) :: CS !< This module's control structure. - type(forcing), pointer :: fluxes !< A structure of surface fluxes that may be updated. - + type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be updated. type(ocean_grid_type), pointer :: G => NULL() ! A pointer to ocean's grid structure. real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa]. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed @@ -957,7 +961,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(ice_shelf_CS), pointer :: CS !< This module's control structure. type(surface), intent(inout) :: sfc_state !< Surface ocean state - type(forcing), pointer :: fluxes !< A structure of surface fluxes that may be used/updated. + type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated. ! local variables real :: frac_shelf !< The fractional area covered by the ice shelf [nondim]. @@ -1305,7 +1309,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, ! Convenience pointers OG => CS%ocn_grid US => CS%US - !CS%diag=>diag ! Are we being called from the solo ice-sheet driver? When called by the ocean ! model solo_ice_sheet_in is not preset. From 8087f9300e4728bd3cb90f883ce69aa870bfbf1f Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 21 Dec 2020 13:56:16 -0500 Subject: [PATCH 03/11] Return early from ice shelf initialization to avoid registering diagnostics and restarts during the call from initialize_MOM --- src/ice_shelf/MOM_ice_shelf.F90 | 244 +++++++++++++++++--------------- 1 file changed, 129 insertions(+), 115 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 3638652bfc..15326b85ec 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1212,7 +1212,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, type(forcing), pointer :: fluxes => NULL() type(surface), pointer :: sfc_state => NULL() type(vardesc) :: u_desc, v_desc - + logical :: complete_initialization ! A flag which is set to true if forces are present + ! This exists for legacy reasons and is a means to avoid some + ! parts of the initilization procedure since the ice shelf + ! is being initialized twice from initialize MOM and from the + ! various driver routines. if (associated(CS)) then call MOM_error(FATAL, "MOM_ice_shelf.F90, initialize_ice_shelf: "// & "called with an associated control structure.") @@ -1220,6 +1224,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif allocate(CS) + complete_initialization=.false. + if (present(forces_in)) complete_initialization = .true. ! Go through all of the infrastructure initialization calls, since this is ! being treated as an independent component that just happens to use the ! MOM's grid and infrastructure. @@ -1245,57 +1251,59 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, !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') -! allocate(CS%Grid_in%HI) call hor_index_init(CS%Grid_in%Domain, CS%Grid_in%HI, param_file, & local_indexing=.not.global_indexing) call MOM_grid_init(CS%Grid_in, param_file, CS%US, CS%Grid_in%HI) - 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 - !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, G%D either analytically or from file - call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file) - call rescale_dyn_horgrid_bathymetry(dG, CS%US%Z_to_m) - call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) - endif + ! 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, G%D either analytically or from file + call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file) + call rescale_dyn_horgrid_bathymetry(dG, CS%US%Z_to_m) + call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) +! endif G=>CS%Grid - allocate(CS%diag) - call diag_mediator_init(G, param_file,CS%diag,component='MOM_IceShelf') - ! This call sets up the diagnostic axes. These are needed, - ! e.g. to generate the target grids below. - call set_axes_info(G, param_file, CS%diag) + if (complete_initialization) then + allocate(CS%diag) + call diag_mediator_init(G, param_file,CS%diag,component='MOM_IceShelf') + ! This call sets up the diagnostic axes. These are needed, + ! e.g. to generate the target grids below. + call set_axes_info(G, param_file, CS%diag) + endif is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -1490,31 +1498,33 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif - call safe_alloc_ptr(CS%utide,isd,ied,jsd,jed) ; CS%utide(:,:) = 0.0 - - if (read_TIDEAMP) then - call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, & + if (complete_initialization) then + call safe_alloc_ptr(CS%utide,isd,ied,jsd,jed) ; CS%utide(:,:) = 0.0 + if (read_TIDEAMP) then + call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, & "The path to the file containing the spatially varying "//& "tidal amplitudes.", & default="tideamp.nc") - call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - TideAmp_file = trim(inputdir) // trim(TideAmp_file) - 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) - deallocate(tmp2d) + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + TideAmp_file = trim(inputdir) // trim(TideAmp_file) + 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) + deallocate(tmp2d) + else + call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) + endif else - call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) - endif - else - call get_param(param_file, mdl, "UTIDE", utide, & + call get_param(param_file, mdl, "UTIDE", utide, & "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & units="m s-1", default=0.0 , scale=US%m_s_to_L_T) - CS%utide(:,:) = utide + CS%utide(:,:) = utide + endif endif + call EOS_init(param_file, CS%eqn_of_state) !! new parameters that need to be in MOM_input @@ -1606,51 +1616,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif endif - ! Set up the restarts. - call restart_init(param_file, CS%restart_CSp, "Shelf.res") - call register_restart_field(ISS%mass_shelf, "shelf_mass", .true., CS%restart_CSp, & - "Ice shelf mass", "kg m-2") - call register_restart_field(ISS%area_shelf_h, "shelf_area", .true., CS%restart_CSp, & - "Ice shelf area in cell", "m2") - call register_restart_field(ISS%h_shelf, "h_shelf", .true., CS%restart_CSp, & - "ice sheet/shelf thickness", "m") - if (PRESENT(sfc_state_in)) then - if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then - u_desc = var_desc("taux_shelf", "Pa", "the zonal stress on the ocean under ice shelves", & - hor_grid='Cu',z_grid='1') - v_desc = var_desc("tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", & - hor_grid='Cv',z_grid='1') - call register_restart_pair(sfc_state%taux_shelf, sfc_state%tauy_shelf, u_desc,v_desc, & - .false., CS%restart_CSp) - endif - endif - - call register_restart_field(ISS%h_shelf, "_shelf", .true., CS%restart_CSp, & - "ice sheet/shelf thickness", "m") - call register_restart_field(US%m_to_Z_restart, "m_to_Z", .false., CS%restart_CSp, & - "Height unit conversion factor", "Z meter-1") - call register_restart_field(US%m_to_L_restart, "m_to_L", .false., CS%restart_CSp, & - "Length unit conversion factor", "L meter-1") - call register_restart_field(US%kg_m3_to_R_restart, "kg_m3_to_R", .false., CS%restart_CSp, & - "Density unit conversion factor", "R m3 kg-1") - if (CS%active_shelf_dynamics) then - call register_restart_field(ISS%hmask, "h_mask", .true., CS%restart_CSp, & - "ice sheet/shelf thickness mask" ,"none") - endif - - if (CS%active_shelf_dynamics) then - ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics - call register_ice_shelf_dyn_restarts(CS%Grid_in, param_file, CS%dCS, CS%restart_CSp) - endif - - !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file - !if (.not. CS%solo_ice_sheet) then - ! call register_restart_field(fluxes%ustar_shelf, "ustar_shelf", .false., CS%restart_CSp, & - ! "Friction velocity under ice shelves", "m s-1") - !endif - - CS%restart_output_dir = dirs%restart_output_dir - new_sim = .false. if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. @@ -1740,12 +1705,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif ! .not. new_sim -! do j=G%jsc,G%jec ; do i=G%isc,G%iec -! ISS%area_shelf_h(i,j) = ISS%area_shelf_h(i,j)*G%mask2dT(i,j) -! enddo; enddo - - CS%Time = Time - id_clock_shelf = cpu_clock_id('Ice shelf', grain=CLOCK_COMPONENT) id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=CLOCK_ROUTINE) @@ -1757,7 +1716,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call pass_var(G%bathyT, G%domain) call cpu_clock_end(id_clock_pass) - do j=jsd,jed ; do i=isd,ied if (ISS%area_shelf_h(i,j) > G%areaT(i,j)) then call MOM_error(WARNING,"Initialize_ice_shelf: area_shelf_h exceeds G%areaT.") @@ -1773,6 +1731,62 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call hchksum(ISS%area_shelf_h, "IS init: area_shelf_h", G%HI, haloshift=0, scale=US%L_to_m*US%L_to_m) endif + + if (.not. complete_initialization) return + + ! Set up the restarts. + + call restart_init(param_file, CS%restart_CSp, "Shelf.res") + call register_restart_field(ISS%mass_shelf, "shelf_mass", .true., CS%restart_CSp, & + "Ice shelf mass", "kg m-2") + call register_restart_field(ISS%area_shelf_h, "shelf_area", .true., CS%restart_CSp, & + "Ice shelf area in cell", "m2") + call register_restart_field(ISS%h_shelf, "h_shelf", .true., CS%restart_CSp, & + "ice sheet/shelf thickness", "m") + if (PRESENT(sfc_state_in)) then + if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then + u_desc = var_desc("taux_shelf", "Pa", "the zonal stress on the ocean under ice shelves", & + hor_grid='Cu',z_grid='1') + v_desc = var_desc("tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", & + hor_grid='Cv',z_grid='1') + call register_restart_pair(sfc_state%taux_shelf, sfc_state%tauy_shelf, u_desc,v_desc, & + .false., CS%restart_CSp) + endif + endif + + call register_restart_field(ISS%h_shelf, "_shelf", .true., CS%restart_CSp, & + "ice sheet/shelf thickness", "m") + call register_restart_field(US%m_to_Z_restart, "m_to_Z", .false., CS%restart_CSp, & + "Height unit conversion factor", "Z meter-1") + call register_restart_field(US%m_to_L_restart, "m_to_L", .false., CS%restart_CSp, & + "Length unit conversion factor", "L meter-1") + call register_restart_field(US%kg_m3_to_R_restart, "kg_m3_to_R", .false., CS%restart_CSp, & + "Density unit conversion factor", "R m3 kg-1") + if (CS%active_shelf_dynamics) then + call register_restart_field(ISS%hmask, "h_mask", .true., CS%restart_CSp, & + "ice sheet/shelf thickness mask" ,"none") + endif + + if (CS%active_shelf_dynamics) then + ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics + call register_ice_shelf_dyn_restarts(CS%Grid_in, param_file, CS%dCS, CS%restart_CSp) + endif + + !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file + !if (.not. CS%solo_ice_sheet) then + ! call register_restart_field(fluxes%ustar_shelf, "ustar_shelf", .false., CS%restart_CSp, & + ! "Friction velocity under ice shelves", "m s-1") + !endif + + CS%restart_output_dir = dirs%restart_output_dir + + +! do j=G%jsc,G%jec ; do i=G%isc,G%iec +! ISS%area_shelf_h(i,j) = ISS%area_shelf_h(i,j)*G%mask2dT(i,j) +! enddo; enddo + + CS%Time = Time + if (present(forces_in)) & call add_shelf_forces(G, US, CS, forces, do_shelf_area=.not.CS%solo_ice_sheet) From 55c6eaa5721e2b18de060a6702276918d1fa86b9 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 21 Dec 2020 15:20:18 -0500 Subject: [PATCH 04/11] Modifications for coupled driver --- config_src/coupled_driver/ocean_model_MOM.F90 | 16 ++++++++-------- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 774201ddb5..b210cd4f81 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -182,13 +182,13 @@ module ocean_model_mod !! processes before time stepping the dynamics. type(directories) :: dirs !< A structure containing several relevant directory paths. - type(mech_forcing), pointer :: forces => NULL() !< A structure with the driving mechanical surface forces - type(forcing), pointer :: fluxes => NULL() !< A structure containing pointers to + type(mech_forcing) :: forces => NULL() !< A structure with the driving mechanical surface forces + type(forcing) :: fluxes => NULL() !< A structure containing pointers to !! the thermodynamic ocean forcing fields. - type(forcing), pointer :: flux_tmp => NULL() !< A secondary structure containing pointers to the + type(forcing) :: flux_tmp => NULL() !< A secondary structure containing pointers to the !! ocean forcing fields for when multiple coupled !! timesteps are taken per thermodynamic step. - type(surface), pointer :: sfc_state => NULL() !< A structure containing pointers to + type(surface) :: sfc_state => NULL() !< A structure containing pointers to !! the ocean surface state fields. type(ocean_grid_type), pointer :: & grid => NULL() !< A pointer to a grid structure containing metrics @@ -273,9 +273,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas endif allocate(OS) - allocate(OS%fluxes) - allocate(OS%forces) - allocate(OS%flux_tmp) +! allocate(OS%fluxes) +! allocate(OS%forces) +! allocate(OS%flux_tmp) OS%is_ocean_pe = Ocean_sfc%is_ocean_pe if (.not.OS%is_ocean_pe) return @@ -379,7 +379,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas if (OS%use_ice_shelf) then call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & - OS%diag_IS, OS%forces, OS%fluxes) + OS%diag, OS%forces, OS%fluxes) endif if (OS%icebergs_alter_ocean) then call marine_ice_init(OS%Time, OS%grid, param_file, OS%diag, OS%marine_ice_CSp) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 15326b85ec..236e0b2b34 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1207,7 +1207,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, real :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting ! does not occur [Z ~> m] real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data - type(mech_forcing), pointer :: forces => NULL() type(forcing), pointer :: fluxes => NULL() type(surface), pointer :: sfc_state => NULL() @@ -1294,6 +1293,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call MOM_initialize_topography(dG%bathyT, CS%Grid%max_depth, dG, param_file) call rescale_dyn_horgrid_bathymetry(dG, CS%US%Z_to_m) call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) + call destroy_dyn_horgrid(dG) ! endif G=>CS%Grid From 19c1a6af2b3afabf32879b5ee6f855092cd9b4cb Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 21 Dec 2020 16:31:29 -0500 Subject: [PATCH 05/11] fix compile issues --- config_src/coupled_driver/ocean_model_MOM.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index b210cd4f81..500c8ba0a2 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -182,13 +182,13 @@ module ocean_model_mod !! processes before time stepping the dynamics. type(directories) :: dirs !< A structure containing several relevant directory paths. - type(mech_forcing) :: forces => NULL() !< A structure with the driving mechanical surface forces - type(forcing) :: fluxes => NULL() !< A structure containing pointers to + type(mech_forcing) :: forces !< A structure with the driving mechanical surface forces + type(forcing) :: fluxes !< A structure containing pointers to !! the thermodynamic ocean forcing fields. - type(forcing) :: flux_tmp => NULL() !< A secondary structure containing pointers to the + type(forcing) :: flux_tmp !< A secondary structure containing pointers to the !! ocean forcing fields for when multiple coupled !! timesteps are taken per thermodynamic step. - type(surface) :: sfc_state => NULL() !< A structure containing pointers to + type(surface) :: sfc_state !< A structure containing pointers to !! the ocean surface state fields. type(ocean_grid_type), pointer :: & grid => NULL() !< A pointer to a grid structure containing metrics @@ -365,7 +365,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas use_melt_pot=.false. endif - allocate(OS%sfc_state) + !allocate(OS%sfc_state) call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., & gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot) From bdcb8588a38dc13276c30357b5941dfd8813a53c Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 22 Dec 2020 12:22:00 -0500 Subject: [PATCH 06/11] Move call to end ice shelf diag mediator into ice_shelf_end - retains original drivers --- config_src/coupled_driver/ocean_model_MOM.F90 | 6 ---- config_src/solo_driver/MOM_driver.F90 | 4 --- src/ice_shelf/MOM_ice_shelf.F90 | 2 ++ src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 | 29 +++++++++++++------ 4 files changed, 22 insertions(+), 19 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 500c8ba0a2..c4bd543bfd 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -48,7 +48,6 @@ 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_IS_diag_mediator, only : diag_IS_ctrl => diag_ctrl, diag_mediator_IS_end=>diag_mediator_end 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 @@ -217,9 +216,6 @@ module ocean_model_mod !! that will be used for MOM restart files. type(diag_ctrl), pointer :: & diag => NULL() !< A pointer to the diagnostic regulatory structure - type(diag_IS_ctrl), pointer :: & - diag_IS => NULL() !< A pointer to the diagnostic regulatory structure - !! for the ice shelf module. end type ocean_state_type contains @@ -728,8 +724,6 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) call ocean_model_save_restart(Ocean_state, Time) call diag_mediator_end(Time, Ocean_state%diag) - if (Ocean_state%use_ice_shelf) & - call diag_mediator_IS_end(Time, Ocean_state%diag_IS) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) end subroutine ocean_model_end diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index d61d6c986a..d36d86c8db 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -28,7 +28,6 @@ program MOM_main use MOM_cpu_clock, only : CLOCK_COMPONENT use MOM_diag_mediator, only : enable_averaging, disable_averaging, diag_mediator_end use MOM_diag_mediator, only : diag_ctrl, diag_mediator_close_registration - use MOM_IS_diag_mediator, only : diag_IS_ctrl=>diag_ctrl, diag_mediator_IS_end=>diag_mediator_end use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only : extract_surface_state, finish_MOM_initialization use MOM, only : get_MOM_state_elements, MOM_state_is_synchronized @@ -199,8 +198,6 @@ program MOM_main !! that will be used for MOM restart files. type(diag_ctrl), pointer :: & diag => NULL() !< A pointer to the diagnostic regulatory structure - type(diag_IS_ctrl), pointer :: & - diag_IS => NULL() !< A pointer to the diagnostic regulatory structure !----------------------------------------------------------------------- character(len=4), parameter :: vers_num = 'v2.0' @@ -665,7 +662,6 @@ program MOM_main call callTree_waypoint("End MOM_main") call diag_mediator_end(Time, diag, end_diag_manager=.true.) - if (use_ice_shelf) call diag_mediator_IS_end(Time, diag_IS) 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/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 236e0b2b34..7d5b7aeedf 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -15,6 +15,7 @@ module MOM_ice_shelf use MOM_IS_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid, diag_ctrl, time_type use MOM_IS_diag_mediator, only : enable_averages, enable_averaging, disable_averaging use MOM_IS_diag_mediator, only : diag_mediator_infrastructure_init, diag_mediator_close_registration +use MOM_IS_diag_mediator, only : diag_mediator_end use MOM_domains, only : MOM_domains_init, clone_MOM_domain use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE, CORNER use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid @@ -2059,6 +2060,7 @@ subroutine ice_shelf_end(CS) if (CS%active_shelf_dynamics) call ice_shelf_dyn_end(CS%dCS) + call diag_mediator_end(CS%diag) deallocate(CS) end subroutine ice_shelf_end diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index ab4245bd83..74d9ed701b 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -304,28 +304,41 @@ subroutine post_data(diag_field_id, field, diag_cs, is_static, mask) used = send_data(fms_diag_id, locfield, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, mask=mask) elseif(i_data .and. associated(diag%mask2d)) then +! used = send_data(fms_diag_id, locfield, & +! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask2d) used = send_data(fms_diag_id, locfield, & - is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask2d) + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev) elseif((.not.i_data) .and. associated(diag%mask2d_comp)) then +! used = send_data(fms_diag_id, locfield, & +! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask2d_comp) used = send_data(fms_diag_id, locfield, & - is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, rmask=diag%mask2d_comp) + is_in=isv, js_in=jsv, ie_in=iev, je_in=jev) else used = send_data(fms_diag_id, locfield, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev) endif elseif (diag_cs%ave_enabled) then if (present(mask)) then +! used = send_data(fms_diag_id, locfield, diag_cs%time_end, & +! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & +! weight=diag_cs%time_int, mask=mask) used = send_data(fms_diag_id, locfield, diag_cs%time_end, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & - weight=diag_cs%time_int, mask=mask) + weight=diag_cs%time_int) elseif(i_data .and. associated(diag%mask2d)) then +! used = send_data(fms_diag_id, locfield, diag_cs%time_end, & +! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & +! weight=diag_cs%time_int, rmask=diag%mask2d) used = send_data(fms_diag_id, locfield, diag_cs%time_end, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & - weight=diag_cs%time_int, rmask=diag%mask2d) + weight=diag_cs%time_int) elseif((.not.i_data) .and. associated(diag%mask2d_comp)) then +! used = send_data(fms_diag_id, locfield, diag_cs%time_end, & +! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & +! weight=diag_cs%time_int, rmask=diag%mask2d_comp) used = send_data(fms_diag_id, locfield, diag_cs%time_end, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & - weight=diag_cs%time_int, rmask=diag%mask2d_comp) + weight=diag_cs%time_int) else used = send_data(fms_diag_id, locfield, diag_cs%time_end, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & @@ -483,7 +496,6 @@ function register_MOM_IS_diag_field(module_name, field_name, axes, init_time, & !Decide what mask to use based on the axes info if (primary_id > 0) then - !3d masks !2d masks if (axes%rank == 2) then diag%mask2d => null() ; diag%mask2d_comp => null() @@ -682,7 +694,7 @@ subroutine diag_masks_set(G, missing_value, diag_cs) type(diag_ctrl), intent(inout) :: diag_cs !< A structure that is used to regulate diagnostic output ! Local variables - integer :: i, j, k, NkIce, CatIce + integer :: i, j diag_cs%mask2dT => G%mask2dT @@ -711,8 +723,7 @@ subroutine diag_mediator_close_registration(diag_CS) end subroutine diag_mediator_close_registration !> Deallocate memory associated with the MOM_IS diag mediator -subroutine diag_mediator_end(time, diag_CS) - type(time_type), intent(in) :: time !< The current model time +subroutine diag_mediator_end(diag_CS) type(diag_ctrl), intent(inout) :: diag_CS !< A structure that is used to regulate diagnostic output if (diag_CS%doc_unit > -1) then From 44e80d59fdce0faa535fc636a8b9faa2fa381237 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 22 Dec 2020 13:03:21 -0500 Subject: [PATCH 07/11] remove rotation related changes in initialize_ice_shelf --- src/ice_shelf/MOM_ice_shelf.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 7d5b7aeedf..8e60d81ba7 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -82,14 +82,14 @@ module MOM_ice_shelf ! Parameters type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart control !! structure for the ice shelves - type(ocean_grid_type) :: Grid_in !< un-rotated input grid metric + 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 + type(ocean_grid_type), pointer :: Grid => NULL() !< Grid for the ice-shelf model type(unit_scale_type), pointer :: & US => NULL() !< A structure containing various unit conversion factors type(ocean_grid_type), pointer :: ocn_grid => NULL() !< A pointer to the ocean model grid @@ -1248,12 +1248,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, ! Set up the ice-shelf domain and grid wd_halos(:)=0 - !allocate(CS%Grid_in) - call MOM_domains_init(CS%Grid_in%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_,& + allocate(CS%Grid) + call MOM_domains_init(CS%Grid%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_,& domain_name='MOM_Ice_Shelf_in') - call hor_index_init(CS%Grid_in%Domain, CS%Grid_in%HI, param_file, & - local_indexing=.not.global_indexing) - call MOM_grid_init(CS%Grid_in, param_file, CS%US, 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) ! if (CS%rotate_index) then ! ! TODO: Index rotation currently only works when index rotation does not @@ -1284,7 +1284,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, ! 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 + !CS%Grid=>CS%Grid_in dG=>NULL() !CS%Grid%HI=>CS%Grid_in%HI call create_dyn_horgrid(dG, CS%Grid%HI) @@ -1296,7 +1296,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) call destroy_dyn_horgrid(dG) ! endif - G=>CS%Grid + G=>CS%Grid;CS%Grid_in=>CS%Grid if (complete_initialization) then allocate(CS%diag) From 88c4102d64614b046754633cf31af52b1ca6fab8 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 22 Dec 2020 13:24:45 -0500 Subject: [PATCH 08/11] move complete_initialization return before diag chksums --- src/ice_shelf/MOM_ice_shelf.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 8e60d81ba7..56461dbc3d 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1727,14 +1727,15 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, if (G%areaT(i,j)>0.) fluxes%frac_shelf_h(i,j) = ISS%area_shelf_h(i,j) / G%areaT(i,j) enddo ; enddo ; endif + if (.not. complete_initialization) return + + if (CS%debug) then call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", G%HI, haloshift=0) call hchksum(ISS%area_shelf_h, "IS init: area_shelf_h", G%HI, haloshift=0, scale=US%L_to_m*US%L_to_m) endif - if (.not. complete_initialization) return - ! Set up the restarts. call restart_init(param_file, CS%restart_CSp, "Shelf.res") From d2d047c5e5665b3a82ed30dd11a4d9512cdd37e7 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 22 Dec 2020 13:25:37 -0500 Subject: [PATCH 09/11] remove masking from ice shelf diagnostics --- src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 index 74d9ed701b..5db9646dae 100644 --- a/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 +++ b/src/ice_shelf/MOM_ice_shelf_diag_mediator.F90 @@ -319,12 +319,12 @@ subroutine post_data(diag_field_id, field, diag_cs, is_static, mask) endif elseif (diag_cs%ave_enabled) then if (present(mask)) then -! used = send_data(fms_diag_id, locfield, diag_cs%time_end, & -! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & -! weight=diag_cs%time_int, mask=mask) used = send_data(fms_diag_id, locfield, diag_cs%time_end, & is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & - weight=diag_cs%time_int) + weight=diag_cs%time_int, mask=mask) +! used = send_data(fms_diag_id, locfield, diag_cs%time_end, & +! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & +! weight=diag_cs%time_int) elseif(i_data .and. associated(diag%mask2d)) then ! used = send_data(fms_diag_id, locfield, diag_cs%time_end, & ! is_in=isv, js_in=jsv, ie_in=iev, je_in=jev, & From a1167461becebf56822d8c7ca08750e789058c4d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 29 Dec 2020 15:15:03 -0500 Subject: [PATCH 10/11] +Added initialize_ice_shelf_fluxes Added the routines initialize_ice_shelf_fluxes and initialize_ice_shelf_forces which call be called from within initialize_ice_shelf or separately, which removes the need to call initialize_ice_shelf multiple times. Also added the new runtime parameters USTAR_SHELF_FROM_VEL and USTAR_SHELF_MAX, which will enable the previous answers for the ISOMIP test cases to be recovered and will facilitate the debugging or control of poorly understood instabilities related to the dynamic or lagged calculation of ustar_shelf. All answers in the usual MOM6-examples test cases are bitwise identical, and the ISOMIP test cases are working once again. --- config_src/solo_driver/MOM_driver.F90 | 20 +- src/ice_shelf/MOM_ice_shelf.F90 | 259 ++++++++++++++------------ 2 files changed, 146 insertions(+), 133 deletions(-) diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index d36d86c8db..9726aa1281 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -70,6 +70,7 @@ program MOM_main use MOM_ice_shelf, only : initialize_ice_shelf, ice_shelf_end, ice_shelf_CS 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_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init use MOM_wave_interface, only: MOM_wave_interface_init_lite, Update_Surface_Waves @@ -308,27 +309,24 @@ program MOM_main if (sum(date) >= 0) then call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, & segment_start_time, offline_tracer_mode=offline_tracer_mode, & - diag_ptr=diag, tracer_flow_CSp=tracer_flow_CSp,ice_shelf_CSp=ice_shelf_CSp) + diag_ptr=diag, tracer_flow_CSp=tracer_flow_CSp, ice_shelf_CSp=ice_shelf_CSp) else call initialize_MOM(Time, Start_time, param_file, dirs, MOM_CSp, restart_CSp, & offline_tracer_mode=offline_tracer_mode, diag_ptr=diag, & - tracer_flow_CSp=tracer_flow_CSp,ice_shelf_CSp=ice_shelf_CSp) + tracer_flow_CSp=tracer_flow_CSp, ice_shelf_CSp=ice_shelf_CSp) endif - call get_param(param_file, mod_name, "ICE_SHELF", use_ice_shelf, & - "If true, enables the ice shelf model.", default=.false.) + call get_MOM_state_elements(MOM_CSp, G=grid, GV=GV, US=US, C_p_scaled=fluxes%C_p) + Master_Time = Time + use_ice_shelf = associated(ice_shelf_CSp) if (use_ice_shelf) then ! These arrays are not initialized in most solo cases, but are needed ! when using an ice shelf - ice_shelf_CSp => NULL() ! Reset the pointer and make another call to reinitialize - ! the ice shelf and associated forcing types - call initialize_ice_shelf(param_file, grid, Time, ice_shelf_CSp, & - diag, forces, fluxes, sfc_state) + call initialize_ice_shelf_fluxes(ice_shelf_CSp, grid, US, fluxes) + call initialize_ice_shelf_forces(ice_shelf_CSp, grid, US, forces) endif - call get_MOM_state_elements(MOM_CSp, G=grid, GV=GV, US=US, C_p_scaled=fluxes%C_p) - Master_Time = Time call callTree_waypoint("done initialize_MOM") @@ -661,6 +659,7 @@ program MOM_main endif call callTree_waypoint("End MOM_main") + if (use_ice_shelf) call ice_shelf_end(ice_shelf_CSp) 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) @@ -668,6 +667,5 @@ program MOM_main call io_infra_end ; call MOM_infra_end call MOM_end(MOM_CSp) - if (use_ice_shelf) call ice_shelf_end(ice_shelf_CSp) end program MOM_main diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 56461dbc3d..9592d17b03 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -11,11 +11,10 @@ module MOM_ice_shelf use MOM_coms, only : num_PEs use MOM_diag_mediator, only : MOM_diag_ctrl=>diag_ctrl use MOM_IS_diag_mediator, only : post_data, register_diag_field=>register_MOM_IS_diag_field, safe_alloc_ptr -use MOM_IS_diag_mediator, only : set_axes_info -use MOM_IS_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid, diag_ctrl, time_type +use MOM_IS_diag_mediator, only : set_axes_info, diag_ctrl, time_type +use MOM_IS_diag_mediator, only : diag_mediator_init, diag_mediator_end, set_diag_mediator_grid use MOM_IS_diag_mediator, only : enable_averages, enable_averaging, disable_averaging use MOM_IS_diag_mediator, only : diag_mediator_infrastructure_init, diag_mediator_close_registration -use MOM_IS_diag_mediator, only : diag_mediator_end use MOM_domains, only : MOM_domains_init, clone_MOM_domain use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE, CORNER use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid @@ -71,6 +70,7 @@ module MOM_ice_shelf public shelf_calc_flux, initialize_ice_shelf, ice_shelf_end, ice_shelf_query public ice_shelf_save_restart, solo_step_ice_shelf, add_shelf_forces +public initialize_ice_shelf_fluxes, initialize_ice_shelf_forces ! 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 @@ -82,14 +82,14 @@ module MOM_ice_shelf ! Parameters 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(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 + type(ocean_grid_type), pointer :: Grid => NULL() !< Grid for the ice-shelf model type(unit_scale_type), pointer :: & US => NULL() !< A structure containing various unit conversion factors type(ocean_grid_type), pointer :: ocn_grid => NULL() !< A pointer to the ocean model grid @@ -105,6 +105,8 @@ module MOM_ice_shelf utide => NULL() !< An unresolved tidal velocity [L T-1 ~> m s-1] real :: ustar_bg !< A minimum value for ustar under ice shelves [Z T-1 ~> m s-1]. + real :: ustar_max !< A maximum value for ustar under ice shelves, or a negative value to + !! have no limit [Z T-1 ~> m s-1]. real :: cdrag !< drag coefficient under ice shelves [nondim]. real :: g_Earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] real :: Cp !< The heat capacity of sea water [Q degC-1 ~> J kg-1 degC-1]. @@ -126,6 +128,8 @@ module MOM_ice_shelf real :: col_mass_melt_threshold !< An ocean column mass below the iceshelf below which melting !! does not occur [R Z ~> kg m-2] logical :: mass_from_file !< Read the ice shelf mass from a file every dt + logical :: ustar_shelf_from_vel !< If true, use the surface velocities, and not the previous + !! values of the stresses to set ustar. !!!! PHYSICAL AND NUMERICAL PARAMETERS FOR ICE DYNAMICS !!!!!! @@ -388,9 +392,14 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS) u2_av = (asu1 * sfc_state%u(I-1,j)**2 + asu2 * sfc_state%u(I,j)**2) * I_au v2_av = (asv1 * sfc_state%v(i,J-1)**2 + asu2 * sfc_state%v(i,J)**2) * I_av - if (taux2 + tauy2 > 0.0) then - fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_to_Z * & - sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%cdrag*CS%utide(i,j)**2)) + if ((taux2 + tauy2 > 0.0) .and. .not.CS%ustar_shelf_from_vel) then + if (CS%ustar_max >= 0.0) then + fluxes%ustar_shelf(i,j) = MIN(CS%ustar_max, MAX(CS%ustar_bg, US%L_to_Z * & + sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%cdrag*CS%utide(i,j)**2))) + else + fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_to_Z * & + sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%cdrag*CS%utide(i,j)**2)) + endif else ! Take care of the cases when taux_shelf is not set or not allocated. fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_TO_Z * & sqrt(CS%cdrag*((u2_av + v2_av) + CS%utide(i,j)**2))) @@ -819,7 +828,7 @@ subroutine add_shelf_forces(Ocn_grid, US, CS, forces, do_shelf_area, external_ca !! is using the input grid metric and needs !! 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), target :: 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. @@ -833,14 +842,10 @@ subroutine add_shelf_forces(Ocn_grid, US, CS, forces, do_shelf_area, external_ca 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.") + call MOM_error(FATAL,"add_shelf_forces: Incompatible Ocean and Ice shelf grids.") if (CS%rotate_index .and. rotate) then call MOM_error(FATAL,"add_shelf_forces: Rotation not implemented for ice shelves.") - ! 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.") - ! allocate(forces) ! call allocate_mech_forcing(forces_in, CS%Grid, forces) ! call rotate_mech_forcing(forces_in, CS%turns, forces) @@ -929,7 +934,8 @@ subroutine add_shelf_pressure(Ocn_grid, US, CS, fluxes) 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), intent(in) :: CS !< This module's control structure. - type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be updated. + type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be updated. + type(ocean_grid_type), pointer :: G => NULL() ! A pointer to ocean's grid structure. real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa]. integer :: i, j, is, ie, js, je, isd, ied, jsd, jed @@ -1162,8 +1168,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure type(time_type), intent(inout) :: Time !< The clock that that will indicate the model time type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure - type(MOM_diag_ctrl), pointer :: diag !< This is a pointer to the MOM diag CS - !! which will be discarded + type(MOM_diag_ctrl), pointer :: diag !< This is a pointer to the MOM diag CS + !! which will be discarded type(mech_forcing), optional, target, intent(inout) :: forces_in !< A structure with the driving mechanical forces type(forcing), optional, target, intent(inout) :: fluxes_in !< A structure containing pointers to any @@ -1208,15 +1214,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, real :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting ! does not occur [Z ~> m] real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data + type(mech_forcing), pointer :: forces => NULL() type(forcing), pointer :: fluxes => NULL() type(surface), pointer :: sfc_state => NULL() type(vardesc) :: u_desc, v_desc - logical :: complete_initialization ! A flag which is set to true if forces are present - ! This exists for legacy reasons and is a means to avoid some - ! parts of the initilization procedure since the ice shelf - ! is being initialized twice from initialize MOM and from the - ! various driver routines. + if (associated(CS)) then call MOM_error(FATAL, "MOM_ice_shelf.F90, initialize_ice_shelf: "// & "called with an associated control structure.") @@ -1224,8 +1227,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif allocate(CS) - complete_initialization=.false. - if (present(forces_in)) complete_initialization = .true. ! Go through all of the infrastructure initialization calls, since this is ! being treated as an independent component that just happens to use the ! MOM's grid and infrastructure. @@ -1251,6 +1252,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, allocate(CS%Grid) call MOM_domains_init(CS%Grid%domain, param_file, min_halo=wd_halos, symmetric=GRID_SYM_,& domain_name='MOM_Ice_Shelf_in') +! 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) @@ -1285,7 +1287,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, ! call copy_dyngrid_to_MOM_grid(dG,CS%Grid,CS%US) ! else !CS%Grid=>CS%Grid_in - dG=>NULL() + 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) @@ -1296,15 +1298,13 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, 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 (complete_initialization) then - allocate(CS%diag) - call diag_mediator_init(G, param_file,CS%diag,component='MOM_IceShelf') - ! This call sets up the diagnostic axes. These are needed, - ! e.g. to generate the target grids below. - call set_axes_info(G, param_file, CS%diag) - endif + G => CS%Grid ; CS%Grid_in => CS%Grid + + allocate(CS%diag) + call diag_mediator_init(G, param_file, CS%diag, component='MOM_IceShelf') + ! This call sets up the diagnostic axes. These are needed, + ! e.g. to generate the target grids below. + call set_axes_info(G, param_file, CS%diag) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -1499,33 +1499,30 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif - if (complete_initialization) then - call safe_alloc_ptr(CS%utide,isd,ied,jsd,jed) ; CS%utide(:,:) = 0.0 - if (read_TIDEAMP) then - call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, & + call safe_alloc_ptr(CS%utide,isd,ied,jsd,jed) ; CS%utide(:,:) = 0.0 + if (read_TIDEAMP) then + call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, & "The path to the file containing the spatially varying "//& "tidal amplitudes.", & default="tideamp.nc") - call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - TideAmp_file = trim(inputdir) // trim(TideAmp_file) - 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) - deallocate(tmp2d) - else - call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) - endif + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + TideAmp_file = trim(inputdir) // trim(TideAmp_file) + 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) + deallocate(tmp2d) else - call get_param(param_file, mdl, "UTIDE", utide, & + call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, CS%Grid%domain, timelevel=1, scale=US%m_s_to_L_T) + endif + else + call get_param(param_file, mdl, "UTIDE", utide, & "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & units="m s-1", default=0.0 , scale=US%m_s_to_L_T) - CS%utide(:,:) = utide - endif + CS%utide(:,:) = utide endif - call EOS_init(param_file, CS%eqn_of_state) !! new parameters that need to be in MOM_input @@ -1565,58 +1562,19 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, "velocity magnitude.", units="m s-1", default=0.0, scale=US%m_to_Z*US%T_to_s) if (CS%cdrag*drag_bg_vel > 0.0) CS%ustar_bg = sqrt(CS%cdrag)*drag_bg_vel endif + call get_param(param_file, mdl, "USTAR_SHELF_FROM_VEL", CS%ustar_shelf_from_vel, & + "If true, use the surface velocities to set the friction velocity under ice "//& + "shelves instead of using the previous values of the stresses.", & + default=.true.) + call get_param(param_file, mdl, "USTAR_SHELF_MAX", CS%ustar_max, & + "The maximum value of ustar under ice shelves, or a negative value for no limit.", & + units="m s-1", default=-1.0, scale=US%m_to_Z*US%T_to_s, & + do_not_log=CS%ustar_shelf_from_vel) ! Allocate and initialize state variables to default values call ice_shelf_state_init(CS%ISS, CS%grid) ISS => CS%ISS - ! Allocate the arrays for passing ice-shelf data through the forcing type. - if (.not. CS%solo_ice_sheet) then - call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes.") - ! GMM: the following assures that water/heat fluxes are just allocated - ! when SHELF_THERMO = True. These fluxes are necessary if one wants to - ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode). - if (present(fluxes_in)) then - call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., & - press=.true., water=CS%isthermo, heat=CS%isthermo) - if (CS%rotate_index) then - allocate(fluxes) - call allocate_forcing_type(fluxes_in, CS%Grid, fluxes) - call rotate_forcing(fluxes_in, fluxes, CS%turns) - else - fluxes=>fluxes_in - endif - endif - if (present(forces_in)) then - call allocate_mech_forcing(CS%Grid_in, forces_in, ustar=.true., shelf=.true., press=.true.) - if (CS%rotate_index) then - allocate(forces) - call allocate_mech_forcing(forces_in, CS%Grid, forces) - call rotate_mech_forcing(forces_in, CS%turns, forces) - else - forces=>forces_in - endif - endif - else - call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.") - if (present(fluxes_in)) then - call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., press=.true.) - if (CS%rotate_index) then - allocate(fluxes) - call allocate_forcing_type(fluxes_in, CS%Grid, fluxes) - call rotate_forcing(fluxes_in, fluxes, CS%turns) - endif - endif - if (present(forces_in)) then - call allocate_mech_forcing(CS%Grid_in, forces_in, ustar=.true., shelf=.true., press=.true.) - if (CS%rotate_index) then - allocate(forces) - call allocate_mech_forcing(forces_in, CS%Grid, forces) - call rotate_mech_forcing(forces_in, CS%turns, forces) - endif - endif - endif - new_sim = .false. if ((dirs%input_filename(1:1) == 'n') .and. & (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. @@ -1706,6 +1664,10 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif ! .not. new_sim +! do j=G%jsc,G%jec ; do i=G%isc,G%iec +! ISS%area_shelf_h(i,j) = ISS%area_shelf_h(i,j)*G%mask2dT(i,j) +! enddo; enddo + id_clock_shelf = cpu_clock_id('Ice shelf', grain=CLOCK_COMPONENT) id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=CLOCK_ROUTINE) @@ -1723,15 +1685,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, ISS%area_shelf_h(i,j) = G%areaT(i,j) endif enddo ; enddo - if (present(fluxes_in)) then ; do j=jsd,jed ; do i=isd,ied - if (G%areaT(i,j)>0.) fluxes%frac_shelf_h(i,j) = ISS%area_shelf_h(i,j) / G%areaT(i,j) - enddo ; enddo ; endif - - if (.not. complete_initialization) return - if (CS%debug) then - call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", G%HI, haloshift=0) call hchksum(ISS%area_shelf_h, "IS init: area_shelf_h", G%HI, haloshift=0, scale=US%L_to_m*US%L_to_m) endif @@ -1782,17 +1737,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, CS%restart_output_dir = dirs%restart_output_dir - -! do j=G%jsc,G%jec ; do i=G%isc,G%iec -! ISS%area_shelf_h(i,j) = ISS%area_shelf_h(i,j)*G%mask2dT(i,j) -! enddo; enddo - CS%Time = Time - if (present(forces_in)) & - call add_shelf_forces(G, US, CS, forces, do_shelf_area=.not.CS%solo_ice_sheet) - - if (present(fluxes_in)) call add_shelf_pressure(ocn_grid, US, CS, fluxes) if (CS%active_shelf_dynamics .and. .not.CS%isthermo) then ISS%water_flux(:,:) = 0.0 @@ -1860,14 +1806,83 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif call diag_mediator_close_registration(CS%diag) - - if (present(fluxes_in) .and. CS%rotate_index) & - call rotate_forcing(fluxes, fluxes_in, -CS%turns) - if (present(forces_in) .and. CS%rotate_index) & - call rotate_mech_forcing(forces, -CS%turns, forces_in) + if (present(fluxes_in)) call initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in) + if (present(forces_in)) call initialize_ice_shelf_forces(CS, ocn_grid, US, forces_in) end subroutine initialize_ice_shelf +subroutine initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_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(forcing), target, intent(inout) :: fluxes_in !< A structure containing pointers to any + !! possible thermodynamic or mass-flux forcing fields. + + ! Local variables + type(ocean_grid_type), pointer :: G => NULL() ! Pointers to grids for convenience. + type(forcing), pointer :: fluxes => NULL() + integer :: i, j, isd, ied, jsd, jed + + G => CS%Grid + isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed + + ! Allocate the arrays for passing ice-shelf data through the forcing type. + if (.not. CS%solo_ice_sheet) then + call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes.") + ! GMM: the following assures that water/heat fluxes are just allocated + ! when SHELF_THERMO = True. These fluxes are necessary if one wants to + ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode). + call allocate_forcing_type(CS%Grid_in, fluxes_in, ustar=.true., shelf=.true., & + press=.true., water=CS%isthermo, heat=CS%isthermo) + 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.) + endif + if (CS%rotate_index) then + allocate(fluxes) + call allocate_forcing_type(fluxes_in, CS%Grid, fluxes) + call rotate_forcing(fluxes_in, fluxes, CS%turns) + else + fluxes=>fluxes_in + endif + + do j=jsd,jed ; do i=isd,ied + if (G%areaT(i,j)>0.) fluxes%frac_shelf_h(i,j) = CS%ISS%area_shelf_h(i,j) / G%areaT(i,j) + enddo ; enddo + if (CS%debug) call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", G%HI, haloshift=0) + call add_shelf_pressure(ocn_grid, US, CS, fluxes) + + if (CS%rotate_index) & + call rotate_forcing(fluxes, fluxes_in, -CS%turns) + +end subroutine initialize_ice_shelf_fluxes + +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(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.") + call allocate_mech_forcing(CS%Grid_in, forces_in, ustar=.true., shelf=.true., press=.true.) + if (CS%rotate_index) then + allocate(forces) + call allocate_mech_forcing(forces_in, CS%Grid, forces) + call rotate_mech_forcing(forces_in, CS%turns, forces) + else + forces=>forces_in + endif + + call add_shelf_forces(ocn_grid, US, CS, forces, do_shelf_area=.not.CS%solo_ice_sheet) + + if (CS%rotate_index) & + call rotate_mech_forcing(forces, -CS%turns, forces_in) + +end subroutine initialize_ice_shelf_forces + !> Initializes shelf mass based on three options (file, zero and user) subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) @@ -1952,7 +1967,7 @@ end subroutine initialize_shelf_mass !> Updates the ice shelf mass using data from a file. subroutine update_shelf_mass(G, US, CS, ISS, Time) type(ocean_grid_type), intent(inout) :: G !< The ocean's 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(ice_shelf_CS), intent(in) :: CS !< A pointer to the ice shelf control structure type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated type(time_type), intent(in) :: Time !< The current model time From 464f39e06edddd2d54dd1af6f6cef67fd617aab1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 3 Jan 2021 06:13:51 -0500 Subject: [PATCH 11/11] Remove extra register_restart call for ISS%h_shelf Removed the duplicate register_restart_field call for ISS%h_shelf under the name "_shelf". This registration call should have never been there in the first place. All answers are bitwise identical, although there are changes to the restart file contents in cases with an active ice shelf. --- src/ice_shelf/MOM_ice_shelf.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 9592d17b03..335db29ccb 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1711,8 +1711,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, endif endif - call register_restart_field(ISS%h_shelf, "_shelf", .true., CS%restart_CSp, & - "ice sheet/shelf thickness", "m") call register_restart_field(US%m_to_Z_restart, "m_to_Z", .false., CS%restart_CSp, & "Height unit conversion factor", "Z meter-1") call register_restart_field(US%m_to_L_restart, "m_to_L", .false., CS%restart_CSp, &