From f14a681079bd563b69b7db5133f23bd0192d63de Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 24 Jan 2024 09:26:50 -0500 Subject: [PATCH 1/3] (*)Use thickness_to_dz with BRYAN_LEWIS_DIFFUSIVITY Use thickness_to_dz to convert the layer thicknesses into the geometric depths used to set the Bryan-Lewis background diffusivities instead of multiplication by GV%H_to_m, thereby avoiding division by the Boussinesq reference density in non-Boussinesq mode. No answers are changed in any Boussinsesq configurations, but answers will change in non-Boussinesq configurations that have BRYAN_LEWIS_DIFFUSIVITY = True. --- src/parameterizations/vertical/MOM_bkgnd_mixing.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 693b9395bd..ee04f3f195 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -13,6 +13,7 @@ module MOM_bkgnd_mixing use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_forcing_type, only : forcing 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_verticalGrid, only : verticalGrid_type use MOM_variables, only : thermo_var_ptrs, vertvisc_type @@ -338,6 +339,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G, real, dimension(SZK_(GV)+1) :: Kv_col !< Viscosities at the interfaces [m2 s-1] real, dimension(SZI_(G)) :: Kd_sfc !< Surface value of the diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1] real, dimension(SZI_(G)) :: depth !< Distance from surface of an interface [H ~> m or kg m-2] + real, dimension(SZI_(G),SZK_(GV)) :: dz !< Height change across layers [Z ~> m] real :: depth_c !< depth of the center of a layer [H ~> m or kg m-2] real :: I_Hmix !< inverse of fixed mixed layer thickness [H-1 ~> m-1 or m2 kg-1] real :: I_2Omega !< 1/(2 Omega) [T ~> s] @@ -364,10 +366,12 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G, ! Set up the background diffusivity. if (CS%Bryan_Lewis_diffusivity) then + call thickness_to_dz(h, tv, dz, j, G, GV) + do i=is,ie depth_int(1) = 0.0 do k=2,nz+1 - depth_int(k) = depth_int(k-1) + GV%H_to_m*h(i,j,k-1) + depth_int(k) = depth_int(k-1) + US%Z_to_m*dz(i,k-1) enddo call CVMix_init_bkgnd(max_nlev=nz, & From 3f7465a57923a75fe23e6b7216ea27264ab7349d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 28 Jan 2024 11:17:45 -0500 Subject: [PATCH 2/3] +tracer_Z_init can rescale tracers Added the new optional scale argument to tracer_Z_init to allow the units of the tracers that it initializes from a z-space file to be rescaled, similarly to what is already done in tracer_z_init_array. This also required the addition of the new missing_scale argument to the private routine read_Z_edges to specify how the missing values are rescaled. Comments were also added or modified describing 28 real variables or their units in this module. By default all answers are bitwise identical, but when this new option is exercised, answers could change at roundoff for some non-Boussinesq cases. --- src/tracer/MOM_tracer_Z_init.F90 | 89 +++++++++++++++++++------------- 1 file changed, 54 insertions(+), 35 deletions(-) diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index c404c560f3..2cf0ba1efe 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -27,50 +27,58 @@ module MOM_tracer_Z_init !> This function initializes a tracer by reading a Z-space file, returning !! .true. if this appears to have been successful, and false otherwise. -function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_val) +function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_val, scale) logical :: tracer_Z_init !< A return code indicating if the initialization has been successful type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure 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(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: tr !< The tracer to initialize + intent(out) :: tr !< The tracer to initialize [CU ~> conc] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] - character(len=*), intent(in) :: filename !< The name of the file to read from - character(len=*), intent(in) :: tr_name !< The name of the tracer in the file - real, optional, intent(in) :: missing_val !< The missing value for the tracer - real, optional, intent(in) :: land_val !< A value to use to fill in land points - -! This include declares and sets the variable "version". -# include "version_variable.h" + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] or other + !! arbitrary units such as [Z ~> m] + character(len=*), intent(in) :: filename !< The name of the file to read from + character(len=*), intent(in) :: tr_name !< The name of the tracer in the file + real, optional, intent(in) :: missing_val !< The missing value for the tracer [CU ~> conc] + real, optional, intent(in) :: land_val !< A value to use to fill in land points [CU ~> conc] + real, optional, intent(in) :: scale !< A factor by which to scale the output tracers from the + !! their units in the file [CU conc-1 ~> 1] + ! Local variables real, allocatable, dimension(:,:,:) :: & - tr_in ! The z-space array of tracer concentrations that is read in. + tr_in ! The z-space array of tracer concentrations that is read in [CU ~> conc] real, allocatable, dimension(:) :: & z_edges, & ! The depths of the cell edges or cell centers (depending on ! the value of has_edges) in the input z* data [Z ~> m]. - tr_1d, & ! A copy of the input tracer concentrations in a column. + tr_1d, & ! A copy of the input tracer concentrations in a column [CU ~> conc] wt, & ! The fractional weight for each layer in the range between - ! k_top and k_bot, nondim. - z1, & ! z1 and z2 are the depths of the top and bottom limits of the part - z2 ! of a z-cell that contributes to a layer, relative to the cell - ! center and normalized by the cell thickness, nondim. + ! k_top and k_bot [nondim] + z1, z2 ! z1 and z2 are the depths of the top and bottom limits of the part + ! of a z-cell that contributes to a layer, relative to the cell + ! center and normalized by the cell thickness [nondim]. ! Note that -1/2 <= z1 <= z2 <= 1/2. real :: e(SZK_(GV)+1) ! The z-star interface heights [Z ~> m]. - real :: landval ! The tracer value to use in land points. + real :: landval ! The tracer value to use in land points [CU ~> conc] real :: sl_tr ! The normalized slope of the tracer - ! within the cell, in tracer units. + ! within the cell, in tracer units [CU ~> conc] real :: htot(SZI_(G)) ! The vertical sum of h [H ~> m or kg m-2]. real :: dilate ! The amount by which the thicknesses are dilated to - ! create a z-star coordinate, nondim or in m3 kg-1. - real :: missing ! The missing value for the tracer. - + ! create a z-star coordinate [Z H-1 ~> nondim or m3 kg-1] + ! or other units reflecting those of h + real :: missing ! The missing value for the tracer [CU ~> conc] + real :: scale_fac ! A factor by which to scale the output tracers from the units in the + ! input file [CU conc-1 ~> 1] + ! This include declares and sets the variable "version". +# include "version_variable.h" logical :: has_edges, use_missing, zero_surface character(len=80) :: loc_msg integer :: k_top, k_bot, k_bot_prev, k_start integer :: i, j, k, kz, is, ie, js, je, nz, nz_in + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + scale_fac = 1.0 ; if (present(scale)) then ; scale_fac = scale ; endif + landval = 0.0 ; if (present(land_val)) landval = land_val zero_surface = .false. ! Make this false for errors to be fatal. @@ -83,7 +91,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va ! Find out the number of input levels and read the depth of the edges, ! also modifying their sign convention to be monotonically decreasing. call read_Z_edges(filename, tr_name, z_edges, nz_in, has_edges, use_missing, & - missing, scale=US%m_to_Z) + missing, scale=US%m_to_Z, missing_scale=scale_fac) if (nz_in < 1) then tracer_Z_init = .false. return @@ -91,7 +99,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va allocate(tr_in(G%isd:G%ied,G%jsd:G%jed,nz_in), source=0.0) allocate(tr_1d(nz_in), source=0.0) - call MOM_read_data(filename, tr_name, tr_in(:,:,:), G%Domain) + call MOM_read_data(filename, tr_name, tr_in(:,:,:), G%Domain, scale=scale_fac) ! Fill missing values from above? Use a "close" test to avoid problems ! from type-conversion rounoff. @@ -297,9 +305,10 @@ subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, n real :: e_1d(nlay+1) ! A 1-d column of interface heights, in the same units as e [Z ~> m] or [m] real :: sl_tr ! The tracer concentration slope times the layer thickness, in tracer units [B] real :: wt(nk_data) ! The fractional weight for each layer in the range between z1 and z2 [nondim] - real :: z1(nk_data) ! z1 and z2 are the fractional depths of the top and bottom - real :: z2(nk_data) ! limits of the part of a z-cell that contributes to a layer, relative - ! to the cell center and normalized by the cell thickness [nondim]. + real :: z1(nk_data) ! The fractional depth of the top limit of the part of a z-cell that contributes to + ! a layer, relative to the cell center and normalized by the cell thickness [nondim]. + real :: z2(nk_data) ! The fractional depth of the bottom limit of the part of a z-cell that contributes to + ! a layer, relative to the cell center and normalized by the cell thickness [nondim]. ! Note that -1/2 <= z1 <= z2 <= 1/2. real :: scale_fac ! A factor by which to scale the output tracers from the input tracers [B A-1 ~> 1] integer :: k_top, k_bot, k_bot_prev, kstart @@ -380,18 +389,21 @@ end subroutine tracer_z_init_array !> This subroutine reads the vertical coordinate data for a field from a NetCDF file. !! It also might read the missing value attribute for that same field. subroutine read_Z_edges(filename, tr_name, z_edges, nz_out, has_edges, & - use_missing, missing, scale) + use_missing, missing, scale, missing_scale) character(len=*), intent(in) :: filename !< The name of the file to read from. character(len=*), intent(in) :: tr_name !< The name of the tracer in the file. real, dimension(:), allocatable, & - intent(out) :: z_edges !< The depths of the vertical edges of the tracer array + intent(out) :: z_edges !< The depths of the vertical edges of the tracer array [Z ~> m] integer, intent(out) :: nz_out !< The number of vertical layers in the tracer array logical, intent(out) :: has_edges !< If true the values in z_edges are the edges of the !! tracer cells, otherwise they are the cell centers logical, intent(inout) :: use_missing !< If false on input, see whether the tracer has a !! missing value, and if so return true - real, intent(inout) :: missing !< The missing value, if one has been found - real, intent(in) :: scale !< A scaling factor for z_edges into new units. + real, intent(inout) :: missing !< The missing value, if one has been found [CU ~> conc] + real, intent(in) :: scale !< A scaling factor for z_edges into new units [Z m-1 ~> 1] + real, intent(in) :: missing_scale !< A scaling factor to use to convert the + !! tracers and their missing value from the units in + !! the file into their internal units [CU conc-1 ~> 1] ! This subroutine reads the vertical coordinate data for a field from a ! NetCDF file. It also might read the missing value attribute for that same field. @@ -419,6 +431,7 @@ subroutine read_Z_edges(filename, tr_name, z_edges, nz_out, has_edges, & if (.not.use_missing) then ! Try to find the missing value from the dataset. call read_attribute(filename, "missing_value", missing, varname=tr_name, found=use_missing, ncid_in=ncid) + if (use_missing) missing = missing * missing_scale endif ! Find out if the Z-axis has an edges attribute call read_attribute(filename, "edges", edge_name, varname=dim_names(3), found=has_edges, ncid_in=ncid) @@ -475,8 +488,12 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z real, dimension(:), intent(out) :: z2 !< Depths of the bottom limit of the part of !! a layer that contributes to a depth level, relative to the cell center and normalized !! by the cell thickness [nondim]. Note that -1/2 <= z1 < z2 <= 1/2. + ! Local variables - real :: Ih, e_c, tot_wt, I_totwt + real :: Ih ! The inverse of the vertical distance across a layer, in the inverse of the units of e [Z-1 ~> m-1] + real :: e_c ! The height of the layer center, in the units of e [Z ~> m] + real :: tot_wt ! The sum of the thicknesses contributing to a layer [Z ~> m] + real :: I_totwt ! The Adcroft reciprocal of tot_wt [Z-1 ~> m-1] integer :: k wt(:) = 0.0 ; z1(:) = 0.0 ; z2(:) = 0.0 ; k_bot = k_max @@ -494,6 +511,7 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z z1(k) = (e_c - MIN(e(K), Z_top)) * Ih z2(k) = (e_c - Z_bot) * Ih else + ! Note that in theis branch, wt temporarily has units of [Z ~> m] wt(k) = MIN(e(K),Z_top) - e(K+1) ; tot_wt = wt(k) ! These are always > 0. if (e(K) /= e(K+1)) then z1(k) = (0.5*(e(K)+e(K+1)) - MIN(e(K), Z_top)) / (e(K)-e(K+1)) @@ -515,6 +533,7 @@ subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z enddo I_totwt = 0.0 ; if (tot_wt > 0.0) I_totwt = 1.0 / tot_wt + ! This loop changes the units of wt from [Z ~> m] to [nondim]. do k=k_top,k_bot ; wt(k) = I_totwt*wt(k) ; enddo endif @@ -523,13 +542,13 @@ end subroutine find_overlap !> This subroutine determines a limited slope for val to be advected with !! a piecewise limited scheme. function find_limited_slope(val, e, k) result(slope) - real, dimension(:), intent(in) :: val !< An column the values that are being interpolated. + real, dimension(:), intent(in) :: val !< A column of the values that are being interpolated, in arbitrary units [A] real, dimension(:), intent(in) :: e !< A column's interface heights [Z ~> m] or other units. integer, intent(in) :: k !< The layer whose slope is being determined. - real :: slope !< The normalized slope in the intracell distribution of val. + real :: slope !< The normalized slope in the intracell distribution of val [A] ! Local variables - real :: amn, cmn - real :: d1, d2 + real :: amn, cmn ! Limited differences and curvatures in the values [A] + real :: d1, d2 ! Layer thicknesses, in the units of e [Z ~> m] if ((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) then slope = 0.0 ! ; curvature = 0.0 From f92e4ac3117bf832d7a6f9b241e3cc084d45b238 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Sun, 25 Feb 2024 15:14:33 -0500 Subject: [PATCH 3/3] Update MOM_mixed_layer_restrat.F90 (#568) * Update MOM_mixed_layer_restrat.F90 Adding and calculating diagnostic front length scale (lf_bodner) to mixed layer restratification code Co-authored-by: Kate Hedstrom Co-authored-by: Marshall Ward * Bodner diag: Fix ustar and dims Fix two issues which emerged from the Bodner length scale diagnostic after the cuberoot ANSWER_DATE was introduced. The split of the loop due to the cuberoot caused ustar to be assigned in a separate loop, leaving it uninitialized for the length scale diagnostics loop. It also had its dimensions removed. However, the whole point of the new cuberoot() function was to preserve its dimensions, so there was no need for the intermediate dimensional manipulation, and all dimensionality should be shifted to the registration. This patch removes the intermediate dimensions and also explicitly defines it as Z2/H rather than L. * Bodner diagnostic: Restore scaling to L The scaling of the Bodner lengthscale diagnostic is correctly restored to L in this patch. The Z2 H-1 constant is factored into the calculation of the lengthscale; I believe this represents a conversion from a [Z2 T-2]-like momentum flux to a [LH T-2]-like momentum flux. A generic `uflux_rescale` term was introduced to be used in both ld_bodner and wpup, to avoid a double division by SpV. The if-block for computation of ld_bodner was also moved outside of the loop, in order to facilitate vectorization. --------- Co-authored-by: Liz Drenkard Co-authored-by: Kate Hedstrom --- .../lateral/MOM_mixed_layer_restrat.F90 | 58 ++++++++++++++++--- 1 file changed, 49 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index c10a55309b..0efe322a1d 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -120,6 +120,7 @@ module MOM_mixed_layer_restrat integer :: id_wpup = -1 integer :: id_ustar = -1 integer :: id_bflux = -1 + integer :: id_lfbod = -1 !>@} end type mixedlayer_restrat_CS @@ -807,6 +808,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d wpup ! Turbulent vertical momentum [L H T-2 ~> m2 s-2 or kg m-1 s-2] real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] + real :: lf_bodner_diag(SZI_(G),SZJ_(G)) ! Front width as in Bodner et al., 2023 (B22), eq 24 [L ~> m] real :: U_star_2d(SZI_(G),SZJ_(G)) ! 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] @@ -829,6 +831,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d real :: u_star3 ! Cube of surface friction velocity [Z3 T-3 ~> m3 s-3] real :: r_wpup ! reciprocal of vertical momentum flux [T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1] real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1] + real :: f_h ! Coriolis parameter at h-points [T-1 ~> s-1] + real :: f2_h ! Coriolis parameter at h-points squared [T-2 ~> s-2] + real :: absurdly_small_freq2 ! Frequency squared used to avoid division by 0 [T-2 ~> s-2] real :: grid_dsd ! combination of grid scales [L2 ~> m2] real :: h_sml ! "Little h", the active mixing depth with diurnal cycle removed [H ~> m or kg m-2] real :: h_big ! "Big H", the mixed layer depth based on a time filtered "little h" [H ~> m or kg m-2] @@ -862,6 +867,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d covTS(:) = 0.0 ! Might be in tv% in the future. Not implemented for the time being. varS(:) = 0.0 ! Ditto. + ! This value is roughly (pi / (the age of the universe) )^2. + absurdly_small_freq2 = 1e-34*US%T_to_s**2 + if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// & "An equation of state must be used with this module.") if (.not.CS%MLE_use_PBL_MLD) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// & @@ -934,8 +942,8 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d ! This expression differs by a factor of 1. / (Rho_0 * SpV_avg) compared with the other ! expressions below, and it is invariant to the value of Rho_0 in non-Boussinesq mode. wpup(i,j) = max((cuberoot( CS%mstar * U_star_2d(i,j)**3 + & - CS%nstar * max(0., -bflux(i,j)) * BLD(i,j) ))**2, CS%min_wstar2) * & - ( US%Z_to_L * GV%RZ_to_H / tv%SpV_avg(i,j,1)) + CS%nstar * max(0., -bflux(i,j)) * BLD(i,j) ))**2, CS%min_wstar2) & + * (US%Z_to_L * GV%RZ_to_H / tv%SpV_avg(i,j,1)) ! The final line above converts from [Z2 T-2 ~> m2 s-2] to [L H T-2 ~> m2 s-2 or Pa]. ! Some rescaling factors and the division by specific volume compensating for other ! factors that are in find_ustar_mech, and others effectively converting the wind @@ -952,14 +960,13 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d w_star3 = max(0., -bflux(i,j)) * BLD(i,j) ! In [Z3 T-3 ~> m3 s-3] u_star3 = U_star_2d(i,j)**3 ! In [Z3 T-3 ~> m3 s-3] wpup(i,j) = max(m2_s2_to_Z2_T2 * (Z3_T3_to_m3_s3 * ( CS%mstar * u_star3 + CS%nstar * w_star3 ) )**two_thirds, & - CS%min_wstar2) * & - ( US%Z_to_L * US%Z_to_m * GV%m_to_H ) ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2] + CS%min_wstar2) * US%Z_to_L * GV%Z_to_H ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2] enddo ; enddo else do j=js-1,je+1 ; do i=is-1,ie+1 w_star3 = max(0., -bflux(i,j)) * BLD(i,j) ! In [Z3 T-3 ~> m3 s-3] - wpup(i,j) = max( (cuberoot(CS%mstar * U_star_2d(i,j)**3 + CS%nstar * w_star3))**2, CS%min_wstar2 ) * & - ( US%Z_to_L * US%Z_to_m * GV%m_to_H ) ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2] + wpup(i,j) = max( (cuberoot(CS%mstar * U_star_2d(i,j)**3 + CS%nstar * w_star3))**2, CS%min_wstar2 ) & + * US%Z_to_L * GV%Z_to_H ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2] enddo ; enddo endif @@ -970,6 +977,35 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d CS%wpup_filtered(i,j) = wpup(i,j) enddo ; enddo + if (CS%id_lfbod > 0) then + do j=js-1,je+1 ; do i=is-1,ie+1 + ! Calculate front length used in B22 formula (eq 24). + w_star3 = max(0., -bflux(i,j)) * BLD(i,j) + u_star3 = U_star_2d(i,j)**3 + + ! Include an absurdly_small_freq2 to prevent division by zero. + f_h = 0.25 * ((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) & + + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) + f2_h = max(f_h**2, absurdly_small_freq2) + + lf_bodner_diag(i,j) = & + 0.25 * cuberoot(CS%mstar * u_star3 + CS%nstar * w_star3)**2 & + / (f2_h * max(little_h(i,j), GV%Angstrom_H)) + enddo ; enddo + + ! Rescale from [Z2 H-1 to L] + if (allocated(tv%SpV_avg) .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then + do j=js-1,je+1 ; do i=is-1,ie+1 + lf_bodner_diag(i,j) = lf_bodner_diag(i,j) & + * (US%Z_to_L * GV%RZ_to_H / tv%SpV_avg(i,j,1)) + enddo ; enddo + else + do j=js-1,je+1 ; do i=is-1,ie+1 + lf_bodner_diag(i,j) = lf_bodner_diag(i,j) * US%Z_to_L * GV%Z_to_H + enddo ; enddo + endif + endif + if (CS%debug) then call hchksum(little_h,'mle_Bodner: little_h', G%HI, haloshift=1, scale=GV%H_to_mks) call hchksum(big_H,'mle_Bodner: big_H', G%HI, haloshift=1, scale=GV%H_to_mks) @@ -1155,6 +1191,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d if (CS%id_vhml > 0) call post_data(CS%id_vhml, vhml, CS%diag) if (CS%id_uDml > 0) call post_data(CS%id_uDml, uDml_diag, CS%diag) if (CS%id_vDml > 0) call post_data(CS%id_vDml, vDml_diag, CS%diag) + if (CS%id_lfbod > 0) call post_data(CS%id_lfbod, lf_bodner_diag, CS%diag) if (CS%id_uml > 0) then do J=js,je ; do i=is-1,ie @@ -1776,14 +1813,17 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, 'm s-1', conversion=US%L_T_to_m_s) if (CS%use_Bodner) then CS%id_wpup = register_diag_field('ocean_model', 'MLE_wpup', diag%axesT1, Time, & - 'Vertical turbulent momentum flux in Bodner mixed layer restratificiation parameterization', & + 'Vertical turbulent momentum flux in Bodner mixed layer restratification parameterization', & 'm2 s-2', conversion=US%L_to_m*GV%H_to_m*US%s_to_T**2) CS%id_ustar = register_diag_field('ocean_model', 'MLE_ustar', diag%axesT1, Time, & - 'Surface turbulent friction velicity, u*, in Bodner mixed layer restratificiation parameterization', & + 'Surface turbulent friction velocity, u*, in Bodner mixed layer restratification parameterization', & 'm s-1', conversion=(US%Z_to_m*US%s_to_T)) CS%id_bflux = register_diag_field('ocean_model', 'MLE_bflux', diag%axesT1, Time, & - 'Surface buoyancy flux, B0, in Bodner mixed layer restratificiation parameterization', & + 'Surface buoyancy flux, B0, in Bodner mixed layer restratification parameterization', & 'm2 s-3', conversion=(US%Z_to_m**2*US%s_to_T**3)) + CS%id_lfbod = register_diag_field('ocean_model', 'lf_bodner', diag%axesT1, Time, & + 'Front length in Bodner mixed layer restratificiation parameterization', & + 'm', conversion=US%L_to_m) endif ! If MLD_filtered is being used, we need to update halo regions after a restart