From 7fb8f55555574a2faa5f9801beb6add833eae8bc Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 30 Mar 2020 14:48:22 -0400 Subject: [PATCH 01/13] Flipped the sign convention for wT_flux Changed the sign conventions of the internal variables wT_flux, wB_flux, dT_ustar and dS_ustar in shelf_calc_flux to follow the vertical flux sign conventions in the rest of the MOM6 code. All answers are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index d05631c621..6fa7aef94e 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -245,19 +245,18 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) real :: hBL_neut !< The neutral boundary layer thickness [Z ~> m]. real :: hBL_neut_h_molec !< The ratio of the neutral boundary layer thickness !! to the molecular boundary layer thickness [nondim]. - !### THESE ARE CURRENTLY POSITIVE UPWARD. - real :: wT_flux !< The upward vertical flux of heat just inside the ocean [degC Z T-1 ~> degC m s-1]. - real :: wB_flux !< The upward vertical flux of buoyancy just inside the ocean [Z2 T-3 ~> m2 s-3]. + real :: wT_flux !< The downward vertical flux of heat just inside the ocean [degC Z T-1 ~> degC m s-1]. + 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 ppt-1 ~> m s-2 ppt-1]. real :: dB_dT !< The derivative of buoyancy with temperature [Z T-2 degC-1 ~> m s-2 degC-1]. real :: I_n_star ! [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] - real :: dT_ustar ! The difference between the ocean boundary layer temperature and the freezing - ! freezing point times the friction velocity [degC Z T-1 ~> degC m s-1] - real :: dS_ustar ! The difference between the ocean boundary layer salinity and the salinity - ! at the ice-ocean interface the friction velocity [ppt Z T-1 ~> ppt m s-1] + real :: dT_ustar ! The difference between the the freezing point and the ocean boundary layer + ! temperature times the friction velocity [degC Z T-1 ~> degC m s-1] + real :: dS_ustar ! The difference between the salinity at the ice-ocean interface and the ocean + ! boundary layer salinity times the friction velocity [ppt 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_mol_t, Gam_mol_s @@ -429,8 +428,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! Determine the potential temperature at the ice-ocean interface. call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state) - dT_ustar = (state%sst(i,j) - ISS%tfreeze(i,j)) * ustar_h - dS_ustar = (state%sss(i,j) - Sbdry(i,j)) * ustar_h + dT_ustar = (ISS%tfreeze(i,j) - state%sst(i,j)) * ustar_h + dS_ustar = (Sbdry(i,j) - 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 @@ -449,7 +448,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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 + if (wB_flux < 0.0) then ! The buoyancy flux is stabilizing and will reduce the tubulent ! fluxes, and iteration is required. n_star_term = (ZETA_N/RC) * (hBL_neut * VK) / (ustar_h)**3 @@ -458,7 +457,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! to the neutral thickness. ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL - I_n_star = sqrt(1.0 + n_star_term * wB_flux) + 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)) + & @@ -484,18 +483,17 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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) < & - 1e-4*(abs(wB_flux_new) + abs(wB_flux))) exit + if (abs(wB_flux_new - wB_flux) < 1e-4*(abs(wB_flux_new) + abs(wB_flux))) exit - dDwB_dwB_in = -dG_dwB * (dB_dS * (dS_ustar * I_Gam_S**2) + & - dB_dT * (dT_ustar * I_Gam_T**2)) - 1.0 + 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 enddo !it3 endif - ISS%tflux_ocn(i,j) = -RhoCp * wT_flux + ISS%tflux_ocn(i,j) = RhoCp * wT_flux exch_vel_t(i,j) = ustar_h * I_Gam_T exch_vel_s(i,j) = ustar_h * I_Gam_S @@ -1109,7 +1107,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl real :: cdrag, drag_bg_vel logical :: new_sim, save_IC, var_force !This include declares and sets the variable "version". -#include "version_variable.h" +# include "version_variable.h" character(len=200) :: config character(len=200) :: IC_file,filename,inputdir character(len=40) :: mdl = "MOM_ice_shelf" ! This module's name. From 63fd8e1e85a9012ab89c09c94c6571da40113ae9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 30 Mar 2020 17:38:25 -0400 Subject: [PATCH 02/13] +Added optional arg area to global_area_integral Added the new optional area argument to global_area_integral, to replace the area in G%areaT. All answers are bitwise identical. --- src/framework/MOM_spatial_means.F90 | 34 +++++++++++++++++++---------- 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90 index 85d5ce452b..d4b687b0a5 100644 --- a/src/framework/MOM_spatial_means.F90 +++ b/src/framework/MOM_spatial_means.F90 @@ -47,14 +47,18 @@ function global_area_mean(var, G, scale) end function global_area_mean -!> Return the global area integral of a variable. This uses reproducing sums. -function global_area_integral(var, G, scale) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure - real, dimension(SZI_(G), SZJ_(G)), intent(in) :: var !< The variable to integrate - real, optional, intent(in) :: scale !< A rescaling factor for the variable - real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming - real :: global_area_integral - +!> Return the global area integral of a variable, by default using the masked area from the +!! grid, but an alternate could be used instead. This uses reproducing sums. +function global_area_integral(var, G, scale, area) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< The variable to integrate + real, optional, intent(in) :: scale !< A rescaling factor for the variable + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: area !< The alternate area to use, including + !! any required masking [L2 ~> m2]. + real :: global_area_integral !< The returned area integral, usually in the units of var times [m2]. + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming real :: scalefac ! An overall scaling factor for the areas and variable. integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -62,9 +66,15 @@ function global_area_integral(var, G, scale) scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale tmpForSumming(:,:) = 0. - do j=js,je ; do i=is, ie - tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j)) - enddo ; enddo + if (present(area)) then + do j=js,je ; do i=is,ie + tmpForSumming(i,j) = var(i,j) * (scalefac * area(i,j)) + enddo ; enddo + else + do j=js,je ; do i=is,ie + tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j)) + enddo ; enddo + endif global_area_integral = reproducing_sum(tmpForSumming) end function global_area_integral @@ -96,7 +106,7 @@ function global_layer_mean(var, h, G, GV, scale) global_temp_scalar = reproducing_sum(tmpForSumming,sums=scalarij) global_weight_scalar = reproducing_sum(weight,sums=weightij) - do k=1, nz + do k=1,nz global_layer_mean(k) = scalarij(k) / weightij(k) enddo From 7121619b29982e37233be66141def8a85d4178c0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 30 Mar 2020 17:55:37 -0400 Subject: [PATCH 03/13] (*)Use global_area_integral in add_shelf_flux Replaced calls to the non-reproducing routines sum_across_PES with calls to global_area_integral that uses the reproducing sums when compensating for the global mean fresh water fluxes in add_shelf_flux. This also includes rescaling the dimensions of mean_melt_flux to [R Z T-1]. This could change answers at roundoff in some cases with an interactive ice shelf and CONST_SEA_LEVEL=True, but all answers in the MOM6-examples test cases are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 44 +++++++++++++++++---------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 6fa7aef94e..f9b397451c 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -45,7 +45,8 @@ module MOM_ice_shelf use MOM_ice_shelf_state, only : ice_shelf_state, ice_shelf_state_end, ice_shelf_state_init use user_shelf_init, only : USER_initialize_shelf_mass, USER_update_shelf_mass use user_shelf_init, only : user_ice_shelf_CS -use MOM_coms, only : reproducing_sum, sum_across_PEs +use MOM_coms, only : reproducing_sum +use MOM_spatial_means, only : global_area_integral use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum use time_interp_external_mod, only : init_external_field, time_interp_external use time_interp_external_mod, only : time_interp_external_init @@ -877,20 +878,21 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) real :: Irho0 !< The inverse of the mean density times unit conversion factors that !! arise because state uses MKS units [Z2 m s2 kg-1 T-2 ~> m3 kg-1]. real :: frac_area !< The fractional area covered by the ice shelf [nondim]. - real :: shelf_mass0 !< Total ice shelf mass at previous time (Time-dt). - real :: shelf_mass1 !< Total ice shelf mass at current time (Time). real :: delta_mass_shelf!< Change in ice shelf mass over one time step [kg s-1] real :: taux2, tauy2 !< The squared surface stresses [Pa]. real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa]. real :: asu1, asu2 !< Ocean areas covered by ice shelves at neighboring u- real :: asv1, asv2 !< and v-points [L2 ~> m2]. real :: fraz !< refreezing rate [kg m-2 s-1] - real :: mean_melt_flux !< spatial mean melt flux [kg s-1] or [kg m-2 s-1] at various points in the code. + real :: mean_melt_flux !< Spatial mean melt flux [R Z T-1 ~> kg m-2 s-1] real :: sponge_area !< total area of sponge region [m2] real :: t0 !< The previous time (Time-dt) [s]. type(time_type) :: Time0!< The previous time (Time-dt) + real, dimension(SZDI_(G),SZDJ_(G)) :: in_sponge !< 1 where the property damping occurs, 0 otherwise [nondim] real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass !! at at previous time (Time-dt) [R Z ~> kg m-2] + real, dimension(SZDI_(G),SZDJ_(G)) :: delta_float_mass !< The change in the floating mass between + !! the two timesteps at (Time) and (Time-dt) [R Z ~> kg m-2]. real, dimension(SZDI_(G),SZDJ_(G)) :: last_h_shelf !< Ice shelf thickness [Z ~> m] !! at at previous time (Time-dt) real, dimension(SZDI_(G),SZDJ_(G)) :: last_hmask !< Ice shelf mask @@ -994,15 +996,13 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je)) fluxes%salt_flux(:,:) = 0.0 ; fluxes%vprec(:,:) = 0.0 - mean_melt_flux = 0.0; sponge_area = 0.0 do j=js,je ; do i=is,ie - frac_area = fluxes%frac_shelf_h(i,j) - if (frac_area > 0.0) & - mean_melt_flux = mean_melt_flux + (ISS%water_flux(i,j)) * US%RZ_T_to_kg_m2s*US%L_to_m**2*ISS%area_shelf_h(i,j) !### These hard-coded limits need to be corrected. They are inappropriate here. if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then - sponge_area = sponge_area + US%L_to_m**2*G%areaT(i,j) + in_sponge(i,j) = 1.0 + else + in_sponge(i,j) = 0.0 endif enddo ; enddo @@ -1027,20 +1027,18 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) last_mass_shelf(:,:) = last_h_shelf(:,:) * CS%density_ice endif - shelf_mass0 = 0.0; shelf_mass1 = 0.0 ! get total ice shelf mass at (Time-dt) and (Time), in kg do j=js,je ; do i=is,ie ! just floating shelf (0.1 is a threshold for min ocean thickness) if (((1.0/CS%density_ocean_avg)*state%ocean_mass(i,j) > 0.1) .and. & (ISS%area_shelf_h(i,j) > 0.0)) then - shelf_mass0 = shelf_mass0 + US%RZ_to_kg_m2*US%L_to_m**2*(last_mass_shelf(i,j) * ISS%area_shelf_h(i,j)) - shelf_mass1 = shelf_mass1 + US%RZ_to_kg_m2*US%L_to_m**2*(ISS%mass_shelf(i,j) * ISS%area_shelf_h(i,j)) + delta_float_mass(i,j) = ISS%mass_shelf(i,j) - last_mass_shelf(i,j) + else + delta_float_mass(i,j) = 0.0 endif enddo ; enddo - call sum_across_PEs(shelf_mass0); call sum_across_PEs(shelf_mass1) - delta_mass_shelf = (shelf_mass1 - shelf_mass0)/CS%time_step -! write(mesg,*) 'delta_mass_shelf = ', delta_mass_shelf -! call MOM_mesg(mesg,5) + delta_mass_shelf = US%kg_m2s_to_RZ_T*(global_area_integral(delta_float_mass, G, scale=US%RZ_to_kg_m2, & + area=ISS%area_shelf_h) / CS%time_step) else! first time step delta_mass_shelf = 0.0 endif @@ -1048,18 +1046,22 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) delta_mass_shelf = 0.0 endif - call sum_across_PEs(mean_melt_flux) - call sum_across_PEs(sponge_area) - ! average total melt flux over sponge area - mean_melt_flux = (mean_melt_flux+delta_mass_shelf) / sponge_area !kg/(m^2 s) + sponge_area = global_area_integral(in_sponge, G) + if (sponge_area > 0.0) then + mean_melt_flux = US%kg_m2s_to_RZ_T*(global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & + area=ISS%area_shelf_h) + & + delta_mass_shelf ) / sponge_area + else + mean_melt_flux = 0.0 + endif ! apply fluxes do j=js,je ; do i=is,ie ! Note the following is hard coded for ISOMIP if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then ! evap is negative, and vprec has units of [R Z T-1 ~> kg m-2 s-1] - fluxes%vprec(i,j) = -US%kg_m2s_to_RZ_T*mean_melt_flux * CS%density_ice / (1000.0*US%kg_m3_to_R) + fluxes%vprec(i,j) = -mean_melt_flux * CS%density_ice / (1000.0*US%kg_m3_to_R) fluxes%sens(i,j) = fluxes%vprec(i,j) * CS%Cp * CS%T0 ! [ Q R Z T-1 ~> W /m^2 ] fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * CS%S0*1.0e-3 ! [kgSalt/kg R Z T-1 ~> kgSalt m-2 s-1] endif From 3c3f72167b8b8e5da306021979bfc4c2697f08ed Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 30 Mar 2020 21:45:40 -0400 Subject: [PATCH 04/13] (*+)Corrected bugs in 3-eqn ice shelf skin salinity Corrected several bugs in the 3-equation ice shelf skin salinity calculation, including renaming variables for greater clarity and using forms for the solutions to a quadratic equation that are accurate without amplifying roundoff errors. In addition, a new runtime parameter, SHELF_3EQ_GAMMA_S, is read and logged when SHELF_3EQ_GAMMA is true. This will change answers and the parameter_doc files with when a thermodynamically active ice shelf is used and SHELF_THREE_EQN is true, but all answers in the MOM6-examples test cases are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 106 +++++++++++++++++--------------- 1 file changed, 57 insertions(+), 49 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index f9b397451c..a4f6169844 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -26,7 +26,7 @@ module MOM_ice_shelf use MOM_io, only : write_field, close_file, SINGLE_FILE, MULTIPLE use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, restore_state, MOM_restart_CS -use MOM_time_manager, only : time_type, time_type_to_real, real_to_time +use MOM_time_manager, only : time_type, time_type_to_real, real_to_time, operator(>), operator(-) use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init, fix_restart_unit_scaling use MOM_variables, only : surface @@ -106,6 +106,7 @@ module MOM_ice_shelf 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 :: col_thick_melt_threshold !< if the mixed layer is below this threshold, melt rate logical :: mass_from_file !< Read the ice shelf mass from a file every dt @@ -150,14 +151,13 @@ module MOM_ice_shelf !! interface. logical :: insulator !< If true, ice shelf is a perfect insulator logical :: const_gamma !< If true, gamma_T is specified by the user. - logical :: find_salt_root !< If true, if true find Sbdry using a quadratic eq. logical :: constant_sea_level !< if true, apply an evaporative, heat and salt !! fluxes. It will avoid large increase in sea level. real :: cutoff_depth !< Depth above which melt is set to zero (>= 0) [Z ~> m]. - ! The following parameters are needed if find_salt_root = true - real :: lambda1 !< liquidus coeff. The freezing point at 0 pressure and 0 salinity [degC] - real :: lambda2 !< Partial derivative of freezing temperature with salinity [degC ppt-1] - real :: lambda3 !< Partial derivative of freezing temperature with pressure [degC Pa-1] + logical :: find_salt_root !< If true, if true find Sbdry using a quadratic eq. + real :: TFr_0_0 !< The freezing point at 0 pressure and 0 salinity [degC] + real :: dTFr_dS !< Partial derivative of freezing temperature with salinity [degC ppt-1] + real :: dTFr_dp !< Partial derivative of freezing temperature with pressure [degC Pa-1] !>@{ Diagnostic handles integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & id_tfreeze = -1, id_tfl_shelf = -1, & @@ -241,7 +241,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: & Sbdry !< Salinities in the ocean at the interface with the ice shelf [ppt]. real :: Sbdry_it - real :: Sbdry1, Sbdry2, S_a, S_b, S_c ! Variables used to find salt roots + real :: Sbdry1, Sbdry2 + real :: S_a, S_b, S_c ! Variables used to find salt roots real :: dS_it !< The interface salinity change during an iteration [ppt]. real :: hBL_neut !< The neutral boundary layer thickness [Z ~> m]. real :: hBL_neut_h_molec !< The ratio of the neutral boundary layer thickness @@ -391,26 +392,31 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec) if (CS%find_salt_root) then - ! read liquidus parameters - - !### This should be CS%lamda2! - S_a = CS%lambda1 * CS%Gamma_T_3EQ * CS%Cp - ! The value of 35.0 here should be a parameter? - !### This should be (CS%lambda1 + CS%lambda3*p_int(i) - state%sst(i,j)) - S_b = CS%Gamma_T_3EQ*CS%Cp*(CS%lambda2 + CS%lambda3*p_int(i)- state%sst(i,j)) - & - CS%Lat_fusion * CS%Gamma_T_3EQ/35.0 - S_c = CS%Lat_fusion * (CS%Gamma_T_3EQ/35.0) * state%sss(i,j) - - !### Depending on the sign of S_b, one of these will be inaccurate! - ! if (S_b >= 0.0) then - Sbdry1 = (-S_b + SQRT(S_b*S_b - 4*S_a*S_c)) / (2*S_a) - ! Sbdry1 = 2*S_c / (S_b + SQRT(S_b*S_b - 4*S_a*S_c)) - Sbdry2 = (-S_b - SQRT(S_b*S_b - 4*S_a*S_c)) / (2*S_a) - ! else - ! Sbdry1 = (-S_b + SQRT(S_b*S_b - 4.*S_a*S_c)) / (2.*S_a) - ! Sbdry2 = -2.*S_c / (-S_b + SQRT(S_b*S_b - 4.*S_a*S_c)) - ! endif - Sbdry(i,j) = MAX(Sbdry1, Sbdry2) + ! Solve for the skin salinity using the linearized liquidus parameters and + ! balancing the turbulent fresh water flux in the near-boundary layer with + ! the net fresh water or salt added by melting: + ! (Cp/Lat_fusion)*Gamma_T_3Eq*(TFr_skin-T_ocn) = Gamma_S_3Eq*(S_skin-S_ocn)/S_skin + + ! S_a is always < 0.0 with a realistic expression for the freezing point. + S_a = CS%dTFr_dS * CS%Gamma_T_3EQ * CS%Cp + S_b = CS%Gamma_T_3EQ*CS%Cp*(CS%TFr_0_0 + CS%dTFr_dp*p_int(i) - state%sst(i,j)) - & + CS%Lat_fusion * CS%Gamma_S_3EQ ! S_b Can take either sign, but is usually negative. + S_c = CS%Lat_fusion * CS%Gamma_S_3EQ * state%sss(i,j) ! Always >= 0 + + if (S_c == 0.0) then ! The solution for fresh water. + Sbdry(i,j) = 0.0 + elseif (S_a < 0.0) then ! This is the usual ocean case + if (S_b < 0.0) then ! This is almost always the case + Sbdry(i,j) = 2.0*S_c / (-S_b + SQRT(S_b*S_b - 4.*S_a*S_c)) + else + Sbdry(i,j) = (S_b + SQRT(S_b*S_b - 4.*S_a*S_c)) / (-2.*S_a) + endif + elseif ((S_a == 0.0) .and. (S_b < 0.0)) then ! It should be the case that S_b < 0. + Sbdry(i,j) = -S_c / S_b + else + call MOM_error(FATAL, "Impossible conditions found in 3-equation skin salinity calculation.") + endif + ! Safety check if (Sbdry(i,j) < 0.) then write(mesg,*) 'state%sss(i,j) = ',state%sss(i,j), 'S_a, S_b, S_c', S_a, S_b, S_c @@ -439,7 +445,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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_T_3EQ/35. + I_Gam_S = CS%Gamma_S_3EQ else Gam_turb = I_VK * (ln_neut + (0.5 * I_ZETA_N - 1.0)) I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb) @@ -474,7 +480,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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_T_3EQ/35. + 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) @@ -883,11 +889,10 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa]. real :: asu1, asu2 !< Ocean areas covered by ice shelves at neighboring u- real :: asv1, asv2 !< and v-points [L2 ~> m2]. - real :: fraz !< refreezing rate [kg m-2 s-1] real :: mean_melt_flux !< Spatial mean melt flux [R Z T-1 ~> kg m-2 s-1] real :: sponge_area !< total area of sponge region [m2] - real :: t0 !< The previous time (Time-dt) [s]. - type(time_type) :: Time0!< The previous time (Time-dt) + type(time_type) :: dTime !< The time step as a time_type + type(time_type) :: Time0 !< The previous time (Time-dt) real, dimension(SZDI_(G),SZDJ_(G)) :: in_sponge !< 1 where the property damping occurs, 0 otherwise [nondim] real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass !! at at previous time (Time-dt) [R Z ~> kg m-2] @@ -989,8 +994,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! This is needed for some of the ISOMIP+ experiments. if (CS%constant_sea_level) then - !### This code has lots of problems with hard coded constants and the use of - !### of non-reproducing sums. It needs to be refactored. -RWH + !### This code has problems with hard coded constants that need to be refactored. -RWH if (.not. associated(fluxes%salt_flux)) allocate(fluxes%salt_flux(ie,je)) if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je)) @@ -1008,11 +1012,11 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! take into account changes in mass (or thickness) when imposing ice shelf mass if (CS%override_shelf_movement .and. CS%mass_from_file) then - t0 = time_type_to_real(CS%Time) - CS%time_step + dTime = real_to_time(CS%time_step) - ! just compute changes in mass after first time step - if (t0>0.0) then - Time0 = real_to_time(t0) + ! Compute changes in mass after at least one full time step + if (CS%Time > dTime) then + Time0 = CS%Time - dTime last_hmask(:,:) = ISS%hmask(:,:) ; last_area_shelf_h(:,:) = ISS%area_shelf_h(:,:) call time_interp_external(CS%id_read_mass, Time0, last_mass_shelf) ! This should only be done if time_interp_external did an update. @@ -1058,9 +1062,9 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! apply fluxes do j=js,je ; do i=is,ie - ! Note the following is hard coded for ISOMIP - if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then + if (in_sponge(i,j) > 0.0) then ! evap is negative, and vprec has units of [R Z T-1 ~> kg m-2 s-1] + !### Why does mean_melt_flux need to be rescaled to get vprec? fluxes%vprec(i,j) = -mean_melt_flux * CS%density_ice / (1000.0*US%kg_m3_to_R) fluxes%sens(i,j) = fluxes%vprec(i,j) * CS%Cp * CS%T0 ! [ Q R Z T-1 ~> W /m^2 ] fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * CS%S0*1.0e-3 ! [kgSalt/kg R Z T-1 ~> kgSalt m-2 s-1] @@ -1073,7 +1077,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) call MOM_forcing_chksum("After constant sea level", fluxes, G, CS%US, haloshift=0) endif - endif !constant_sea_level + endif ! constant_sea_level end subroutine add_shelf_flux @@ -1238,9 +1242,14 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "If true, user specifies a constant nondimensional heat-transfer coefficient "//& "(GAMMA_T_3EQ), from which the salt-transfer coefficient is then computed "//& " as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.) - if (CS%const_gamma) call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", CS%Gamma_T_3EQ, & - "Nondimensional heat-transfer coefficient.",default=2.2E-2, & - units="nondim.", fail_if_missing=.true.) + if (CS%const_gamma) then + call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", CS%Gamma_T_3EQ, & + "Nondimensional heat-transfer coefficient.", & + units="nondim", default=2.2e-2) + call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_S", CS%Gamma_S_3EQ, & + "Nondimensional salt-transfer coefficient.", & + default=CS%Gamma_T_3EQ/35.0, units="nondim") + endif call get_param(param_file, mdl, "ICE_SHELF_MASS_FROM_FILE", & CS%mass_from_file, "Read the mass of the "//& @@ -1252,14 +1261,13 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "is computed from a quadratic equation. Otherwise, the previous "//& "interactive method to estimate Sbdry is used.", default=.false.) if (CS%find_salt_root) then ! read liquidus coeffs. - call get_param(param_file, mdl, "TFREEZE_S0_P0", CS%lambda1, & + call get_param(param_file, mdl, "TFREEZE_S0_P0", CS%TFr_0_0, & "this is the freezing potential temperature at "//& "S=0, P=0.", units="degC", default=0.0, do_not_log=.true.) - call get_param(param_file, mdl, "DTFREEZE_DS", CS%lambda1, & !### This should be CS%lambda2! + call get_param(param_file, mdl, "DTFREEZE_DS", CS%dTFr_dS, & "this is the derivative of the freezing potential "//& - "temperature with salinity.", & - units="degC psu-1", default=-0.054, do_not_log=.true.) - call get_param(param_file, mdl, "DTFREEZE_DP", CS%lambda3, & + "temperature with salinity.", units="degC psu-1", default=-0.054, do_not_log=.true.) + call get_param(param_file, mdl, "DTFREEZE_DP", CS%dTFr_dp, & "this is the derivative of the freezing potential "//& "temperature with pressure.", & units="degC Pa-1", default=0.0, do_not_log=.true.) From 880da802441d95b4a78c11f6ad89a2cafee56b11 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 31 Mar 2020 12:04:34 -0400 Subject: [PATCH 05/13] (*+)Corrected a bug in setting ustar_shelf Corrected a bug in setting fluxes%ustar_shelf in shelf_calc_flux, that will change answers with an active ice shelf when UTIDE is nonzero. Also rescaled the units of utide in MOM_ice_shelf.F90 to [L T-1] and added a units argument to get_param calls for 5 ISOMIP or ice-shelf related variables. This commit can change answers and the parameter_doc files in some cases when a thermodynamically active ice shelf is used, but all answers in the MOM6-examples test cases are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 25 ++++++++++++------------- src/user/ISOMIP_initialization.F90 | 10 +++++----- 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index a4f6169844..6eb82f15d5 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -87,7 +87,7 @@ module MOM_ice_shelf type(ice_shelf_dyn_CS), pointer :: dCS => NULL() !< The control structure for the ice-shelf dynamics. real, pointer, dimension(:,:) :: & - utide => NULL() !< tidal velocity [m s-1] + utide => NULL() !< An unresolved tidal velocity [L T-1 ~> m s-1] real :: ustar_bg !< A minimum value for ustar under ice shelves [Z T-1 ~> m s-1]. real :: cdrag !< drag coefficient under ice shelves [nondim]. @@ -360,13 +360,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! Iteratively determine a self-consistent set of fluxes, with the ocean ! salinity just below the ice-shelf as the variable that is being ! iterated for. - ! ### SHOULD USTAR_SHELF BE SET YET? - !### I think that CS%utide**1 should be CS%utide**2 - ! Also I think that if taux_shelf and tauy_shelf have been calculated by the + ! ### SHOULD USTAR_SHELF BE SET YET, or should it be set from taux_shelf & tauy_shelf? + ! I think that if taux_shelf and tauy_shelf have been calculated by the ! ocean stress calculation, they should be used here or later to set ustar_shelf. - RWH - fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%m_to_Z*US%T_to_s * & - sqrt(CS%cdrag*((state%u(i,j)**2 + state%v(i,j)**2) + CS%utide(i,j)**1))) + fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_TO_Z * & + sqrt(CS%cdrag*(US%m_s_to_L_T**2*(state%u(i,j)**2 + state%v(i,j)**2) + CS%utide(i,j)**2))) ustar_h = fluxes%ustar_shelf(i,j) @@ -495,7 +494,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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? + ! ### SHOULD BOUNDS BE NEEDED IN THIS NEWTONS METHOD SOLVER? wB_flux_new = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB_in enddo !it3 endif @@ -1121,7 +1120,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl integer :: wd_halos(2) logical :: read_TideAmp, shelf_mass_is_dynamic, debug character(len=240) :: Tideamp_file - real :: utide + real :: utide ! A tidal velocity [L T-1 ~> m s-1] if (associated(CS)) then call MOM_error(FATAL, "MOM_ice_shelf.F90, initialize_ice_shelf: "// & "called with an associated control structure.") @@ -1218,7 +1217,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "(no conduction).", default=.false.) call get_param(param_file, mdl, "MELTING_CUTOFF_DEPTH", CS%cutoff_depth, & "Depth above which the melt is set to zero (it must be >= 0) "//& - "Default value won't affect the solution.", default=0.0, scale=US%m_to_Z) !###, units="m" + "Default value won't affect the solution.", units="m", default=0.0, scale=US%m_to_Z) if (CS%cutoff_depth < 0.) & call MOM_error(WARNING,"Initialize_ice_shelf: MELTING_CUTOFF_DEPTH must be >= 0.") @@ -1342,11 +1341,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) TideAmp_file = trim(inputdir) // trim(TideAmp_file) - call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, G%domain, timelevel=1) + call MOM_read_data(TideAmp_file, 'tideamp', CS%utide, G%domain, timelevel=1, scale=US%m_s_to_L_T) else call get_param(param_file, mdl, "UTIDE", utide, & "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & - units="m s-1", default=0.0) + units="m s-1", default=0.0 , scale=US%m_s_to_L_T) CS%utide(:,:) = utide endif @@ -1443,10 +1442,10 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "ice sheet/shelf thickness mask" ,"none") endif - ! if (CS%active_shelf_dynamics) then !### Consider adding an ice shelf dynamics switch. + if (CS%active_shelf_dynamics) then ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics call register_ice_shelf_dyn_restarts(G, param_file, CS%dCS, CS%restart_CSp) - ! endif + endif !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file !if (.not. CS%solo_ice_sheet) then diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index 1d14ff9cc5..aa7de04dac 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -172,14 +172,14 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read select case ( coordinateMode(verticalCoordinate) ) case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Initial thicknesses for isopycnal coordinates - call get_param(param_file, mdl, "ISOMIP_T_SUR",t_sur, & - 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read) + call get_param(param_file, mdl, "ISOMIP_T_SUR", t_sur, & + 'Temperature at the surface (interface)', units="degC", default=-1.9, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, & - 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read) + 'Salinity at the surface (interface)', units="ppt", default=33.8, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, & - 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read) + 'Temperature at the bottom (interface)', units="degC", default=-1.9, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot,& - 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read) + 'Salinity at the bottom (interface)', units="ppt", default=34.55, do_not_log=just_read) if (just_read) return ! All run-time parameters have been read, so return. From 1a8c75aa66807801409cbe2130abd50c4f46458b Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 31 Mar 2020 20:21:08 -0400 Subject: [PATCH 06/13] (*+)Set ice shelf latent heat consistently Set the latent heat of fusion and the heat capacity of water used by the ice shelf code consistently with the rest of MOM6, including their default values. This changes answers in all cases with active ice shelf thermodynamics. Also corrected a scaling factor for Rho0 and added several new chksum calls. Also added units for 6 ISOMIP-related input variables, and reordered the calls for several ice shelf parameters to make sure they are all being set when needed. This changes the solutions and the MOM_parameter_doc files in an updated ISOMIP test case but all answers in the MOM6-examples test cases are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 68 +++++++++++-------- .../vertical/MOM_diabatic_aux.F90 | 4 +- src/user/ISOMIP_initialization.F90 | 24 +++---- 3 files changed, 54 insertions(+), 42 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 6eb82f15d5..ab3d52c6ae 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -5,6 +5,7 @@ module MOM_ice_shelf ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_constants, only : hlf use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr @@ -351,7 +352,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! DNG - to allow this everywhere Hml>0.0 allows for melting under grounded cells ! propose instead to allow where Hml > [some threshold] - + !### I do not know what the Hml flag adds; consider removing it. if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. & (CS%isthermo) .and. (state%Hml(i,j) > 0.0) ) then @@ -961,6 +962,10 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) enddo ; enddo endif + if (CS%debug) then + call MOM_forcing_chksum("Before adding shelf fluxes", fluxes, G, CS%US, haloshift=0) + endif + do j=js,je ; do i=is,ie ; if (ISS%area_shelf_h(i,j) > 0.0) then frac_area = fluxes%frac_shelf_h(i,j) !### Should this be 1-frac_shelf_h? if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0 @@ -986,6 +991,12 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) fluxes%salt_flux(i,j) = frac_area * ISS%salt_flux(i,j)*CS%flux_factor endif ; enddo ; enddo + if (CS%debug) then + call hchksum(ISS%water_flux, "water_flux add shelf fluxes", G%HI, haloshift=0, scale=US%RZ_T_to_kg_m2s) + call hchksum(ISS%tflux_ocn, "tflux_ocn add shelf fluxes", G%HI, haloshift=0, scale=US%QRZ_T_to_W_m2) + call MOM_forcing_chksum("After adding shelf fluxes", fluxes, G, CS%US, haloshift=0) + endif + ! keep sea level constant by removing mass in the sponge ! region (via virtual precip, vprec). Apply additional ! salt/heat fluxes so that the resultant surface buoyancy @@ -1071,7 +1082,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) enddo ; enddo if (CS%debug) then - write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux, CS%time_step + write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux*US%RZ_T_to_kg_m2s, CS%time_step call MOM_mesg(mesg) call MOM_forcing_chksum("After constant sea level", fluxes, G, CS%US, haloshift=0) endif @@ -1177,8 +1188,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl isd = G%isd ; jsd = G%jsd ; ied = G%ied ; jed = G%jed Isdq = G%IsdB ; Iedq = G%IedB ; Jsdq = G%JsdB ; Jedq = G%JedB - !### This should be a run-time parameter that is read in consistently with MOM6 and SIS2. - CS%Lat_fusion = 3.34e5*US%J_kg_to_Q CS%override_shelf_movement = .false. ; CS%active_shelf_dynamics = .false. call log_version(param_file, mdl, version, "") @@ -1208,6 +1217,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "SHELF_THERMO", CS%isthermo, & "If true, use a thermodynamically interactive ice shelf.", & default=.false.) + call get_param(param_file, mdl, "LATENT_HEAT_FUSION", CS%Lat_fusion, & + "The latent heat of fusion.", units="J/kg", default=hlf, scale=US%J_kg_to_Q) call get_param(param_file, mdl, "SHELF_THREE_EQN", CS%threeeq, & "If true, use the three equation expression of "//& "consistency to calculate the fluxes at the ice-ocean "//& @@ -1223,25 +1234,36 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "CONST_SEA_LEVEL", CS%constant_sea_level, & "If true, apply evaporative, heat and salt fluxes in "//& - "the sponge region. This will avoid a large increase "//& + "the sponge region. This will avoid a large increase "//& "in sea level. This option is needed for some of the "//& "ISOMIP+ experiments (Ocean3 and Ocean4). "//& "IMPORTANT: it is not currently possible to do "//& "prefect restarts using this flag.", default=.false.) - call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", & - CS%S0, "Surface salinity in the resoring region.", & - default=33.8, do_not_log=.true.) + call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", CS%S0, & + "Surface salinity in the restoring region.", & + default=33.8, units='ppt', do_not_log=.true.) - call get_param(param_file, mdl, "ISOMIP_T_SUR_SPONGE", & - CS%T0, "Surface temperature in the resoring region.", & - default=-1.9, do_not_log=.true.) + call get_param(param_file, mdl, "ISOMIP_T_SUR_SPONGE", CS%T0, & + "Surface temperature in the restoring region.", & + default=-1.9, units='degC', do_not_log=.true.) call get_param(param_file, mdl, "SHELF_3EQ_GAMMA", CS%const_gamma, & "If true, user specifies a constant nondimensional heat-transfer coefficient "//& - "(GAMMA_T_3EQ), from which the salt-transfer coefficient is then computed "//& - " as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.) - if (CS%const_gamma) then + "(GAMMA_T_3EQ), from which the default salt-transfer coefficient is set "//& + "as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.) + if (CS%threeeq) then + call get_param(param_file, mdl, "SHELF_S_ROOT", CS%find_salt_root, & + "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) "//& + "is computed from a quadratic equation. Otherwise, the previous "//& + "interactive method to estimate Sbdry is used.", default=.false.) + else + call get_param(param_file, mdl, "SHELF_2EQ_GAMMA_T", CS%gamma_t, & + "If SHELF_THREE_EQN is false, this the fixed turbulent "//& + "exchange velocity at the ice-ocean interface.", & + units="m s-1", scale=US%m_to_Z*US%T_to_s, fail_if_missing=.true.) + endif + if (CS%const_gamma .or. CS%find_salt_root) then call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", CS%Gamma_T_3EQ, & "Nondimensional heat-transfer coefficient.", & units="nondim", default=2.2e-2) @@ -1254,11 +1276,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl CS%mass_from_file, "Read the mass of the "//& "ice shelf (every time step) from a file.", default=.false.) - if (CS%threeeq) & - call get_param(param_file, mdl, "SHELF_S_ROOT", CS%find_salt_root, & - "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) "//& - "is computed from a quadratic equation. Otherwise, the previous "//& - "interactive method to estimate Sbdry is used.", default=.false.) if (CS%find_salt_root) then ! read liquidus coeffs. call get_param(param_file, mdl, "TFREEZE_S0_P0", CS%TFr_0_0, & "this is the freezing potential temperature at "//& @@ -1272,24 +1289,19 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl units="degC Pa-1", default=0.0, do_not_log=.true.) endif - if (.not.CS%threeeq) & - call get_param(param_file, mdl, "SHELF_2EQ_GAMMA_T", CS%gamma_t, & - "If SHELF_THREE_EQN is false, this the fixed turbulent "//& - "exchange velocity at the ice-ocean interface.", & - units="m s-1", scale=US%m_to_Z*US%T_to_s, fail_if_missing=.true.) - call get_param(param_file, mdl, "G_EARTH", CS%g_Earth, & "The gravitational acceleration of the Earth.", & units="m s-2", default = 9.80, scale=US%m_to_Z*US%T_to_s**2) call get_param(param_file, mdl, "C_P", CS%Cp, & - "The heat capacity of sea water.", units="J kg-1 K-1", scale=US%J_kg_to_Q, & - fail_if_missing=.true.) + "The heat capacity of sea water, approximated as a constant. "//& + "The default value is from the TEOS-10 definition of conservative temperature.", & + units="J kg-1 K-1", default=3991.86795711963, scale=US%J_kg_to_Q) call get_param(param_file, mdl, "RHO_0", CS%Rho0, & "The mean ocean density used with BOUSSINESQ true to "//& "calculate accelerations and the mass for conservation "//& "properties, or with BOUSSINSEQ false to convert some "//& "parameters from vertical units of m to kg m-2.", & - units="kg m-3", default=1035.0, scale=US%R_to_kg_m3) !### MAKE THIS A SEPARATE PARAMETER. + units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) !### MAKE THIS A SEPARATE PARAMETER. call get_param(param_file, mdl, "C_P_ICE", CS%Cp_ice, & "The heat capacity of ice.", units="J kg-1 K-1", scale=US%J_kg_to_Q, & default=2.10e3) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 343423a221..b17f6e4323 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -931,7 +931,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t real, dimension(max(nsw,1),SZI_(G),SZK_(G)) :: & opacityBand ! The opacity (inverse of the exponential absorption length) of each frequency ! band of shortwave radation in each layer [H-1 ~> m-1 or m2 kg-1] - real, dimension(maxGroundings) :: hGrounding + real, dimension(maxGroundings) :: hGrounding ! Thickness added by each grounding event [H ~> m or kg m-2] real :: Temp_in, Salin_in real :: g_Hconv2 ! A conversion factor for use in the TKE calculation ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2]. @@ -1380,7 +1380,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t do i = 1, min(numberOfGroundings, maxGroundings) call forcing_SinglePointPrint(fluxes,G,iGround(i),jGround(i),'applyBoundaryFluxesInOut (grounding)') write(mesg(1:45),'(3es15.3)') G%geoLonT( iGround(i), jGround(i) ), & - G%geoLatT( iGround(i), jGround(i)) , hGrounding(i) + G%geoLatT( iGround(i), jGround(i)), hGrounding(i)*GV%H_to_m call MOM_error(WARNING, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//& "Mass created. x,y,dh= "//trim(mesg), all_print=.true.) enddo diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index aa7de04dac..ba8dc1162f 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -173,13 +173,13 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read case ( REGRIDDING_LAYER, REGRIDDING_RHO ) ! Initial thicknesses for isopycnal coordinates call get_param(param_file, mdl, "ISOMIP_T_SUR", t_sur, & - 'Temperature at the surface (interface)', units="degC", default=-1.9, do_not_log=just_read) + "Temperature at the surface (interface)", units="degC", default=-1.9, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, & - 'Salinity at the surface (interface)', units="ppt", default=33.8, do_not_log=just_read) + "Salinity at the surface (interface)", units="ppt", default=33.8, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, & - 'Temperature at the bottom (interface)', units="degC", default=-1.9, do_not_log=just_read) + "Temperature at the bottom (interface)", units="degC", default=-1.9, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot,& - 'Salinity at the bottom (interface)', units="ppt", default=34.55, do_not_log=just_read) + "Salinity at the bottom (interface)", units="ppt", default=34.55, do_not_log=just_read) if (just_read) return ! All run-time parameters have been read, so return. @@ -293,13 +293,13 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalCoordinate, & default=DEFAULT_COORDINATE_MODE, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_T_SUR",t_sur, & - 'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read) + "Temperature at the surface (interface)", units="degC", default=-1.9, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, & - 'Salinity at the surface (interface)', default=33.8, do_not_log=just_read) + "Salinity at the surface (interface)", units="ppt", default=33.8, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, & - 'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read) + "Temperature at the bottom (interface)", units="degC", default=-1.9, do_not_log=just_read) call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot, & - 'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read) + "Salinity at the bottom (interface)", units="ppt", default=34.55, do_not_log=just_read) call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state, scale=US%kg_m3_to_R) ! write(mesg,*) 'Density in the surface layer:', rho_sur @@ -481,16 +481,16 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) do_not_log=.true.) call get_param(PF, mdl, "ISOMIP_S_SUR_SPONGE", s_sur, & - 'Surface salinity in sponge layer.', default=s_ref) ! units="ppt") + "Surface salinity in sponge layer.", units="ppt", default=s_ref) ! units="ppt") call get_param(PF, mdl, "ISOMIP_S_BOT_SPONGE", s_bot, & - 'Bottom salinity in sponge layer.', default=s_ref) ! units="ppt") + "Bottom salinity in sponge layer.", units="ppt", default=s_ref) ! units="ppt") call get_param(PF, mdl, "ISOMIP_T_SUR_SPONGE", t_sur, & - 'Surface temperature in sponge layer.', default=t_ref) ! units="degC") + "Surface temperature in sponge layer.", units="degC", default=t_ref) ! units="degC") call get_param(PF, mdl, "ISOMIP_T_BOT_SPONGE", t_bot, & - 'Bottom temperature in sponge layer.', default=t_ref) ! units="degC") + "Bottom temperature in sponge layer.", units="degC", default=t_ref) ! units="degC") T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 !; RHO(:,:,:) = 0.0 From e88732c3eae3ee31a04bda1464e137651955f0e0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 1 Apr 2020 08:32:30 -0400 Subject: [PATCH 07/13] (*)Do not use Hml>0 as an ice shelf melt filter Avoid using the criterion that Hml>0 as a filter of when ice shelf melt can occur. I am not aware of a good justification for this filter, and it seems to be an historical artefact. Removing it causes melting to start one time step earlier at the start of a run, giving apparently better answers. This changes the answers in cases with a thermodynamically active ice shelf, including an ISOMIP test case, but all answers in the MOM6-examples test suite are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index ab3d52c6ae..13dd8940a3 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -350,12 +350,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! but it won't make a difference otherwise. fluxes%ustar_shelf(i,j)= 0.0 - ! DNG - to allow this everywhere Hml>0.0 allows for melting under grounded cells - ! propose instead to allow where Hml > [some threshold] - !### I do not know what the Hml flag adds; consider removing it. if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & - (ISS%area_shelf_h(i,j) > 0.0) .and. & - (CS%isthermo) .and. (state%Hml(i,j) > 0.0) ) then + (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo ) then if (CS%threeeq) then ! Iteratively determine a self-consistent set of fluxes, with the ocean @@ -602,8 +598,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) do j=js,je ; do i=is,ie if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & - (ISS%area_shelf_h(i,j) > 0.0) .and. & - (CS%isthermo) .and. (state%Hml(i,j) > 0.0) ) then + (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then ! Set melt to zero above a cutoff pressure (CS%Rho0*CS%cutoff_depth*CS%g_Earth). ! This is needed for the ISOMIP test case. From e87c66427f66dd9d944c44e35ccdf09e174f2a64 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 1 Apr 2020 08:50:23 -0400 Subject: [PATCH 08/13] (*)Removed rescaling from mean_melt_flux to vprec Removed inappropriate density ratio rescaling from the conversion of mean_melt_flux to fluxes%vprec when CONST_SEA_LEVEL is true. Both are cast as mass fluxes, not thicknesses fluxes, so this ratio is not needed. This will change answers in some regional cases with thermodynamically active ice shelves, but all answers in the MOM6-examples test suite are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 50 +++++++++++++++++---------------- 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 13dd8940a3..aa20a7ccb4 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -351,7 +351,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) fluxes%ustar_shelf(i,j)= 0.0 if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & - (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo ) then + (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then if (CS%threeeq) then ! Iteratively determine a self-consistent set of fluxes, with the ocean @@ -485,13 +485,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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) < 1e-4*(abs(wB_flux_new) + abs(wB_flux))) exit + ! Find the root where wB_flux_new = wB_flux. Make the 1.0e-4 below into a parameter? + if (abs(wB_flux_new - wB_flux) < 1.0e-4*(abs(wB_flux_new) + abs(wB_flux))) exit 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 IN THIS NEWTONS METHOD SOLVER? + ! 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 enddo !it3 endif @@ -500,13 +499,12 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) exch_vel_t(i,j) = ustar_h * I_Gam_T exch_vel_s(i,j) = ustar_h * I_Gam_S - !Calculate the heat flux inside the ice shelf. - - !vertical adv/diff as in H+J 1999, eqns (26) & approx from (31). - ! Q_ice = density_ice * CS%Cp_ice * K_ice * dT/dz (at interface) - !vertical adv/diff as in H+J 199, eqs (31) & (26)... - ! dT/dz ~= min( (lprec/(density_ice*K_ice))*(CS%Temp_Ice-T_freeze) , 0.0 ) - !If this approximation is not made, iterations are required... See H+J Fig 3. + ! Calculate the heat flux inside the ice shelf. + ! Vertical adv/diff as in H+J 1999, eqns (26) & approx from (31). + ! Q_ice = density_ice * CS%Cp_ice * K_ice * dT/dz (at interface) + ! vertical adv/diff as in H+J 1999, eqs (31) & (26)... + ! dT/dz ~= min( (lprec/(density_ice*K_ice))*(CS%Temp_Ice-T_freeze) , 0.0 ) + ! If this approximation is not made, iterations are required... See H+J Fig 3. if (ISS%tflux_ocn(i,j) >= 0.0) then ! Freezing occurs due to downward ocean heat flux, so zero iout ce heat flux. @@ -548,19 +546,17 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) if (dS_it < 0.0) then ! Sbdry is now the upper bound. if (Sb_max_set .and. (Sbdry(i,j) > Sb_max)) & - call MOM_error(FATAL,"shelf_calc_flux: Irregular iteration for Sbdry (max).") + call MOM_error(FATAL,"shelf_calc_flux: Irregular iteration for Sbdry (max).") Sb_max = Sbdry(i,j) ; dS_max = dS_it ; Sb_max_set = .true. else ! Sbdry is now the lower bound. if (Sb_min_set .and. (Sbdry(i,j) < Sb_min)) & - call MOM_error(FATAL, & - "shelf_calc_flux: Irregular iteration for Sbdry (min).") - Sb_min = Sbdry(i,j) ; dS_min = dS_it ; Sb_min_set = .true. + call MOM_error(FATAL, "shelf_calc_flux: Irregular iteration for Sbdry (min).") + Sb_min = Sbdry(i,j) ; dS_min = dS_it ; Sb_min_set = .true. endif ! dS_it < 0.0 if (Sb_min_set .and. Sb_max_set) then ! Use the false position method for the next iteration. - Sbdry(i,j) = Sb_min + (Sb_max-Sb_min) * & - (dS_min / (dS_min - dS_max)) + Sbdry(i,j) = Sb_min + (Sb_max-Sb_min) * (dS_min / (dS_min - dS_max)) else Sbdry(i,j) = Sbdry_it endif ! Sb_min_set @@ -569,7 +565,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) endif ! CS%find_salt_root enddo !it1 - ! Check for non-convergence and/or non-boundedness? + ! Check for non-convergence and/or non-boundedness? else ! In the 2-equation form, the mixed layer turbulent exchange velocity @@ -584,7 +580,9 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ISS%water_flux(i,j) = -I_LF * ISS%tflux_ocn(i,j) Sbdry(i,j) = 0.0 endif - else !not shelf + elseif (ISS%area_shelf_h(i,j) > 0.0) then ! This is an ice-sheet, not a floating shelf. + ISS%tflux_ocn(i,j) = 0.0 + else ! There is no ice shelf or sheet here. ISS%tflux_ocn(i,j) = 0.0 endif @@ -598,7 +596,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) do j=js,je ; do i=is,ie if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & - (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then + (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then ! Set melt to zero above a cutoff pressure (CS%Rho0*CS%cutoff_depth*CS%g_Earth). ! This is needed for the ISOMIP test case. @@ -627,8 +625,13 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) write(mesg,*) "|melt| = ",fluxes%iceshelf_melt(i,j)," > 0 and ustar_shelf = 0. at i,j", i, j call MOM_error(FATAL, "shelf_calc_flux: "//trim(mesg)) endif - endif ! area_shelf_h !!!!!!!!!!!!!!!!!!!!!!!!!!!!End of safety checks !!!!!!!!!!!!!!!!!!! + elseif (ISS%area_shelf_h(i,j) > 0.0) then + ! This is grounded ice, that could be modified to melt if a geothermal heat flux were used. + haline_driving(i,j) = 0.0 + ISS%water_flux(i,j) = 0.0 + fluxes%iceshelf_melt(i,j) = 0.0 + endif ! area_shelf_h enddo ; enddo ! i- and j-loops ! mass flux [kg s-1], part of ISOMIP diags. @@ -1069,8 +1072,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) do j=js,je ; do i=is,ie if (in_sponge(i,j) > 0.0) then ! evap is negative, and vprec has units of [R Z T-1 ~> kg m-2 s-1] - !### Why does mean_melt_flux need to be rescaled to get vprec? - fluxes%vprec(i,j) = -mean_melt_flux * CS%density_ice / (1000.0*US%kg_m3_to_R) + fluxes%vprec(i,j) = -mean_melt_flux fluxes%sens(i,j) = fluxes%vprec(i,j) * CS%Cp * CS%T0 ! [ Q R Z T-1 ~> W /m^2 ] fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * CS%S0*1.0e-3 ! [kgSalt/kg R Z T-1 ~> kgSalt m-2 s-1] endif From fcb1e92a12eedd7315448a9a8c605e9da81e0b30 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 1 Apr 2020 12:09:04 -0400 Subject: [PATCH 09/13] +Add optional arg to ice_shelf_min_thickness_calve Added a new optional argument, halo, to ice_shelf_min_thickness_calve to specify the range of indices over which to work. All answers are bitwise identical, but there is a new argument in a public interface. --- src/ice_shelf/MOM_ice_shelf_dynamics.F90 | 31 ++++++++++++++---------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index ca8faf55f3..be3ae1ecde 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -1650,7 +1650,7 @@ subroutine shelf_advance_front(CS, ISS, G, hmask, uh_ice, vh_ice) end subroutine shelf_advance_front !> Apply a very simple calving law using a minimum thickness rule -subroutine ice_shelf_min_thickness_calve(G, h_shelf, area_shelf_h, hmask, thickness_calve) +subroutine ice_shelf_min_thickness_calve(G, h_shelf, area_shelf_h, hmask, thickness_calve, halo) type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: h_shelf !< The ice shelf thickness [Z ~> m]. real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: area_shelf_h !< The area per cell covered by @@ -1658,20 +1658,25 @@ subroutine ice_shelf_min_thickness_calve(G, h_shelf, area_shelf_h, hmask, thickn real, dimension(SZDI_(G),SZDJ_(G)), intent(inout) :: hmask !< A mask indicating which tracer points are !! partly or fully covered by an ice-shelf real, intent(in) :: thickness_calve !< The thickness at which to trigger calving [Z ~> m]. + integer, optional, intent(in) :: halo !< The number of halo points to use. If not present, + !! work on the entire data domain. + integer :: i, j, is, ie, js, je - integer :: i,j + if (present(halo)) then + is = G%isc - halo ; ie = G%iec + halo ; js = G%jsc - halo ; je = G%jec + halo + else + is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed + endif - do j=G%jsd,G%jed - do i=G%isd,G%ied -! if ((h_shelf(i,j) < CS%thickness_calve) .and. (hmask(i,j) == 1) .and. & -! (CS%ground_frac(i,j) == 0.0)) then - if ((h_shelf(i,j) < thickness_calve) .and. (area_shelf_h(i,j) > 0.)) then - h_shelf(i,j) = 0.0 - area_shelf_h(i,j) = 0.0 - hmask(i,j) = 0.0 - endif - enddo - enddo + do j=js,je ; do i=is,ie +! if ((h_shelf(i,j) < CS%thickness_calve) .and. (hmask(i,j) == 1) .and. & +! (CS%ground_frac(i,j) == 0.0)) then + if ((h_shelf(i,j) < thickness_calve) .and. (area_shelf_h(i,j) > 0.)) then + h_shelf(i,j) = 0.0 + area_shelf_h(i,j) = 0.0 + hmask(i,j) = 0.0 + endif + enddo ; enddo end subroutine ice_shelf_min_thickness_calve From 04d5a83dcb93dce3bc22d7ceeeddaf303d203870 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 1 Apr 2020 12:10:23 -0400 Subject: [PATCH 10/13] (*)Use a blend of ice-shelf and open water fluxes Use a blend of ice-shelf and open water fluxes when there is partial cover by an ice shelf. This will change answers in some cases with temporally evolving ice shelves, but all answers in the MOM6-examples test suite are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 74 +++++++++++++++++++-------------- 1 file changed, 42 insertions(+), 32 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index aa20a7ccb4..ea9162afc5 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -881,7 +881,8 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! local variables real :: Irho0 !< The inverse of the mean density times unit conversion factors that !! arise because state uses MKS units [Z2 m s2 kg-1 T-2 ~> m3 kg-1]. - real :: frac_area !< The fractional area covered by the ice shelf [nondim]. + 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]. real :: delta_mass_shelf!< Change in ice shelf mass over one time step [kg s-1] real :: taux2, tauy2 !< The squared surface stresses [Pa]. real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa]. @@ -956,7 +957,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) if (CS%active_shelf_dynamics .or. CS%override_shelf_movement) then do j=jsd,jed ; do i=isd,ied if (G%areaT(i,j) > 0.0) & - fluxes%frac_shelf_h(i,j) = ISS%area_shelf_h(i,j) * G%IareaT(i,j) + fluxes%frac_shelf_h(i,j) = min(1.0, ISS%area_shelf_h(i,j) * G%IareaT(i,j)) enddo ; enddo endif @@ -965,28 +966,32 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) endif do j=js,je ; do i=is,ie ; if (ISS%area_shelf_h(i,j) > 0.0) then - frac_area = fluxes%frac_shelf_h(i,j) !### Should this be 1-frac_shelf_h? - if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0 - if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = 0.0 - if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = 0.0 - if (associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir(i,j) = 0.0 - if (associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif(i,j) = 0.0 - if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0 - if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0 - if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0 + ! Replace fluxes intercepted by the ice shelf with fluxes from the ice shelf + frac_shelf = min(1.0, ISS%area_shelf_h(i,j) * G%IareaT(i,j)) + frac_open = max(0.0, 1.0 - frac_shelf) + + if (associated(fluxes%sw)) fluxes%sw(i,j) = frac_open * fluxes%sw(i,j) + if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = frac_open * fluxes%sw_vis_dir(i,j) + if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = frac_open * fluxes%sw_vis_dif(i,j) + if (associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir(i,j) = frac_open * fluxes%sw_nir_dir(i,j) + if (associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif(i,j) = frac_open * fluxes%sw_nir_dif(i,j) + if (associated(fluxes%lw)) fluxes%lw(i,j) = frac_open * fluxes%lw(i,j) + if (associated(fluxes%latent)) fluxes%latent(i,j) = frac_open * fluxes%latent(i,j) + if (associated(fluxes%evap)) fluxes%evap(i,j) = frac_open * fluxes%evap(i,j) if (associated(fluxes%lprec)) then if (ISS%water_flux(i,j) > 0.0) then - fluxes%lprec(i,j) = frac_area*ISS%water_flux(i,j)*CS%flux_factor + fluxes%lprec(i,j) = frac_shelf*ISS%water_flux(i,j)*CS%flux_factor + frac_open * fluxes%lprec(i,j) else - fluxes%lprec(i,j) = 0.0 - fluxes%evap(i,j) = frac_area*ISS%water_flux(i,j)*CS%flux_factor + fluxes%lprec(i,j) = frac_open * fluxes%lprec(i,j) + fluxes%evap(i,j) = fluxes%evap(i,j) + frac_shelf*ISS%water_flux(i,j)*CS%flux_factor endif endif if (associated(fluxes%sens)) & - fluxes%sens(i,j) = frac_area*ISS%tflux_ocn(i,j)*CS%flux_factor + fluxes%sens(i,j) = frac_shelf*ISS%tflux_ocn(i,j)*CS%flux_factor + frac_open * fluxes%sens(i,j) + ! The salt flux should be mostly from sea ice, so perhaps none should be intercepted and this should be changed. if (associated(fluxes%salt_flux)) & - fluxes%salt_flux(i,j) = frac_area * ISS%salt_flux(i,j)*CS%flux_factor + fluxes%salt_flux(i,j) = frac_shelf * ISS%salt_flux(i,j)*CS%flux_factor + frac_open * fluxes%salt_flux(i,j) endif ; enddo ; enddo if (CS%debug) then @@ -1008,16 +1013,6 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je)) fluxes%salt_flux(:,:) = 0.0 ; fluxes%vprec(:,:) = 0.0 - do j=js,je ; do i=is,ie - - !### These hard-coded limits need to be corrected. They are inappropriate here. - if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then - in_sponge(i,j) = 1.0 - else - in_sponge(i,j) = 0.0 - endif - enddo ; enddo - ! take into account changes in mass (or thickness) when imposing ice shelf mass if (CS%override_shelf_movement .and. CS%mass_from_file) then dTime = real_to_time(CS%time_step) @@ -1025,18 +1020,24 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! Compute changes in mass after at least one full time step if (CS%Time > dTime) then Time0 = CS%Time - dTime - last_hmask(:,:) = ISS%hmask(:,:) ; last_area_shelf_h(:,:) = ISS%area_shelf_h(:,:) + do j=js,je ; do i=is,ie + last_hmask(i,j) = ISS%hmask(i,j) ; last_area_shelf_h(i,j) = ISS%area_shelf_h(i,j) + enddo ; enddo call time_interp_external(CS%id_read_mass, Time0, last_mass_shelf) + do j=js,je ; do i=is,ie ! This should only be done if time_interp_external did an update. - last_mass_shelf(:,:) = US%kg_m3_to_R*US%m_to_Z * last_mass_shelf(:,:) ! Rescale after time_interp - last_h_shelf(:,:) = last_mass_shelf(:,:) / CS%density_ice + last_mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * last_mass_shelf(i,j) ! Rescale after time_interp + last_h_shelf(i,j) = last_mass_shelf(i,j) / CS%density_ice + enddo ; enddo ! apply calving if (CS%min_thickness_simple_calve > 0.0) then call ice_shelf_min_thickness_calve(G, last_h_shelf, last_area_shelf_h, last_hmask, & - CS%min_thickness_simple_calve) + CS%min_thickness_simple_calve, halo=0) ! convert to mass again - last_mass_shelf(:,:) = last_h_shelf(:,:) * CS%density_ice + do j=js,je ; do i=is,ie + last_mass_shelf(i,j) = last_h_shelf(i,j) * CS%density_ice + enddo ; enddo endif ! get total ice shelf mass at (Time-dt) and (Time), in kg @@ -1059,6 +1060,15 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) endif ! average total melt flux over sponge area + do j=js,je ; do i=is,ie + !### These hard-coded limits need to be corrected. They are inappropriate here. + if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then + in_sponge(i,j) = 1.0 + else + in_sponge(i,j) = 0.0 + endif + enddo ; enddo + sponge_area = global_area_integral(in_sponge, G) if (sponge_area > 0.0) then mean_melt_flux = US%kg_m2s_to_RZ_T*(global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & @@ -1754,7 +1764,7 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) if (CS%min_thickness_simple_calve > 0.0) then call ice_shelf_min_thickness_calve(G, ISS%h_shelf, ISS%area_shelf_h, ISS%hmask, & - CS%min_thickness_simple_calve) + CS%min_thickness_simple_calve, halo=0) endif call pass_var(ISS%area_shelf_h, G%domain) From 684681b2b34403b7b293767beba7504ab7184ffe Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 1 Apr 2020 14:58:18 -0400 Subject: [PATCH 11/13] +Ice shelf code cleanup Added a new run-time parameter to specify how much water has to be under an ice shelf for it to float; this is only logged when the CONST_SEA_LEVEL option is true. The description of another parameter is corrected, which changes the MOM_parameter_doc files with ice shelf thermodynamics. In addition, merged the mass_shelf into the same loop as h_shelf in change_thickness_using_melt and reduced the loop extents for setting forces%p_surf. Use column masses instead of thicknesses for thresholds in MOM_ice_shelf. Also combined CS%Rho0 and CS%density_ocean_avg as CS%Rho_ocn in ice_shelf_CS. All answers are bitwise identical, but there are changes in output files in some cases. --- src/ice_shelf/MOM_ice_shelf.F90 | 73 ++++++++++++++++----------------- 1 file changed, 35 insertions(+), 38 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index ea9162afc5..cd9903bc82 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -94,7 +94,7 @@ module MOM_ice_shelf real :: cdrag !< drag coefficient under ice shelves [nondim]. real :: g_Earth !< The gravitational acceleration [Z T-2 ~> m s-2] real :: Cp !< The heat capacity of sea water [Q degC-1 ~> J kg-1 degC-1]. - real :: Rho0 !< A reference ocean density [R ~> kg m-3]. + real :: Rho_ocn !< A reference ocean density [R ~> kg m-3]. real :: Cp_ice !< The heat capacity of fresh ice [Q degC-1 ~> J kg-1 degC-1]. real :: gamma_t !< The (fixed) turbulent exchange velocity in the !< 2-equation formulation [Z T-1 ~> m s-1]. @@ -109,7 +109,8 @@ module MOM_ice_shelf 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 :: col_thick_melt_threshold !< if the mixed layer is below this threshold, melt rate + real :: col_mass_melt_threshold !< An ocean column mass below the iceshelf below which melting + !! does not occur [kg m-2] logical :: mass_from_file !< Read the ice shelf mass from a file every dt !!!! PHYSICAL AND NUMERICAL PARAMETERS FOR ICE DYNAMICS !!!!!! @@ -127,10 +128,6 @@ module MOM_ice_shelf !!determined by ocean column thickness means update_OD_ffrac !! will be called (note: GL_regularize and GL_couple !! should be exclusive) - real :: density_ocean_avg !< this does not affect ocean circulation OR thermodynamics - !! it is to estimate the gravitational driving force at the - !! shelf front (until we think of a better way to do it, - !! but any difference will be negligible) logical :: calve_to_mask !< If true, calve any ice that passes outside of a masked area real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m]. real :: T0 !< temperature at ocean surface in the restoring region [degC] @@ -154,6 +151,9 @@ module MOM_ice_shelf logical :: const_gamma !< If true, gamma_T is specified by the user. logical :: constant_sea_level !< if true, apply an evaporative, heat and salt !! fluxes. It will avoid large increase in sea level. + real :: min_ocean_mass_float !< The minimum ocean mass per unit area before the ice + !! shelf is considered to float when constant_sea_level + !! is used [kg m-2] real :: cutoff_depth !< Depth above which melt is set to zero (>= 0) [Z ~> m]. logical :: find_salt_root !< If true, if true find Sbdry using a quadratic eq. real :: TFr_0_0 !< The freezing point at 0 pressure and 0 salinity [degC] @@ -270,7 +270,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) real :: dS_min, dS_max ! Variables used in iterating for wB_flux. real :: wB_flux_new, dDwB_dwB_in - real :: I_Gam_T, I_Gam_S, iDens + real :: I_Gam_T, I_Gam_S real :: dG_dwB ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2] real :: Isqrt2 logical :: Sb_min_set, Sb_max_set @@ -296,15 +296,13 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) SC = CS%kv_molec/CS%kd_molec_salt PR = CS%kv_molec/CS%kd_molec_temp I_VK = 1.0/VK - RhoCp = CS%Rho0 * CS%Cp + RhoCp = CS%Rho_ocn * CS%Cp Isqrt2 = 1.0/sqrt(2.0) !first calculate molecular component Gam_mol_t = 12.5 * (PR**c2_3) - 6 Gam_mol_s = 12.5 * (SC**c2_3) - 6 - iDens = 1.0/CS%density_ocean_avg - ! GMM, zero some fields of the ice shelf structure (ice_shelf_CS) ! these fields are already set to zero during initialization ! However, they seem to be changed somewhere and, for diagnostic @@ -350,7 +348,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! but it won't make a difference otherwise. fluxes%ustar_shelf(i,j)= 0.0 - if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & + if ((state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then if (CS%threeeq) then @@ -370,7 +368,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then ! ! These arrays are supposed to be stress components at C-grid points, which is ! ! inconsistent with what is coded up here. - ! state%taux_shelf(i,j) = US%RZ_T_to_kg_m2s*US%Z_to_m*US%s_to_T * ustar_h**2 * CS%Rho0*Isqrt2 + ! state%taux_shelf(i,j) = US%RZ_T_to_kg_m2s*US%Z_to_m*US%s_to_T * ustar_h**2 * CS%Rho_ocn*Isqrt2 ! state%tauy_shelf(i,j) = state%taux_shelf(i,j) ! endif @@ -537,7 +535,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) exit ! no need to do interaction, so exit loop else - mass_exch = exch_vel_s(i,j) * CS%Rho0 + mass_exch = exch_vel_s(i,j) * CS%Rho_ocn Sbdry_it = (state%sss(i,j) * mass_exch + CS%Salin_ice * ISS%water_flux(i,j)) / & (mass_exch + ISS%water_flux(i,j)) dS_it = Sbdry_it - Sbdry(i,j) @@ -595,17 +593,17 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) fluxes%iceshelf_melt(:,:) = ISS%water_flux(:,:) * CS%flux_factor do j=js,je ; do i=is,ie - if ((iDens*state%ocean_mass(i,j) > CS%col_thick_melt_threshold) .and. & + if ((state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then - ! Set melt to zero above a cutoff pressure (CS%Rho0*CS%cutoff_depth*CS%g_Earth). + ! Set melt to zero above a cutoff pressure (CS%Rho_ocn*CS%cutoff_depth*CS%g_Earth). ! This is needed for the ISOMIP test case. - if (ISS%mass_shelf(i,j) < CS%Rho0*CS%cutoff_depth) then + if (ISS%mass_shelf(i,j) < CS%Rho_ocn*CS%cutoff_depth) then ISS%water_flux(i,j) = 0.0 fluxes%iceshelf_melt(i,j) = 0.0 endif ! Compute haline driving, which is one of the diags. used in ISOMIP - haline_driving(i,j) = (ISS%water_flux(i,j) * Sbdry(i,j)) / (CS%Rho0 * exch_vel_s(i,j)) + haline_driving(i,j) = (ISS%water_flux(i,j) * Sbdry(i,j)) / (CS%Rho_ocn * exch_vel_s(i,j)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!! !1)Check if haline_driving computed above is consistent with @@ -739,20 +737,13 @@ subroutine change_thickness_using_melt(ISS, G, US, time_step, fluxes, density_ic ISS%hmask(i,j) = 0.0 ISS%area_shelf_h(i,j) = 0.0 endif + ISS%mass_shelf(i,j) = ISS%h_shelf(i,j) * density_ice endif enddo ; enddo call pass_var(ISS%area_shelf_h, G%domain) call pass_var(ISS%h_shelf, G%domain) call pass_var(ISS%hmask, G%domain) - - !### combine this with the loops above. - do j=G%jsd,G%jed ; do i=G%isd,G%ied - if ((ISS%hmask(i,j) == 1) .or. (ISS%hmask(i,j) == 2)) then - ISS%mass_shelf(i,j) = ISS%h_shelf(i,j) * density_ice - endif - enddo ; enddo - call pass_var(ISS%mass_shelf, G%domain) end subroutine change_thickness_using_melt @@ -801,8 +792,7 @@ subroutine add_shelf_forces(G, US, CS, forces, do_shelf_area) call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, G%domain, TO_ALL, CGRID_NE) endif - !### Consider working over a smaller array range. - do j=jsd,jed ; do i=isd,ied + do j=js,je ; do i=is,ie press_ice = (ISS%area_shelf_h(i,j) * G%IareaT(i,j)) * & US%RZ_to_kg_m2*US%Z_to_m*US%s_to_T**2*(CS%g_Earth * ISS%mass_shelf(i,j)) if (associated(forces%p_surf)) then @@ -935,7 +925,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) call pass_vector(state%taux_shelf, state%tauy_shelf, G%domain, TO_ALL, CGRID_NE) ! GMM: melting is computed using ustar_shelf (and not ustar), which has already ! been passed, I so believe we do not need to update fluxes%ustar. -! Irho0 = US%m_to_Z*US%T_to_s*US%kg_m2s_to_RZ_T / CS%Rho0 +! Irho0 = US%m_to_Z*US%T_to_s*US%kg_m2s_to_RZ_T / CS%Rho_ocn ! do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then ! ### THIS SHOULD BE AN AREA WEIGHTED AVERAGE OF THE ustar_shelf POINTS. ! taux2 = 0.0 ; tauy2 = 0.0 @@ -969,7 +959,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! Replace fluxes intercepted by the ice shelf with fluxes from the ice shelf frac_shelf = min(1.0, ISS%area_shelf_h(i,j) * G%IareaT(i,j)) frac_open = max(0.0, 1.0 - frac_shelf) - + if (associated(fluxes%sw)) fluxes%sw(i,j) = frac_open * fluxes%sw(i,j) if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = frac_open * fluxes%sw_vis_dir(i,j) if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = frac_open * fluxes%sw_vis_dif(i,j) @@ -1042,8 +1032,8 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! get total ice shelf mass at (Time-dt) and (Time), in kg do j=js,je ; do i=is,ie - ! just floating shelf (0.1 is a threshold for min ocean thickness) - if (((1.0/CS%density_ocean_avg)*state%ocean_mass(i,j) > 0.1) .and. & + ! Just consider the change in the mass of the floating shelf. + if ((state%ocean_mass(i,j) > CS%min_ocean_mass_float) .and. & (ISS%area_shelf_h(i,j) > 0.0)) then delta_float_mass(i,j) = ISS%mass_shelf(i,j) - last_mass_shelf(i,j) else @@ -1127,6 +1117,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl real :: L_rescale ! A rescaling factor for horizontal lengths from the representation in ! a restart file to the internal representation in this run. real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic. + real :: dz_ocean_min_float ! The minimum ocean thickness above which the ice shelf is considered + ! to be floating when CONST_SEA_LEVEL = True [m]. real :: cdrag, drag_bg_vel logical :: new_sim, save_IC, var_force !This include declares and sets the variable "version". @@ -1139,6 +1131,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl logical :: read_TideAmp, shelf_mass_is_dynamic, debug character(len=240) :: Tideamp_file 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 [m] if (associated(CS)) then call MOM_error(FATAL, "MOM_ice_shelf.F90, initialize_ice_shelf: "// & "called with an associated control structure.") @@ -1246,6 +1240,10 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "ISOMIP+ experiments (Ocean3 and Ocean4). "//& "IMPORTANT: it is not currently possible to do "//& "prefect restarts using this flag.", default=.false.) + call get_param(param_file, mdl, "MIN_OCEAN_FLOAT_THICK", dz_ocean_min_float, & + "The minimum ocean thickness above which the ice shelf is considered to be "//& + "floating when CONST_SEA_LEVEL = True.", & + default=0.1, units="m", do_not_log=.not.CS%constant_sea_level) call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", CS%S0, & "Surface salinity in the restoring region.", & @@ -1303,15 +1301,16 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "The heat capacity of sea water, approximated as a constant. "//& "The default value is from the TEOS-10 definition of conservative temperature.", & units="J kg-1 K-1", default=3991.86795711963, scale=US%J_kg_to_Q) - call get_param(param_file, mdl, "RHO_0", CS%Rho0, & + call get_param(param_file, mdl, "RHO_0", CS%Rho_ocn, & "The mean ocean density used with BOUSSINESQ true to "//& "calculate accelerations and the mass for conservation "//& "properties, or with BOUSSINSEQ false to convert some "//& "parameters from vertical units of m to kg m-2.", & - units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) !### MAKE THIS A SEPARATE PARAMETER. + units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "C_P_ICE", CS%Cp_ice, & "The heat capacity of ice.", units="J kg-1 K-1", scale=US%J_kg_to_Q, & default=2.10e3) + if (CS%constant_sea_level) CS%min_ocean_mass_float = dz_ocean_min_float*CS%Rho_ocn call get_param(param_file, mdl, "ICE_SHELF_FLUX_FACTOR", CS%flux_factor, & "Non-dimensional factor applied to shelf thermodynamic "//& @@ -1334,17 +1333,15 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "KD_TEMP_MOLECULAR", CS%kd_molec_temp, & "The molecular diffusivity of heat in sea water at the "//& "freezing point.", units="m2 s-1", default=1.41e-7, scale=US%m2_s_to_Z2_T) - call get_param(param_file, mdl, "RHO_0", CS%density_ocean_avg, & - "avg ocean density used in floatation cond", & - units="kg m-3", default=1035.) call get_param(param_file, mdl, "DT_FORCING", CS%time_step, & "The time step for changing forcing, coupling with other "//& "components, or potentially writing certain diagnostics. "//& "The default value is given by DT.", units="s", default=0.0) - call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", CS%col_thick_melt_threshold, & - "The minimum ML thickness where melting is allowed.", units="m", & + call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", col_thick_melt_thresh, & + "The minimum ocean column thickness where melting is allowed.", units="m", & default=0.0) + CS%col_mass_melt_threshold = CS%Rho_ocn * col_thick_melt_thresh call get_param(param_file, mdl, "READ_TIDEAMP", read_TIDEAMP, & "If true, read a file (given by TIDEAMP_FILE) containing "//& From 7cdffb103b4344899025498569f5b5e0ac23efe1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 1 Apr 2020 18:12:25 -0400 Subject: [PATCH 12/13] (*)Use state%taux_shelf to set state%ustar_shelf If possible, use state%taux_shelf and %tauy_shelf to set state%ustar_shelf, or if these are not available use the (now appropriately staggered) state%u and state%v to set ustar%shelf. In addition, ustar%shelf is set in all cases with a thermodynamically interactive ice shelf, and not just those that use the 3-equation expressions for the skin salinity. In addition, when CONST_SEA_LEVEL is true, the balancing flux occurs over all open ocean area, although the previous mode of using a hard-coded region is still there in commented out code. These code changes alter answers in all cases with a thermodynamically interactive ice shelf, but the solutions in the MOM6-examples test cases are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 147 +++++++++++++++----------------- 1 file changed, 68 insertions(+), 79 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index cd9903bc82..733245b1ce 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -110,7 +110,7 @@ module MOM_ice_shelf real :: Gamma_S_3EQ !< Nondimensional salt-transfer coefficient, used in the 3Eq. formulation !< 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 [kg m-2] + !! does not occur [R Z ~> kg m-2] logical :: mass_from_file !< Read the ice shelf mass from a file every dt !!!! PHYSICAL AND NUMERICAL PARAMETERS FOR ICE DYNAMICS !!!!!! @@ -153,7 +153,7 @@ module MOM_ice_shelf !! fluxes. It will avoid large increase in sea level. real :: min_ocean_mass_float !< The minimum ocean mass per unit area before the ice !! shelf is considered to float when constant_sea_level - !! is used [kg m-2] + !! is used [R Z ~> kg m-2] real :: cutoff_depth !< Depth above which melt is set to zero (>= 0) [Z ~> m]. logical :: find_salt_root !< If true, if true find Sbdry using a quadratic eq. real :: TFr_0_0 !< The freezing point at 0 pressure and 0 salinity [degC] @@ -272,7 +272,13 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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 :: Isqrt2 + real :: taux2, tauy2 ! The squared surface stresses [Pa]. + 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 :: 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 unit conversion factors that + ! arise because state uses MKS units [L2 m s2 kg-1 T-2 ~> m3 kg-1]. logical :: Sb_min_set, Sb_max_set logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities. logical :: coupled_GL ! If true, the grouding line position is determined based on @@ -297,7 +303,6 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) PR = CS%kv_molec/CS%kd_molec_temp I_VK = 1.0/VK RhoCp = CS%Rho_ocn * CS%Cp - Isqrt2 = 1.0/sqrt(2.0) !first calculate molecular component Gam_mol_t = 12.5 * (PR**c2_3) - 6 @@ -332,6 +337,37 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) call hchksum(state%ocean_mass, "ocean_mass before apply melting", G%HI, haloshift=0) endif + ! Calculate the friction velocity under ice shelves, using taux_shelf and tauy_shelf if possible. + if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then + call pass_vector(state%taux_shelf, state%tauy_shelf, G%domain, TO_ALL, CGRID_NE) + endif + Irho0 = US%m_s_to_L_T**2*US%kg_m3_to_R / CS%Rho_ocn + do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then + taux2 = 0.0 ; tauy2 = 0.0 ; u2_av = 0.0 ; v2_av = 0.0 + asu1 = (ISS%area_shelf_h(i-1,j) + ISS%area_shelf_h(i,j)) + asu2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i+1,j)) + asv1 = (ISS%area_shelf_h(i,j-1) + ISS%area_shelf_h(i,j)) + asv2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i,j+1)) + I_au = 0.0 ; if (asu1 + asu2 > 0.0) I_au = 1.0 / (asu1 + asu2) + I_av = 0.0 ; if (asv1 + asv2 > 0.0) I_av = 1.0 / (asv1 + asv2) + if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then + taux2 = (asu1 * state%taux_shelf(I-1,j)**2 + asu2 * state%taux_shelf(I,j)**2 ) * I_au + tauy2 = (asv1 * state%tauy_shelf(i,J-1)**2 + asv2 * state%tauy_shelf(i,J)**2 ) * I_av + endif + u2_av = US%m_s_to_L_T**2*(asu1 * state%u(I-1,j)**2 + asu2 * state%u(I,j)**2) * I_au + v2_av = US%m_s_to_L_T**2*(asv1 * state%v(i,J-1)**2 + asu2 * state%v(i,J)**2) * I_av + + if (taux2 + tauy2 > 0.0) then + fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_to_Z * & + sqrt(Irho0 * sqrt(taux2 + tauy2) + CS%cdrag*CS%utide(i,j)**2)) + 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))) + endif + else ! There is no shelf here. + fluxes%ustar_shelf(i,j) = 0.0 + endif ; enddo ; enddo + do j=js,je ! Find the pressure at the ice-ocean interface, averaged only over the ! part of the cell covered by ice shelf. @@ -344,11 +380,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) dR0_dT, dR0_dS, is, ie-is+1, CS%eqn_of_state) do i=is,ie - ! set ustar_shelf to zero. This is necessary if shelf_mass_is_dynamic - ! but it won't make a difference otherwise. - fluxes%ustar_shelf(i,j)= 0.0 - - if ((state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & + if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. CS%isthermo) then if (CS%threeeq) then @@ -356,22 +388,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! salinity just below the ice-shelf as the variable that is being ! iterated for. - ! ### SHOULD USTAR_SHELF BE SET YET, or should it be set from taux_shelf & tauy_shelf? - ! I think that if taux_shelf and tauy_shelf have been calculated by the - ! ocean stress calculation, they should be used here or later to set ustar_shelf. - RWH - fluxes%ustar_shelf(i,j) = MAX(CS%ustar_bg, US%L_TO_Z * & - sqrt(CS%cdrag*(US%m_s_to_L_T**2*(state%u(i,j)**2 + state%v(i,j)**2) + CS%utide(i,j)**2))) - ustar_h = fluxes%ustar_shelf(i,j) - ! I think that the following can be deleted without causing any problems. - ! if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then - ! ! These arrays are supposed to be stress components at C-grid points, which is - ! ! inconsistent with what is coded up here. - ! state%taux_shelf(i,j) = US%RZ_T_to_kg_m2s*US%Z_to_m*US%s_to_T * ustar_h**2 * CS%Rho_ocn*Isqrt2 - ! state%tauy_shelf(i,j) = state%taux_shelf(i,j) - ! endif - ! Estimate the neutral ocean boundary layer thickness as the minimum of the ! reported ocean mixed layer thickness and the neutral Ekman depth. absf = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & @@ -593,7 +611,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) fluxes%iceshelf_melt(:,:) = ISS%water_flux(:,:) * CS%flux_factor do j=js,je ; do i=is,ie - if ((state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & + if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then ! Set melt to zero above a cutoff pressure (CS%Rho_ocn*CS%cutoff_depth*CS%g_Earth). @@ -869,20 +887,16 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated. ! local variables - real :: Irho0 !< The inverse of the mean density times unit conversion factors that - !! arise because state uses MKS units [Z2 m s2 kg-1 T-2 ~> m3 kg-1]. 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]. real :: delta_mass_shelf!< Change in ice shelf mass over one time step [kg s-1] - real :: taux2, tauy2 !< The squared surface stresses [Pa]. real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa]. - real :: asu1, asu2 !< Ocean areas covered by ice shelves at neighboring u- - real :: asv1, asv2 !< and v-points [L2 ~> m2]. - real :: mean_melt_flux !< Spatial mean melt flux [R Z T-1 ~> kg m-2 s-1] - real :: sponge_area !< total area of sponge region [m2] + real :: balancing_flux !< The fresh water flux that balances the integrated melt flux [R Z T-1 ~> kg m-2 s-1] + real :: balancing_area !< total area where the balancing flux is applied [m2] type(time_type) :: dTime !< The time step as a time_type type(time_type) :: Time0 !< The previous time (Time-dt) - real, dimension(SZDI_(G),SZDJ_(G)) :: in_sponge !< 1 where the property damping occurs, 0 otherwise [nondim] + real, dimension(SZDI_(G),SZDJ_(G)) :: bal_frac !< Fraction of the cel1 where the mass flux + !! balancing the net melt flux occurs, 0 to 1 [nondim] real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass !! at at previous time (Time-dt) [R Z ~> kg m-2] real, dimension(SZDI_(G),SZDJ_(G)) :: delta_float_mass !< The change in the floating mass between @@ -921,29 +935,6 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) endif endif - if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then - call pass_vector(state%taux_shelf, state%tauy_shelf, G%domain, TO_ALL, CGRID_NE) - ! GMM: melting is computed using ustar_shelf (and not ustar), which has already - ! been passed, I so believe we do not need to update fluxes%ustar. -! Irho0 = US%m_to_Z*US%T_to_s*US%kg_m2s_to_RZ_T / CS%Rho_ocn -! do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then - ! ### THIS SHOULD BE AN AREA WEIGHTED AVERAGE OF THE ustar_shelf POINTS. - ! taux2 = 0.0 ; tauy2 = 0.0 - ! asu1 = (ISS%area_shelf_h(i-1,j) + ISS%area_shelf_h(i,j)) - ! asu2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i+1,j)) - ! asv1 = (ISS%area_shelf_h(i,j-1) + ISS%area_shelf_h(i,j)) - ! asv2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i,j+1)) - ! if ((asu1 + asu2 > 0.0) .and. associated(state%taux_shelf)) & - ! taux2 = (asu1 * state%taux_shelf(I-1,j)**2 + & - ! asu2 * state%taux_shelf(I,j)**2 ) / (asu1 + asu2) - ! if ((asv1 + asv2 > 0.0) .and. associated(state%tauy_shelf)) & - ! tauy2 = (asv1 * state%tauy_shelf(i,J-1)**2 + & - ! asv2 * state%tauy_shelf(i,J)**2 ) / (asv1 + asv2) - - ! fluxes%ustar(i,j) = MAX(CS%ustar_bg, sqrt(Irho0 * sqrt(taux2 + tauy2))) -! endif ; enddo ; enddo - endif - if (CS%active_shelf_dynamics .or. CS%override_shelf_movement) then do j=jsd,jed ; do i=isd,ied if (G%areaT(i,j) > 0.0) & @@ -990,15 +981,12 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) call MOM_forcing_chksum("After adding shelf fluxes", fluxes, G, CS%US, haloshift=0) endif - ! keep sea level constant by removing mass in the sponge - ! region (via virtual precip, vprec). Apply additional - ! salt/heat fluxes so that the resultant surface buoyancy - ! forcing is ~ 0. + ! Keep sea level constant by removing mass via a balancing flux that might be applied + ! in the open ocean or the sponge region (via virtual precip, vprec). Apply additional + ! salt/heat fluxes so that the resultant surface buoyancy forcing is ~ 0. ! This is needed for some of the ISOMIP+ experiments. if (CS%constant_sea_level) then - !### This code has problems with hard coded constants that need to be refactored. -RWH - if (.not. associated(fluxes%salt_flux)) allocate(fluxes%salt_flux(ie,je)) if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je)) fluxes%salt_flux(:,:) = 0.0 ; fluxes%vprec(:,:) = 0.0 @@ -1033,7 +1021,7 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! get total ice shelf mass at (Time-dt) and (Time), in kg do j=js,je ; do i=is,ie ! Just consider the change in the mass of the floating shelf. - if ((state%ocean_mass(i,j) > CS%min_ocean_mass_float) .and. & + if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%min_ocean_mass_float) .and. & (ISS%area_shelf_h(i,j) > 0.0)) then delta_float_mass(i,j) = ISS%mass_shelf(i,j) - last_mass_shelf(i,j) else @@ -1051,35 +1039,36 @@ subroutine add_shelf_flux(G, US, CS, state, fluxes) ! average total melt flux over sponge area do j=js,je ; do i=is,ie - !### These hard-coded limits need to be corrected. They are inappropriate here. - if (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then - in_sponge(i,j) = 1.0 + if ((G%mask2dT(i,j) > 0.0) .AND. (ISS%area_shelf_h(i,j) * G%IareaT(i,j) < 1.0)) then + ! Uncomment this for some ISOMIP cases: + ! .AND. (G%geoLonT(i,j) >= 790.0) .AND. (G%geoLonT(i,j) <= 800.0)) then + bal_frac(i,j) = max(1.0 - ISS%area_shelf_h(i,j) * G%IareaT(i,j), 0.0) else - in_sponge(i,j) = 0.0 + bal_frac(i,j) = 0.0 endif enddo ; enddo - sponge_area = global_area_integral(in_sponge, G) - if (sponge_area > 0.0) then - mean_melt_flux = US%kg_m2s_to_RZ_T*(global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & + balancing_area = global_area_integral(bal_frac, G) + if (balancing_area > 0.0) then + balancing_flux = US%kg_m2s_to_RZ_T*(global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & area=ISS%area_shelf_h) + & - delta_mass_shelf ) / sponge_area + delta_mass_shelf ) / balancing_area else - mean_melt_flux = 0.0 + balancing_flux = 0.0 endif ! apply fluxes do j=js,je ; do i=is,ie - if (in_sponge(i,j) > 0.0) then + if (bal_frac(i,j) > 0.0) then ! evap is negative, and vprec has units of [R Z T-1 ~> kg m-2 s-1] - fluxes%vprec(i,j) = -mean_melt_flux + fluxes%vprec(i,j) = -balancing_flux fluxes%sens(i,j) = fluxes%vprec(i,j) * CS%Cp * CS%T0 ! [ Q R Z T-1 ~> W /m^2 ] fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * CS%S0*1.0e-3 ! [kgSalt/kg R Z T-1 ~> kgSalt m-2 s-1] endif enddo ; enddo if (CS%debug) then - write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux*US%RZ_T_to_kg_m2s, CS%time_step + write(mesg,*) 'Balancing flux (kg/(m^2 s)), dt = ', balancing_flux*US%RZ_T_to_kg_m2s, CS%time_step call MOM_mesg(mesg) call MOM_forcing_chksum("After constant sea level", fluxes, G, CS%US, haloshift=0) endif @@ -1118,7 +1107,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl ! a restart file to the internal representation in this run. real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic. real :: dz_ocean_min_float ! The minimum ocean thickness above which the ice shelf is considered - ! to be floating when CONST_SEA_LEVEL = True [m]. + ! to be floating when CONST_SEA_LEVEL = True [Z ~> m]. real :: cdrag, drag_bg_vel logical :: new_sim, save_IC, var_force !This include declares and sets the variable "version". @@ -1243,7 +1232,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl call get_param(param_file, mdl, "MIN_OCEAN_FLOAT_THICK", dz_ocean_min_float, & "The minimum ocean thickness above which the ice shelf is considered to be "//& "floating when CONST_SEA_LEVEL = True.", & - default=0.1, units="m", do_not_log=.not.CS%constant_sea_level) + default=0.1, units="m", scale=US%m_to_Z, do_not_log=.not.CS%constant_sea_level) call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", CS%S0, & "Surface salinity in the restoring region.", & @@ -1339,8 +1328,8 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "The default value is given by DT.", units="s", default=0.0) call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", col_thick_melt_thresh, & - "The minimum ocean column thickness where melting is allowed.", units="m", & - default=0.0) + "The minimum ocean column thickness where melting is allowed.", & + units="m", scale=US%m_to_Z, default=0.0) CS%col_mass_melt_threshold = CS%Rho_ocn * col_thick_melt_thresh call get_param(param_file, mdl, "READ_TIDEAMP", read_TIDEAMP, & @@ -1389,7 +1378,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl units="m", default=0.0, scale=US%m_to_Z) call get_param(param_file, mdl, "USTAR_SHELF_BG", CS%ustar_bg, & - "The minimum value of ustar under ice sheves.", & + "The minimum value of ustar under ice shelves.", & units="m s-1", default=0.0, scale=US%m_to_Z*US%T_to_s) call get_param(param_file, mdl, "CDRAG_SHELF", cdrag, & "CDRAG is the drag coefficient relating the magnitude of "//& From 43ebd34bb1f79042f787a06b13ade68f4f375217 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 2 Apr 2020 14:55:00 -0400 Subject: [PATCH 13/13] Removed spaces from a blank line --- src/framework/MOM_spatial_means.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90 index d4b687b0a5..2423a19433 100644 --- a/src/framework/MOM_spatial_means.F90 +++ b/src/framework/MOM_spatial_means.F90 @@ -56,7 +56,7 @@ function global_area_integral(var, G, scale, area) real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: area !< The alternate area to use, including !! any required masking [L2 ~> m2]. real :: global_area_integral !< The returned area integral, usually in the units of var times [m2]. - + ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming real :: scalefac ! An overall scaling factor for the areas and variable.