Merge branch 'dev/gfdl' into better_Stokes_PGF_calcs

.github/workflows/macos-regression.yml

@@ -11,7 +11,6 @@ jobs:
       CC: gcc
       FC: gfortran
       FMS_COMMIT: 2019.01.03
-      FRAMEWORK: fms1
 
     defaults:
       run:

.github/workflows/macos-stencil.yml

@@ -11,7 +11,6 @@ jobs:
       CC: gcc
       FC: gfortran
       FMS_COMMIT: 2019.01.03
-      FRAMEWORK: fms1
 
     defaults:
       run:

.testing/Makefile

@@ -20,7 +20,6 @@
 #
 # General test configuration:
 #   MPIRUN                  MPI job launcher (mpirun, srun, etc)
-#   FRAMEWORK               Model framework (fms1 or fms2)
 #   DO_REPRO_TESTS          Enable production ("repro") testing equivalence
 #   DO_REGRESSION_TESTS     Enable regression tests (usually dev/gfdl)
 #   DO_COVERAGE             Enable code coverage and generate .gcov reports
@@ -74,8 +73,11 @@ AC_SRCDIR := $(dir $(abspath $(lastword $(MAKEFILE_LIST))))../ac
 # User-defined configuration
 -include config.mk
 
-# Set the infra framework
-FRAMEWORK ?= fms2
+# Set the FMS library
+FMS_COMMIT ?= 2023.03
+FMS_URL ?= https://github.com/NOAA-GFDL/FMS.git
+export FMS_COMMIT
+export FMS_URL
 
 # Set the MPI launcher here
 # TODO: This needs more automated configuration

.testing/README.rst

@@ -47,9 +47,11 @@ Several of the following may require configuration for particular systems.
    Name of the MPI launcher.  Often this is ``mpirun`` or ``mpiexec`` but may
    all need to run through a scheduler, e.g. ``srun`` if using Slurm.
 
-``FRAMEWORK`` (*default:* ``fms1``)
-   Select either the legacy FMS framework (``fms1``) or an FMS2 I/O compatible
-   version (``fms2``).
+``FMS_COMMIT`` (*default:* ``2023.03``)
+   Set the FMS version, either by tag or commit (as defined in ``FMS_URL``).
+
+``FMS_URL`` (*default*: ``https://github.com/NOAA-GFDL/FMS.git``)
+   Set the URL of the FMS repository.
 
 ``DO_REPRO_TESTS`` (*default:* *none*)
    Set to ``true`` to test the REPRO build and confirm equivalence of DEBUG and

config_src/drivers/FMS_cap/ocean_model_MOM.F90

@@ -55,6 +55,7 @@ module ocean_model_mod
 use MOM_ice_shelf, only : initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS
 use MOM_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces
 use MOM_ice_shelf, only : add_shelf_forces, ice_shelf_end, ice_shelf_save_restart
+use MOM_ice_shelf, only : ice_sheet_calving_to_ocean_sfc
 use MOM_wave_interface, only: wave_parameters_CS, MOM_wave_interface_init
 use MOM_wave_interface, only: Update_Surface_Waves
 use iso_fortran_env, only : int64
@@ -121,7 +122,10 @@ module ocean_model_mod
                         !! formation in the ocean.
     melt_potential => NULL(), & !< Instantaneous heat used to melt sea ice [J m-2].
     OBLD => NULL(),   & !< Ocean boundary layer depth [m].
-    area => NULL()      !< cell area of the ocean surface [m2].
+    area => NULL(),   & !< cell area of the ocean surface [m2].
+    calving => NULL(), &!< The mass per unit area of the ice shelf to convert to
+                        !! bergs [kg m-2].
+    calving_hflx => NULL() !< Calving heat flux [W m-2].
   type(coupler_2d_bc_type) :: fields    !< A structure that may contain named
                                         !! arrays of tracer-related surface fields.
   integer                  :: avg_kount !< A count of contributions to running
@@ -157,6 +161,8 @@ module ocean_model_mod
                               !! ocean dynamics and forcing fluxes.
   real :: press_to_z          !< A conversion factor between pressure and ocean depth,
                               !! usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1].
+  logical :: calve_ice_shelf_bergs = .false. !< If true, bergs are initialized according to
+                              !! ice shelf flux through the ice front
   real :: C_p                 !< The heat capacity of seawater [J degC-1 kg-1].
   logical :: offline_tracer_mode = .false. !< If false, use the model in prognostic mode
                               !! with the barotropic and baroclinic dynamics, thermodynamics,
