From 1d89598c042029a90d223aa3828385cc54da3888 Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Mon, 21 Dec 2020 09:25:52 -0700 Subject: [PATCH 01/18] correct CFC index check in src/tracer/MOM_OCMIP2_CFC.F90 --- src/tracer/MOM_OCMIP2_CFC.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90 index 9aad84a6dd..c568f4cacc 100644 --- a/src/tracer/MOM_OCMIP2_CFC.F90 +++ b/src/tracer/MOM_OCMIP2_CFC.F90 @@ -126,7 +126,7 @@ function register_OCMIP2_CFC(HI, GV, param_file, CS, tr_Reg, restart_CS) ! This call sets default properties for the air-sea CFC fluxes and obtains the ! indicies for the CFC11 and CFC12 flux coupling. call flux_init_OCMIP2_CFC(CS, verbosity=3) - if ((CS%ind_cfc_11_flux < 0) .or. (CS%ind_cfc_11_flux < 0)) then + if ((CS%ind_cfc_11_flux < 0) .or. (CS%ind_cfc_12_flux < 0)) then ! This is most likely to happen with the dummy version of aof_set_coupler_flux ! used in ocean-only runs. call MOM_ERROR(WARNING, "CFCs are currently only set up to be run in " // & From 2c93933a971e9054058ff3512c79622ff3a9cd92 Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Mon, 21 Dec 2020 12:43:08 -0500 Subject: [PATCH 02/18] 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 03/18] 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 04/18] 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 05/18] 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 06/18] 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 07/18] 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 08/18] 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 09/18] 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 10/18] 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 11/18] +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 12/18] 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, & From 1e6bdd6a1788bdafa7842d66e47d7326e41764f1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 3 Jan 2021 06:44:25 -0500 Subject: [PATCH 13/18] Remove unused module use statements Removed unused module use statements in various modules, to help eliminate apparent but inaccurate module dependencies, and to facilitate the migration to FMS2. All answers are bitwise identical. --- src/core/MOM_dynamics_split_RK2.F90 | 2 +- src/core/MOM_dynamics_unsplit.F90 | 1 - src/core/MOM_dynamics_unsplit_RK2.F90 | 1 - src/framework/MOM_diag_remap.F90 | 1 - src/framework/MOM_horizontal_regridding.F90 | 11 ++++----- src/framework/MOM_transform_FMS.F90 | 2 +- .../MOM_tracer_initialization_from_Z.F90 | 24 ++++++------------- src/tracer/advection_test_tracer.F90 | 2 +- src/tracer/boundary_impulse_tracer.F90 | 5 ++-- src/tracer/dye_example.F90 | 2 +- src/tracer/pseudo_salt_tracer.F90 | 5 ++-- src/user/BFB_surface_forcing.F90 | 1 - src/user/ISOMIP_initialization.F90 | 4 +--- src/user/dumbbell_surface_forcing.F90 | 1 - src/user/user_revise_forcing.F90 | 2 +- 15 files changed, 21 insertions(+), 43 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 519f510239..4ea6734511 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -27,7 +27,7 @@ module MOM_dynamics_split_RK2 use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories -use MOM_io, only : MOM_io_init, vardesc, var_desc +use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, register_restart_pair use MOM_restart, only : query_initialized, save_restart use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 6b9aa8e759..30544b0193 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -68,7 +68,6 @@ module MOM_dynamics_unsplit use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories -use MOM_io, only : MOM_io_init use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, MOM_restart_CS use MOM_time_manager, only : time_type, real_to_time, operator(+) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 4181ab519d..2f93561c3f 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -67,7 +67,6 @@ module MOM_dynamics_unsplit_RK2 use MOM_error_handler, only : MOM_set_verbosity use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories -use MOM_io, only : MOM_io_init use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, MOM_restart_CS use MOM_time_manager, only : time_type, time_type_to_real, operator(+) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 4e12abaa5b..6d1fa7b6fa 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -80,7 +80,6 @@ module MOM_diag_remap use coord_sigma, only : build_sigma_column use coord_rho, only : build_rho_column -use diag_axis_mod, only : get_diag_axis_name use diag_manager_mod, only : diag_axis_init use MOM_debugging, only : check_column_integrals diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 66f58b5b9d..8af6129812 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -24,15 +24,12 @@ module MOM_horizontal_regridding use MOM_time_manager, only : get_external_field_axes, get_external_field_missing use MOM_transform_FMS, only : time_interp_external => rotated_time_interp_external use MOM_variables, only : thermo_var_ptrs -use mpp_io_mod, only : axistype -use mpp_domains_mod, only : mpp_global_field, mpp_get_compute_domain -use mpp_mod, only : mpp_broadcast,mpp_root_pe,mpp_sync,mpp_sync_self -use mpp_mod, only : mpp_max -use horiz_interp_mod, only : horiz_interp_new, horiz_interp,horiz_interp_type + +use mpp_io_mod, only : axistype, mpp_get_axis_data +use mpp_mod, only : mpp_broadcast, mpp_sync, mpp_sync_self, mpp_max +use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_type use horiz_interp_mod, only : horiz_interp_init, horiz_interp_del -use mpp_io_mod, only : mpp_get_axis_data -use mpp_io_mod, only : MPP_SINGLE use netcdf implicit none ; private diff --git a/src/framework/MOM_transform_FMS.F90 b/src/framework/MOM_transform_FMS.F90 index 97e0be85f6..572a9717dc 100644 --- a/src/framework/MOM_transform_FMS.F90 +++ b/src/framework/MOM_transform_FMS.F90 @@ -7,7 +7,7 @@ module MOM_transform_FMS use MOM_error_handler, only : MOM_error, FATAL use MOM_io, only : fieldtype, write_field use mpp_domains_mod, only : domain2D -use fms_mod, only : mpp_chksum +use mpp_mod, only : mpp_chksum use time_manager_mod, only : time_type use time_interp_external_mod, only : time_interp_external diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 1a4c5bd011..12235ddd87 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -4,27 +4,17 @@ module MOM_tracer_initialization_from_Z ! This file is part of MOM6. See LICENSE.md for the license. use MOM_debugging, only : hchksum -use MOM_coms, only : max_across_PEs, min_across_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP -use MOM_density_integrals, only : int_specific_vol_dp -use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast -use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID -use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe +use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP +use MOM_domains, only : pass_var +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint -use MOM_file_parser, only : get_param, read_param, log_param, param_file_type -use MOM_file_parser, only : log_version -use MOM_get_input, only : directories -use MOM_grid, only : ocean_grid_type, isPointInCell +use MOM_file_parser, only : get_param, param_file_type, log_version +use MOM_grid, only : ocean_grid_type use MOM_horizontal_regridding, only : myStats, horiz_interp_and_extrap_tracer -use MOM_regridding, only : regridding_CS use MOM_remapping, only : remapping_CS, initialize_remapping -use MOM_remapping, only : remapping_core_h -use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type, setVerticalGridAxes -use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type +use MOM_verticalGrid, only : verticalGrid_type use MOM_ALE, only : ALE_remap_scalar implicit none ; private @@ -42,7 +32,7 @@ module MOM_tracer_initialization_from_Z contains -!> Initializes a tracer from a z-space data file. +!> Initializes a tracer from a z-space data file, including any lateral regridding that is needed. subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_nam, & src_var_unit_conversion, src_var_record, homogenize, & useALEremapping, remappingScheme, src_var_gridspec ) diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90 index b1d657d6e2..3aa65e8b3c 100644 --- a/src/tracer/advection_test_tracer.F90 +++ b/src/tracer/advection_test_tracer.F90 @@ -9,7 +9,7 @@ module advection_test_tracer use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_io, only : slasher, vardesc, var_desc, query_vardesc use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS diff --git a/src/tracer/boundary_impulse_tracer.F90 b/src/tracer/boundary_impulse_tracer.F90 index da76cb3026..fc85b5c3ec 100644 --- a/src/tracer/boundary_impulse_tracer.F90 +++ b/src/tracer/boundary_impulse_tracer.F90 @@ -9,7 +9,7 @@ module boundary_impulse_tracer use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_io, only : vardesc, var_desc, query_vardesc use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS @@ -18,8 +18,7 @@ module boundary_impulse_tracer use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : surface -use MOM_variables, only : thermo_var_ptrs +use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use coupler_types_mod, only : coupler_type_set_data, ind_csurf diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index cd17415b21..8a970fa9ca 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -9,7 +9,7 @@ module regional_dyes use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_io, only : vardesc, var_desc, query_vardesc use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS diff --git a/src/tracer/pseudo_salt_tracer.F90 b/src/tracer/pseudo_salt_tracer.F90 index 95396a3b58..11238fee89 100644 --- a/src/tracer/pseudo_salt_tracer.F90 +++ b/src/tracer/pseudo_salt_tracer.F90 @@ -11,7 +11,7 @@ module pseudo_salt_tracer use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type -use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_io, only : vardesc, var_desc, query_vardesc use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS @@ -20,8 +20,7 @@ module pseudo_salt_tracer use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : surface -use MOM_variables, only : thermo_var_ptrs +use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use coupler_types_mod, only : coupler_type_set_data, ind_csurf diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index d06262b7cf..3963d4d90d 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -10,7 +10,6 @@ module BFB_surface_forcing use MOM_file_parser, only : get_param, param_file_type, log_version use MOM_forcing_type, only : forcing, allocate_forcing_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists, read_data use MOM_safe_alloc, only : safe_alloc_ptr use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_tracer_flow_control, only : call_tracer_set_forcing diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index d125495d7f..a0b8990e62 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -10,9 +10,7 @@ module ISOMIP_initialization use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists -use MOM_io, only : MOM_read_data -use MOM_io, only : slasher +use MOM_io, only : file_exists, MOM_read_data, slasher use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index 4b5bf5a2fb..ea27d01cdc 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -10,7 +10,6 @@ module dumbbell_surface_forcing use MOM_file_parser, only : get_param, param_file_type, log_version use MOM_forcing_type, only : forcing, allocate_forcing_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists, read_data use MOM_safe_alloc, only : safe_alloc_ptr use MOM_time_manager, only : time_type, operator(+), operator(/), get_time use MOM_tracer_flow_control, only : call_tracer_set_forcing diff --git a/src/user/user_revise_forcing.F90 b/src/user/user_revise_forcing.F90 index c53451f4e8..bf31ca02f8 100644 --- a/src/user/user_revise_forcing.F90 +++ b/src/user/user_revise_forcing.F90 @@ -8,7 +8,7 @@ module user_revise_forcing use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists, read_data +use MOM_io, only : file_exists, MOM_read_data use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_tracer_flow_control, only : call_tracer_set_forcing From 9b0b8db88ad92a62b36afbd8005cab182e60fb77 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 3 Jan 2021 09:45:39 -0500 Subject: [PATCH 14/18] Use MOM framework routines in MOM_open_boundary Use MOM framework interfaces in MOM_open_boundary in place of direct calls to mpp routines, to facilitate the migration to FMS2. All answers are bitwise identical. --- src/core/MOM_open_boundary.F90 | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 46d144a8c6..0232ff91ff 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -5,13 +5,12 @@ module MOM_open_boundary use MOM_array_transform, only : rotate_array, rotate_array_pair use MOM_array_transform, only : allocate_rotated_array -use MOM_coms, only : sum_across_PEs +use MOM_coms, only : sum_across_PEs, Set_PElist, Get_PElist, PE_here, num_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_domains, only : pass_var, pass_vector use MOM_domains, only : To_All, SCALAR_PAIR, CGRID_NE, CORNER -use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe -use MOM_error_handler, only : NOTE +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, NOTE, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type, log_param use MOM_grid, only : ocean_grid_type, hor_index_type use MOM_dyn_horgrid, only : dyn_horgrid_type @@ -651,13 +650,11 @@ end subroutine open_boundary_config !> Allocate space for reading OBC data from files. It sets up the required vertical !! remapping. In the process, it does funky stuff with the MPI processes. subroutine initialize_segment_data(G, OBC, PF) - use mpp_mod, only : mpp_pe, mpp_set_current_pelist, mpp_get_current_pelist,mpp_npes - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure - type(param_file_type), intent(in) :: PF !< Parameter file handle + type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure + type(param_file_type), intent(in) :: PF !< Parameter file handle - integer :: n,m,num_fields + integer :: n, m, num_fields character(len=1024) :: segstr character(len=256) :: filename character(len=20) :: segnam, suffix @@ -697,11 +694,11 @@ subroutine initialize_segment_data(G, OBC, PF) !< temporarily disable communication in order to read segment data independently - allocate(saved_pelist(0:mpp_npes()-1)) - call mpp_get_current_pelist(saved_pelist) - current_pe = mpp_pe() + allocate(saved_pelist(0:num_PEs()-1)) + call Get_PElist(saved_pelist) + current_pe = PE_here() single_pelist(1) = current_pe - call mpp_set_current_pelist(single_pelist) + call Set_PElist(single_pelist) do n=1, OBC%number_of_segments segment => OBC%segment(n) @@ -955,7 +952,7 @@ subroutine initialize_segment_data(G, OBC, PF) endif enddo - call mpp_set_current_pelist(saved_pelist) + call Set_PElist(saved_pelist) end subroutine initialize_segment_data From 611132731dd731b7f2f4d26bcde2188b218fafef Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 3 Jan 2021 09:46:01 -0500 Subject: [PATCH 15/18] Use MOM_read_data in RGC_initialization Use MOM_read_data in place of read_data in RGC_initialization to match the routines used in other modules and facilitate migration to FMS2. All answers are bitwise identical. --- src/user/RGC_initialization.F90 | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/user/RGC_initialization.F90 b/src/user/RGC_initialization.F90 index 70b9fcd4dc..1600aca5bd 100644 --- a/src/user/RGC_initialization.F90 +++ b/src/user/RGC_initialization.F90 @@ -28,8 +28,7 @@ module RGC_initialization use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type -use MOM_io, only : file_exists, read_data -use MOM_io, only : slasher +use MOM_io, only : file_exists, MOM_read_data, slasher use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge use MOM_sponge, only : set_up_sponge_ML_density use MOM_unit_scaling, only : unit_scale_type @@ -173,12 +172,12 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) filename = trim(inputdir)//trim(state_file) if (.not.file_exists(filename, G%Domain)) & call MOM_error(FATAL, " RGC_initialize_sponges: Unable to open "//trim(filename)) - call read_data(filename,temp_var,T(:,:,:), domain=G%Domain%mpp_domain) - call read_data(filename,salt_var,S(:,:,:), domain=G%Domain%mpp_domain) + call MOM_read_data(filename, temp_var, T(:,:,:), G%Domain) + call MOM_read_data(filename, salt_var, S(:,:,:), G%Domain) if (use_ALE) then - call read_data(filename,h_var,h(:,:,:), domain=G%Domain%mpp_domain) + call MOM_read_data(filename, h_var, h(:,:,:), G%Domain) call pass_var(h, G%domain) !call initialize_ALE_sponge(Idamp, h, nz, G, PF, ACSp) @@ -201,7 +200,7 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) else ! layer mode !read eta - call read_data(filename,eta_var,eta(:,:,:), domain=G%Domain%mpp_domain) + call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain) ! Set the inverse damping rates so that the model will know where to ! apply the sponges, along with the interface heights. From 0b019b66773b6b8589b68fe565193cab10e6307f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 3 Jan 2021 09:58:00 -0500 Subject: [PATCH 16/18] Avoid using memory macros in MOM_random.F90 Expanded the SZI_ and SZJ_ macros in random_2d_ routines to eliminate any dependence on MOM_memory.h and facilitate the future compilation of MOM_random as a part of a MOM framework library. All answers are bitwise identical. --- src/framework/MOM_random.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/framework/MOM_random.F90 b/src/framework/MOM_random.F90 index 161236572c..21e3223a03 100644 --- a/src/framework/MOM_random.F90 +++ b/src/framework/MOM_random.F90 @@ -11,7 +11,7 @@ module MOM_random use MersenneTwister_mod, only : getRandomReal ! Generates a random number use MersenneTwister_mod, only : getRandomPositiveInt ! Generates a random positive integer -use MOM_io, only : stdout, stderr +use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit implicit none ; private @@ -23,8 +23,6 @@ module MOM_random public :: random_2d_norm public :: random_unit_tests -#include - !> Container for pseudo-random number generators type, public :: PRNG ; private @@ -63,7 +61,7 @@ end function random_norm subroutine random_2d_01(CS, HI, rand) type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators type(hor_index_type), intent(in) :: HI !< Horizontal index structure - real, dimension(SZI_(HI),SZJ_(HI)), intent(out) :: rand !< Random numbers between 0 and 1 + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(out) :: rand !< Random numbers between 0 and 1 ! Local variables integer :: i,j @@ -80,7 +78,7 @@ end subroutine random_2d_01 subroutine random_2d_norm(CS, HI, rand) type(PRNG), intent(inout) :: CS !< Container for pseudo-random number generators type(hor_index_type), intent(in) :: HI !< Horizontal index structure - real, dimension(SZI_(HI),SZJ_(HI)), intent(out) :: rand !< Random numbers between 0 and 1 + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), intent(out) :: rand !< Random numbers between 0 and 1 ! Local variables integer :: i,j,n From aed5a680babcf4123c79df1fb4dd4758dbd79ef8 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 3 Jan 2021 10:28:46 -0500 Subject: [PATCH 17/18] +Add the new routine read_field_chksum to MOM_io Added the new routine read_field_chksum to MOM_io.F90, so that all calls to the FMS i/o layer can be directed via MOM_io.F90, in order to facilitate the painless and compartmentalized migration to FMS2. Also added a 0-d variant for MOM_read_data, and standardized the control-flag and subroutine aliases used in MOM_io.F90. All answers are bitwise identical, but there are new public interfaces. --- src/framework/MOM_io.F90 | 78 +++++++++++++++++++++++++++++----------- 1 file changed, 57 insertions(+), 21 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index d13dddc3c7..529c725274 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -16,19 +16,19 @@ module MOM_io use ensemble_manager_mod, only : get_ensemble_id use fms_mod, only : write_version_number, open_namelist_file, check_nml_error use fms_io_mod, only : file_exist, field_size, read_data -use fms_io_mod, only : field_exists => field_exist, io_infra_end=>fms_io_exit -use fms_io_mod, only : get_filename_appendix => get_filename_appendix +use fms_io_mod, only : field_exists=>field_exist, io_infra_end=>fms_io_exit +use fms_io_mod, only : get_filename_appendix=>get_filename_appendix use mpp_domains_mod, only : domain1d, domain2d, mpp_get_domain_components use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST use mpp_io_mod, only : open_file => mpp_open, close_file => mpp_close -use mpp_io_mod, only : mpp_write_meta, write_field => mpp_write, mpp_get_info -use mpp_io_mod, only : mpp_get_atts, mpp_get_axes, get_axis_data=>mpp_get_axis_data, axistype -use mpp_io_mod, only : mpp_get_fields, fieldtype, axistype, flush_file => mpp_flush +use mpp_io_mod, only : mpp_write_meta, write_field => mpp_write +use mpp_io_mod, only : mpp_get_atts, mpp_attribute_exist +use mpp_io_mod, only : mpp_get_axes, axistype, get_axis_data=>mpp_get_axis_data +use mpp_io_mod, only : mpp_get_fields, fieldtype, flush_file=>mpp_flush use mpp_io_mod, only : APPEND_FILE=>MPP_APPEND, ASCII_FILE=>MPP_ASCII use mpp_io_mod, only : MULTIPLE=>MPP_MULTI, NETCDF_FILE=>MPP_NETCDF use mpp_io_mod, only : OVERWRITE_FILE=>MPP_OVERWR, READONLY_FILE=>MPP_RDONLY use mpp_io_mod, only : SINGLE_FILE=>MPP_SINGLE, WRITEONLY_FILE=>MPP_WRONLY -use mpp_io_mod, only : MPP_APPEND, MPP_MULTI, MPP_OVERWR, MPP_NETCDF, MPP_RDONLY use mpp_io_mod, only : get_file_info=>mpp_get_info, get_file_atts=>mpp_get_atts use mpp_io_mod, only : get_file_fields=>mpp_get_fields, get_file_times=>mpp_get_times use mpp_io_mod, only : io_infra_init=>mpp_io_init @@ -40,7 +40,7 @@ module MOM_io public :: close_file, create_file, field_exists, field_size, fieldtype, get_filename_appendix public :: file_exists, flush_file, get_file_info, get_file_atts, get_file_fields -public :: get_file_times, open_file, read_axis_data, read_data +public :: get_file_times, open_file, read_axis_data, read_data, read_field_chksum public :: num_timelevels, MOM_read_data, MOM_read_vector, ensembler public :: reopen_file, slasher, write_field, write_version_number, MOM_io_init public :: open_namelist_file, check_nml_error, io_infra_init, io_infra_end @@ -77,6 +77,7 @@ module MOM_io module procedure MOM_read_data_3d module procedure MOM_read_data_2d module procedure MOM_read_data_1d + module procedure MOM_read_data_0d end interface !> Read a pair of data fields representing the two components of a vector from a file @@ -162,9 +163,9 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then - call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, threading=thread) + call open_file(unit, filename, OVERWRITE_FILE, NETCDF_FILE, threading=thread) else - call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, domain=Domain%mpp_domain) + call open_file(unit, filename, OVERWRITE_FILE, NETCDF_FILE, domain=Domain%mpp_domain) endif ! Define the coordinates. @@ -404,13 +405,13 @@ subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then - call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, threading=thread) + call open_file(unit, filename, APPEND_FILE, NETCDF_FILE, threading=thread) else - call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, domain=Domain%mpp_domain) + call open_file(unit, filename, APPEND_FILE, NETCDF_FILE, domain=Domain%mpp_domain) endif if (unit < 0) return - call mpp_get_info(unit, ndim, nvar, natt, ntime) + call get_file_info(unit, ndim, nvar, natt, ntime) if (nvar == -1) then write (mesg,*) "Reopening file ",trim(filename)," apparently had ",nvar,& @@ -449,11 +450,11 @@ subroutine read_axis_data(filename, axis_name, var) type(axistype) :: time_axis character(len=32) :: name, units - call open_file(unit, trim(filename), action=MPP_RDONLY, form=MPP_NETCDF, & - threading=MPP_MULTI, fileset=SINGLE_FILE) + call open_file(unit, trim(filename), action=READONLY_FILE, form=NETCDF_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) !Find the number of variables (nvar) in this file - call mpp_get_info(unit, ndim, nvar, natt, ntime) + call get_file_info(unit, ndim, nvar, natt, ntime) ! ------------------------------------------------------------------- ! Allocate space for the number of axes in the data file. ! ------------------------------------------------------------------- @@ -462,7 +463,7 @@ subroutine read_axis_data(filename, axis_name, var) axis_found = .false. do i = 1, ndim - call mpp_get_atts(axes(i), name=name,len=len,units=units) + call get_file_atts(axes(i), name=name, len=len, units=units) if (name == axis_name) then axis_found = .true. call get_axis_data(axes(i),var) @@ -477,6 +478,23 @@ subroutine read_axis_data(filename, axis_name, var) end subroutine read_axis_data +subroutine read_field_chksum(field, chksum, valid_chksum) + type(fieldtype), intent(in) :: field !< The field whose checksum attribute is to be read. + integer(kind=8), intent(out) :: chksum !< The checksum for the field. + logical, intent(out) :: valid_chksum !< If true, chksum has been successfully read. + ! Local variables + integer(kind=8), dimension(3) :: checksum_file + + checksum_file(:) = -1 + valid_chksum = mpp_attribute_exist(field, "checksum") + if (valid_chksum) then + call mpp_get_atts(field, checksum=checksum_file) + chksum = checksum_file(1) + else + chksum = -1 + endif +end subroutine read_field_chksum + !> This function determines how many time levels a variable has. function num_timelevels(filename, varname, min_dims) result(n_time) character(len=*), intent(in) :: filename !< name of the file to read @@ -519,7 +537,6 @@ function num_timelevels(filename, varname, min_dims) result(n_time) return endif - allocate(varids(nvars)) status = nf90_inq_varids(ncid, nvars, varids) @@ -848,7 +865,26 @@ function FMS_file_exists(filename, domain, no_domain) end function FMS_file_exists -!> This function uses the fms_io function read_data to read 1-D + +!> This function uses the fms_io function read_data to read a scalar +!! data field named "fieldname" from file "filename". +subroutine MOM_read_data_0d(filename, fieldname, data, timelevel, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, intent(inout) :: data !< The 1-dimensional array into which the data + integer, optional, intent(in) :: timelevel !< The time level in the file to read + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) + + if (present(scale)) then ; if (scale /= 1.0) then + data = scale*data + endif ; endif + +end subroutine MOM_read_data_0d + +!> This function uses the fms_io function read_data to read a 1-D !! data field named "fieldname" from file "filename". subroutine MOM_read_data_1d(filename, fieldname, data, timelevel, scale) character(len=*), intent(in) :: filename !< The name of the file to read @@ -879,7 +915,7 @@ subroutine MOM_read_data_2d(filename, fieldname, data, MOM_Domain, & integer, optional, intent(in) :: timelevel !< The time level in the file to read integer, optional, intent(in) :: position !< A flag indicating where this data is located real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before they are returned. + !! by before it is returned. integer :: is, ie, js, je @@ -907,7 +943,7 @@ subroutine MOM_read_data_3d(filename, fieldname, data, MOM_Domain, & integer, optional, intent(in) :: timelevel !< The time level in the file to read integer, optional, intent(in) :: position !< A flag indicating where this data is located real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before they are returned. + !! by before it is returned. integer :: is, ie, js, je @@ -935,7 +971,7 @@ subroutine MOM_read_data_4d(filename, fieldname, data, MOM_Domain, & integer, optional, intent(in) :: timelevel !< The time level in the file to read integer, optional, intent(in) :: position !< A flag indicating where this data is located real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before they are returned. + !! by before it is returned. integer :: is, ie, js, je From aea16f7302d52fed890b8c4e79876029a4ebbb5c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 3 Jan 2021 11:43:23 -0500 Subject: [PATCH 18/18] Use read_field_chksum in MOM_restart Use read_field_chksum and MOM_read_data in MOM_restart. Also internally renamed mpp_chksum to just chksum in MOM_restart to aid in identifying unfiltered dependencies on FMS in MOM_restart.F90. All answers are bitwise identical. --- src/framework/MOM_restart.F90 | 71 +++++++++++++++-------------------- 1 file changed, 31 insertions(+), 40 deletions(-) diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 7181a1f1b9..6e4e3d745b 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -3,24 +3,22 @@ module MOM_restart ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_domains, only : pe_here, num_PEs +use MOM_domains, only : PE_here, num_PEs use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_string_functions, only : lowercase use MOM_grid, only : ocean_grid_type use MOM_io, only : create_file, fieldtype, file_exists, open_file, close_file -use MOM_io, only : MOM_read_data, read_data, get_filename_appendix +use MOM_io, only : MOM_read_data, read_data, get_filename_appendix, read_field_chksum use MOM_io, only : get_file_info, get_file_atts, get_file_fields, get_file_times use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc use MOM_io, only : MULTIPLE, NETCDF_FILE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE -use MOM_time_manager, only : time_type, time_type_to_real, real_to_time -use MOM_time_manager, only : days_in_month, get_date, set_date -use MOM_transform_FMS, only : mpp_chksum => rotated_mpp_chksum +use MOM_time_manager, only : time_type, time_type_to_real, real_to_time +use MOM_time_manager, only : days_in_month, get_date, set_date +use MOM_transform_FMS, only : chksum => rotated_mpp_chksum use MOM_transform_FMS, only : write_field => rotated_write_field -use MOM_verticalGrid, only : verticalGrid_type -use mpp_io_mod, only : mpp_attribute_exist, mpp_get_atts -use mpp_mod, only : mpp_pe +use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -116,7 +114,7 @@ module MOM_restart module procedure register_restart_field_ptr0d, register_restart_field_0d end interface -!> Register a pair of restart fieilds whose rotations map onto each other +!> Register a pair of restart fields whose rotations map onto each other interface register_restart_pair module procedure register_restart_pair_ptr2d module procedure register_restart_pair_ptr3d @@ -1010,18 +1008,15 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ endif do m=start_var,next_var-1 if (associated(CS%var_ptr3d(m)%p)) then - check_val(m-start_var+1,1) = & - mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:), turns=-turns) + check_val(m-start_var+1,1) = chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:), turns=-turns) elseif (associated(CS%var_ptr2d(m)%p)) then - check_val(m-start_var+1,1) = & - mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL), turns=-turns) + check_val(m-start_var+1,1) = chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL), turns=-turns) elseif (associated(CS%var_ptr4d(m)%p)) then - check_val(m-start_var+1,1) = & - mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:), turns=-turns) + check_val(m-start_var+1,1) = chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:), turns=-turns) elseif (associated(CS%var_ptr1d(m)%p)) then - check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr1d(m)%p) + check_val(m-start_var+1,1) = chksum(CS%var_ptr1d(m)%p) elseif (associated(CS%var_ptr0d(m)%p)) then - check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr0d(m)%p,pelist=(/mpp_pe()/)) + check_val(m-start_var+1,1) = chksum(CS%var_ptr0d(m)%p, pelist=(/PE_here()/)) endif enddo @@ -1100,9 +1095,9 @@ subroutine restore_state(filename, directory, day, G, CS) real :: t1, t2 ! Two times. real, allocatable :: time_vals(:) type(fieldtype), allocatable :: fields(:) - logical :: check_exist, is_there_a_checksum - integer(kind=8),dimension(3) :: checksum_file - integer(kind=8) :: checksum_data + logical :: is_there_a_checksum ! Is there a valid checksum that should be checked. + integer(kind=8) :: checksum_file ! The checksum value recorded in the input file. + integer(kind=8) :: checksum_data ! The checksum value for the data that was read in. if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "restore_state: Module must be initialized before it is used.") @@ -1202,25 +1197,21 @@ subroutine restore_state(filename, directory, day, G, CS) do i=1, nvar call get_file_atts(fields(i),name=varname) if (lowercase(trim(varname)) == lowercase(trim(CS%restart_field(m)%var_name))) then - check_exist = mpp_attribute_exist(fields(i),"checksum") - checksum_file(:) = -1 checksum_data = -1 - is_there_a_checksum = .false. - if ( check_exist ) then - call mpp_get_atts(fields(i),checksum=checksum_file) - is_there_a_checksum = .true. + if (CS%checksum_required) then + call read_field_chksum(fields(i), checksum_file, is_there_a_checksum) + else + checksum_file = -1 + is_there_a_checksum = .false. ! Do not need to do data checksumming. endif - if (.NOT. CS%checksum_required) is_there_a_checksum = .false. ! Do not need to do data checksumming. if (associated(CS%var_ptr1d(m)%p)) then ! Read a 1d array, which should be invariant to domain decomposition. - call read_data(unit_path(n), varname, CS%var_ptr1d(m)%p, & - G%Domain%mpp_domain, timelevel=1) - if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr1d(m)%p) + call MOM_read_data(unit_path(n), varname, CS%var_ptr1d(m)%p, timelevel=1) + if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr1d(m)%p) elseif (associated(CS%var_ptr0d(m)%p)) then ! Read a scalar... - call read_data(unit_path(n), varname, CS%var_ptr0d(m)%p, & - G%Domain%mpp_domain, timelevel=1) - if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr0d(m)%p,pelist=(/mpp_pe()/)) + call MOM_read_data(unit_path(n), varname, CS%var_ptr0d(m)%p, timelevel=1) + if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr0d(m)%p, pelist=(/PE_here()/)) elseif (associated(CS%var_ptr2d(m)%p)) then ! Read a 2d array. if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr2d(m)%p, & @@ -1229,7 +1220,7 @@ subroutine restore_state(filename, directory, day, G, CS) call read_data(unit_path(n), varname, CS%var_ptr2d(m)%p, & no_domain=.true., timelevel=1) endif - if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) + if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) elseif (associated(CS%var_ptr3d(m)%p)) then ! Read a 3d array. if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, & @@ -1238,7 +1229,7 @@ subroutine restore_state(filename, directory, day, G, CS) call read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, & no_domain=.true., timelevel=1) endif - if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) + if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) elseif (associated(CS%var_ptr4d(m)%p)) then ! Read a 4d array. if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & @@ -1247,14 +1238,14 @@ subroutine restore_state(filename, directory, day, G, CS) call read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & no_domain=.true., timelevel=1) endif - if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) + if (is_there_a_checksum) checksum_data = chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) else call MOM_error(FATAL, "MOM_restart restore_state: No pointers set for "//trim(varname)) endif - if (is_root_pe() .and. is_there_a_checksum .and. (checksum_file(1) /= checksum_data)) then + if (is_root_pe() .and. is_there_a_checksum .and. (checksum_file /= checksum_data)) then write (mesg,'(a,Z16,a,Z16,a)') "Checksum of input field "// trim(varname)//" ",checksum_data,& - " does not match value ", checksum_file(1), & + " does not match value ", checksum_file, & " stored in "//trim(unit_path(n)//"." ) call MOM_error(FATAL, "MOM_restart(restore_state): "//trim(mesg) ) endif @@ -1455,7 +1446,7 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & if (fexists) then if (present(units)) & call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & - threading = MULTIPLE, fileset = SINGLE_FILE) + threading=MULTIPLE, fileset=SINGLE_FILE) if (present(global_files)) global_files(n) = .true. elseif (CS%parallel_restartfiles) then ! Look for decomposed files using the I/O Layout. @@ -1484,7 +1475,7 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & if (fexists) then if (present(units)) & call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & - threading = MULTIPLE, fileset = SINGLE_FILE) + threading=MULTIPLE, fileset=SINGLE_FILE) if (present(global_files)) global_files(n) = .true. if (present(file_paths)) file_paths(n) = filepath n = n + 1