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_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..1441a98724 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,9 @@ 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_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 @@ -134,7 +137,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. @@ -210,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. @@ -388,14 +392,20 @@ 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) 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_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) @@ -740,7 +750,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 @@ -749,7 +759,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 @@ -758,7 +768,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 @@ -819,7 +829,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..8f94a640c9 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 @@ -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 @@ -197,12 +198,14 @@ 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 :: 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] + 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, & @@ -297,11 +300,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]. @@ -321,7 +324,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] @@ -330,34 +334,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.") @@ -397,8 +409,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 @@ -457,19 +469,26 @@ 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 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 @@ -505,11 +524,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 @@ -559,68 +579,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 @@ -691,7 +795,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 @@ -1141,7 +1245,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]. @@ -1199,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 @@ -1380,10 +1492,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" @@ -1399,7 +1513,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 @@ -1543,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.", & @@ -1677,7 +1795,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 "//& @@ -1689,11 +1809,14 @@ 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_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, & + "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) @@ -1712,6 +1835,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.", & @@ -1732,7 +1858,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 @@ -1906,7 +2033,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 @@ -2359,6 +2487,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 @@ -2375,14 +2504,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) @@ -2462,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 ! 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) + logical, optional :: data_override_melt !< If true, basal melt overridden integer :: i, j @@ -2491,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 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 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', & 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 "// &