@@ -221,7 +227,7 @@ module ocean_model_mod
 !!   This subroutine initializes both the ocean state and the ocean surface type.
 !! Because of the way that indices and domains are handled, Ocean_sfc must have
 !! been used in a previous call to initialize_ocean_type.
-subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas_fields_ocn)
+subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas_fields_ocn, calve_ice_shelf_bergs)
   type(ocean_public_type), target, &
                        intent(inout) :: Ocean_sfc !< A structure containing various publicly
                                 !! visible ocean surface properties after initialization,
@@ -239,6 +245,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
                                               !! in the calculation of additional gas or other
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
+  logical, optional,   intent(in)    :: calve_ice_shelf_bergs !< If true, track ice shelf flux through a
+                                              !! static ice shelf, so that it can be converted into icebergs
   ! Local variables
   real :: Rho0        ! The Boussinesq ocean density [R ~> kg m-3]
   real :: G_Earth     ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
@@ -247,6 +255,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
                       !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
                       !! If HFrz <= 0 (default), melt potential will not be computed.
   logical :: use_melt_pot !< If true, allocate melt_potential array
+  logical :: point_calving ! Equals calve_ice_shelf_bergs if calve_ice_shelf_bergs is present
 
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
@@ -274,11 +283,11 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
   OS%Time = Time_in ; OS%Time_dyn = Time_in
   ! Call initialize MOM with an optional Ice Shelf CS which, if present triggers
   ! initialization of ice shelf parameters and arrays.
-
+  point_calving=.false.; if (present(calve_ice_shelf_bergs)) point_calving=calve_ice_shelf_bergs
   call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, &
                       Time_in, offline_tracer_mode=OS%offline_tracer_mode, &
                       diag_ptr=OS%diag, count_calls=.true., ice_shelf_CSp=OS%ice_shelf_CSp, &
-                      waves_CSp=OS%Waves)
+                      waves_CSp=OS%Waves, calve_ice_shelf_bergs=point_calving)
   call get_MOM_state_elements(OS%MOM_CSp, G=OS%grid, GV=OS%GV, US=OS%US, C_p=OS%C_p, &
                               C_p_scaled=OS%fluxes%C_p, use_temp=use_temperature)
 
@@ -406,6 +415,13 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
 
   endif
 
+  if (present(calve_ice_shelf_bergs)) then
+    if (calve_ice_shelf_bergs) then
+      call convert_shelf_state_to_ocean_type(Ocean_sfc, OS%Ice_shelf_CSp, OS%US)
+      OS%calve_ice_shelf_bergs=.true.
+    endif
+  endif
+
   call close_param_file(param_file)
   call diag_mediator_close_registration(OS%diag)
 
@@ -668,6 +684,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
 !  call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US, &
 !                                   OS%fluxes%p_surf_full, OS%press_to_z)
   call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
+  if (OS%calve_ice_shelf_bergs) call convert_shelf_state_to_ocean_type(Ocean_sfc,OS%Ice_shelf_CSp, OS%US)
   Time1 = OS%Time ; if (do_dyn) Time1 = OS%Time_dyn
   call coupler_type_send_data(Ocean_sfc%fields, Time1)
 
@@ -789,6 +806,8 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, gas_field
              Ocean_sfc%u_surf (isc:iec,jsc:jec), &
              Ocean_sfc%v_surf (isc:iec,jsc:jec), &
              Ocean_sfc%sea_lev(isc:iec,jsc:jec), &
+             Ocean_sfc%calving(isc:iec,jsc:jec), &
+             Ocean_sfc%calving_hflx(isc:iec,jsc:jec), &
              Ocean_sfc%area   (isc:iec,jsc:jec), &
              Ocean_sfc%melt_potential(isc:iec,jsc:jec), &
              Ocean_sfc%OBLD   (isc:iec,jsc:jec), &
@@ -799,6 +818,8 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, gas_field
   Ocean_sfc%u_surf(:,:)  = 0.0  ! time averaged u-current (m/sec) passed to atmosphere/ice models
   Ocean_sfc%v_surf(:,:)  = 0.0  ! time averaged v-current (m/sec)  passed to atmosphere/ice models
   Ocean_sfc%sea_lev(:,:) = 0.0  ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav
+  Ocean_sfc%calving(:,:)  = 0.0  ! time accumulated ice sheet calving (kg m-2) passed to ice model
+  Ocean_sfc%calving_hflx(:,:) = 0.0 ! time accumulated ice sheet calving heat flux (W m-2) passed to ice model
   Ocean_sfc%frazil(:,:)  = 0.0  ! time accumulated frazil (J/m^2) passed to ice model
   Ocean_sfc%melt_potential(:,:)  = 0.0  ! time accumulated melt potential (J/m^2) passed to ice model
   Ocean_sfc%OBLD(:,:)    = 0.0  ! ocean boundary layer depth (m)
