From 4226800cbfe8d4fbc3103030aeaf16d3e2e53638 Mon Sep 17 00:00:00 2001 From: Theresa Morrison Date: Fri, 2 Aug 2024 13:16:44 -0400 Subject: [PATCH] Add option to use non-surface density in MLD_003 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A new parameter is added to diagnoseMLDbyDensityDifference that allows a user to specify a reference for the surface density that is not the top model layer. This feature should make the MLD_003 diagnostic more consistent with the de Boyer Montégut MLD climatology. In addition, new diagnostics have been added to save the actual depth of the density used and the "surface" density used in the MLD calculation. These options have only been added for the MLD_003 option. --- .../vertical/MOM_diabatic_aux.F90 | 89 +++++++++++++++---- .../vertical/MOM_diabatic_driver.F90 | 19 +++- 2 files changed, 91 insertions(+), 17 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index aa31024b24..59b37860ab 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -679,7 +679,7 @@ end subroutine set_pen_shortwave !> Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface. !> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping. subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, diagPtr, & - id_N2subML, id_MLDsq, dz_subML) + ref_h_mld, id_ref_z, id_ref_rho, id_N2subML, id_MLDsq, dz_subML) type(ocean_grid_type), intent(in) :: G !< Grid type type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -690,6 +690,9 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, !! available thermodynamic fields. real, intent(in) :: densityDiff !< Density difference to determine MLD [R ~> kg m-3] type(diag_ctrl), pointer :: diagPtr !< Diagnostics structure + real, intent(in) :: ref_h_mld !< Depth of the calculated "surface" densisty [Z ~> m] + integer, intent(in) :: id_ref_z !< Handle (ID) of reference depth diagnostic + integer, intent(in) :: id_ref_rho !< Handle (ID) of reference density diagnostic integer, optional, intent(in) :: id_N2subML !< Optional handle (ID) of subML stratification integer, optional, intent(in) :: id_MLDsq !< Optional handle (ID) of squared MLD real, optional, intent(in) :: dz_subML !< The distance over which to calculate N2subML @@ -698,6 +701,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, ! Local variables real, dimension(SZI_(G)) :: deltaRhoAtKm1, deltaRhoAtK ! Density differences [R ~> kg m-3]. real, dimension(SZI_(G)) :: pRef_MLD, pRef_N2 ! Reference pressures [R L2 T-2 ~> Pa]. + real, dimension(SZI_(G)) :: hRef_MLD ! Reference depth [Z ~> m]. real, dimension(SZI_(G)) :: H_subML, dH_N2 ! Summed thicknesses used in N2 calculation [H ~> m or kg m-2] real, dimension(SZI_(G)) :: dZ_N2 ! Summed vertical distance used in N2 calculation [Z ~> m] real, dimension(SZI_(G)) :: T_subML, T_deeper ! Temperatures used in the N2 calculation [C ~> degC]. @@ -706,6 +710,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, real, dimension(SZI_(G),SZK_(GV)) :: dZ_2d ! Layer thicknesses in depth units [Z ~> m] real, dimension(SZI_(G)) :: dZ, dZm1 ! Layer thicknesses associated with interfaces [Z ~> m] real, dimension(SZI_(G)) :: rhoSurf ! Density used in finding the mixed layer depth [R ~> kg m-3]. + real, dimension(SZI_(G), SZJ_(G)) :: z_ref_diag ! The actual depth of the reference density [Z ~> m]. real, dimension(SZI_(G), SZJ_(G)) :: MLD ! Diagnosed mixed layer depth [Z ~> m]. real, dimension(SZI_(G), SZJ_(G)) :: subMLN2 ! Diagnosed stratification below ML [T-2 ~> s-2]. real, dimension(SZI_(G), SZJ_(G)) :: MLD2 ! Diagnosed MLD^2 [Z2 ~> m2]. @@ -716,6 +721,10 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, real :: dZ_sub_ML ! Depth below ML over which to diagnose stratification [Z ~> m] real :: aFac ! A nondimensional factor [nondim] real :: ddRho ! A density difference [R ~> kg m-3] + real :: dddpth ! A depth difference [Z ~> m] + real :: rhoSurf_k, rhoSurf_km1 ! Desisty in the layers below and above the target reference depth [R ~> kg m-3]. + real, dimension(SZI_(G), SZJ_(G)) :: rhoSurf_2d ! The density that is considered the "surface" when calculating + ! the MLD. It can be saved as a diagnostic [R ~> kg m-3]. integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, k, nz, id_N2, id_SQ @@ -737,24 +746,72 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - pRef_MLD(:) = 0.0 + hRef_MLD(:) = ref_h_mld + pRef_MLD(:) = GV%H_to_RZ*GV%g_Earth*ref_h_mld + z_ref_diag(:,:) = 0. + EOSdom(:) = EOS_domain(G%HI) do j=js,je ! Find the vertical distances across layers. call thickness_to_dz(h, tv, dZ_2d, j, G, GV) - do i=is,ie ; dZ(i) = 0.5 * dZ_2d(i,1) ; enddo ! Depth of center of surface layer - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom) - do i=is,ie - deltaRhoAtK(i) = 0. - MLD(i,j) = 0. - if (id_N2>0) then - subMLN2(i,j) = 0.0 - H_subML(i) = h(i,j,1) ; dH_N2(i) = 0.0 ; dZ_N2(i) = 0.0 - T_subML(i) = 0.0 ; S_subML(i) = 0.0 ; T_deeper(i) = 0.0 ; S_deeper(i) = 0.0 - N2_region_set(i) = (G%mask2dT(i,j)<0.5) ! Only need to work on ocean points. - endif - enddo + if (pRef_MLD(is) /= 0.0) then + rhoSurf(:) = 0.0 + do i=is,ie + dZ(i) = 0.5 * dZ_2d(i,1) ! Depth of center of surface layer + if (dZ(i) >= hRef_MLD(i)) then + call calculate_density(tv%T(i,j,1), tv%S(i,j,1), pRef_MLD(i), rhoSurf_k, tv%eqn_of_state) + rhoSurf(i) = rhoSurf_k + endif + enddo + do k=2,nz + do i=is,ie + dZm1(i) = dZ(i) ! Depth of center of layer K-1 + dZ(i) = dZ(i) + 0.5 * ( dZ_2d(i,k) + dZ_2d(i,k-1) ) ! Depth of center of layer K + dddpth = dZ(i) - dZm1(i) + if ((rhoSurf(i) == 0.) .and. & + (dZm1(i) < hRef_MLD(i)) .and. (dZ(i) >= hRef_MLD(i))) then + aFac = ( hRef_MLD(i) - dZm1(i) ) / dddpth + z_ref_diag(i,j) = (dZ(i) * aFac + dZm1(i) * (1. - aFac)) + call calculate_density(tv%T(i,j,k) , tv%S(i,j,k) , pRef_MLD(i), rhoSurf_k, tv%eqn_of_state) + call calculate_density(tv%T(i,j,k-1), tv%S(i,j,k-1), pRef_MLD(i), rhoSurf_km1, tv%eqn_of_state) + rhoSurf(i) = (rhoSurf_k * aFac + rhoSurf_km1 * (1. - aFac)) + H_subML(i) = h(i,j,k) + elseif ((rhoSurf(i) == 0.) .and. (k >= nz)) then + call calculate_density(tv%T(i,j,1), tv%S(i,j,1), pRef_MLD(i), rhoSurf_k, tv%eqn_of_state) + rhoSurf(i) = rhoSurf_k + endif + enddo + enddo + do i=is,ie + dZ(i) = 0.5 * dZ_2d(i,1) ! reset dZ to surface depth + rhoSurf_2d(i,j) = rhoSurf(i) + deltaRhoAtK(i) = 0. + MLD(i,j) = 0. + if (id_N2>0) then + subMLN2(i,j) = 0.0 + dH_N2(i) = 0.0 ; dZ_N2(i) = 0.0 + T_subML(i) = 0.0 ; S_subML(i) = 0.0 ; T_deeper(i) = 0.0 ; S_deeper(i) = 0.0 + N2_region_set(i) = (G%mask2dT(i,j)<0.5) ! Only need to work on ocean points. + endif + enddo + elseif (pRef_MLD(is) == 0.0) then + rhoSurf(:) = 0.0 + do i=is,ie ; dZ(i) = 0.5 * dZ_2d(i,1) ; enddo ! Depth of center of surface layer + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom) + do i=is,ie + rhoSurf_2d(i,j) = rhoSurf(i) + deltaRhoAtK(i) = 0. + MLD(i,j) = 0. + if (id_N2>0) then + subMLN2(i,j) = 0.0 + H_subML(i) = h(i,j,1) ; dH_N2(i) = 0.0 ; dZ_N2(i) = 0.0 + T_subML(i) = 0.0 ; S_subML(i) = 0.0 ; T_deeper(i) = 0.0 ; S_deeper(i) = 0.0 + N2_region_set(i) = (G%mask2dT(i,j)<0.5) ! Only need to work on ocean points. + endif + enddo + endif + do k=2,nz do i=is,ie dZm1(i) = dZ(i) ! Depth of center of layer K-1 @@ -823,6 +880,9 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, if (id_N2 > 0) call post_data(id_N2, subMLN2, diagPtr) if (id_SQ > 0) call post_data(id_SQ, MLD2, diagPtr) + if ((id_ref_z > 0) .and. (pRef_MLD(is)/=0.)) call post_data(id_ref_z, z_ref_diag , diagPtr) + if (id_ref_rho > 0) call post_data(id_ref_rho, rhoSurf_2d , diagPtr) + end subroutine diagnoseMLDbyDensityDifference !> Diagnose a mixed layer depth (MLD) determined by the depth a given energy value would mix. @@ -1964,5 +2024,4 @@ end subroutine diabatic_aux_end !! salinities due to the application of the surface forcing. It may also calculate the implied !! turbulent kinetic energy requirements for this forcing to be mixed over the model's finite !! vertical resolution in the surface layers. - end module MOM_diabatic_aux diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 18bbc45e13..9386b599a7 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -175,6 +175,8 @@ module MOM_diabatic_driver real :: dz_subML_N2 !< The distance over which to calculate a diagnostic of the !! average stratification at the base of the mixed layer [Z ~> m]. real :: MLD_En_vals(3) !< Energy values for energy mixed layer diagnostics [R Z3 T-2 ~> J m-2] + real :: ref_h_mld = 0.0 !< The depth of the "surface" density used in a difference mixed based + !! MLD calculation [Z ~> m]. !>@{ Diagnostic IDs integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic @@ -183,6 +185,7 @@ module MOM_diabatic_driver integer :: id_Tdif = -1, id_Sdif = -1, id_Tadv = -1, id_Sadv = -1 ! These are handles to diagnostics related to the mixed layer properties. integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1 + integer :: id_MLD_003_zr = -1, id_MLD_003_rr = -1 integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3 = -1, id_subMLN2 = -1 ! These are handles to diagnostics that are only available in non-ALE layered mode. @@ -479,13 +482,16 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, & call enable_averages(dt, Time_end, CS%diag) if (CS%id_MLD_003 > 0 .or. CS%id_subMLN2 > 0 .or. CS%id_mlotstsq > 0) then call diagnoseMLDbyDensityDifference(CS%id_MLD_003, h, tv, 0.03*US%kg_m3_to_R, G, GV, US, CS%diag, & + CS%ref_h_mld, CS%id_MLD_003_zr, CS%id_MLD_003_rr, & id_N2subML=CS%id_subMLN2, id_MLDsq=CS%id_mlotstsq, dz_subML=CS%dz_subML_N2) endif if (CS%id_MLD_0125 > 0) then - call diagnoseMLDbyDensityDifference(CS%id_MLD_0125, h, tv, 0.125*US%kg_m3_to_R, G, GV, US, CS%diag) + call diagnoseMLDbyDensityDifference(CS%id_MLD_0125, h, tv, 0.125*US%kg_m3_to_R, G, GV, US, CS%diag, & + ref_H_MLD=0.0, id_ref_z=-1, id_ref_rho=-1) endif if (CS%id_MLD_user > 0) then - call diagnoseMLDbyDensityDifference(CS%id_MLD_user, h, tv, CS%MLDdensityDifference, G, GV, US, CS%diag) + call diagnoseMLDbyDensityDifference(CS%id_MLD_user, h, tv, CS%MLDdensityDifference, G, GV, US, CS%diag, & + ref_H_MLD=0.0, id_ref_z=-1, id_ref_rho=-1) endif if ((CS%id_MLD_EN1 > 0) .or. (CS%id_MLD_EN2 > 0) .or. (CS%id_MLD_EN3 > 0)) then call diagnoseMLDbyEnergy((/CS%id_MLD_EN1, CS%id_MLD_EN2, CS%id_MLD_EN3/),& @@ -3273,6 +3279,15 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Squared buoyancy frequency below mixed layer', units='s-2', conversion=US%s_to_T**2) CS%id_MLD_user = register_diag_field('ocean_model', 'MLD_user', diag%axesT1, Time, & 'Mixed layer depth (used defined)', units='m', conversion=US%Z_to_m) + if (CS%id_MLD_003 > 0) then + call get_param(param_file, mdl, "HREF_FOR_MLD", CS%ref_h_mld, & + "Reference depth used to calculate the potential density used to find the mixed layer depth "//& + "based on a delta rho = 0.03 kg/m3.", units='m', default=0.0, scale=US%m_to_Z) + CS%id_MLD_003_zr = register_diag_field('ocean_model', 'MLD_003_refZ', diag%axesT1, Time, & + 'Depth of reference density for MLD (delta rho = 0.03)', units='m', conversion=US%Z_to_m) + CS%id_MLD_003_rr = register_diag_field('ocean_model', 'MLD_003_refRho', diag%axesT1, Time, & + 'Reference density for MLD (delta rho = 0.03)', units='kg/m3', conversion=US%R_to_kg_m3) + endif endif call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", CS%MLDdensityDifference, & "The density difference used to determine a diagnostic mixed "//&