From 45cd5c644b6450d921b07e687e69f14f8879c0bd Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 20 Jul 2023 22:40:44 -0400 Subject: [PATCH 1/8] Move file parser inquire calls to root PE MOM_file_parser's open_param_file() contains explicit inquire() calls when assessing the correctness of opening such a file. As written, these could be called by any rank, and are not thread safe. In rare cases (usually related to testing), this would cause a race condition and raise an error. Even ignoring the errors, it is probably better if only one rank makes these calls, rather than all of them. The following patch modifies the function so that only root PE invokes inquire(). There is not much to celebrate about this patch; it does not try to clean up the intrinsic weirdness of the IO handling. But it does appear to fix some of the most apparent problems. --- src/framework/MOM_file_parser.F90 | 40 ++++++++++++++++++------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index 35d75cff7f..88fabf0c74 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -152,28 +152,34 @@ subroutine open_param_file(filename, CS, checkable, component, doc_file_dir) ! Check that this file has not already been opened if (CS%nfiles > 0) then reopened_file = .false. - inquire(file=trim(filename), number=iounit) - if (iounit /= -1) then - do i = 1, CS%nfiles - if (CS%iounit(i) == iounit) then - call assert(trim(CS%filename(1)) == trim(filename), & - "open_param_file: internal inconsistency! "//trim(filename)// & - " is registered as open but has the wrong unit number!") - call MOM_error(WARNING, & - "open_param_file: file "//trim(filename)// & - " has already been opened. This should NOT happen!"// & - " Did you specify the same file twice in a namelist?") - reopened_file = .true. - endif ! unit numbers - enddo ! i + + if (is_root_pe()) then + inquire(file=trim(filename), number=iounit) + if (iounit /= -1) then + do i = 1, CS%nfiles + if (CS%iounit(i) == iounit) then + call assert(trim(CS%filename(1)) == trim(filename), & + "open_param_file: internal inconsistency! "//trim(filename)// & + " is registered as open but has the wrong unit number!") + call MOM_error(WARNING, & + "open_param_file: file "//trim(filename)// & + " has already been opened. This should NOT happen!"// & + " Did you specify the same file twice in a namelist?") + reopened_file = .true. + endif ! unit numbers + enddo ! i + endif endif + if (any_across_PEs(reopened_file)) return endif ! Check that the file exists to readstdlog - inquire(file=trim(filename), exist=file_exists) - if (.not.file_exists) call MOM_error(FATAL, & - "open_param_file: Input file '"// trim(filename)//"' does not exist.") + if (is_root_pe()) then + inquire(file=trim(filename), exist=file_exists) + if (.not.file_exists) call MOM_error(FATAL, & + "open_param_file: Input file '"// trim(filename)//"' does not exist.") + endif Netcdf_file = .false. if (strlen > 3) then From f01d256ccce458ce16ca31544b6726bca5645002 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 3 Aug 2023 11:22:19 -0400 Subject: [PATCH 2/8] +Non-Boussinesq revisions to wave_interface The surface gravity waves invariably work with depths in meters to match the units of the horizontal wavenumbers, even in non-Boussinesq mode, so several of the variables that are passed to or used in the MOM_wave_interface module have been revised to use units of [Z ~> m] and not [H ~> m or kg m-2]. There are changes to the names or units of a total of 8 arguments to Stokes_PGF, Update_Stokes_Drift, Get_StokesSL_LiFoxKemper, Get_SL_Average_Prof and get_Langmuir_Number, and a new argument to StokesMixing. To accommodate these changes, there are now calls to thickness_to_dz that precede the Update_Stokes_Drift or Stokes_PGF calls and a new 3d array in step_MOM. Additionally, four hard-coded dimensional constants in Stokes_PGF were given the needed scaling factors to pass dimensional consistency testing. This change also required the addition of a unit_scale_type argument to Stokes_PGF, and corresponding changes to the calls to this routine from step_MOM_dyn_split_RK2. Three comments pointing out probable bugs or instances of inaccurate algorithms were also added. Incorrect units were also fixed in the comments describing one internal variable in Update_Stokes_Drift. Also added the new runtime parameter RHO_SFC_WAVES to set the surface seawater density that is used in comparison with the typical density of air set by RHO_AIR to estimate the 10-meter wind speed from the ocean friction velocity and the COARE 3.5 stress relationships inside of get_StokesSL_LiFoxKemper() when the Li and Fox-Kemper (2017) wave model (i.e. LF17) is used. As a part of this change, there is a new verticalGrid_type argument to the internal routine set_LF17_wave_params(). By default, RHO_SFC_WAVES is set to RHO_0. Other specific changes include: - Changed the units of the publicly visible KvS element of the wave_parameters_CS from [Z2 T-1] to [H Z T-1] - Replaced the layer thickness argument (h with units of [H ~> m or kg m-2]) to Update_Stokes_Drift(), Get_SL_Average_Prof(), get_Langmuir_Number() and Stokes_PGF() with a vertical layer extent argument (dz in units of [Z ~> m]) - Add a vertical layer extent argument (dz) to StokesMixing() - Add a unit_scale_type argument to Stokes_PGF() - Multiply a hard-coded length scale in Stokes_PGF by the correct unit scaling factor and added comments to highlight these hard-coded lengths - Corrected the unit description of the DecayScale internal variable - Changed the units of 3 internal variables and added 1 new internal variable in StokesMixing - Added two new arrays to step_MOM for the layer vertical extent and friction velocity being set via thickness_to_dz and find_ustar and being passed to Update_Stokes_Drift - Added new arrays for the layer vertical extent and calls to thickness_to_dz in KPP_compute_BLD and energetic_PBL, as well as a new dz argument to ePBL_column to provide an altered argument to get_Langmuir_Number - Use find_ustar to set the friction velocity passed to Update_Stokes_Drift - Update fluxes%ustar or fluxes%tau_mag halos as necessary when the waves are in use A total of 28 unit conversion factors were eliminated with these changes. All answers are bitwise identical in Boussinesq mode, but they are changed in non-Boussinesq mode by the use of the layer specific volume to convert between thicknesses and vertical extents. The units of several arguments to publicly visible routines are altered as is a publicly visible element in wave_parameters_CS and there is a new runtime parameter that appears in some MOM_parameter_doc files. --- src/core/MOM.F90 | 26 +++- src/core/MOM_dynamics_split_RK2.F90 | 6 +- .../vertical/MOM_CVMix_KPP.F90 | 19 ++- .../vertical/MOM_energetic_PBL.F90 | 26 ++-- src/user/MOM_wave_interface.F90 | 143 ++++++++++-------- 5 files changed, 132 insertions(+), 88 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2f001305d1..9cbb744560 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -34,7 +34,7 @@ module MOM use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : read_param, get_param, log_version, param_file_type -use MOM_forcing_type, only : forcing, mech_forcing +use MOM_forcing_type, only : forcing, mech_forcing, find_ustar use MOM_forcing_type, only : MOM_forcing_chksum, MOM_mech_forcing_chksum use MOM_get_input, only : Get_MOM_Input, directories use MOM_io, only : MOM_io_init, vardesc, var_desc @@ -91,7 +91,7 @@ module MOM use MOM_grid, only : set_first_direction, rescale_grid_bathymetry use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_hor_index, only : rotate_hor_index -use MOM_interface_heights, only : find_eta, calc_derived_thermo +use MOM_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz use MOM_interface_filter, only : interface_filter, interface_filter_init, interface_filter_end use MOM_interface_filter, only : interface_filter_CS use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end @@ -544,8 +544,13 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! the end of a stepping cycle (whatever that may mean). logical :: therm_reset ! If true, reset running sums of thermodynamic quantities. real :: cycle_time ! The length of the coupled time-stepping cycle [T ~> s]. + real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: & + U_star ! The wind friction velocity, calculated using the Boussinesq reference density or + ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1] real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: & ssh ! sea surface height, which may be based on eta_av [Z ~> m] + real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%GV)) :: & + dz ! Vertical distance across layers [Z ~> m] real, dimension(:,:,:), pointer :: & u => NULL(), & ! u : zonal velocity component [L T-1 ~> m s-1] @@ -672,13 +677,18 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS dt = time_interval / real(n_max) dt_therm = dt ; ntstep = 1 + + if (CS%UseWaves .and. associated(fluxes%ustar)) & + call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass, halo=1) + if (CS%UseWaves .and. associated(fluxes%tau_mag)) & + call pass_var(fluxes%tau_mag, G%Domain, clock=id_clock_pass, halo=1) + if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf CS%tv%p_surf => NULL() if (CS%use_p_surf_in_EOS .and. associated(fluxes%p_surf)) then CS%tv%p_surf => fluxes%p_surf if (allocated(CS%tv%SpV_avg)) call pass_var(fluxes%p_surf, G%Domain, clock=id_clock_pass) endif - if (CS%UseWaves) call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass) endif if (therm_reset) then @@ -722,12 +732,16 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS if (CS%UseWaves) then ! Update wave information, which is presently kept static over each call to step_mom call enable_averages(time_interval, Time_start + real_to_time(US%T_to_s*time_interval), CS%diag) - call Update_Stokes_Drift(G, GV, US, Waves, h, forces%ustar, time_interval, do_dyn) + call find_ustar(forces, CS%tv, U_star, G, GV, US, halo=1) + call thickness_to_dz(h, CS%tv, dz, G, GV, US, halo_size=1) + call Update_Stokes_Drift(G, GV, US, Waves, dz, U_star, time_interval, do_dyn) call disable_averaging(CS%diag) endif else ! not do_dyn. if (CS%UseWaves) then ! Diagnostics are not enabled in this call. - call Update_Stokes_Drift(G, GV, US, Waves, h, fluxes%ustar, time_interval, do_dyn) + call find_ustar(fluxes, CS%tv, U_star, G, GV, US, halo=1) + call thickness_to_dz(h, CS%tv, dz, G, GV, US, halo_size=1) + call Update_Stokes_Drift(G, GV, US, Waves, dz, U_star, time_interval, do_dyn) endif endif @@ -3261,7 +3275,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS) ! Write initial conditions if (CS%write_IC) then allocate(restart_CSp_tmp) - restart_CSP_tmp = CS%restart_CS + restart_CSp_tmp = CS%restart_CS call restart_registry_lock(restart_CSp_tmp, unlocked=.true.) allocate(z_interface(SZI_(G),SZJ_(G),SZK_(GV)+1)) call find_eta(CS%h, CS%tv, G, GV, US, z_interface, dZref=G%Z_ref) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 5ce9ec8962..eebb7d6b8a 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -481,7 +481,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s Use_Stokes_PGF = associated(Waves) if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF if (Use_Stokes_PGF) then - call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) + call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) ! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv ! will therefore report the sum total PGF and we avoid other @@ -748,7 +749,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s Use_Stokes_PGF = associated(Waves) if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF if (Use_Stokes_PGF) then - call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) + call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1) + call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves) if (.not.Waves%Passive_Stokes_PGF) then do k=1,nz do j=js,je ; do I=Isq,Ieq diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 5e56098c98..d24c3e2954 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -12,6 +12,7 @@ module MOM_CVMix_KPP use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_grid, only : ocean_grid_type, isPointInCell +use MOM_interface_heights, only : thickness_to_dz use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -604,7 +605,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, & type(ocean_grid_type), intent(in) :: G !< Ocean grid type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3] @@ -863,14 +864,14 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, & Kt(i,j,k) = Kt(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,1) Ks(i,j,k) = Ks(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,2) Kv(i,j,k) = Kv(i,j,k) + GV%m2_s_to_HZ_T * Kviscosity(k) - if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = GV%H_to_Z*Kv(i,j,k) + if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = Kv(i,j,k) enddo else ! KPP replaces prior diffusivity when former is non-zero do k=1, GV%ke+1 if (Kdiffusivity(k,1) /= 0.) Kt(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,1) if (Kdiffusivity(k,2) /= 0.) Ks(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,2) if (Kviscosity(k) /= 0.) Kv(i,j,k) = GV%m2_s_to_HZ_T * Kviscosity(k) - if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = GV%H_to_Z*Kv(i,j,k) + if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = Kv(i,j,k) enddo endif endif @@ -912,7 +913,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl type(ocean_grid_type), intent(inout) :: G !< Ocean grid type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Temp !< potential/cons temp [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Salt !< Salinity [S ~> ppt] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Velocity i-component [L T-1 ~> m s-1] @@ -924,6 +925,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement factor [nondim] ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m] + ! Variables for passing to CVMix routines, often in MKS units real, dimension( GV%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars in MKS units [m s-1] real, dimension( GV%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3] @@ -940,7 +943,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real :: Coriolis ! Coriolis parameter at tracer points in MKS units [s-1] real :: KPP_OBL_depth ! Boundary layer depth calculated by CVMix_kpp_compute_OBL_depth in MKS units [m] - ! Variables for EOS calculations real, dimension( 3*GV%ke ) :: rho_1D ! A column of densities [R ~> kg m-3] real, dimension( 3*GV%ke ) :: pres_1D ! A column of pressures [R L2 T-2 ~> Pa] @@ -996,6 +998,9 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl GoRho = US%Z_to_m*US%s_to_T**2 * GoRho_Z_L2 buoy_scale = US%L_to_m**2*US%s_to_T**3 + ! Find the vertical distances across layers. + call thickness_to_dz(h, tv, dz, G, GV, US) + ! loop over horizontal points on processor !$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, & !$OMP surfBuoyFlux, U_H, V_H, Coriolis, pRef, SLdepth_0d, vt2_1d, & @@ -1005,7 +1010,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl !$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_guess, LA, rho_1D, & !$OMP deltarho, N2_1d, ws_1d, LangEnhVT2,KPP_OBL_depth, z_cell, & !$OMP z_inter, OBL_depth, BulkRi_1d, zBottomMinusOffset) & - !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, & + !$OMP shared(G, GV, CS, US, uStar, h, dz, buoy_scale, buoyFlux, & !$OMP Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult) do j = G%jsc, G%jec do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then @@ -1127,7 +1132,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if ( (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) .and. .not. present(lamult)) then MLD_guess = max( CS%MLD_guess_min, abs(CS%OBLdepthprev(i,j) ) ) call get_Langmuir_Number(LA, G, GV, US, MLD_guess, uStar(i,j), i, j, & - H=H(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES) + dz=dz(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES) CS%La_SL(i,j) = LA endif diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 1556955424..06c3915d84 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -12,6 +12,7 @@ module MOM_energetic_PBL use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : thickness_to_dz use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs @@ -309,6 +310,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Local variables real, dimension(SZI_(G),SZK_(GV)) :: & h_2d, & ! A 2-d slice of the layer thickness [H ~> m or kg m-2]. + dz_2d, & ! A 2-d slice of the vertical distance across layers [Z ~> m]. T_2d, & ! A 2-d slice of the layer temperatures [C ~> degC]. S_2d, & ! A 2-d slice of the layer salinities [S ~> ppt]. TKE_forced_2d, & ! A 2-d slice of TKE_forced [R Z3 T-2 ~> J m-2]. @@ -320,6 +322,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS Kd_2d ! A 2-d version of the diapycnal diffusivity [Z2 T-1 ~> m2 s-1]. real, dimension(SZK_(GV)) :: & h, & ! The layer thickness [H ~> m or kg m-2]. + dz, & ! The vertical distance across layers [Z ~> m]. T0, & ! The initial layer temperatures [C ~> degC]. S0, & ! The initial layer salinities [S ~> ppt]. dSV_dT_1d, & ! The partial derivatives of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]. @@ -362,7 +365,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Zero out diagnostics before accumulation. if (CS%TKE_diagnostics) then -!!OMP parallel do default(none) shared(is,ie,js,je,CS) + !!OMP parallel do default(none) shared(is,ie,js,je,CS) do j=js,je ; do i=is,ie CS%diag_TKE_wind(i,j) = 0.0 ; CS%diag_TKE_MKE(i,j) = 0.0 CS%diag_TKE_conv(i,j) = 0.0 ; CS%diag_TKE_forcing(i,j) = 0.0 @@ -373,8 +376,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! if (CS%id_Mixing_Length>0) CS%Mixing_Length(:,:,:) = 0.0 ! if (CS%id_Velocity_Scale>0) CS%Velocity_Scale(:,:,:) = 0.0 -!!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt, & -!!OMP CS,G,GV,US,fluxes,TKE_forced,dSV_dT,dSV_dS,Kd_int) + !!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt, & + !!OMP CS,G,GV,US,fluxes,TKE_forced,dSV_dT,dSV_dS,Kd_int) do j=js,je ! Copy the thicknesses and other fields to 2-d arrays. do k=1,nz ; do i=is,ie @@ -383,6 +386,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS TKE_forced_2d(i,k) = TKE_forced(i,j,k) dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; dSV_dS_2d(i,k) = dSV_dS(i,j,k) enddo ; enddo + call thickness_to_dz(h_3d, tv, dz_2d, j, G, GV) ! Determine the initial mech_TKE and conv_PErel, including the energy required ! to mix surface heating through the topmost cell, the energy released by mixing @@ -394,7 +398,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Copy the thicknesses and other fields to 1-d arrays. do k=1,nz - h(k) = h_2d(i,k) + GV%H_subroundoff ; u(k) = u_2d(i,k) ; v(k) = v_2d(i,k) + h(k) = h_2d(i,k) + GV%H_subroundoff ; dz(k) = dz_2d(i,k) + GV%dZ_subroundoff + u(k) = u_2d(i,k) ; v(k) = v_2d(i,k) T0(k) = T_2d(i,k) ; S0(k) = S_2d(i,k) ; TKE_forcing(k) = TKE_forced_2d(i,k) dSV_dT_1d(k) = dSV_dT_2d(i,k) ; dSV_dS_1d(k) = dSV_dS_2d(i,k) enddo @@ -421,15 +426,15 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Perhaps provide a first guess for MLD based on a stored previous value. MLD_io = -1.0 - if (CS%MLD_iteration_guess .and. (CS%ML_Depth(i,j) > 0.0)) MLD_io = CS%ML_Depth(i,j) + if (CS%MLD_iteration_guess .and. (CS%ML_depth(i,j) > 0.0)) MLD_io = CS%ML_depth(i,j) if (stoch_CS%pert_epbl) then ! stochastics are active - call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, & US, CS, eCD, Waves, G, i, j, & TKE_gen_stoch=stoch_CS%epbl1_wts(i,j), TKE_diss_stoch=stoch_CS%epbl2_wts(i,j)) else - call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, & u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, & US, CS, eCD, Waves, G, i, j) endif @@ -499,12 +504,13 @@ end subroutine energetic_PBL !> This subroutine determines the diffusivities from the integrated energetics !! mixed layer model for a single column of water. -subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, & +subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, & u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, & Waves, G, i, j, TKE_gen_stoch, TKE_diss_stoch) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + real, dimension(SZK_(GV)), intent(in) :: dz !< The vertical distance across layers [Z ~> m]. real, dimension(SZK_(GV)), intent(in) :: u !< Zonal velocities interpolated to h points !! [L T-1 ~> m s-1]. real, dimension(SZK_(GV)), intent(in) :: v !< Zonal velocities interpolated to h points @@ -828,7 +834,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !/ Here we get MStar, which is the ratio of convective TKE driven mixing to UStar**3 MLD_guess_z = GV%H_to_Z*MLD_guess ! Convert MLD from thickness to height coordinates for these calls if (CS%Use_LT) then - call get_Langmuir_Number(LA, G, GV, US, abs(MLD_guess_z), u_star_mean, i, j, h, Waves, & + call get_Langmuir_Number(LA, G, GV, US, abs(MLD_guess_z), u_star_mean, i, j, dz, Waves, & U_H=u, V_H=v) call find_mstar(CS, US, B_flux, u_star, u_star_Mean, MLD_guess_z, absf, & MStar_total, Langmuir_Number=La, Convect_Langmuir_Number=LAmod,& @@ -1931,7 +1937,7 @@ subroutine energetic_PBL_get_MLD(CS, MLD, G, US, m_to_MLD_units) scale = 1.0 ; if (present(m_to_MLD_units)) scale = US%Z_to_m * m_to_MLD_units do j=G%jsc,G%jec ; do i=G%isc,G%iec - MLD(i,j) = scale*CS%ML_Depth(i,j) + MLD(i,j) = scale*CS%ML_depth(i,j) enddo ; enddo end subroutine energetic_PBL_get_MLD diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 321528b739..580e293f4f 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -97,7 +97,7 @@ module MOM_wave_interface !! Horizontal -> V points !! Vertical -> Mid-points real, allocatable, dimension(:,:,:), public :: & - KvS !< Viscosity for Stokes Drift shear [Z2 T-1 ~> m2 s-1] + KvS !< Viscosity for Stokes Drift shear [H Z T-1 ~> m2 s-1 or Pa s] ! The remainder of this control structure is private integer :: WaveMethod = -99 !< Options for including wave information @@ -197,6 +197,8 @@ module MOM_wave_interface real :: VonKar = -1.0 !< The von Karman coefficient as used in the MOM_wave_interface module [nondim] real :: rho_air !< A typical density of air at sea level, as used in wave calculations [R ~> kg m-3] real :: nu_air !< The viscosity of air, as used in wave calculations [Z2 T-1 ~> m2 s-1] + real :: rho_ocn !< A typical surface density of seawater, as used in wave calculations in + !! comparison with the density of air [R ~> kg m-3]. The default is RHO_0. real :: SWH_from_u10sq !< A factor for converting the square of the 10 m wind speed to the !! significant wave height [Z T2 L-2 ~> s2 m-1] real :: Charnock_min !< The minimum value of the Charnock coefficient, which relates the square of @@ -334,7 +336,7 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) if (StatisticalWaves) then CS%WaveMethod = LF17 - call set_LF17_wave_params(param_file, mdl, US, CS) + call set_LF17_wave_params(param_file, mdl, GV, US, CS) if (.not.use_waves) return else CS%WaveMethod = NULL_WaveMethod @@ -500,7 +502,7 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) "Flag to disable updating DHH85 Stokes drift.", default=.false.) case (LF17_STRING) !Li and Fox-Kemper 17 wind-sea Langmuir number CS%WaveMethod = LF17 - call set_LF17_wave_params(param_file, mdl, US, CS) + call set_LF17_wave_params(param_file, mdl, GV, US, CS) case (EFACTOR_STRING) !Li and Fox-Kemper 16 CS%WaveMethod = EFACTOR case default @@ -578,9 +580,10 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) end subroutine MOM_wave_interface_init !> Set the parameters that are used to determine the averaged Stokes drift and Langmuir numbers -subroutine set_LF17_wave_params(param_file, mdl, US, CS) +subroutine set_LF17_wave_params(param_file, mdl, GV, US, CS) type(param_file_type), intent(in) :: param_file !< Input parameter structure character(len=*), intent(in) :: mdl !< A module name to use in the get_param calls + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(wave_parameters_CS), pointer :: CS !< Wave parameter control structure @@ -596,6 +599,10 @@ subroutine set_LF17_wave_params(param_file, mdl, US, CS) call get_param(param_file, mdl, "RHO_AIR", CS%rho_air, & "A typical density of air at sea level, as used in wave calculations", & units="kg m-3", default=1.225, scale=US%kg_m3_to_R) + call get_param(param_file, mdl, "RHO_SFC_WAVES", CS%Rho_ocn, & + "A typical surface density of seawater, as used in wave calculations in "//& + "comparison with the density of air. The default is RHO_0.", & + units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "WAVE_HEIGHT_SCALE_FACTOR", CS%SWH_from_u10sq, & "A factor relating the square of the 10 m wind speed to the significant "//& "wave height, with a default value based on the Pierson-Moskowitz spectrum.", & @@ -713,13 +720,13 @@ end subroutine Update_Surface_Waves !> Constructs the Stokes Drift profile on the model grid based on !! desired coupling options -subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) +subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step) type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure type(ocean_grid_type), intent(inout) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Thickness [H ~> m or kg m-2] + intent(in) :: dz !< Thickness in height units [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: ustar !< Wind friction velocity [Z T-1 ~> m s-1]. real, intent(in) :: dt !< Time-step for computing Stokes-tendency [T ~> s] @@ -728,7 +735,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) ! Local Variables real :: Top, MidPoint, Bottom ! Positions within the layer [Z ~> m] real :: level_thick ! The thickness of each layer [Z ~> m] - real :: DecayScale ! A vertical decay scale in the test profile [Z ~> m] + real :: DecayScale ! A vertical decay scale in the test profile [Z-1 ~> m-1] real :: CMN_FAC ! A nondimensional factor [nondim] real :: WN ! Model wavenumber [Z-1 ~> m-1] real :: UStokes ! A Stokes drift velocity [L T-1 ~> m s-1] @@ -755,8 +762,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) MidPoint = 0.0 do kk = 1,GV%ke Top = Bottom - MidPoint = Bottom - GV%H_to_Z*0.25*(h(II,jj,kk)+h(IIm1,jj,kk)) - Bottom = Bottom - GV%H_to_Z*0.5*(h(II,jj,kk)+h(IIm1,jj,kk)) + MidPoint = Bottom - 0.25*(dz(II,jj,kk)+dz(IIm1,jj,kk)) + Bottom = Bottom - 0.5*(dz(II,jj,kk)+dz(IIm1,jj,kk)) CS%Us_x(II,jj,kk) = CS%TP_STKX0*exp(MidPoint*DecayScale) enddo enddo @@ -768,8 +775,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) MidPoint = 0.0 do kk = 1,GV%ke Top = Bottom - MidPoint = Bottom - GV%H_to_Z*0.25*(h(ii,JJ,kk)+h(ii,JJm1,kk)) - Bottom = Bottom - GV%H_to_Z*0.5*(h(ii,JJ,kk)+h(ii,JJm1,kk)) + MidPoint = Bottom - 0.25*(dz(ii,JJ,kk)+dz(ii,JJm1,kk)) + Bottom = Bottom - 0.5*(dz(ii,JJ,kk)+dz(ii,JJm1,kk)) CS%Us_y(ii,JJ,kk) = CS%TP_STKY0*exp(MidPoint*DecayScale) enddo enddo @@ -796,7 +803,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) do kk = 1,GV%ke Top = Bottom IIm1 = max(II-1,1) - level_thick = 0.5*GV%H_to_Z*(h(II,jj,kk)+h(IIm1,jj,kk)) + level_thick = 0.5*(dz(II,jj,kk)+dz(IIm1,jj,kk)) MidPoint = Top - 0.5*level_thick Bottom = Top - level_thick @@ -854,7 +861,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) do kk = 1,GV%ke Top = Bottom JJm1 = max(JJ-1,1) - level_thick = 0.5*GV%H_to_Z*(h(ii,JJ,kk)+h(ii,JJm1,kk)) + level_thick = 0.5*(dz(ii,JJ,kk)+dz(ii,JJm1,kk)) MidPoint = Top - 0.5*level_thick Bottom = Top - level_thick @@ -908,8 +915,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) do kk = 1,GV%ke Top = Bottom IIm1 = max(II-1,1) - MidPoint = Top - GV%H_to_Z*0.25*(h(II,jj,kk)+h(IIm1,jj,kk)) - Bottom = Top - GV%H_to_Z*0.5*(h(II,jj,kk)+h(IIm1,jj,kk)) + MidPoint = Top - 0.25*(dz(II,jj,kk)+dz(IIm1,jj,kk)) + Bottom = Top - 0.5*(dz(II,jj,kk)+dz(IIm1,jj,kk)) !bgr note that this is using a u-point ii on h-point ustar ! this code has only been previous used for uniform ! grid cases. This needs fixed if DHH85 is used for non @@ -926,8 +933,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) do kk=1, GV%ke Top = Bottom JJm1 = max(JJ-1,1) - MidPoint = Bottom - GV%H_to_Z*0.25*(h(ii,JJ,kk)+h(ii,JJm1,kk)) - Bottom = Bottom - GV%H_to_Z*0.5*(h(ii,JJ,kk)+h(ii,JJm1,kk)) + MidPoint = Bottom - 0.25*(dz(ii,JJ,kk)+dz(ii,JJm1,kk)) + Bottom = Bottom - 0.5*(dz(ii,JJ,kk)+dz(ii,JJm1,kk)) !bgr note that this is using a v-point jj on h-point ustar ! this code has only been previous used for uniform ! grid cases. This needs fixed if DHH85 is used for non @@ -965,9 +972,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar, dt, dynamics_step) ! in the routine it is needed by (e.g. KPP or ePBL). do jj = G%jsc, G%jec do ii = G%isc,G%iec - Top = h(ii,jj,1)*GV%H_to_Z - call get_Langmuir_Number( La, G, GV, US, Top, ustar(ii,jj), ii, jj, & - h(ii,jj,:), CS, Override_MA=.false.) + call get_Langmuir_Number( La, G, GV, US, dz(ii,jj,1), ustar(ii,jj), ii, jj, & + dz(ii,jj,:), CS, Override_MA=.false.) CS%La_turb(ii,jj) = La enddo enddo @@ -1138,7 +1144,7 @@ end subroutine Surface_Bands_by_data_override !! Note this can be called with an unallocated Waves pointer, which is okay if we !! want the wind-speed only dependent Langmuir number. Therefore, we need to be !! careful about what we try to access here. -subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, h, Waves, & +subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, dz, Waves, & U_H, V_H, Override_MA ) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -1148,7 +1154,7 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, h, Waves, & real, intent(in) :: ustar !< Friction velocity [Z T-1 ~> m s-1] integer, intent(in) :: i !< Meridional index of h-point integer, intent(in) :: j !< Zonal index of h-point - real, dimension(SZK_(GV)), intent(in) :: h !< Grid layer thickness [H ~> m or kg m-2] + real, dimension(SZK_(GV)), intent(in) :: dz !< Grid layer thickness [Z ~> m] type(Wave_parameters_CS), pointer :: Waves !< Surface wave control structure. real, dimension(SZK_(GV)), & optional, intent(in) :: U_H !< Zonal velocity at H point [L T-1 ~> m s-1] or [m s-1] @@ -1161,7 +1167,7 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, h, Waves, & !Local Variables - real :: Top, bottom, midpoint ! Positions within each layer [Z ~> m] + real :: Top, Bottom, MidPoint ! Positions within each layer [Z ~> m] real :: Dpt_LASL ! Averaging depth for Stokes drift [Z ~> m] real :: ShearDirection ! Shear angular direction from atan2 [radians] real :: WaveDirection ! Wave angular direction from atan2 [radians] @@ -1185,8 +1191,11 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, h, Waves, & bottom = 0.0 do kk = 1,GV%ke Top = Bottom - MidPoint = Bottom + GV%H_to_Z*0.5*h(kk) - Bottom = Bottom + GV%H_to_Z*h(kk) + MidPoint = Bottom + 0.5*dz(kk) + Bottom = Bottom + dz(kk) + !### Given the sign convention that Dpt_LASL is negative, the next line seems to have a bug. + ! To correct this bug, this line should be changed to: + ! if (MidPoint > abs(Dpt_LASL) .and. (kk > 1) .and. ContinueLoop) then if (MidPoint > Dpt_LASL .and. kk > 1 .and. ContinueLoop) then ShearDirection = atan2(V_H(1)-V_H(kk),U_H(1)-U_H(kk)) ContinueLoop = .false. @@ -1199,8 +1208,8 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, h, Waves, & US_H(kk) = 0.5*(Waves%US_X(I,j,kk)+Waves%US_X(I-1,j,kk)) VS_H(kk) = 0.5*(Waves%US_Y(i,J,kk)+Waves%US_Y(i,J-1,kk)) enddo - call Get_SL_Average_Prof( GV, Dpt_LASL, h, US_H, LA_STKx) - call Get_SL_Average_Prof( GV, Dpt_LASL, h, VS_H, LA_STKy) + call Get_SL_Average_Prof( GV, Dpt_LASL, dz, US_H, LA_STKx) + call Get_SL_Average_Prof( GV, Dpt_LASL, dz, VS_H, LA_STKy) LA_STK = sqrt(LA_STKX*LA_STKX+LA_STKY*LA_STKY) elseif (Waves%WaveMethod==SURFBANDS) then allocate(StkBand_X(Waves%NumBands), StkBand_Y(Waves%NumBands)) @@ -1218,11 +1227,11 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, h, Waves, & US_H(kk) = 0.5*(Waves%US_X(I,j,kk)+Waves%US_X(I-1,j,kk)) VS_H(kk) = 0.5*(Waves%US_Y(i,J,kk)+Waves%US_Y(i,J-1,kk)) enddo - call Get_SL_Average_Prof( GV, Dpt_LASL, h, US_H, LA_STKx) - call Get_SL_Average_Prof( GV, Dpt_LASL, h, VS_H, LA_STKy) + call Get_SL_Average_Prof( GV, Dpt_LASL, dz, US_H, LA_STKx) + call Get_SL_Average_Prof( GV, Dpt_LASL, dz, VS_H, LA_STKy) LA_STK = sqrt(LA_STKX**2 + LA_STKY**2) elseif (Waves%WaveMethod==LF17) then - call get_StokesSL_LiFoxKemper(ustar, hbl*Waves%LA_FracHBL, GV, US, Waves, LA_STK, LA) + call get_StokesSL_LiFoxKemper(ustar, HBL*Waves%LA_FracHBL, GV, US, Waves, LA_STK, LA) elseif (Waves%WaveMethod==Null_WaveMethod) then call MOM_error(FATAL, "Get_Langmuir_number called without defining a WaveMethod. "//& "Suggest to make sure USE_LT is set/overridden to False or choose "//& @@ -1322,7 +1331,7 @@ subroutine get_StokesSL_LiFoxKemper(ustar, hbl, GV, US, CS, UStokes_SL, LA) ! This code should be revised to minimize the number of divisions and cancel out common factors. ! Computing u10 based on u_star and COARE 3.5 relationships - call ust_2_u10_coare3p5(ustar*sqrt(GV%Rho0/CS%rho_air), u10, GV, US, CS) + call ust_2_u10_coare3p5(ustar*sqrt(CS%rho_ocn/CS%rho_air), u10, GV, US, CS) ! surface Stokes drift UStokes = us_to_u10*u10 ! @@ -1406,19 +1415,19 @@ subroutine get_StokesSL_LiFoxKemper(ustar, hbl, GV, US, CS, UStokes_SL, LA) end subroutine Get_StokesSL_LiFoxKemper !> Get SL Averaged Stokes drift from a Stokes drift Profile -subroutine Get_SL_Average_Prof( GV, AvgDepth, H, Profile, Average ) +subroutine Get_SL_Average_Prof( GV, AvgDepth, dz, Profile, Average ) type(verticalGrid_type), & intent(in) :: GV !< Ocean vertical grid structure - real, intent(in) :: AvgDepth !< Depth to average over (negative) [Z ~> m]. + real, intent(in) :: AvgDepth !< Depth to average over (negative) [Z ~> m] real, dimension(SZK_(GV)), & - intent(in) :: H !< Grid thickness [H ~> m or kg m-2] + intent(in) :: dz !< Grid thickness [Z ~> m] real, dimension(SZK_(GV)), & intent(in) :: Profile !< Profile of quantity to be averaged in arbitrary units [A] !! (used here for Stokes drift) real, intent(out) :: Average !< Output quantity averaged over depth AvgDepth [A] !! (used here for Stokes drift) !Local variables - real :: top, midpoint, bottom ! Depths, negative downward [Z ~> m]. + real :: Top, Bottom ! Depths, negative downward [Z ~> m] real :: Sum ! The depth weighted vertical sum of a quantity [A Z ~> A m] integer :: kk @@ -1429,10 +1438,9 @@ subroutine Get_SL_Average_Prof( GV, AvgDepth, H, Profile, Average ) bottom = 0.0 do kk = 1, GV%ke Top = Bottom - MidPoint = Bottom - GV%H_to_Z * 0.5*h(kk) - Bottom = Bottom - GV%H_to_Z * h(kk) + Bottom = Bottom - dz(kk) if (AvgDepth < Bottom) then ! The whole cell is within H_LA - Sum = Sum + Profile(kk) * (GV%H_to_Z * H(kk)) + Sum = Sum + Profile(kk) * dz(kk) elseif (AvgDepth < Top) then ! A partial cell is within H_LA Sum = Sum + Profile(kk) * (Top-AvgDepth) exit @@ -1546,7 +1554,7 @@ end subroutine DHH85_mid !> Explicit solver for Stokes mixing. !! Still in development do not use. -subroutine StokesMixing(G, GV, dt, h, u, v, Waves ) +subroutine StokesMixing(G, GV, dt, h, dz, u, v, Waves ) type(ocean_grid_type), & intent(in) :: G !< Ocean grid type(verticalGrid_type), & @@ -1554,6 +1562,8 @@ subroutine StokesMixing(G, GV, dt, h, u, v, Waves ) real, intent(in) :: dt !< Time step of MOM6 [T ~> s] for explicit solver real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: dz !< Vertical distance between interfaces around a layer [Z ~> m] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: u !< Velocity i-component [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -1561,8 +1571,9 @@ subroutine StokesMixing(G, GV, dt, h, u, v, Waves ) type(Wave_parameters_CS), & pointer :: Waves !< Surface wave related control structure. ! Local variables - real :: dTauUp, dTauDn ! Vertical momentum fluxes [Z L T-2 ~> m2 s-2] - real :: h_Lay ! The layer thickness at a velocity point [Z ~> m]. + real :: dTauUp, dTauDn ! Vertical momentum fluxes [H L T-2 ~> m2 s-2 or Pa] + real :: h_lay ! The layer thickness at a velocity point [H ~> m or kg m-2] + real :: dz_lay ! The distance between interfaces at a velocity point [Z ~> m] integer :: i, j, k ! This is a template to think about down-Stokes mixing. @@ -1571,18 +1582,19 @@ subroutine StokesMixing(G, GV, dt, h, u, v, Waves ) do k = 1, GV%ke do j = G%jsc, G%jec do I = G%iscB, G%iecB - h_lay = GV%H_to_Z*0.5*(h(i,j,k)+h(i+1,j,k)) + h_lay = 0.5*(h(i,j,k)+h(i+1,j,k)) + dz_lay = 0.5*(dz(i,j,k)+dz(i+1,j,k)) dTauUp = 0.0 if (k > 1) & - dTauUp = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i+1,j,k)) * & + dTauUp = (0.5*(waves%Kvs(i,j,k)+waves%Kvs(i+1,j,k))) * & (waves%us_x(i,j,k-1)-waves%us_x(i,j,k)) / & - (0.5*(h_lay + GV%H_to_Z*0.5*(h(i,j,k-1)+h(i+1,j,k-1)) )) + (0.5*(dz_lay + 0.5*(dz(i,j,k-1)+dz(i+1,j,k-1)) )) dTauDn = 0.0 if (k < GV%ke-1) & - dTauDn = 0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i+1,j,k+1)) * & + dTauDn = (0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i+1,j,k+1))) * & (waves%us_x(i,j,k)-waves%us_x(i,j,k+1)) / & - (0.5*(h_lay + GV%H_to_Z*0.5*(h(i,j,k+1)+h(i+1,j,k+1)) )) - u(i,j,k) = u(i,j,k) + dt * (dTauUp-dTauDn) / h_Lay + (0.5*(dz_lay + 0.5*(dz(i,j,k+1)+dz(i+1,j,k+1)) )) + u(i,j,k) = u(i,j,k) + dt * (dTauUp-dTauDn) / h_lay enddo enddo enddo @@ -1590,18 +1602,19 @@ subroutine StokesMixing(G, GV, dt, h, u, v, Waves ) do k = 1, GV%ke do J = G%jscB, G%jecB do i = G%isc, G%iec - h_Lay = GV%H_to_Z*0.5*(h(i,j,k)+h(i,j+1,k)) + h_lay = 0.5*(h(i,j,k)+h(i,j+1,k)) + dz_lay = 0.5*(dz(i,j,k)+dz(i,j+1,k)) dTauUp = 0. if (k > 1) & - dTauUp = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i,j+1,k)) * & + dTauUp = (0.5*(waves%Kvs(i,j,k)+waves%Kvs(i,j+1,k))) * & (waves%us_y(i,j,k-1)-waves%us_y(i,j,k)) / & - (0.5*(h_lay + GV%H_to_Z*0.5*(h(i,j,k-1)+h(i,j+1,k-1)) )) + (0.5*(dz_lay + 0.5*(dz(i,j,k-1)+dz(i,j+1,k-1)) )) dTauDn = 0.0 if (k < GV%ke-1) & - dTauDn =0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i,j+1,k+1)) * & + dTauDn = (0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i,j+1,k+1))) * & (waves%us_y(i,j,k)-waves%us_y(i,j,k+1)) / & - (0.5*(h_lay + GV%H_to_Z*0.5*(h(i,j,k+1)+h(i,j+1,k+1)) )) - v(i,J,k) = v(i,J,k) + dt * (dTauUp-dTauDn) / h_Lay + (0.5*(dz_lay + 0.5*(dz(i,j,k+1)+dz(i,j+1,k+1)) )) + v(i,J,k) = v(i,J,k) + dt * (dTauUp-dTauDn) / h_lay enddo enddo enddo @@ -1658,13 +1671,15 @@ end subroutine CoriolisStokes !! including analytical integration of Stokes shear using multiple-exponential decay !! Stokes drift profile and vertical integration of the resulting pressure !! anomaly to the total pressure gradient force -subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) +subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS ) type(ocean_grid_type), & intent(in) :: G !< Ocean grid type(verticalGrid_type), & intent(in) :: GV !< Ocean vertical grid + type(unit_scale_type), & + intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)),& - intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + intent(in) :: dz !< Layer thicknesses in height units [Z ~> m] real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & intent(in) :: u !< Lagrangian Velocity i-component [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & @@ -1737,12 +1752,13 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) zi_l(1) = 0.0 zi_r(1) = 0.0 do k = 1, G%ke - h_l = h(i,j,k)*GV%H_to_Z - h_r = h(i+1,j,k)*GV%H_to_Z + h_l = dz(i,j,k) + h_r = dz(i+1,j,k) zi_l(k+1) = zi_l(k) - h_l zi_r(k+1) = zi_r(k) - h_r - Idz_l(k) = 1./max(0.1,h_l) - Idz_r(k) = 1./max(0.1,h_r) + !### If the code were properly refactored, the following hard-coded constants would be unnecessary. + Idz_l(k) = 1./max(0.1*US%m_to_Z, h_l) + Idz_r(k) = 1./max(0.1*US%m_to_Z, h_r) enddo do k = 1,G%ke ! Computing (left/right) Eulerian velocities assuming the velocity passed to this routine is the @@ -1830,12 +1846,13 @@ subroutine Stokes_PGF(G, GV, h, u, v, PFu_Stokes, PFv_Stokes, CS ) zi_l(1) = 0.0 zi_r(1) = 0.0 do k = 1, G%ke - h_l = h(i,j,k)*GV%H_to_Z - h_r = h(i,j+1,k)*GV%H_to_Z + h_l = dz(i,j,k) + h_r = dz(i,j+1,k) zi_l(k+1) = zi_l(k) - h_l zi_r(k+1) = zi_r(k) - h_r - Idz_l(k) = 1./max(0.1,h_l) - Idz_r(k) = 1./max(0.1,h_r) + !### If the code were properly refactored, the following hard-coded constants would be unnecessary. + Idz_l(k) = 1. / max(0.1*US%m_to_Z, h_l) + Idz_r(k) = 1. / max(0.1*US%m_to_Z, h_r) enddo do k = 1,G%ke ! Computing (left/right) Eulerian velocities assuming the velocity passed to this routine is the From 8f7cc0e19022ba1bb961e1e874aff5474b60d2fc Mon Sep 17 00:00:00 2001 From: claireyung <61528379+claireyung@users.noreply.github.com> Date: Mon, 7 Aug 2023 19:07:47 +1000 Subject: [PATCH 3/8] Ice shelf melt parameterization fixes (#395) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two fixes to the ice shelf melt parameterization in MOM_ice_shelf.F90: (1) the removal of an extra von Kármán constant which differs from the Holland and Jenkins (1999) formulation. My working is detailed in this document. (2) a modification to the it3 loop in shelf_calc_flux subroutine, which currently does not iterate because wB_flux does not get updated to its new value via the Newton solver method. Instead, the loop either runs 30 times with the same value, or is below the threshold and exits on the first loop. I added a line to update the value of wB_flux. Specific changes include: * Remove KV in n_star_term definition * Add line to ice shelf param it3 to correct iteration * Reinstate von Karman constant and rearrange to remove double division * Remove unneeded salt iteration line that overrides option of using false position method * Change value of ZETA_N to be consistent with Holland & Jenkins 99 * Fix typo, ZETA_N should be 0.13 not 0.013 * Add parameter for ZETA_N * Make RC and VK runtime parameters * Add runtime parameters to fix buoyancy iteration bug and salt overwriting false position method bug * Add runtime parameter for buoyancy iteration Newton's method convergence criteria --- src/ice_shelf/MOM_ice_shelf.F90 | 49 +++++++++++++++++++++++++++------ 1 file changed, 40 insertions(+), 9 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 8e0e58c1b6..d0faeb3aae 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -185,6 +185,14 @@ module MOM_ice_shelf !! salinity [C S-1 ~> degC ppt-1] real :: dTFr_dp !< Partial derivative of freezing temperature with !! pressure [C T2 R-1 L-2 ~> degC Pa-1] + real :: Zeta_N !< The stability constant xi_N = 0.052 from Holland & Jenkins '99 + !! divided by the von Karman constant VK. Was 1/8. + real :: Vk !< Von Karman's constant - dimensionless + real :: Rc !< critical flux Richardson number. + logical :: buoy_flux_itt_bug !< If true, fixes buoyancy iteration bug + logical :: salt_flux_itt_bug !< If true, fixes salt iteration bug + real :: buoy_flux_itt_threshold !< Buoyancy iteration threshold for convergence + !>@{ Diagnostic handles integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & id_tfreeze = -1, id_tfl_shelf = -1, & @@ -261,10 +269,10 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) !! interface, positive for melting and negative for freezing [S ~> ppt]. !! This is computed as part of the ISOMIP diagnostics. real :: time_step !< Length of time over which these fluxes will be applied [T ~> s]. - real, parameter :: VK = 0.40 !< Von Karman's constant - dimensionless - real :: ZETA_N = 0.052 !> The fraction of the boundary layer over which the - !! viscosity is linearly increasing [nondim]. (Was 1/8. Why?) - real, parameter :: RC = 0.20 ! critical flux Richardson number. + real :: VK !< Von Karman's constant - dimensionless + real :: ZETA_N !< This is the stability constant xi_N = 0.052 from Holland & Jenkins '99 + !! divided by the von Karman constant VK. Was 1/8. [nondim] + real :: RC !< critical flux Richardson number. real :: I_ZETA_N !< The inverse of ZETA_N [nondim]. real :: I_LF !< The inverse of the latent heat of fusion [Q-1 ~> kg J-1]. real :: I_VK !< The inverse of the Von Karman constant [nondim]. @@ -346,6 +354,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) endif ! useful parameters is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; ied = G%ied ; jed = G%jed + ZETA_N = CS%Zeta_N + VK = CS%Vk + RC = CS%Rc I_ZETA_N = 1.0 / ZETA_N I_LF = 1.0 / CS%Lat_fusion SC = CS%kv_molec/CS%kd_molec_salt @@ -527,7 +538,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) if (wB_flux < 0.0) then ! The buoyancy flux is stabilizing and will reduce the turbulent ! fluxes, and iteration is required. - n_star_term = (ZETA_N/RC) * (hBL_neut * VK) / (ustar_h)**3 + n_star_term = (ZETA_N * hBL_neut * VK) / (RC * ustar_h**3) do it3 = 1,30 ! n_star <= 1.0 is the ratio of working boundary layer thickness ! to the neutral thickness. @@ -558,13 +569,15 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) 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. 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 + ! Find the root where wB_flux_new = wB_flux. + if (abs(wB_flux_new - wB_flux) < CS%buoy_flux_itt_threshold*(abs(wB_flux_new) + abs(wB_flux))) exit dDwB_dwB_in = dG_dwB * (dB_dS * (dS_ustar * I_Gam_S**2) + & dB_dT * (dT_ustar * I_Gam_T**2)) - 1.0 ! This is Newton's method without any bounds. Should bounds be needed? wB_flux_new = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB_in + ! Update wB_flux + if (CS%buoy_flux_itt_bug) wB_flux = wB_flux_new enddo !it3 endif @@ -637,7 +650,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) Sbdry(i,j) = Sbdry_it endif ! Sb_min_set - Sbdry(i,j) = Sbdry_it + if (.not.CS%salt_flux_itt_bug) Sbdry(i,j) = Sbdry_it + endif ! CS%find_salt_root enddo !it1 @@ -1514,7 +1528,24 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in, call get_param(param_file, mdl, "READ_TIDEAMP", read_TIDEAMP, & "If true, read a file (given by TIDEAMP_FILE) containing "//& "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) - + call get_param(param_file, mdl, "ICE_SHELF_LINEAR_SHELF_FRAC", CS%Zeta_N, & + "Ratio of HJ99 stability constant xi_N (ratio of maximum "//& + "mixing length to planetary boundary layer depth in "//& + "neutrally stable conditions) to the von Karman constant", & + units="nondim", default=0.13) + call get_param(param_file, mdl, "ICE_SHELF_VK_CNST", CS%Vk, & + "Von Karman constant.", & + units="nondim", default=0.40) + call get_param(param_file, mdl, "ICE_SHELF_RC", CS%Rc, & + "Critical flux Richardson number for ice melt ", & + units="nondim", default=0.20) + call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUG", CS%buoy_flux_itt_bug, & + "Bug fix of buoyancy iteration", default=.true.) + call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUG", CS%salt_flux_itt_bug, & + "Bug fix of salt iteration", default=.true.) + call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_THRESHOLD", CS%buoy_flux_itt_threshold, & + "Convergence criterion of Newton's method for ice shelf "//& + "buoyancy iteration.", units="nondim", default=1.0e-4) if (PRESENT(sfc_state_in)) then allocate(sfc_state) From 46c52622cf1d7a39ea9919279bed0b652aad96b7 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 20 Jun 2023 07:36:26 -0400 Subject: [PATCH 4/8] *Use tv%SpV_avg in non-Boussinesq regridding Use SpV_avg in both hybgen_unmix and regridding_main to estimate the nominal ocean bottom depth in thickness units when in fully non-Boussinesq mode. This change eliminates certain dependencies on the Boussinesq reference density via GV%Z_to_H. Also added a call to calc_derived_thermo in ALE_regrid_accelerated that is necessary for the specific volume to be updated, along with tests for the validity of the specific volume with the 1-point halos that are currently used with ALE regridding and remapping. Also, use the total thickness place of the nominal depth in build_grid_rho and build_grid_sigma when in fully non-Boussinesq mode. This is mathematically equivalent but changes answers at roundoff. All answers are bitwise identical in Boussinesq mode, but ALE-based solutions change in fully non-Boussinesq mode with this commit. --- src/ALE/MOM_ALE.F90 | 5 ++- src/ALE/MOM_hybgen_unmix.F90 | 42 ++++++++++++++++++++----- src/ALE/MOM_regridding.F90 | 61 ++++++++++++++++++++++++++---------- 3 files changed, 83 insertions(+), 25 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 2e16933525..b8f60b2830 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -20,7 +20,7 @@ module MOM_ALE use MOM_hybgen_unmix, only : hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix, hybgen_unmix_CS use MOM_hybgen_regrid, only : hybgen_regrid_CS use MOM_file_parser, only : get_param, param_file_type, log_param -use MOM_interface_heights,only : find_eta +use MOM_interface_heights,only : find_eta, calc_derived_thermo use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_regridding, only : initialize_regridding, regridding_main, end_regridding @@ -659,6 +659,9 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d ! generate new grid if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local) + ! Update the layer specific volumes if necessary + if (allocated(tv_local%SpV_avg)) call calc_derived_thermo(tv_local, h, G, GV, US, halo=1) + call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface) dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:) diff --git a/src/ALE/MOM_hybgen_unmix.F90 b/src/ALE/MOM_hybgen_unmix.F90 index 024a9baffa..6ddb828abe 100644 --- a/src/ALE/MOM_hybgen_unmix.F90 +++ b/src/ALE/MOM_hybgen_unmix.F90 @@ -9,6 +9,7 @@ module MOM_hybgen_unmix use MOM_file_parser, only : get_param, param_file_type, log_param use MOM_hybgen_regrid, only : hybgen_column_init use MOM_hybgen_regrid, only : hybgen_regrid_CS, get_hybgen_regrid_params +use MOM_interface_heights, only : calc_derived_thermo use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chkinv use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : ocean_grid_type, thermo_var_ptrs @@ -146,7 +147,8 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h) real :: p_col(GV%ke) ! A column of reference pressures [R L2 T-2 ~> Pa] real :: tracer(GV%ke,max(ntr,1)) ! Columns of each tracer [Conc] real :: h_tot ! Total thickness of the water column [H ~> m or kg m-2] - real :: nominalDepth ! Depth of ocean bottom (positive downward) [H ~> m or kg m-2] + real :: dz_tot ! Vertical distance between the top and bottom of the water column [Z ~> m] + real :: nominalDepth ! Depth of ocean bottom in thickness units (positive downward) [H ~> m or kg m-2] real :: h_thin ! A negligibly small thickness to identify essentially ! vanished layers [H ~> m or kg m-2] real :: dilate ! A factor by which to dilate the target positions from z to z* [nondim] @@ -169,6 +171,15 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h) h_thin = 1e-6*GV%m_to_H debug_conservation = .false. ! Set this to true for debugging + if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then + if (tv%valid_SpV_halo < 0) then + mesg = "invalid values of SpV_avg." + else + mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed." + endif + call MOM_error(FATAL, "hybgen_unmix called in fully non-Boussinesq mode with "//trim(mesg)) + endif + p_col(:) = CS%ref_pressure do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 ; if (G%mask2dT(i,j)>0.) then @@ -203,13 +214,27 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h) endif ! The following block of code is used to trigger z* stretching of the targets heights. - nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H - if (h_tot <= CS%min_dilate*nominalDepth) then - dilate = CS%min_dilate - elseif (h_tot >= CS%max_dilate*nominalDepth) then - dilate = CS%max_dilate + if (allocated(tv%SpV_avg)) then ! This is the fully non-Boussiesq version + dz_tot = 0.0 + do k=1,nk + dz_tot = dz_tot + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h_col(k) + enddo + if (dz_tot <= CS%min_dilate*(G%bathyT(i,j)+G%Z_ref)) then + dilate = CS%min_dilate + elseif (dz_tot >= CS%max_dilate*(G%bathyT(i,j)+G%Z_ref)) then + dilate = CS%max_dilate + else + dilate = dz_tot / (G%bathyT(i,j)+G%Z_ref) + endif else - dilate = h_tot / nominalDepth + nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H + if (h_tot <= CS%min_dilate*nominalDepth) then + dilate = CS%min_dilate + elseif (h_tot >= CS%max_dilate*nominalDepth) then + dilate = CS%max_dilate + else + dilate = h_tot / nominalDepth + endif endif terrain_following = (h_tot < dilate*CS%dpns) .and. (CS%dpns >= CS%dsns) @@ -268,6 +293,9 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h) endif endif ; enddo ; enddo !i & j. + ! Update the layer properties + if (allocated(tv%SpV_avg)) call calc_derived_thermo(tv, h, G, GV, US, halo=1) + end subroutine hybgen_unmix diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 5c4a76c7e5..c238c2aa61 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -814,20 +814,47 @@ subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, & ! Local variables real :: nom_depth_H(SZI_(G),SZJ_(G)) !< The nominal ocean depth at each point in thickness units [H ~> m or kg m-2] + real :: tot_h(SZI_(G),SZJ_(G)) !< The total thickness of the water column [H ~> m or kg m-2] + real :: tot_dz(SZI_(G),SZJ_(G)) !< The total distance between the top and bottom of the water column [Z ~> m] real :: Z_to_H ! A conversion factor used by some routines to convert coordinate ! parameters to depth units [H Z-1 ~> nondim or kg m-3] real :: trickGnuCompiler - integer :: i, j + character(len=128) :: mesg ! A string for error messages + integer :: i, j, k if (present(PCM_cell)) PCM_cell(:,:,:) = .false. Z_to_H = US%Z_to_m * GV%m_to_H ! Often this is equivalent to GV%Z_to_H. - do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 - nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * Z_to_H - ! Consider using the following instead: - ! nom_depth_H(i,j) = max( (G%bathyT(i,j)+G%Z_ref) * Z_to_H , CS%min_nom_depth ) - ! if (G%mask2dT(i,j)==0.) nom_depth_H(i,j) = 0.0 - enddo ; enddo + + if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then + if (tv%valid_SpV_halo < 0) then + mesg = "invalid values of SpV_avg." + else + mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed." + endif + call MOM_error(FATAL, "Regridding_main called in fully non-Boussinesq mode with "//trim(mesg)) + endif + + if (allocated(tv%SpV_avg)) then ! This is the fully non-Boussinesq case + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + tot_h(i,j) = 0.0 ; tot_dz(i,j) = 0.0 + enddo ; enddo + do k=1,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + tot_h(i,j) = tot_h(i,j) + h(i,j,k) + tot_dz(i,j) = tot_dz(i,j) + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h(i,j,k) + enddo ; enddo ; enddo + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + if ((tot_dz(i,j) > 0.0) .and. (G%bathyT(i,j)+G%Z_ref > 0.0)) then + nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * (tot_h(i,j) / tot_dz(i,j)) + else + nom_depth_H(i,j) = 0.0 + endif + enddo ; enddo + else + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + nom_depth_H(i,j) = max((G%bathyT(i,j)+G%Z_ref) * Z_to_H, 0.0) + enddo ; enddo + endif select case ( CS%regridding_scheme ) @@ -1308,12 +1335,12 @@ subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface ) ! In sigma coordinates, the bathymetric depth is only used as an arbitrary offset that ! cancels out when determining coordinate motion, so referencing the column postions to ! the surface is perfectly acceptable, but for preservation of previous answers the - ! referencing is done relative to the bottom when in Boussinesq mode. - ! if (GV%Boussinesq) then + ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode. + if (GV%Boussinesq .or. GV%semi_Boussinesq) then nominalDepth = nom_depth_H(i,j) - ! else - ! nominalDepth = totalThickness - ! endif + else + nominalDepth = totalThickness + endif call build_sigma_column(CS%sigma_CS, nominalDepth, totalThickness, zNew) @@ -1436,12 +1463,12 @@ subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, ! In rho coordinates, the bathymetric depth is only used as an arbitrary offset that ! cancels out when determining coordinate motion, so referencing the column postions to ! the surface is perfectly acceptable, but for preservation of previous answers the - ! referencing is done relative to the bottom when in Boussinesq mode. - ! if (GV%Boussinesq) then + ! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode. + if (GV%Boussinesq .or. GV%semi_Boussinesq) then nominalDepth = nom_depth_H(i,j) - ! else - ! nominalDepth = totalThickness - ! endif + else + nominalDepth = totalThickness + endif ! Determine absolute interface positions zOld(nz+1) = - nominalDepth From ba70663e3e2d688e6a11bedd891fe3b5bee7b3fc Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 2 Aug 2023 08:38:42 -0400 Subject: [PATCH 5/8] *Non-Boussinesq revision of diabatic_driver This commit completes the non-Boussinesq revision of the diabatic_driver code. The revised code uses thickness_to_dz and works with internal variables in vertical distances in the denominator of diffusive flux calculations in diabatic_ALE_legacy, diabatic_ALE and layered_diabatic. The code now uses find_ustar to extract the friction velocities passed to KPP_compute. The (tiny) boundary layer tracer diffusivity is also rescaled to [H2 T-1]. With this set of changes, all implicit references to Boussinesq reference density are eliminated from the calculations in diabatic_driver when in non-Boussinesq mode. A total of 14 thickness rescaling factors were cancelled out, and there are 15 new or renamed variables in the diabatic routines. 4 new checksum calls are also added to diabatic_ALE_legacy to assist in debugging. Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are bitwise identical in that mode, but in non-Boussinesq mode the use of the layer averaged specific volume to convert thicknesses to vertical distances leads to changing answers. --- .../vertical/MOM_diabatic_driver.F90 | 170 +++++++++++------- 1 file changed, 107 insertions(+), 63 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index dae52592e9..b2b8527819 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -12,8 +12,9 @@ module MOM_diabatic_driver use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS use MOM_diabatic_aux, only : make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS -use MOM_diabatic_aux, only : triDiagTS_Eulerian, find_uv_at_h, diagnoseMLDbyDensityDifference -use MOM_diabatic_aux, only : applyBoundaryFluxesInOut, diagnoseMLDbyEnergy, set_pen_shortwave +use MOM_diabatic_aux, only : triDiagTS_Eulerian, find_uv_at_h +use MOM_diabatic_aux, only : applyBoundaryFluxesInOut, set_pen_shortwave +use MOM_diabatic_aux, only : diagnoseMLDbyDensityDifference, diagnoseMLDbyEnergy use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : post_product_sum_u, post_product_sum_v use MOM_diag_mediator, only : diag_ctrl, time_type, diag_update_remap_grids @@ -36,14 +37,14 @@ module MOM_diabatic_driver use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery,MOM_mesg use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_version, param_file_type, read_param -use MOM_forcing_type, only : forcing, MOM_forcing_chksum +use MOM_forcing_type, only : forcing, MOM_forcing_chksum, find_ustar use MOM_forcing_type, only : calculateBuoyancyFlux2d, forcing_SinglePointPrint use MOM_geothermal, only : geothermal_entraining, geothermal_in_place use MOM_geothermal, only : geothermal_init, geothermal_end, geothermal_CS use MOM_grid, only : ocean_grid_type use MOM_int_tide_input, only : set_int_tide_input, int_tide_input_init use MOM_int_tide_input, only : int_tide_input_end, int_tide_input_CS, int_tide_input_type -use MOM_interface_heights, only : find_eta, calc_derived_thermo +use MOM_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz use MOM_internal_tides, only : propagate_int_tide use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used @@ -145,8 +146,8 @@ module MOM_diabatic_driver !! diffusivity of Kd_min_tr (see below) were operating. real :: Kd_BBL_tr !< A bottom boundary layer tracer diffusivity that !! will allow for explicitly specified bottom fluxes - !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. The entrainment at the bottom is at - !! least sqrt(Kd_BBL_tr*dt) over the same distance. + !! [H2 T-1 ~> m2 s-1 or kg2 m-4 s-2]. The entrainment at the + !! bottom is at least sqrt(Kd_BBL_tr*dt) over the same distance. real :: Kd_min_tr !< A minimal diffusivity that should always be !! applied to tracers, especially in massless layers !! near the bottom [H Z T-1 ~> m2 s-1 or kg m-1 s-1] @@ -170,8 +171,6 @@ module MOM_diabatic_driver real :: MLD_En_vals(3) !< Energy values for energy mixed layer diagnostics [R Z3 T-2 ~> J m-2] !>@{ Diagnostic IDs - integer :: id_cg1 = -1 ! diagnostic handle for mode-1 speed - integer, allocatable, dimension(:) :: id_cn ! diagnostic handle for all mode speeds integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic integer :: id_ea_t = -1, id_eb_t = -1, id_ea_s = -1, id_eb_s = -1 integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_int = -1, id_Kd_ePBL = -1 @@ -530,6 +529,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & h_orig, & ! Initial layer thicknesses [H ~> m or kg m-2] + dz, & ! The vertical distance between interfaces around a layer [Z ~> m] dSV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] dSV_dS, & ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. cTKE, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2]. @@ -555,6 +555,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim Sdif_flx ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] real, dimension(SZI_(G),SZJ_(G)) :: & + U_star, & ! The friction velocity [Z T-1 ~> m s-1]. SkinBuoyFlux ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL logical, dimension(SZI_(G)) :: & @@ -562,11 +563,11 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! sufficiently thick that the no-flux boundary conditions have not restricted ! the entrainment - usually sqrt(Kd*dt). - real :: h_neglect ! A thickness that is so small it is usually lost - ! in roundoff and can be neglected [H ~> m or kg m-2] - real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4] + real :: dz_neglect ! A vertical distance that is so small it is usually lost + ! in roundoff and can be neglected [Z ~> m] + real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2] real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2] - real :: I_hval ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1] + real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1] real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is ! coupled to the bottom within a timestep [H ~> m or kg m-2] real :: Kd_add_here ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. @@ -580,7 +581,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect + dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect + Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 showCallTree = callTree_showQuery() @@ -674,19 +676,22 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, & CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) + ! Determine the friction velocity, perhaps using the evovling surface density. + call find_ustar(fluxes, tv, U_star, G, GV, US) + ! The KPP scheme calculates boundary layer diffusivities and non-local transport. if ( associated(fluxes%lamult) ) then call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) + U_star, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) else call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + U_star, CS%KPP_buoy_flux, Waves=Waves) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) endif if (associated(Hml)) then @@ -775,15 +780,18 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call calculate_CVMix_conv(h, tv, G, GV, US, CS%CVMix_conv, Hml, Kd_int, visc%Kv_shear) endif + ! Find the vertical distances across layers. + call thickness_to_dz(h, tv, dz, G, GV, US) + ! This block sets ent_t and ent_s from h and Kd_int. do j=js,je ; do i=is,ie ent_s(i,j,1) = 0.0 ; ent_s(i,j,nz+1) = 0.0 ent_t(i,j,1) = 0.0 ; ent_t(i,j,nz+1) = 0.0 enddo ; enddo - !$OMP parallel do default(shared) private(I_hval) + !$OMP parallel do default(shared) private(I_dzval) do K=2,nz ; do j=js,je ; do i=is,ie - I_hval = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - ent_s(i,j,K) = GV%Z_to_H * dt * I_hval * Kd_int(i,j,K) + I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k))) + ent_s(i,j,K) = dt * I_dzval * Kd_int(i,j,K) ent_t(i,j,K) = ent_s(i,j,K) enddo ; enddo ; enddo if (showCallTree) call callTree_waypoint("done setting ent_s and ent_t from Kd_int (diabatic)") @@ -826,6 +834,11 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim scale=US%kg_m3_to_R*US%degC_to_C) call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS", G%HI, haloshift=0, & scale=US%kg_m3_to_R*US%ppt_to_S) + call hchksum(h, "after applyBoundaryFluxes h", G%HI, haloshift=0, scale=GV%H_to_mks) + call hchksum(tv%T, "after applyBoundaryFluxes tv%T", G%HI, haloshift=0, scale=US%C_to_degC) + call hchksum(tv%S, "after applyBoundaryFluxes tv%S", G%HI, haloshift=0, scale=US%S_to_ppt) + call hchksum(SkinBuoyFlux, "after applyBdryFlux SkinBuoyFlux", G%HI, haloshift=0, & + scale=US%Z_to_m**2*US%s_to_T**3) endif call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) @@ -846,6 +859,9 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call pass_var(visc%sfc_buoy_flx, G%domain, halo=1) endif + ! Find the vertical distances across layers, which may have been modified by the net surface flux + call thickness_to_dz(h, tv, dz, G, GV, US) + ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL. do K=2,nz ; do j=js,je ; do i=is,ie if (CS%ePBL_is_additive) then @@ -856,7 +872,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), CS%ePBL_Prandtl*Kd_ePBL(i,j,K)) endif - Ent_int = Kd_add_here * (GV%Z_to_H * dt) / (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) + Ent_int = Kd_add_here * dt / (0.5*(dz(i,j,k-1) + dz(i,j,k)) + dz_neglect) ent_s(i,j,K) = ent_s(i,j,K) + Ent_int Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here @@ -877,6 +893,9 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, & CS%evap_CFL_limit, CS%minimum_forcing_depth, MLD=visc%MLD) + ! Find the vertical distances across layers, which may have been modified by the net surface flux + call thickness_to_dz(h, tv, dz, G, GV, US) + endif ! endif for CS%use_energetic_PBL ! diagnose the tendencies due to boundary forcing @@ -1002,7 +1021,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim call cpu_clock_begin(id_clock_tracers) if (CS%mix_boundary_tracer_ALE) then - Tr_ea_BBL = sqrt(dt * GV%Z_to_H * CS%Kd_BBL_tr) + Tr_ea_BBL = sqrt(dt * CS%Kd_BBL_tr) !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) do j=js,je @@ -1021,8 +1040,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim ! in the calculation of the fluxes in the first place. Kd_min_tr ! should be much less than the values that have been set in Kd_int, ! perhaps a molecular diffusivity. - add_ent = ((dt * CS%Kd_min_tr) * GV%Z_to_H) * & - ((h(i,j,k-1)+h(i,j,k)+h_neglect) / (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - & + add_ent = ((dt * CS%Kd_min_tr)) * & + ((dz(i,j,k-1)+dz(i,j,k)+dz_neglect) / (dz(i,j,k-1)*dz(i,j,k)+dz_neglect2)) - & 0.5*(ent_s(i,j,K) + ent_s(i,j,K)) if (htot(i) < Tr_ea_BBL) then add_ent = max(0.0, add_ent, (Tr_ea_BBL - htot(i)) - ent_s(i,j,K)) @@ -1034,8 +1053,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim endif if (CS%double_diffuse) then ; if (Kd_extra_S(i,j,k) > 0.0) then - add_ent = ((dt * Kd_extra_S(i,j,k)) * GV%Z_to_H) / & - (0.5 * (h(i,j,k-1) + h(i,j,k)) + h_neglect) + add_ent = (dt * Kd_extra_S(i,j,k)) / & + (0.5 * (dz(i,j,k-1) + dz(i,j,k)) + dz_neglect) ent_s(i,j,K) = ent_s(i,j,K) + add_ent endif ; endif enddo ; enddo @@ -1045,8 +1064,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim !$OMP parallel do default(shared) private(add_ent) do k=nz,2,-1 ; do j=js,je ; do i=is,ie if (Kd_extra_S(i,j,k) > 0.0) then - add_ent = ((dt * Kd_extra_S(i,j,k)) * GV%Z_to_H) / & - (0.5 * (h(i,j,k-1) + h(i,j,k)) + h_neglect) + add_ent = (dt * Kd_extra_S(i,j,k)) / & + (0.5 * (dz(i,j,k-1) + dz(i,j,k)) + dz_neglect) else add_ent = 0.0 endif @@ -1126,6 +1145,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & h_orig, & ! Initial layer thicknesses [H ~> m or kg m-2] + dz, & ! The vertical distance between interfaces around a layer [Z ~> m] dSV_dT, & ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] dSV_dS, & ! The partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. cTKE, & ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2]. @@ -1151,18 +1171,18 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, Sdif_flx ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] real, dimension(SZI_(G),SZJ_(G)) :: & + U_star, & ! The friction velocity [Z T-1 ~> m s-1]. SkinBuoyFlux ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL logical, dimension(SZI_(G)) :: & in_boundary ! True if there are no massive layers below, where massive is defined as ! sufficiently thick that the no-flux boundary conditions have not restricted ! the entrainment - usually sqrt(Kd*dt). - - real :: h_neglect ! A thickness that is so small it is usually lost - ! in roundoff and can be neglected [H ~> m or kg m-2] - real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4] + real :: dz_neglect ! A vertical distance that is so small it is usually lost + ! in roundoff and can be neglected [Z ~> m] + real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2] real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2] - real :: I_hval ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1] + real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1] real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is ! coupled to the bottom within a timestep [H ~> m or kg m-2] real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2]. @@ -1174,7 +1194,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect + dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect + Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 ent_s(:,:,:) = 0.0 ; ent_t(:,:,:) = 0.0 @@ -1276,18 +1297,21 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call calculateBuoyancyFlux2d(G, GV, US, fluxes, CS%optics, h, tv%T, tv%S, tv, & CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) + ! Determine the friction velocity, perhaps using the evovling surface density. + call find_ustar(fluxes, tv, U_star, G, GV, US) + ! The KPP scheme calculates boundary layer diffusivities and non-local transport. if ( associated(fluxes%lamult) ) then call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) + U_star, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) else call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + U_star, CS%KPP_buoy_flux, Waves=Waves) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) endif @@ -1464,17 +1488,20 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, enddo ; enddo ; enddo endif + ! Find the vertical distances across layers, which may have been modified by the net surface flux + call thickness_to_dz(h, tv, dz, G, GV, US) + ! set ent_t=dt*Kd_heat/h_int and est_s=dt*Kd_salt/h_int on interfaces for use in the tridiagonal solver. do j=js,je ; do i=is,ie ent_t(i,j,1) = 0. ; ent_t(i,j,nz+1) = 0. ent_s(i,j,1) = 0. ; ent_s(i,j,nz+1) = 0. enddo ; enddo - !$OMP parallel do default(shared) private(I_hval) + !$OMP parallel do default(shared) private(I_dzval) do K=2,nz ; do j=js,je ; do i=is,ie - I_hval = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - ent_t(i,j,K) = GV%Z_to_H * dt * I_hval * Kd_heat(i,j,k) - ent_s(i,j,K) = GV%Z_to_H * dt * I_hval * Kd_salt(i,j,k) + I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k))) + ent_t(i,j,K) = dt * I_dzval * Kd_heat(i,j,k) + ent_s(i,j,K) = dt * I_dzval * Kd_salt(i,j,k) enddo ; enddo ; enddo if (showCallTree) call callTree_waypoint("done setting ent_t and ent_t from Kd_heat and " //& "Kd_salt (diabatic_ALE)") @@ -1540,7 +1567,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, call cpu_clock_begin(id_clock_tracers) if (CS%mix_boundary_tracer_ALE) then - Tr_ea_BBL = sqrt(dt * GV%Z_to_H * CS%Kd_BBL_tr) + Tr_ea_BBL = sqrt(dt * CS%Kd_BBL_tr) !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) do j=js,je do i=is,ie @@ -1554,8 +1581,8 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, ! bottom, add some mixing of tracers between these layers. This flux is based on the ! harmonic mean of the two thicknesses, following what is done in layered mode. Kd_min_tr ! should be much less than the values in Kd_salt, perhaps a molecular diffusivity. - add_ent = ((dt * CS%Kd_min_tr) * GV%Z_to_H) * & - ((h(i,j,k-1)+h(i,j,k) + h_neglect) / (h(i,j,k-1)*h(i,j,k) + h_neglect2)) - & + add_ent = (dt * CS%Kd_min_tr) * & + ((dz(i,j,k-1)+dz(i,j,k) + dz_neglect) / (dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - & ent_s(i,j,K) if (htot(i) < Tr_ea_BBL) then add_ent = max(0.0, add_ent, (Tr_ea_BBL - htot(i)) - ent_s(i,j,K)) @@ -1648,13 +1675,17 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! one time step [H ~> m or kg m-2] Kd_lay, & ! diapycnal diffusivity of layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] h_orig, & ! initial layer thicknesses [H ~> m or kg m-2] + dz, & ! The vertical distance between interfaces around a layer [Z ~> m] hold, & ! layer thickness before diapycnal entrainment, and later the initial ! layer thicknesses (if a mixed layer is used) [H ~> m or kg m-2] + dz_old, & ! The initial vertical distance between interfaces around a layer + ! or the distance before entrainment [Z ~> m] u_h, & ! Zonal velocities at thickness points after entrainment [L T-1 ~> m s-1] v_h, & ! Meridional velocities at thickness points after entrainment [L T-1 ~> m s-1] temp_diag, & ! Diagnostic array of previous temperatures [C ~> degC] saln_diag ! Diagnostic array of previous salinity [S ~> ppt] real, dimension(SZI_(G),SZJ_(G)) :: & + U_star, & ! The friction velocity [Z T-1 ~> m s-1]. Rcv_ml ! Coordinate density of mixed layer [R ~> kg m-3], used for applying sponges real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: & @@ -1697,7 +1728,9 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2] - real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4] + real :: dz_neglect ! A vertical distance that is so small it is usually lost + ! in roundoff and can be neglected [Z ~> m] + real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2] real :: net_ent ! The net of ea-eb at an interface [H ~> m or kg m-2] real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2] real :: eaval ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2] @@ -1724,7 +1757,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB nkmb = GV%nk_rho_varies - h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect + h_neglect = GV%H_subroundoff + dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 @@ -1885,17 +1919,20 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e enddo ; enddo ; enddo endif + ! Determine the friction velocity, perhaps using the evovling surface density. + call find_ustar(fluxes, tv, U_star, G, GV, US) + if ( associated(fluxes%lamult) ) then call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) + U_star, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) else call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + U_star, CS%KPP_buoy_flux, Waves=Waves) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, tv, U_star, CS%KPP_buoy_flux, Kd_heat, & Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) endif @@ -2300,8 +2337,15 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! mixing of passive tracers from massless boundary layers to interior call cpu_clock_begin(id_clock_tracers) + + ! Find the vertical distances across layers. + if (CS%mix_boundary_tracers .or. CS%double_diffuse) & + call thickness_to_dz(h, tv, dz, G, GV, US) + if (CS%double_diffuse) & + call thickness_to_dz(hold, tv, dz_old, G, GV, US) + if (CS%mix_boundary_tracers) then - Tr_ea_BBL = sqrt(dt * GV%Z_to_H * CS%Kd_BBL_tr) + Tr_ea_BBL = sqrt(dt * CS%Kd_BBL_tr) !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) do j=js,je do i=is,ie @@ -2320,9 +2364,9 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! in the calculation of the fluxes in the first place. Kd_min_tr ! should be much less than the values that have been set in Kd_lay, ! perhaps a molecular diffusivity. - add_ent = ((dt * CS%Kd_min_tr) * GV%Z_to_H) * & - ((h(i,j,k-1)+h(i,j,k)+h_neglect) / & - (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - & + add_ent = (dt * CS%Kd_min_tr) * & + ((dz(i,j,k-1) + dz(i,j,k) + dz_neglect) / & + (dz(i,j,k-1)*dz(i,j,k) + dz_neglect2)) - & 0.5*(ea(i,j,k) + eb(i,j,k-1)) if (htot(i) < Tr_ea_BBL) then add_ent = max(0.0, add_ent, & @@ -2337,9 +2381,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k) endif if (CS%double_diffuse) then ; if (Kd_extra_S(i,j,K) > 0.0) then - add_ent = ((dt * Kd_extra_S(i,j,K)) * GV%Z_to_H) / & - (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + & - h_neglect) + add_ent = (dt * Kd_extra_S(i,j,K)) / & + (0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect) ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent eatr(i,j,k) = eatr(i,j,k) + add_ent endif ; endif @@ -2361,9 +2404,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e !$OMP parallel do default(shared) private(add_ent) do k=nz,2,-1 ; do j=js,je ; do i=is,ie if (Kd_extra_S(i,j,K) > 0.0) then - add_ent = ((dt * Kd_extra_S(i,j,K)) * GV%Z_to_H) / & - (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + & - h_neglect) + add_ent = (dt * Kd_extra_S(i,j,K)) / & + (0.25 * ((dz(i,j,k-1) + dz(i,j,k)) + (dz_old(i,j,k-1) + dz_old(i,j,k))) + dz_neglect) else add_ent = 0.0 endif @@ -3095,7 +3137,9 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di "A bottom boundary layer tracer diffusivity that will "//& "allow for explicitly specified bottom fluxes. The "//& "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//& - "over the same distance.", units="m2 s-1", default=0., scale=GV%m2_s_to_HZ_T) + "over the same distance.", & + units="m2 s-1", default=0., scale=GV%m2_s_to_HZ_T*(US%Z_to_m*GV%m_to_H)) + ! The scaling factor here is usually equivalent to GV%m2_s_to_HZ_T*GV%Z_to_H. endif call get_param(param_file, mdl, "TRACER_TRIDIAG", CS%tracer_tridiag, & From 648012edff219ab8f65d6376d6524b60ac12ed29 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 4 Aug 2023 11:39:28 -0800 Subject: [PATCH 6/8] Adding a knob for strength of brine plume mixing. --- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index f176b0d726..6a5e454d19 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -71,6 +71,7 @@ module MOM_diabatic_aux logical :: do_brine_plume !< If true, insert salt flux below the surface according to !! a parameterization by \cite Nguyen2009. integer :: brine_plume_n !< The exponent in the brine plume parameterization. + real :: plume_strength !< Fraction of the available brine to take to the bottom of the mixed layer. type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output @@ -1425,7 +1426,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! Place forcing into this layer by depth for brine plume parameterization. if (k == 1) then dK(i) = 0.5 * h(i,j,k) * GV%H_to_Z ! Depth of center of layer K - plume_flux = - (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) * GV%RZ_to_H + plume_flux = - (1000.0*US%ppt_to_S * (CS%plume_strength * fluxes%salt_left_behind(i,j))) * GV%RZ_to_H plume_fraction = 1.0 else dK(i) = dK(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * GV%H_to_Z ! Depth of center of layer K @@ -1440,8 +1441,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t plume_fraction = min(fraction_left_brine, h2d(i,k)*IforcingDepthScale) endif fraction_left_brine = fraction_left_brine - plume_fraction - plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * fluxes%salt_left_behind(i,j)) & - * GV%RZ_to_H + plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * (CS%plume_strength * & + fluxes%salt_left_behind(i,j))) * GV%RZ_to_H else plume_flux = 0.0 endif @@ -1767,6 +1768,9 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori call get_param(param_file, mdl, "BRINE_PLUME_EXPONENT", CS%brine_plume_n, & "If using the brine plume parameterization, set the integer exponent.", & default=5, do_not_log=.not.CS%do_brine_plume) + call get_param(param_file, mdl, "BRINE_PLUME_FRACTION", CS%plume_strength, & + "Fraction of the available brine to mix down using the brine plume parameterization.", & + units="nondim", default=1.0, do_not_log=.not.CS%do_brine_plume) if (useALEalgorithm) then CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, & From 07713af190222538c08aa7fafa99c46ab226d09e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 30 Jun 2023 17:50:21 -0400 Subject: [PATCH 7/8] +*Ignore SURFACE_ANSWER_DATE when non-Boussinesq Ignore SURFACE_ANSWER_DATE in non-Boussinesq mode and always allocate tv%SpV_avg in fully non-Boussinesq mode, setting it to the inverse of the layer densities in calc_derived_thermo() if there is no equation of state. Also when in fully non-Boussinesq mode cancelled out some rescaling factors in dz_to_thickness_EOS and dz_to_thickness_tv. Also revised dz_to_thickness_simple in non-Boussinesq and non-layered mode to use RHO_KV_CONVERT instead of RHO_0 to rescale vertical distances to thicknesses. Also set the default value of CALC_RHO_FOR_SEA_LEVEL to true when fully non-Boussinesq. All Boussinesq answers are bitwise identical, but some non-Boussinesq answers do change and become less dependent on the Boussinesq reference density. Because SURFACE_ANSWER_DATE is no longer being used in non-Boussinesq mode, is is no longer being logged in the MOM_parameter_doc files for these experiments. --- src/core/MOM.F90 | 6 ++-- src/core/MOM_interface_heights.F90 | 44 ++++++++++++++++++++++-------- 2 files changed, 36 insertions(+), 14 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9cbb744560..c9b8ea42c0 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2103,7 +2103,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call get_param(param_file, "MOM", "CALC_RHO_FOR_SEA_LEVEL", CS%calc_rho_for_sea_lev, & "If true, the in-situ density is used to calculate the "//& "effective sea level that is returned to the coupler. If false, "//& - "the Boussinesq parameter RHO_0 is used.", default=.false.) + "the Boussinesq parameter RHO_0 is used.", default=non_Bous) call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, & "If true, Temperature and salinity are used as state "//& "variables.", default=.true.) @@ -2892,8 +2892,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & endif endif - ! Allocate any derived equation of state fields. - if (use_temperature .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then + ! Allocate any derived densities or other equation of state derived fields. + if (.not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then allocate(CS%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0) CS%tv%valid_SpV_halo = -1 ! This array does not yet have any valid data. endif diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 1893859fe7..0a579db299 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -279,6 +279,8 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug) ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: p_t ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa] real, dimension(SZI_(G),SZJ_(G)) :: dp ! Pressure change across a layer [R L2 T-2 ~> Pa] + real, dimension(SZK_(GV)) :: SpV_lay ! The specific volume of each layer when no equation of + ! state is used [R-1 ~> m3 kg-1] logical :: do_debug ! If true, write checksums for debugging. integer :: i, j, k, is, ie, js, je, halos, nz @@ -310,6 +312,12 @@ subroutine calc_derived_thermo(tv, h, G, GV, US, halo, debug) call hchksum(tv%T, "derived_thermo T", G%HI, haloshift=halos, scale=US%C_to_degC) call hchksum(tv%S, "derived_thermo S", G%HI, haloshift=halos, scale=US%S_to_ppt) endif + elseif (allocated(tv%Spv_avg)) then + do k=1,nz ; SpV_lay(k) = 1.0 / GV%Rlay(k) ; enddo + do k=1,nz ; do j=js,je ; do i=is,ie + tv%SpV_avg(i,j,k) = SpV_lay(k) + enddo ; enddo ; enddo + tv%valid_SpV_halo = halos endif end subroutine calc_derived_thermo @@ -481,9 +489,7 @@ subroutine dz_to_thickness_tv(dz, tv, h, G, GV, US, halo_size) endif else do k=1,nz ; do j=js,je ; do i=is,ie - h(i,j,k) = (GV%Z_to_H*dz(i,j,k)) * (GV%Rlay(k) / GV%Rho0) - ! Consider revising this to the mathematically equivalent expression: - ! h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k) + h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k) enddo ; enddo ; enddo endif endif @@ -551,10 +557,16 @@ subroutine dz_to_thickness_EOS(dz, Temp, Saln, EoS, h, G, GV, US, halo_size, p_s do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo call calculate_density(Temp(:,j,k), Saln(:,j,k), p_top(:,j), rho, & EoS, EOSdom) - do i=is,ie - ! This could be simplified, but it would change answers at roundoff. - p_bot(i,j) = p_top(i,j) + (GV%g_Earth*GV%H_to_Z) * ((GV%Z_to_H*dz(i,j,k)) * rho(i)) - enddo + ! The following two expressions are mathematically equivalent. + if (GV%semi_Boussinesq) then + do i=is,ie + p_bot(i,j) = p_top(i,j) + (GV%g_Earth*GV%H_to_Z) * ((GV%Z_to_H*dz(i,j,k)) * rho(i)) + enddo + else + do i=is,ie + p_bot(i,j) = p_top(i,j) + rho(i) * (GV%g_Earth * dz(i,j,k)) + enddo + endif enddo do itt=1,max_itt @@ -565,9 +577,15 @@ subroutine dz_to_thickness_EOS(dz, Temp, Saln, EoS, h, G, GV, US, halo_size, p_s EoS, EOSdom) ! Use Newton's method to correct the bottom value. ! The hydrostatic equation is sufficiently linear that no bounds-checking is needed. - do i=is,ie - p_bot(i,j) = p_bot(i,j) + rho(i) * ((GV%g_Earth*GV%H_to_Z)*(GV%Z_to_H*dz(i,j,k)) - dz_geo(i,j)) - enddo + if (GV%semi_Boussinesq) then + do i=is,ie + p_bot(i,j) = p_bot(i,j) + rho(i) * ((GV%g_Earth*GV%H_to_Z)*(GV%Z_to_H*dz(i,j,k)) - dz_geo(i,j)) + enddo + else + do i=is,ie + p_bot(i,j) = p_bot(i,j) + rho(i) * (GV%g_Earth*dz(i,j,k) - dz_geo(i,j)) + enddo + endif enddo ; endif enddo @@ -608,7 +626,7 @@ subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode) layered = .false. ; if (present(layer_mode)) layered = layer_mode is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke - if (GV%Boussinesq .or. (.not.layered)) then + if (GV%Boussinesq) then do k=1,nz ; do j=js,je ; do i=is,ie h(i,j,k) = GV%Z_to_H * dz(i,j,k) enddo ; enddo ; enddo @@ -616,6 +634,10 @@ subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode) do k=1,nz ; do j=js,je ; do i=is,ie h(i,j,k) = (GV%RZ_to_H * GV%Rlay(k)) * dz(i,j,k) enddo ; enddo ; enddo + else + do k=1,nz ; do j=js,je ; do i=is,ie + h(i,j,k) = (US%Z_to_m * GV%m_to_H) * dz(i,j,k) + enddo ; enddo ; enddo endif end subroutine dz_to_thickness_simple From 597bbf119e83405c005b98f280c82224a9f535bb Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 5 Aug 2023 04:59:06 -0400 Subject: [PATCH 8/8] +*Non-Boussinesq revision of full_convection This commit revises full_convection to avoid any dependency on the Boussinesq reference density when in non-Boussinesq mode. Specifically, it changes the units of the Kddt_smooth argument to full_convection and its counterpart in smoothed_dRdT_dRdS to use units of [H Z ~> m2 or kg s-1], which becomes a mass based diffusivity when in non-Boussinesq mode. This change also uses vertical distances for internal calculations in smoothed_dRdT_dRdS, and includes a call to thickness_to_dz in the full_convection routine. This commit also revises the unit conversion of the (absurdly large) mixing length in unstable points in full_convection and of the (tiny) added smoothing distance used in the denominators of some expressions in smoothed_dRdT_dRdS to avoid division by zero with massless layers so that they are both based on the density used in rescaling input parameters (RHO_KV_CONVERT), and do not depend directly on the Boussinesq reference density (RHO_0). These parameters would have been set via get_param calls, but there is no control structure for full_convection parameters. All Boussinesq answers are bitwise identical, but non-Boussinesq answers are altered by the use of the layer specific volumes, rather than the Boussinesq reference density, to convert layer thicknesses into vertical distances. This commit includes a change to the units of an argument (Kddt_smooth) to a publicly visible interface (full_convection). --- .../vertical/MOM_full_convection.F90 | 44 +++++++++++-------- .../vertical/MOM_set_diffusivity.F90 | 2 +- 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/src/parameterizations/vertical/MOM_full_convection.F90 b/src/parameterizations/vertical/MOM_full_convection.F90 index 344511bf29..a5fba3adc6 100644 --- a/src/parameterizations/vertical/MOM_full_convection.F90 +++ b/src/parameterizations/vertical/MOM_full_convection.F90 @@ -3,11 +3,12 @@ module MOM_full_convection ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_grid, only : ocean_grid_type -use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density_derivs, EOS_domain +use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : thickness_to_dz +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type +use MOM_EOS, only : calculate_density_derivs, EOS_domain implicit none ; private @@ -31,15 +32,16 @@ subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: S_adj !< Adjusted salinity [S ~> ppt]. real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL). - real, intent(in) :: Kddt_smooth !< A smoothing vertical - !! diffusivity times a timestep [H2 ~> m2 or kg2 m-4]. + real, intent(in) :: Kddt_smooth !< A smoothing vertical diffusivity + !! times a timestep [H Z ~> m2 or kg m-1]. integer, intent(in) :: halo !< Halo width over which to compute ! Local variables real, dimension(SZI_(G),SZK_(GV)+1) :: & dRho_dT, & ! The derivative of density with temperature [R C-1 ~> kg m-3 degC-1] dRho_dS ! The derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]. - real :: h_neglect, h0 ! A thickness that is so small it is usually lost + real :: dz(SZI_(G),SZK_(GV)) ! Height change across layers [Z ~> m] + real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. ! logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. real, dimension(SZI_(G),SZK0_(G)) :: & @@ -90,15 +92,17 @@ subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, if (.not.associated(tv%eqn_of_state)) return h_neglect = GV%H_subroundoff - mix_len = (1.0e20 * nz) * (G%max_depth * GV%Z_to_H) - h0 = 1.0e-16*sqrt(Kddt_smooth) + h_neglect + mix_len = (1.0e20 * nz) * (G%max_depth * US%Z_to_m * GV%m_to_H) do j=js,je mix(:,:) = 0.0 ; d_b(:,:) = 1.0 ! These would be Te_b(:,:) = tv%T(:,j,:), etc., but the values are not used Te_b(:,:) = 0.0 ; Se_b(:,:) = 0.0 - call smoothed_dRdT_dRdS(h, tv, Kddt_smooth, dRho_dT, dRho_dS, G, GV, US, j, p_surf, halo) + ! Find the vertical distances across layers. + call thickness_to_dz(h, tv, dz, j, G, GV, halo_size=halo) + + call smoothed_dRdT_dRdS(h, dz, tv, Kddt_smooth, dRho_dT, dRho_dS, G, GV, US, j, p_surf, halo) do i=is,ie do_i(i) = (G%mask2dT(i,j) > 0.0) @@ -306,14 +310,16 @@ end function is_unstable !> Returns the partial derivatives of locally referenced potential density with !! temperature and salinity after the properties have been smoothed with a small !! constant diffusivity. -subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, halo) +subroutine smoothed_dRdT_dRdS(h, dz, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, halo) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZK_(GV)), & + intent(in) :: dz !< Height change across layers [Z ~> m] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables - real, intent(in) :: Kddt !< A diffusivity times a time increment [H2 ~> m2 or kg2 m-4]. + real, intent(in) :: Kddt !< A diffusivity times a time increment [H Z ~> m2 or kg m-1]. real, dimension(SZI_(G),SZK_(GV)+1), & intent(out) :: dR_dT !< Derivative of locally referenced !! potential density with temperature [R C-1 ~> kg m-3 degC-1] @@ -336,8 +342,9 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h real :: pres(SZI_(G)) ! Interface pressures [R L2 T-2 ~> Pa]. real :: T_EOS(SZI_(G)) ! Filtered and vertically averaged temperatures [C ~> degC] real :: S_EOS(SZI_(G)) ! Filtered and vertically averaged salinities [S ~> ppt] - real :: kap_dt_x2 ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4]. - real :: h_neglect, h0 ! Negligible thicknesses to allow for zero thicknesses, + real :: kap_dt_x2 ! The product of 2*kappa*dt [H Z ~> m2 or kg m-1]. + real :: dz_neglect, h0 ! A negligible vertical distances [Z ~> m] + real :: h_neglect ! A negligible thickness to allow for zero thicknesses ! [H ~> m or kg m-2]. real :: h_tr ! The thickness at tracer points, plus h_neglect [H ~> m or kg m-2]. integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state @@ -347,6 +354,7 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h nz = GV%ke h_neglect = GV%H_subroundoff + dz_neglect = GV%dz_subroundoff kap_dt_x2 = 2.0*Kddt if (Kddt <= 0.0) then @@ -354,9 +362,9 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h T_f(i,k) = tv%T(i,j,k) ; S_f(i,k) = tv%S(i,j,k) enddo ; enddo else - h0 = 1.0e-16*sqrt(Kddt) + h_neglect + h0 = 1.0e-16*sqrt(GV%H_to_m*US%m_to_Z*Kddt) + dz_neglect do i=is,ie - mix(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0) + mix(i,2) = kap_dt_x2 / ((dz(i,1)+dz(i,2)) + h0) h_tr = h(i,j,1) + h_neglect b1(i) = 1.0 / (h_tr + mix(i,2)) @@ -365,7 +373,7 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h S_f(i,1) = (b1(i)*h_tr)*tv%S(i,j,1) enddo do k=2,nz-1 ; do i=is,ie - mix(i,K+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0) + mix(i,K+1) = kap_dt_x2 / ((dz(i,k)+dz(i,k+1)) + h0) c1(i,k) = mix(i,K) * b1(i) h_tr = h(i,j,k) + h_neglect diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index b3b49e0772..32553de3d1 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -344,7 +344,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i call cpu_clock_begin(id_clock_kappaShear) if (CS%Vertex_shear) then call full_convection(G, GV, US, h, tv, T_f, S_f, fluxes%p_surf, & - GV%Z_to_H*kappa_dt_fill, halo=1) + kappa_dt_fill, halo=1) call calc_kappa_shear_vertex(u, v, h, T_f, S_f, tv, fluxes%p_surf, visc%Kd_shear, & visc%TKE_turb, visc%Kv_shear_Bu, dt, G, GV, US, CS%kappaShear_CSp)