@@ -932,6 +953,24 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_
 
 end subroutine convert_state_to_ocean_type
 
+!> Converts the ice-shelf-to-ocean calving and calving_hflx variables from the ice-shelf state (ISS) type
+!! to the ocean public type
+subroutine convert_shelf_state_to_ocean_type(Ocean_sfc, CS, US)
+  type(ocean_public_type), &
+               target, intent(inout) :: Ocean_sfc !< A structure containing various publicly
+                                                  !! visible ocean surface fields, whose elements
+                                                  !! have their data set here.
+  type(ice_shelf_CS),      pointer :: CS        !< A pointer to the ice shelf control structure
+  type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
+  integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd, i, j
+
+  call get_domain_extent(Ocean_sfc%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
+
+  call ice_sheet_calving_to_ocean_sfc(CS,US,Ocean_sfc%calving(isc_bnd:iec_bnd,jsc_bnd:jec_bnd),&
+    Ocean_sfc%calving_hflx(isc_bnd:iec_bnd,jsc_bnd:jec_bnd))
+
+end subroutine convert_shelf_state_to_ocean_type
+
 !>   This subroutine extracts the surface properties from the ocean's internal
 !! state and stores them in the ocean type returned to the calling ice model.
 !! It has to be separate from the ocean_initialization call because the coupler

config_src/drivers/ice_solo_driver/ice_shelf_driver.F90

@@ -26,6 +26,7 @@ program Shelf_main
   use MOM_debugging,       only : MOM_debugging_init
   use MOM_diag_mediator,   only : diag_mediator_init, diag_mediator_infrastructure_init, set_axes_info
   use MOM_diag_mediator,   only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration
+  use MOM_diag_manager_infra, only : diag_manager_set_time_end_infra
   use MOM_domains,         only : MOM_infra_init, MOM_infra_end
   use MOM_domains,         only : MOM_domains_init, clone_MOM_domain, pass_var
   use MOM_dyn_horgrid,     only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
@@ -325,6 +326,8 @@ program Shelf_main
     Time_end = daymax
   endif
 
+  call diag_manager_set_time_end_infra (Time_end)
+
   if (Time >= Time_end) call MOM_error(FATAL, &
     "Shelf_driver: The run has been started at or after the end time of the run.")
 

config_src/drivers/nuopc_cap/mom_cap.F90

@@ -136,6 +136,7 @@ module MOM_cap_mod
 logical              :: grid_attach_area = .false.
 logical              :: use_coldstart = .true.
 logical              :: use_mommesh = .true.
+logical              :: restart_eor = .false.
 character(len=128)   :: scalar_field_name = ''
 integer              :: scalar_field_count = 0
 integer              :: scalar_field_idx_grid_nx = 0
@@ -381,6 +382,13 @@ subroutine InitializeP0(gcomp, importState, exportState, clock, rc)
     geomtype = ESMF_GEOMTYPE_GRID
   endif
 
+  ! Read end of run restart config option
+  call NUOPC_CompAttributeGet(gcomp, name="write_restart_at_endofrun", value=value, &
+                              isPresent=isPresent, isSet=isSet, rc=rc)
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  if (isPresent .and. isSet) then
+     if (trim(value) .eq. '.true.') restart_eor = .true.
+  end if
 
 end subroutine
 
@@ -1637,6 +1645,8 @@ subroutine ModelAdvance(gcomp, rc)
   real(8)                                :: MPI_Wtime, timers
   logical                                :: write_restart
   logical                                :: write_restartfh
+  logical                                :: write_restart_eor
+
 
   rc = ESMF_SUCCESS
   if(profile_memory) call ESMF_VMLogMemInfo("Entering MOM Model_ADVANCE: ")
@@ -1776,7 +1786,6 @@ subroutine ModelAdvance(gcomp, rc)
   !---------------
   ! Get the stop alarm
   !---------------
-
   call ESMF_ClockGetAlarm(clock, alarmname='stop_alarm', alarm=stop_alarm, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
@@ -1807,7 +1816,18 @@ subroutine ModelAdvance(gcomp, rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
     end if
 
-    if (write_restart .or. write_restartfh) then
+    write_restart_eor = .false.
+    if (restart_eor) then
+      if (ESMF_AlarmIsRinging(stop_alarm, rc=rc)) then
+         if (ChkErr(rc,__LINE__,u_FILE_u)) return
+         write_restart_eor = .true.
+         ! turn off the alarm
+         call ESMF_AlarmRingerOff(stop_alarm, rc=rc )
+         if (ChkErr(rc,__LINE__,u_FILE_u)) return
+       end if
+    end if
+
+    if (write_restart .or. write_restartfh .or. write_restart_eor) then
       ! determine restart filename
       call ESMF_ClockGetNextTime(clock, MyTime, rc=rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return

config_src/drivers/solo_driver/MOM_driver.F90

@@ -28,6 +28,7 @@ program MOM6
   use MOM_cpu_clock,       only : CLOCK_COMPONENT
   use MOM_data_override,   only : data_override_init
   use MOM_diag_mediator,   only : diag_mediator_end, diag_ctrl, diag_mediator_close_registration
+  use MOM_diag_manager_infra, only : diag_manager_set_time_end_infra
   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
@@ -375,6 +376,8 @@ program MOM6
     Time_end = daymax
   endif
 
+  call diag_manager_set_time_end_infra(Time_end)
+
   call get_param(param_file, mod_name, "SINGLE_STEPPING_CALL", single_step_call, &
                  "If true, advance the state of MOM with a single step "//&
                  "including both dynamics and thermodynamics.  If false "//&

config_src/infra/FMS1/MOM_diag_manager_infra.F90

@@ -57,6 +57,8 @@ module MOM_diag_manager_infra
 public MOM_diag_manager_init
 public MOM_diag_manager_end
 public send_data_infra
+public diag_send_complete_infra
+public diag_manager_set_time_end_infra
 public MOM_diag_field_add_attribute
 public register_diag_field_infra
 public register_static_field_infra
@@ -451,4 +453,13 @@ subroutine MOM_diag_field_add_attribute_i1d(diag_field_id, att_name, att_value)
 
 end subroutine MOM_diag_field_add_attribute_i1d
 
+!> Needed for backwards compatibility, does nothing
+subroutine diag_send_complete_infra ()
+end subroutine diag_send_complete_infra
+
+!> Needed for backwards compatibility, does nothing
+subroutine diag_manager_set_time_end_infra(time)
+  type(time_type), intent(in) :: time !< The model time that simulation ends
+end subroutine diag_manager_set_time_end_infra
+
 end module MOM_diag_manager_infra

config_src/infra/FMS2/MOM_diag_manager_infra.F90

@@ -14,13 +14,15 @@ module MOM_diag_manager_infra
 use diag_data_mod,    only : null_axis_id
 use diag_manager_mod, only : fms_diag_manager_init => diag_manager_init
 use diag_manager_mod, only : fms_diag_manager_end => diag_manager_end
+use diag_manager_mod, only : diag_send_complete
+use diag_manager_mod, only : diag_manager_set_time_end
 use diag_manager_mod, only : send_data_fms => send_data
 use diag_manager_mod, only : fms_diag_field_add_attribute => diag_field_add_attribute
 use diag_manager_mod, only : DIAG_FIELD_NOT_FOUND
 use diag_manager_mod, only : register_diag_field_fms => register_diag_field
 use diag_manager_mod, only : register_static_field_fms => register_static_field
 use diag_manager_mod, only : get_diag_field_id_fms => get_diag_field_id
-use MOM_time_manager, only : time_type
+use MOM_time_manager, only : time_type, set_time
 use MOM_domain_infra, only : MOM_domain_type
 use MOM_error_infra,  only : MOM_error => MOM_err, FATAL, WARNING
 
@@ -57,6 +59,8 @@ module MOM_diag_manager_infra
 public MOM_diag_manager_init
 public MOM_diag_manager_end
 public send_data_infra
+public diag_send_complete_infra
+public diag_manager_set_time_end_infra
 public MOM_diag_field_add_attribute
 public register_diag_field_infra
 public register_static_field_infra
@@ -451,4 +455,19 @@ subroutine MOM_diag_field_add_attribute_i1d(diag_field_id, att_name, att_value)
 
 end subroutine MOM_diag_field_add_attribute_i1d
 
+!> Finishes the diag manager reduction methods as needed for the time_step
+subroutine diag_send_complete_infra ()
+  !! The time_step in the diag_send_complete call is a dummy argument, needed for backwards compatibility
+  !! It won't be used at all when diag_manager_nml::use_modern_diag=.true.
+  !! It won't have any impact when diag_manager_nml::use_modern_diag=.false.
+  call diag_send_complete (set_time(0))
+end subroutine diag_send_complete_infra
+
+!> Sets the time that the simulation ends in the diag manager
+subroutine diag_manager_set_time_end_infra(time)
+  type(time_type),           optional, intent(in) :: time  !< The time the simulation ends
+
+  call diag_manager_set_time_end(time)
+end subroutine diag_manager_set_time_end_infra
+
 end module MOM_diag_manager_infra