From 5c477003fad3955721662db13b961336c085b512 Mon Sep 17 00:00:00 2001 From: Claire Yung Date: Mon, 12 Jan 2026 15:16:46 +1100 Subject: [PATCH 1/9] Allow MOM6 ice shelves to be initialised with NUOPC - turn FATAL error when ESMF mask and MOM6 mask differ intol WARNING, so that we can set up a MOM6 domain including ice shelves but avoid ESMF and the mediator transferrring atmospheric data into the ice shelves - use ocean_state%fluxes%frac_shelf_h when computing the So_omask for the mediator from MOM6 data, so that the coupler does not get sent fields in regions it does not expect. A ceiling is used for frac_shelf_h (fraction of ice shelf values) - partial ice shelf coverage may need to be better thought out. - for above to work, the ocean_state_type is made public. - Also initialise ice shelf fluxes and forces separately. Otherwise, shelf_sfc_mass_flux field is initialised which should only be used with dynamic ice shelf (and results in crash when it can't find the file if you're running from a restart.) --- config_src/drivers/nuopc_cap/mom_cap.F90 | 4 ++-- config_src/drivers/nuopc_cap/mom_cap_methods.F90 | 2 +- config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 | 7 +++++-- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 8d04ba19c7..0d20b1b29c 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -22,7 +22,7 @@ module MOM_cap_mod use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file use MOM_get_input, only: get_MOM_input, directories use MOM_domains, only: pass_var, pe_here -use MOM_error_handler, only: MOM_error, FATAL, is_root_pe +use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING use MOM_grid, only: ocean_grid_type, get_global_grid_size use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type @@ -1383,7 +1383,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) frmt = "('ERROR: ESMF mesh and MOM6 domain masks are inconsistent! - "//& "MOM n, maskMesh(n), mask(n) = ',3(i8,2x))" write(err_msg, frmt)n,maskMesh(n),mask(n) - call MOM_error(FATAL, err_msg) + call MOM_error(WARNING, err_msg) endif end do diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 1accf7472a..11d1f88390 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -696,7 +696,7 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, jg = j + ocean_grid%jsc - jsc do i = isc, iec ig = i + ocean_grid%isc - isc - omask(i,j) = nint(ocean_grid%mask2dT(ig,jg)) + omask(i,j) = nint(ocean_grid%mask2dT(ig,jg)) - ceiling(ocean_state%fluxes%frac_shelf_h(ig,jg)) enddo enddo diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index b8bb87bff3..bd25fdc5c0 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -50,6 +50,7 @@ module MOM_ocean_model_nuopc 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_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces use MOM_coupler_types, only : coupler_1d_bc_type, coupler_2d_bc_type use MOM_coupler_types, only : coupler_type_spawn, coupler_type_write_chksums use MOM_coupler_types, only : coupler_type_initialized, coupler_type_copy_data @@ -134,7 +135,7 @@ module MOM_ocean_model_nuopc !> The ocean_state_type contains all information about the state of the ocean, !! with a format that is private so it can be readily changed without disrupting !! other coupled components. -type, public :: ocean_state_type ; private +type, public :: ocean_state_type ! This type is private, and can therefore vary between different ocean models. logical :: is_ocean_PE = .false. !< True if this is an ocean PE. type(time_type) :: Time !< The ocean model's time and master clock. @@ -395,7 +396,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i if (OS%use_ice_shelf) then call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & - OS%diag, Time_init, OS%dirs%output_directory, OS%forces, OS%fluxes) + OS%diag, Time_init, OS%dirs%output_directory) + call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%fluxes) + call initialize_ice_shelf_forces(OS%ice_shelf_CSp, OS%grid, OS%US, OS%forces) endif if (OS%icebergs_alter_ocean) then call marine_ice_init(OS%Time, OS%grid, param_file, OS%diag, OS%marine_ice_CSp) From a38ebbe94dd8f02fbd1fda32fc66542ea3a7fc32 Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Fri, 14 Nov 2025 10:22:48 +1100 Subject: [PATCH 2/9] Initialise OS%flux_tmp in nuopc cap Avoids segmentation fault https://github.com/ACCESS-NRI/access-om3-configs/issues/863 --- config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index bd25fdc5c0..5d6009e231 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -398,6 +398,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & OS%diag, Time_init, OS%dirs%output_directory) call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%fluxes) + call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%flux_tmp) call initialize_ice_shelf_forces(OS%ice_shelf_CSp, OS%grid, OS%US, OS%forces) endif if (OS%icebergs_alter_ocean) then From f753ce6bf475772a485f6f4df0d391c6ca4caa8b Mon Sep 17 00:00:00 2001 From: Claire Yung Date: Wed, 17 Sep 2025 17:17:22 +1000 Subject: [PATCH 3/9] Hardcoded name changes to allow ice shelf to restart with NUOPC. See also https://github.com/ACCESS-NRI/MOM6/issues/29 Hardcode ice shelf restart to be in 'RESTART/' in nuopc cap. Also hardcodes the name of the ice shelf restart Shelf.res.nc in initialize_ice_shelf. --- config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 | 8 ++++---- src/ice_shelf/MOM_ice_shelf.F90 | 3 ++- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 5d6009e231..878848260d 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -744,7 +744,7 @@ subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, nu OS%dirs%restart_output_dir) ! Is this needed? if (OS%use_ice_shelf) then call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, & - OS%dirs%restart_output_dir) + './RESTART/') !OS%dirs%restart_output_dir) endif else if (BTEST(OS%Restart_control,1)) then @@ -753,7 +753,7 @@ subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, nu call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir, time_stamped=.true.) if (OS%use_ice_shelf) then - call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir, .true.) + call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, './RESTART/', .true.) ! OS%dirs%restart_output_dir, .true.) endif endif if (BTEST(OS%Restart_control,0)) then @@ -762,7 +762,7 @@ subroutine ocean_model_restart(OS, timestamp, restartname, stoch_restartname, nu call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & OS%dirs%restart_output_dir) if (OS%use_ice_shelf) then - call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir) + call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, './RESTART/') !OS%dirs%restart_output_dir) endif endif endif @@ -823,7 +823,7 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix) call forcing_save_restart(OS%forcing_CSp, OS%grid, Time, restart_dir) if (OS%use_ice_shelf) then - call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir) + call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, './RESTART/') ! OS%dirs%restart_output_dir) endif end subroutine ocean_model_save_restart diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index bdde9ba6e4..aa99722167 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1906,7 +1906,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, elseif (.not.new_sim) then ! This line calls a subroutine that reads the initial conditions from a restart file. call MOM_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.") - call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, CS%restart_CSp) + ! call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, CS%restart_CSp) + call restore_state('Shelf.res.nc', dirs%restart_input_dir, Time, G, CS%restart_CSp) endif ! .not. new_sim From da84b0feecbbe4e2e8a3d3d892989d04c6abdff0 Mon Sep 17 00:00:00 2001 From: Claire Yung Date: Sat, 11 Oct 2025 09:20:30 +1100 Subject: [PATCH 4/9] Add parameter ICE_SHELF_TIDEAMP_SCALING_FACTOR This parameter scales the friction velocity computed due to prescribed tides either through UTIDE = constant value or through TIDEAMP = True and associated TIDEAMP_FILE. Jourdain et al. 2019 suggest a value of 0.66 to estimate tidal ice shelf boundary layer velocity from total barotropic tide amplitude. Default is set to 1 so should not change answers --- src/ice_shelf/MOM_ice_shelf.F90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index aa99722167..6c7b3308a6 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -203,6 +203,7 @@ module MOM_ice_shelf logical :: buoy_flux_itt_bug !< If true, fixes buoyancy iteration bug logical :: salt_flux_itt_bug !< If true, fixes salt iteration bug real :: buoy_flux_itt_threshold !< Buoyancy iteration threshold for convergence + real :: tideamp_scaling_factor !< Tideamp scaling factor in friction velocity calculation !>@{ Diagnostic handles integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & @@ -462,14 +463,17 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) 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))) + sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%tideamp_scaling_factor * & + (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)) + sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%tideamp_scaling_factor * & + (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))) + sqrt(CS%cdrag*((u2_av + v2_av) + CS%tideamp_scaling_factor * & + (CS%utide(i,j)**2)))) endif else ! There is no shelf here. fluxes%ustar_shelf(i,j) = 0.0 @@ -1712,6 +1716,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, call safe_alloc_ptr(CS%utide,isd,ied,jsd,jed) ; CS%utide(:,:) = 0.0 + call get_param(param_file, mdl, "ICE_SHELF_TIDEAMP_SCALING_FACTOR", CS%tideamp_scaling_factor, & + "Scale TIDEAMP_FILE or UTIDE by number in melt parameterisation "//& + "friction velocity calculation.", units = "nondim", default=1.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.", & From 3daa613c16db2b1bd6204f345a83e464d7235caf Mon Sep 17 00:00:00 2001 From: Claire Yung Date: Mon, 12 Jan 2026 15:40:36 +1100 Subject: [PATCH 5/9] Add more detail to namelist documentation of tidal amplitude Parameters READ_TIDEAMP and UTIDE not only affect INT_TIDE_DISSIPATION but also ice shelf melt rates via the ice shelf basal melt rate friction velocity calculation. In this commit the effect on ice shelves is added to the description of the parameters to make this clearer. --- src/ice_shelf/MOM_ice_shelf.F90 | 7 +++++-- .../vertical/MOM_internal_tide_input.F90 | 9 ++++++--- src/parameterizations/vertical/MOM_tidal_mixing.F90 | 9 ++++++--- 3 files changed, 17 insertions(+), 8 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 6c7b3308a6..163bf4c390 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1681,7 +1681,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, call get_param(param_file, mdl, "READ_TIDEAMP", read_TIDEAMP, & "If true, read a file (given by TIDEAMP_FILE) containing "//& - "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) + "the tidal amplitude with INT_TIDE_DISSIPATION. If true, "//& + "also used as the tidal amplitude in the ice shelf melt "//& + "parameterisation.", default=.false.) call get_param(param_file, mdl, "ICE_SHELF_LINEAR_SHELF_FRAC", CS%Zeta_N, & "Ratio of HJ99 stability constant xi_N (ratio of maximum "//& "mixing length to planetary boundary layer depth in "//& @@ -1739,7 +1741,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, endif else call get_param(param_file, mdl, "UTIDE", utide, & - "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & + "The constant tidal amplitude used with INT_TIDE_DISSIPATION "//& + "and the ice shelf melt parameterisation.", & units="m s-1", default=0.0 , scale=US%m_s_to_L_T) CS%utide(:,:) = utide endif diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index a03dca73a8..35711feda2 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -457,7 +457,8 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) call get_param(param_file, mdl, "UTIDE", utide, & - "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & + "The constant tidal amplitude used with INT_TIDE_DISSIPATION "//& + "and the ice shelf melt parameterisation.", & units="m s-1", default=0.0, scale=US%m_s_to_L_T) call read_param(param_file, "INTERNAL_TIDE_FREQS", num_freq) @@ -483,9 +484,11 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) "above the bottom boundary layer with INT_TIDE_DISSIPATION.", & units="W m-2", default=1.0e3, scale=W_m2_to_HZ2_T3) - call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, & + call get_param(param_file, mdl, "READ_TIDEAMP", read_TIDEAMP, & "If true, read a file (given by TIDEAMP_FILE) containing "//& - "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) + "the tidal amplitude with INT_TIDE_DISSIPATION. If true, "//& + "also used as the tidal amplitude in the ice shelf melt "//& + "parameterisation.", default=.false.) if (read_tideamp) then call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, & "The path to the file containing the spatially varying "//& diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index e596a9af2f..236f297420 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -455,7 +455,8 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di units="m-1", default=8.e-4*atan(1.0), scale=US%Z_to_m) call get_param(param_file, mdl, "UTIDE", CS%utide, & - "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & + "The constant tidal amplitude used with INT_TIDE_DISSIPATION "//& + "and the ice shelf melt parameterisation.", & units="m s-1", default=0.0, scale=US%m_to_Z*US%T_to_s) allocate(CS%tideamp(is:ie,js:je), source=CS%utide) @@ -467,9 +468,11 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di "above the bottom boundary layer with INT_TIDE_DISSIPATION.", & units="W m-2", default=1.0e3, scale=US%W_m2_to_RZ3_T3) - call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, & + call get_param(param_file, mdl, "READ_TIDEAMP", read_TIDEAMP, & "If true, read a file (given by TIDEAMP_FILE) containing "//& - "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) + "the tidal amplitude with INT_TIDE_DISSIPATION. If true, "//& + "also used as the tidal amplitude in the ice shelf melt "//& + "parameterisation.", default=.false.) if (read_tideamp) then if (CS%use_CVMix_tidal) then call MOM_error(FATAL, "tidal_mixing_init: Tidal amplitude files are "// & From 9573ede7c7647c811d7e1d5bacc15a265e6b77f1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 1 Oct 2025 10:48:32 -0400 Subject: [PATCH 6/9] +*Fix 3-equation ice-ocean flux iteration (#972) Fix the 3-equation iteration for the buoyancy flux between the ocean and an overlying ice-shelf when ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX is true and SHELF_3EQ_GAMMA it false. This code now uses proper bounding of the self-consistent solution, avoiding further amplifying the fluxes in the cases when the differences between the diffusivities of heat and salt to make the buoyancy flux destabilizing for finite turbulent mixing. Both the false-position iterations and the (appropriately chosen) Newton's method iterations have been extensively examined and determined to be working correctly via print statements that have subsequently been removed for efficiency. Previously, the code to determine the 3-equation solution for the buoyancy flux between the ocean and an ice shelf had been skipping iteration altogether or doing un-bounded Newton's method iterations with a sign error in part of the derivative, including taking the square root of negative numbers, leading to the issue described at https://github.com/NOAA-GFDL/MOM6/issues/945. That issue has now been corrected and can be closed once this commit has been merged into the dev/gfdl branch of MOM6. This commit also changes the names of the runtime parameters to correct the ice shelf flux iteration bugs from ICE_SHELF_BUOYANCY_FLUX_ITT_BUG and ICE_SHELF_SALT_FLUX_ITT_BUG to ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX and ICE_SHELF_SALT_FLUX_ITT_BUGFIX to avoid confusion with other ..._BUG parameters where `true` is to turn the bugs on, whereas here `true` fixes them. The old names are retained via `old_name` arguments to the `get_param()` calls, so no existing configurations will be disrupted by these changes. Additionally, an expression to determine a scaling factor to limit ice-shelf bottom slopes in `calc_shelf_driving_stress()` was refactored to avoid the possibility of division by zero. This commit will change (and correct) answers for cases with ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX set to true, but as these would often fail with a NaN from taking the square root of a negative value, it is very unlikely that any such configurations are actively being used, and there seems little point in retaining the previous answers. No answers are changed in cases that do not use an active ice shelf. Co-authored-by: Alistair Adcroft --- src/ice_shelf/MOM_ice_shelf.F90 | 274 +++++++++++++++-------- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 2 +- 2 files changed, 186 insertions(+), 90 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 163bf4c390..971fc5005e 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -129,9 +129,9 @@ module MOM_ice_shelf real :: kd_molec_salt!< The molecular diffusivity of salt [Z2 T-1 ~> m2 s-1]. real :: kd_molec_temp!< The molecular diffusivity of heat [Z2 T-1 ~> m2 s-1]. real :: Lat_fusion !< The latent heat of fusion [Q ~> J kg-1]. - real :: Gamma_T_3EQ !< Nondimensional heat-transfer coefficient, used in the 3Eq. formulation - real :: Gamma_S_3EQ !< Nondimensional salt-transfer coefficient, used in the 3Eq. formulation - !< This number should be specified by the user. + real :: Gamma_T_3EQ !< Nondimensional heat-transfer coefficient, used in the 3Eq. formulation [nondim] + real :: Gamma_S_3EQ !< Nondimensional salt-transfer coefficient, used in the 3Eq. formulation [nondim] + !< This number should be specified by the user. 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 @@ -197,12 +197,12 @@ module MOM_ice_shelf real :: dTFr_dp !< Partial derivative of freezing temperature with !! pressure [C T2 R-1 L-2 ~> degC Pa-1] real :: Zeta_N !< The stability constant xi_N = 0.052 from Holland & Jenkins '99 - !! divided by the von Karman constant VK. Was 1/8. - real :: Vk !< Von Karman's constant - dimensionless - real :: Rc !< critical flux Richardson number. - logical :: buoy_flux_itt_bug !< If true, fixes buoyancy iteration bug - logical :: salt_flux_itt_bug !< If true, fixes salt iteration bug - real :: buoy_flux_itt_threshold !< Buoyancy iteration threshold for convergence + !! divided by the von Karman constant VK [nondim]. Was 1/8. + real :: Vk !< Von Karman's constant [nondim] + real :: Rc !< critical flux Richardson number [nondim] + logical :: buoy_flux_itt_bugfix !< If true, fixes buoyancy iteration bug + logical :: salt_flux_itt_bugfix !< If true, fixes salt iteration bug + real :: buoy_flux_tol !< Fractional buoyancy iteration tolerance for convergence [nondim] real :: tideamp_scaling_factor !< Tideamp scaling factor in friction velocity calculation !>@{ Diagnostic handles @@ -298,11 +298,11 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) !! This is computed as part of the ISOMIP diagnostics. real :: time_step !< Length of time over which these fluxes will be applied [T ~> s]. real :: Itime_step !< Inverse of the length of time over which these fluxes will be applied [T-1 ~> s-1] - real :: VK !< Von Karman's constant - dimensionless + real :: VK !< Von Karman's constant [nondim] real :: ZETA_N !< This is the stability constant xi_N = 0.052 from Holland & Jenkins '99 !! divided by the von Karman constant VK. Was 1/8. [nondim] - real :: RC !< critical flux Richardson number. - real :: I_ZETA_N !< The inverse of ZETA_N [nondim]. + real :: Rf_crit !< critical flux Richardson number [nondim] + real :: I_2Zeta_N !< Half the inverse of Zeta_N [nondim]. real :: I_LF !< The inverse of the latent heat of fusion [Q-1 ~> kg J-1]. real :: I_VK !< The inverse of the Von Karman constant [nondim]. real :: PR, SC !< The Prandtl number and Schmidt number [nondim]. @@ -322,7 +322,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) real :: wB_flux !< The downward vertical flux of buoyancy just inside the ocean [Z2 T-3 ~> m2 s-3]. real :: dB_dS !< The derivative of buoyancy with salinity [Z T-2 S-1 ~> m s-2 ppt-1]. real :: dB_dT !< The derivative of buoyancy with temperature [Z T-2 C-1 ~> m s-2 degC-1]. - real :: I_n_star ! [nondim] + real :: I_n_star ! The inverse of the ratio of working boundary layer thickness + ! to the neutral thickness [nondim] real :: n_star_term ! A term in the expression for nstar [T3 Z-2 ~> s3 m-2] real :: absf ! The absolute value of the Coriolis parameter [T-1 ~> s-1] real :: dIns_dwB !< The partial derivative of I_n_star with wB_flux, in [T3 Z-2 ~> s3 m-2] @@ -331,34 +332,42 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) real :: dS_ustar ! The difference between the salinity at the ice-ocean interface and the ocean ! boundary layer salinity times the friction velocity [S Z T-1 ~> ppt m s-1] real :: ustar_h ! The friction velocity in the water below the ice shelf [Z T-1 ~> m s-1] - real :: Gam_turb ! [nondim] + real :: Gam_turb ! A relative turbluent diffusivity [nondim] real :: Gam_mol_t, Gam_mol_s ! Relative coefficients of molecular diffusivities [nondim] real :: RhoCp ! A typical ocean density times the heat capacity of water [Q R C-1 ~> J m-3 degC-1] - real :: ln_neut + real :: ln_neut ! The log of the ratio of the neutral boundary layer thickness to the molecular + ! boundary layer thickness if it is greater than 1 or 0 otherwise [nondim] real :: mass_exch ! A mass exchange rate [R Z T-1 ~> kg m-2 s-1] real :: Sb_min, Sb_max ! Minimum and maximum boundary salinities [S ~> ppt] real :: dS_min, dS_max ! Minimum and maximum salinity changes [S ~> ppt] ! Variables used in iterating for wB_flux. - real :: wB_flux_new, dDwB_dwB_in - real :: I_Gam_T, I_Gam_S - real :: dG_dwB ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2] + real :: wB_flux_next ! The next interation's guess for wB_flux [Z2 T-3 ~> m2 s-2] + real :: wB_flux_new ! An updated value of wB_flux when Gam_turb is based on wB_flux [Z2 T-3 ~> m2 s-2] + real :: wB_flux_max ! The upper bound on wB_flux [Z2 T-3 ~> m2 s-2] + real :: wB_flux_min ! The lower bound on wB_flux [Z2 T-3 ~> m2 s-2] + real :: dDwB_dwB ! The slope of the change in wB_flux between iterations with wB_flux [nondim] + real :: DwB_max ! The change in wB_flux when it is wB_flux_max [Z2 T-3 ~> m2 s-2] + real :: DwB_min ! The change in wB_flux when it is wB_flux_min [Z2 T-3 ~> m2 s-2] + real :: I_Gam_T, I_Gam_S ! Terms that vary inversely with Gam_mol_T or Gam_mol_S and Gam_turb [nondim] + real :: dG_dwB ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2] real :: taux2, tauy2 ! The squared surface stresses [R2 L2 Z2 T-4 ~> Pa2]. real :: u2_av, v2_av ! The ice-area weighted average squared ocean velocities [L2 T-2 ~> m2 s-2] - real :: asu1, asu2 ! Ocean areas covered by ice shelves at neighboring u- - real :: asv1, asv2 ! and v-points [L2 ~> m2]. + real :: asu1, asu2 ! Ocean areas covered by ice shelves at neighboring u-points [L2 ~> m2] + real :: asv1, asv2 ! Ocean areas covered by ice shelves at neighboring v-points [L2 ~> m2] real :: I_au, I_av ! The Adcroft reciprocals of the ice shelf areas at adjacent points [L-2 ~> m-2] real :: Irho0 ! The inverse of the mean density times a unit conversion factor [R-1 L Z-1 ~> m3 kg-1] logical :: Sb_min_set, Sb_max_set + logical :: root_found logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities. logical :: coupled_GL ! If true, the grounding line position is determined based on ! coupled ice-ocean dynamics. - real, parameter :: c2_3 = 2.0/3.0 - character(len=160) :: mesg ! The text of an error message + real, parameter :: c2_3 = 2.0/3.0 ! Two thirds [nondim] + character(len=320) :: mesg ! The text of an error message integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, ied, jed, it1, it3 - real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation [Z L2 ~> m3] - !for all ice sheets, Antarctica only, or Greenland only [Z L2 ~> m3] + real :: vaf0, vaf0_A, vaf0_G ! The previous volumes above floatation [Z L2 ~> m3] + ! for all ice sheets, Antarctica only, or Greenland only if (.not. associated(CS)) call MOM_error(FATAL, "shelf_calc_flux: "// & "initialize_ice_shelf must be called before shelf_calc_flux.") @@ -398,8 +407,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) ! useful parameters ZETA_N = CS%Zeta_N VK = CS%Vk - RC = CS%Rc - I_ZETA_N = 1.0 / ZETA_N + Rf_crit = CS%Rc + I_2Zeta_N = 0.5 / CS%Zeta_N I_LF = 1.0 / CS%Lat_fusion SC = CS%kv_molec/CS%kd_molec_salt PR = CS%kv_molec/CS%kd_molec_temp @@ -509,11 +518,12 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) if (absf*sfc_state%Hml(i,j) <= VK*ustar_h) then ; hBL_neut = sfc_state%Hml(i,j) else ; hBL_neut = (VK*ustar_h) / absf ; endif hBL_neut_h_molec = ZETA_N * ((hBL_neut * ustar_h) / (5.0 * CS%kv_molec)) + ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec) + n_star_term = (ZETA_N * hBL_neut * VK) / (Rf_crit * ustar_h**3) ! Determine the mixed layer buoyancy flux, wB_flux. dB_dS = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * dR0_dS(i) dB_dT = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * dR0_dT(i) - ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec) if (CS%find_salt_root) then ! Solve for the skin salinity using the linearized liquidus parameters and @@ -563,68 +573,152 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) dT_ustar = (ISS%tfreeze(i,j) - sfc_state%sst(i,j)) * ustar_h dS_ustar = (Sbdry(i,j) - sfc_state%sss(i,j)) * ustar_h - ! First, determine the buoyancy flux assuming no effects of stability - ! on the turbulence. Following H & J '99, this limit also applies - ! when the buoyancy flux is destabilizing. - - if (CS%const_gamma) then ! if using a constant gamma_T - ! note the different form, here I_Gam_T is NOT 1/Gam_T! + if (CS%const_gamma) then + ! If using a constant gamma_T, there are no effects of the buoyancy flux on the turbulence. I_Gam_T = CS%Gamma_T_3EQ I_Gam_S = CS%Gamma_S_3EQ - else - Gam_turb = I_VK * (ln_neut + (0.5 * I_ZETA_N - 1.0)) + wT_flux = dT_ustar * CS%Gamma_T_3EQ + wB_flux = dB_dS * (dS_ustar * CS%Gamma_S_3EQ) + dB_dT * wT_flux + elseif (.not.CS%buoy_flux_itt_bugfix) then + ! Gamma_T and gamma_S are a function of the buoyancy flux, and there should have been + ! iteration to find the root where wB_flux is consistent with the values of gamma with + ! that flux, but it was omitted. + Gam_turb = I_VK * (ln_neut + (I_2Zeta_N - 1.0)) I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) - endif + wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * (dT_ustar * I_Gam_T) - wT_flux = dT_ustar * I_Gam_T - wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux + if (wB_flux < 0.0) then ! The stabilising buoyancy flux reduces the turbulent fluxes. + I_n_star = sqrt(1.0 - n_star_term * wB_flux) + if (hBL_neut_h_molec > I_n_star**2) then + Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + else ! The layer dominated by molecular viscosity is smaller than the boundary layer. + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + endif + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + endif + wT_flux = dT_ustar * I_Gam_T + else ! gamma_T and gamma_S are a function of the buoyancy flux with proper iteration. + ! Find the root where wB_flux is consistent with the values of gamma with that flux. + + ! First, determine the buoyancy flux assuming no effects of stability + ! on the turbulence. Following H & J '99, this limit also applies + ! when the buoyancy flux is destabilizing. + Gam_turb = I_VK * (ln_neut + (I_2Zeta_N - 1.0)) + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + wB_flux = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T - if (wB_flux < 0.0) then - ! The buoyancy flux is stabilizing and will reduce the turbulent - ! fluxes, and iteration is required. - n_star_term = (ZETA_N * hBL_neut * VK) / (RC * ustar_h**3) - do it3 = 1,30 - ! n_star <= 1.0 is the ratio of working boundary layer thickness - ! to the neutral thickness. - ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL + if (wB_flux < 0.0) then + ! The buoyancy flux is stabilizing and will reduce the turbulent + ! fluxes, and iteration is required. + ! n_star <= 1.0 is the ratio of working boundary layer thickness + ! to the neutral thickness. I_n_star is its inverse. I_n_star = sqrt(1.0 - n_star_term * wB_flux) - dIns_dwB = 0.5 * n_star_term / I_n_star if (hBL_neut_h_molec > I_n_star**2) then - Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + & - (0.5*I_ZETA_N*I_n_star - 1.0)) - dG_dwB = I_VK * ( -2.0 / I_n_star + (0.5 * I_ZETA_N)) * dIns_dwB - else - ! The layer dominated by molecular viscosity is smaller than - ! the assumed boundary layer. This should be rare! - Gam_turb = I_VK * (0.5 * I_ZETA_N*I_n_star - 1.0) - dG_dwB = I_VK * (0.5 * I_ZETA_N) * dIns_dwB + Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + else ! The layer dominated by molecular viscosity is smaller than the boundary layer. + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) endif - - if (CS%const_gamma) then ! if using a constant gamma_T - ! note the different form, here I_Gam_T is NOT 1/Gam_T! - I_Gam_T = CS%Gamma_T_3EQ - I_Gam_S = CS%Gamma_S_3EQ - else - I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) - I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + + wB_flux_new = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T + root_found = (abs(wB_flux_new - wB_flux) < CS%buoy_flux_tol*(abs(wB_flux_new) + abs(wB_flux))) + ! Do not update the flux if its maagnitude would be increased by the otherwise + ! stabilizing buoyancy fluxes. This can happen when the buoyancy flux + ! is stabilizing when one of the heat or salt fluxes are destabilizing due + ! to their different molecular properties. + if (wB_flux_new <= wB_flux) root_found = .true. + + if (.not.root_found) then + wB_flux_max = 0.0 ; DwB_max = wB_flux + wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux + + if ((wB_flux_min*n_star_term < (1.0 - hBL_neut_h_molec)) .and. & + ((1.0 - hBL_neut_h_molec) < wB_flux_max*n_star_term)) then + ! The derivative of Gam_turb with wB_flux has a discontinuous change within the + ! bracketed range of values. Take this discontinous slope value for a first + ! guess, because Newton's method and the false position method may not converge + ! quickly when this discontinuity is between a guess and the solution. + wB_flux = (1.0 - hBL_neut_h_molec) / n_star_term + I_n_star = sqrt(hBL_neut_h_molec) + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + wB_flux_new = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T + + if (abs(wB_flux_new - wB_flux) <= CS%buoy_flux_tol*(abs(wB_flux_new) + abs(wB_flux))) then + ! The root has been found to within the tolerance at the kink. This should be very rare. + root_found = .true. + elseif (wB_flux_new > wB_flux) then + ! The solution is in the limit where abs(wB_flux) is small and + ! Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux + else + ! The solution is in the limt where abs(wB_flux) is large and + ! Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + wB_flux_max = wB_flux ; DwB_max = wB_flux_new - wB_flux + endif + endif endif - wT_flux = dT_ustar * I_Gam_T - wB_flux_new = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux - - ! Find the root where wB_flux_new = wB_flux. - if (abs(wB_flux_new - wB_flux) < CS%buoy_flux_itt_threshold*(abs(wB_flux_new) + abs(wB_flux))) exit + if (.not.root_found) then + ! Use the false position for the next guess. + wB_flux = wB_flux_min + (wB_flux_max-wB_flux_min) * (DwB_min / (DwB_min - DwB_max)) + + do it3 = 1,30 + ! Iterate using Newton's method with bounds or the false position method to find the root. + + I_n_star = sqrt(1.0 - n_star_term * wB_flux) + dIns_dwB = -0.5 * n_star_term / I_n_star + if (hBL_neut_h_molec > I_n_star**2) then + Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0)) + dG_dwB = I_VK * (( -2.0 / I_n_star + I_2Zeta_N) * dIns_dwB) + else + ! The layer dominated by molecular viscosity is smaller than the boundary layer. + Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0) + dG_dwB = I_VK * (I_2Zeta_N * dIns_dwB) + endif + I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) + I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb) + wB_flux_new = (dB_dS * dS_ustar) * I_Gam_S + (dB_dT * dT_ustar) * I_Gam_T + + ! Test for convergence to within tolerance at the point where wB_flux_new = wB_flux. + if (abs(wB_flux_new - wB_flux) <= CS%buoy_flux_tol*(abs(wB_flux_new) + abs(wB_flux))) & + root_found = .true. + if (root_found) exit + + dDwB_dwB = -dG_dwB * ((dB_dS * dS_ustar) * I_Gam_S**2 + & + (dB_dT * dT_ustar) * I_Gam_T**2) - 1.0 + if ((dDwB_dwB >= 0.0) .or. & + ( wB_flux - wB_flux_new >= abs(dDwB_dwB)*(wB_flux_max - wB_flux)) .or. & + ( wB_flux - wB_flux_new <= abs(dDwB_dwB)*(wB_flux_min - wB_flux)) ) then + ! Use the False position method to determine the guess for the next iteration when + ! Newton's method would go out of bounds + wB_flux_next = wB_flux_min + (wB_flux_max-wB_flux_min) * (DwB_min / (DwB_min - DwB_max)) + else + ! Use Newton's method for the next guess. + wB_flux_next = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB + endif + + ! Reset one of the bounds inward. + if (wB_flux_new - wB_flux > 0) then + wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux + else + wB_flux_max = wB_flux ; DwB_max = wB_flux_new - wB_flux + endif + + ! Update wB_flux + wB_flux = wB_flux_next + enddo ! it3 + endif - dDwB_dwB_in = dG_dwB * (dB_dS * (dS_ustar * I_Gam_S**2) + & - dB_dT * (dT_ustar * I_Gam_T**2)) - 1.0 - ! This is Newton's method without any bounds. Should bounds be needed? - wB_flux_new = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB_in - ! Update wB_flux - if (CS%buoy_flux_itt_bug) wB_flux = wB_flux_new - enddo !it3 - endif + endif ! End of test for first guess of wB_flux < 0. + wT_flux = dT_ustar * I_Gam_T + endif ! End of test for CS%const_gamma ISS%tflux_ocn(i,j) = RhoCp * wT_flux exch_vel_t(i,j) = ustar_h * I_Gam_T @@ -695,7 +789,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) Sbdry(i,j) = Sbdry_it endif ! Sb_min_set - if (.not.CS%salt_flux_itt_bug) Sbdry(i,j) = Sbdry_it + if (.not.CS%salt_flux_itt_bugfix) Sbdry(i,j) = Sbdry_it endif ! CS%find_salt_root @@ -1145,7 +1239,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step) type(ice_shelf_CS), pointer :: CS !< This module's control structure. type(surface), intent(inout) :: sfc_state !< Surface ocean state type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated. - real, intent(in) :: time_step !< Time step over which fluxes are applied + real, intent(in) :: time_step !< Time step over which fluxes are applied [T ~> s] ! local variables real :: frac_shelf !< The fractional area covered by the ice shelf [nondim]. real :: frac_open !< The fractional area of the ocean that is not covered by the ice shelf [nondim]. @@ -1384,10 +1478,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, type(directories) :: dirs type(dyn_horgrid_type), pointer :: dG => NULL() type(dyn_horgrid_type), pointer :: dG_in => NULL() - real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic. + real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic + ! [T kg R-1 Z-1 m-2 s-1 ~> nondim] real :: dz_ocean_min_float ! The minimum ocean thickness above which the ice shelf is considered ! to be floating when CONST_SEA_LEVEL = True [Z ~> m]. - real :: cdrag, drag_bg_vel + real :: cdrag ! The drag coefficient at the ice-ocean interface [nondim] + real :: drag_bg_vel ! A background velocity used in the quadratic drag [Z T-1 ~> m s-1] logical :: new_sim, save_IC !This include declares and sets the variable "version". # include "version_variable.h" @@ -1403,7 +1499,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, real :: utide ! A tidal velocity [L T-1 ~> m s-1] 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 + real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for ice shelf input data [L T-1 ~> m s-1] type(surface), pointer :: sfc_state => NULL() type(vardesc) :: u_desc, v_desc @@ -1695,11 +1791,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, call get_param(param_file, mdl, "ICE_SHELF_RC", CS%Rc, & "Critical flux Richardson number for ice melt ", & units="nondim", default=0.20) - call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUG", CS%buoy_flux_itt_bug, & - "Bug fix of buoyancy iteration", default=.true.) - call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUG", CS%salt_flux_itt_bug, & - "Bug fix of salt iteration", default=.true.) - call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_THRESHOLD", CS%buoy_flux_itt_threshold, & + call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX", CS%buoy_flux_itt_bugfix, & + "Bug fix of buoyancy iteration", default=.true., old_name="ICE_SHELF_BUOYANCY_FLUX_ITT_BUG") + call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUGFIX", CS%salt_flux_itt_bugfix, & + "Bug fix of salt iteration", default=.true., old_name="ICE_SHELF_SALT_FLUX_ITT_BUG") + call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_THRESHOLD", CS%buoy_flux_tol, & "Convergence criterion of Newton's method for ice shelf "//& "buoyancy iteration.", units="nondim", default=1.0e-4) @@ -2370,6 +2466,7 @@ subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim) end select end subroutine initialize_shelf_mass + !> This subroutine applies net accumulation/ablation at the top surface to the dynamic ice shelf. !! acc_rate[m-s]=surf_mass_flux/density_ice is ablation/accumulation rate !! positive for accumulation negative for ablation @@ -2386,14 +2483,13 @@ subroutine change_thickness_using_precip(CS, ISS, G, US, fluxes, time_step, Time ! locals integer :: i, j - real ::I_rho_ice + real :: I_rho_ice ! The specific volume of ice [R-1 ~> m3 kg-1] I_rho_ice = 1.0 / CS%density_ice !update time ! CS%Time = Time - ! CS%time_step = time_step ! update surface mass flux rate ! if (CS%surf_mass_flux_from_file) call update_surf_mass_flux(G, US, CS, ISS, Time) @@ -2477,7 +2573,7 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_ type(ice_shelf_CS), pointer :: CS !< ice shelf control structure type(ocean_grid_type), intent(in) :: G !< A pointer to an ocean grid control structure. real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: frac_shelf_h !< Ice shelf area fraction [nondim]. - real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf ! kg m-2] + real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf !< Ice shelf mass [R Z ~> kg m-2] logical, optional :: data_override_shelf_fluxes !< If true, shelf fluxes can be written using !! the data_override capability (only for MOSAIC grids) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 9ad992bb97..2fd03e30bc 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -2540,7 +2540,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif if (CS%max_surface_slope>0) then - scale = min(CS%max_surface_slope/sqrt((sx**2)+(sy**2)),1.0) + scale = CS%max_surface_slope / max( sqrt((sx**2) + (sy**2)), CS%max_surface_slope ) sx = scale*sx; sy = scale*sy endif From 44e9bd57232ef9d12e07849462560f71f7fe383c Mon Sep 17 00:00:00 2001 From: Alex Huth Date: Tue, 9 Dec 2025 16:34:45 -0500 Subject: [PATCH 7/9] Fix for ice-shelf friction velocity bugs (#995) * Fix for ice-shelf friction velocity bugs Fixed an incorrect area used to calculate cell-centered ocean surface velocity under the ice_shelf, which can impact the calculation of ice-shelf friction velocity. Added missing flags to some allocate_surface_state calls so that sfc_state%taux_shelf and sfc_state%tauy_shelf are allocated. This is required for the surface-stress-based (rather than surface-velocity-based) calculation of ice-shelf friction velocity. Also added taux_shelf and tauy_shelf as diagnostics for the surface stress under the ice shelf. * Removed unneeded taux_shelf and tauy_shelf diagnostics * Added ustar_from_vel_bugfix flag, which if true, fixes the ustar from ocean velocity bug --- config_src/drivers/FMS_cap/ocean_model_MOM.F90 | 3 ++- .../drivers/STALE_mct_cap/mom_ocean_model_mct.F90 | 3 ++- config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 | 3 ++- src/ice_shelf/MOM_ice_shelf.F90 | 10 +++++++++- 4 files changed, 15 insertions(+), 4 deletions(-) diff --git a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index 110415c6e0..13ab3eecef 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -374,7 +374,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas !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) + gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot, & + use_iceshelves=OS%use_ice_shelf) if (present(wind_stagger)) then call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & diff --git a/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 b/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 index 6e1545efe1..990a52c427 100644 --- a/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 +++ b/config_src/drivers/STALE_mct_cap/mom_ocean_model_mct.F90 @@ -364,7 +364,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! Consider using a run-time flag to determine whether to do the diagnostic ! vertical integrals, since the related 3-d sums are not negligible in cost. 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) + do_integrals=.true., gas_fields_ocn=gas_fields_ocn, & + use_meltpot=use_melt_pot, use_iceshelves=OS%use_ice_shelf) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 878848260d..4dc68d7b0a 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -389,7 +389,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i ! vertical integrals, since the related 3-d sums are not negligible in cost. 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, use_marbl_tracers=OS%use_MARBL) + use_meltpot=use_melt_pot, use_iceshelves=OS%use_ice_shelf, & + use_marbl_tracers=OS%use_MARBL) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 971fc5005e..5be4e8c6f0 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -200,6 +200,7 @@ module MOM_ice_shelf !! divided by the von Karman constant VK [nondim]. Was 1/8. real :: Vk !< Von Karman's constant [nondim] real :: Rc !< critical flux Richardson number [nondim] + logical :: ustar_from_vel_bugfix !< If true, fixes ustar from ocean velocity bug logical :: buoy_flux_itt_bugfix !< If true, fixes buoyancy iteration bug logical :: salt_flux_itt_bugfix !< If true, fixes salt iteration bug real :: buoy_flux_tol !< Fractional buoyancy iteration tolerance for convergence [nondim] @@ -467,7 +468,11 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) tauy2 = (((asv1 * (sfc_state%tauy_shelf(i,J-1)**2)) + (asv2 * (sfc_state%tauy_shelf(i,J)**2)) ) * I_av) endif 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 (CS%ustar_from_vel_bugfix) then + v2_av = (((asv1 * (sfc_state%v(i,J-1)**2)) + (asv2 * sfc_state%v(i,J)**2)) * I_av) + else + v2_av = (((asv1 * (sfc_state%v(i,J-1)**2)) + (asu2 * sfc_state%v(i,J)**2)) * I_av) + endif if ((taux2 + tauy2 > 0.0) .and. .not.CS%ustar_shelf_from_vel) then if (CS%ustar_max >= 0.0) then @@ -1791,6 +1796,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, call get_param(param_file, mdl, "ICE_SHELF_RC", CS%Rc, & "Critical flux Richardson number for ice melt ", & units="nondim", default=0.20) + call get_param(param_file, mdl, "ICE_SHELF_USTAR_FROM_VEL_BUGFIX", CS%ustar_from_vel_bugfix, & + "Bug fix for ice-area weighting of squared ocean velocities "//& + "used to calculate friction velocity under ice shelves", default=.false.) call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX", CS%buoy_flux_itt_bugfix, & "Bug fix of buoyancy iteration", default=.true., old_name="ICE_SHELF_BUOYANCY_FLUX_ITT_BUG") call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUGFIX", CS%salt_flux_itt_bugfix, & From 42fd96d104e72a1ab2de53e20cca52ec73f925a7 Mon Sep 17 00:00:00 2001 From: Claire Yung Date: Fri, 10 Oct 2025 14:09:08 +1100 Subject: [PATCH 8/9] Add FRAZIL_NOT_UNDER_ICESHELF parameter to turn off frazil scheme under ice shelves If True, frazil scheme is not applied where frac_shelf_h > 0. If False, frazil calculated everywhere as usual if FRAZIL = True. Defaults to False so shouldn't change answers. --- .../vertical/MOM_diabatic_aux.F90 | 38 ++++++++++++++++++- .../vertical/MOM_diabatic_driver.F90 | 20 +++++++++- 2 files changed, 55 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index bcde4feb34..516c3bec26 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -109,7 +109,7 @@ module MOM_diabatic_aux !! This subroutine warms any water that is colder than the (currently !! surface) freezing point up to the freezing point and accumulates !! the required heat (in [Q R Z ~> J m-2]) in tv%frazil. -subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) +subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo, frac_shelf_h) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -122,6 +122,9 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]. integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage + ! Local variables real, dimension(SZI_(G)) :: & fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2]. @@ -189,7 +192,39 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) endif endif ; enddo endif + if (PRESENT(frac_shelf_h)) then + do k=nz,1,-1 + T_fr_set = .false. + do i=is,ie + if (((G%mask2dT(i,j) > 0.0) .and. (.not. (frac_shelf_h(i,j) > 0.0))) .and. & + ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then + if (.not.T_fr_set) then + call calculate_TFreeze(tv%S(i:ie,j,k), pressure(i:ie,k), T_freeze(i:ie), & + tv%eqn_of_state) + T_fr_set = .true. + endif + hc = (tv%C_p*GV%H_to_RZ) * h(i,j,k) + if (h(i,j,k) <= 10.0*(GV%Angstrom_H + GV%H_subroundoff)) then + ! Very thin layers should not be cooled by the frazil flux. + if (tv%T(i,j,k) < T_freeze(i)) then + fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) + tv%T(i,j,k) = T_freeze(i) + endif + elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < T_freeze(i))) then + if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) < 0.0) then + tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc + fraz_col(i) = 0.0 + else + fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) + tv%T(i,j,k) = T_freeze(i) + endif + endif + endif + enddo + enddo + + else do k=nz,1,-1 T_fr_set = .false. do i=is,ie @@ -220,6 +255,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) endif enddo enddo + endif do i=is,ie tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i) enddo diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 5389de3e1d..06b5ca1b6a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -188,6 +188,8 @@ module MOM_diabatic_driver logical :: Use_KdWork_diag = .false. !< Logical flag to indicate if any Kd_work diagnostics are on. logical :: Use_N2_diag = .false. !< Logical flag to indicate if any N2 diagnostics are on. + logical :: frazil_not_under_iceshelves !< If true, turn off frazil under ice shelves + !>@{ Diagnostic IDs integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1 @@ -380,7 +382,12 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (associated(fluxes%p_surf_full)) then - call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) + if (CS%frazil_not_under_iceshelves) then + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff, & + frac_shelf_h=fluxes%frac_shelf_h) + else + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) + endif else call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, halo=CS%halo_TS_diff) endif @@ -440,7 +447,12 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (associated(fluxes%p_surf_full)) then - call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full) + if (CS%frazil_not_under_iceshelves) then + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff, & + frac_shelf_h=fluxes%frac_shelf_h) + else + call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) + endif else call make_frazil(h, tv, G, GV, US, CS%diabatic_aux_CSp) endif @@ -3721,6 +3733,10 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di endif endif + call get_param(param_file, mdl, "FRAZIL_NOT_UNDER_ICESHELF", CS%frazil_not_under_iceshelves, & + "If true, do not use frazil scheme underneath ice shelves "//& + "defined by frac_shelf_h greater than 0.", default=.false.) + ! diagnostics for tendencies of temp and heat due to frazil CS%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, Time, & long_name='Cell Thickness', standard_name='cell_thickness', & From ef4fc4547345cae8ef50957145803b930f68f4f9 Mon Sep 17 00:00:00 2001 From: Claire Yung Date: Wed, 10 Dec 2025 12:11:37 +1100 Subject: [PATCH 9/9] Add option to override ice shelf melt mass and heat fluxes This commit modifies the nuopc cap and ice shelf code to add an option to override the prognostically calculated ice shelf basal melt with an input file, defined in data_table, using the data override functionality. To use this, set DATA_OVERRIDE_MELT = True. If False, should not change answers. --- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 5 +++++ src/ice_shelf/MOM_ice_shelf.F90 | 21 ++++++++++++++++++- 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 4dc68d7b0a..1441a98724 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -51,6 +51,8 @@ module MOM_ocean_model_nuopc 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_ice_shelf, only : initialize_ice_shelf_fluxes, initialize_ice_shelf_forces +use MOM_ice_shelf, only : ice_shelf_query +use MOM_data_override, only : data_override_init use MOM_coupler_types, only : coupler_1d_bc_type, coupler_2d_bc_type use MOM_coupler_types, only : coupler_type_spawn, coupler_type_write_chksums use MOM_coupler_types, only : coupler_type_initialized, coupler_type_copy_data @@ -211,6 +213,7 @@ module MOM_ocean_model_nuopc Ice_shelf_CSp => NULL() !< A pointer to the control structure for the !! ice shelf model that couples with MOM6. This !! is null if there is no ice shelf. + logical :: override_melt !< If true, override melt using data_override type(marine_ice_CS), pointer :: & marine_ice_CSp => NULL() !< A pointer to the control structure for the !! marine ice effects module. @@ -401,6 +404,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%fluxes) call initialize_ice_shelf_fluxes(OS%ice_shelf_CSp, OS%grid, OS%US, OS%flux_tmp) call initialize_ice_shelf_forces(OS%ice_shelf_CSp, OS%grid, OS%US, OS%forces) + call ice_shelf_query(OS%ice_shelf_CSp, OS%grid, data_override_melt=OS%override_melt) + if (OS%override_melt) call data_override_init(Ocean_Domain_in=OS%grid%domain%mpp_domain) 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 5be4e8c6f0..8f94a640c9 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -171,6 +171,7 @@ module MOM_ice_shelf !! shelf dynamics will be initialized logical :: data_override_shelf_fluxes !< True if the ice shelf surface mass fluxes can be !! written using the data_override feature (only for MOSAIC grids) + logical :: data_override_melt !< True if overriding melt rate and heat flux logical :: override_shelf_movement !< If true, user code specifies the shelf movement !! instead of using the dynamic ice-shelf mode. logical :: isthermo !< True if the ice shelf can exchange heat and @@ -1302,6 +1303,14 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step) enddo ; enddo endif + !CLAIRE + if (CS%data_override_melt) then + call data_override(G%Domain, 'water_flux', ISS%water_flux(is:ie,js:je), CS%Time, & + scale=US%kg_m2s_to_RZ_T) + call data_override(G%Domain, 'tflux_ocn', ISS%tflux_ocn(is:ie,js:je), CS%Time, & + scale=US%W_m2_to_QRZ_T) + endif + if (CS%debug) then call MOM_forcing_chksum("Before adding shelf fluxes", fluxes, G, CS%US, haloshift=0) endif @@ -1648,6 +1657,10 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, if (CS%GL_regularize) CS%GL_couple = .false. if (CS%solo_ice_sheet) CS%GL_couple = .false. endif + call get_param(param_file, mdl, "DATA_OVERRIDE_MELT", & + CS%data_override_melt, & + "If true, the data override feature is used for melt rate "//& + "and heat flux from ice shelves.", default=.false.) call get_param(param_file, mdl, "SHELF_THERMO", CS%isthermo, & "If true, use a thermodynamically interactive ice shelf.", & @@ -2577,13 +2590,14 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) end subroutine update_shelf_mass !> Save the ice shelf restart file -subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_fluxes) +subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_fluxes, data_override_melt) type(ice_shelf_CS), pointer :: CS !< ice shelf control structure type(ocean_grid_type), intent(in) :: G !< A pointer to an ocean grid control structure. real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: frac_shelf_h !< Ice shelf area fraction [nondim]. real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf !< Ice shelf mass [R Z ~> kg m-2] logical, optional :: data_override_shelf_fluxes !< If true, shelf fluxes can be written using !! the data_override capability (only for MOSAIC grids) + logical, optional :: data_override_melt !< If true, basal melt overridden integer :: i, j @@ -2606,6 +2620,11 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_ if (CS%active_shelf_dynamics) data_override_shelf_fluxes = CS%data_override_shelf_fluxes endif + if (present(data_override_melt)) then + data_override_melt = .false. + if (CS%data_override_melt) data_override_melt = CS%data_override_melt + endif + end subroutine ice_shelf_query !> Save the ice shelf restart file