Skip to content

Commit 2a45abf

Browse files
Theresa MorrisonTheresa Morrison
Theresa Morrison
authored and
Theresa Morrison
committed
Add new method for finding target depth
Add new method for finding target depth that should be more efficient. The conversion from h to p has been corrected and should be correct, although I am not sure about the factor of 0.5?
1 parent 9737d2c commit 2a45abf

File tree

2 files changed

+89
-37
lines changed

2 files changed

+89
-37
lines changed

‎src/parameterizations/vertical/MOM_diabatic_aux.F90

+76-25
Original file line numberDiff line numberDiff line change
@@ -679,7 +679,7 @@ end subroutine set_pen_shortwave
679679
!> Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface.
680680
!> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping.
681681
subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, diagPtr, &
682-
id_N2subML, id_MLDsq, dz_subML, ref_p_mld, id_ref_z)
682+
id_N2subML, id_MLDsq, dz_subML, ref_h_mld, id_ref_z, id_ref_rho)
683683
type(ocean_grid_type), intent(in) :: G !< Grid type
684684
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
685685
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
@@ -694,9 +694,10 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
694694
integer, optional, intent(in) :: id_MLDsq !< Optional handle (ID) of squared MLD
695695
real, optional, intent(in) :: dz_subML !< The distance over which to calculate N2subML
696696
!! or 50 m if missing [Z ~> m]
697-
real, optional, intent(in) :: ref_p_mld !< Optional reference pressure used to calculate the
698-
!! densisty, defults to 0.0 is not present [R L2 T-2 ~> Pa].
697+
real, optional, intent(in) :: ref_h_mld !< Optional reference depth used to calculate the
698+
!! densisty, defults to 0.0 is not present [Z ~> m].
699699
integer, optional, intent(in) :: id_ref_z !< Handle (ID) of reference depth diagnostic
700+
integer, optional, intent(in) :: id_ref_rho !< Handle (ID) of reference depth diagnostic
700701

701702
! Local variables
702703
real, dimension(SZI_(G)) :: deltaRhoAtKm1, deltaRhoAtK ! Density differences [R ~> kg m-3].
@@ -710,7 +711,6 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
710711
real, dimension(SZI_(G),SZK_(GV)) :: dZ_2d ! Layer thicknesses in depth units [Z ~> m]
711712
real, dimension(SZI_(G)) :: dZ, dZm1 ! Layer thicknesses associated with interfaces [Z ~> m]
712713
real, dimension(SZI_(G)) :: rhoSurf ! Density used in finding the mixed layer depth [R ~> kg m-3].
713-
real :: rhoSurf_ref ! Density used in finding the mixed layer depth [R ~> kg m-3].
714714
real, dimension(SZI_(G), SZJ_(G)) :: z_ref_diag ! the actual depth of the k interface [Z ~> m].
715715
real, dimension(SZI_(G), SZJ_(G)) :: MLD ! Diagnosed mixed layer depth [Z ~> m].
716716
real, dimension(SZI_(G), SZJ_(G)) :: subMLN2 ! Diagnosed stratification below ML [T-2 ~> s-2].
@@ -719,11 +719,14 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
719719
! have been stored already.
720720
real :: gE_Rho0 ! The gravitational acceleration, sometimes divided by the Boussinesq
721721
! reference density [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2].
722+
real :: H_to_RL2_T2 ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
722723
real :: dZ_sub_ML ! Depth below ML over which to diagnose stratification [Z ~> m]
723724
real :: aFac ! A nondimensional factor [nondim]
724725
real :: ddRho ! A density difference [R ~> kg m-3]
725-
real, dimension(SZI_(G)) :: ks
726-
real :: errZ, errZm1
726+
real :: dddpth !, thresh
727+
real :: rhoSurf_k, rhoSurf_km1
728+
real, dimension(SZI_(G), SZJ_(G)) :: rhoSurf_2d
729+
!real, dimension(SZI_(G)) :: ks
727730
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
728731
integer :: i, j, is, ie, js, je, k, nz, id_N2, id_SQ
729732

@@ -746,41 +749,51 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
746749
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
747750

748751
pRef_MLD(:) = 0.0; hRef_MLD(:) = 0.0
749-
if (present(ref_p_mld)) pRef_MLD(is:ie) = ref_p_mld
750-
if (present(ref_p_mld)) hRef_MLD(is:ie) = ref_p_mld ! horrible pressure to depth calc
752+
if (present(ref_h_mld)) hRef_MLD(:) = ref_h_mld
753+
if (present(ref_h_mld)) then
754+
H_to_RL2_T2 = GV%H_to_RZ * GV%g_Earth
755+
pRef_MLD(:) = (0.5*H_to_RL2_T2)*ref_h_mld
756+
endif
751757
z_ref_diag(:,:) = 0.
758+
!thresh = 0.0011
752759

753760
EOSdom(:) = EOS_domain(G%HI)
754761
do j=js,je
755762
! Find the vertical distances across layers.
756763
call thickness_to_dz(h, tv, dZ_2d, j, G, GV)
757764

758765
if (pRef_MLD(is) .ne. 0.0) then
759-
ks(is:ie) = 0
766+
!ks(is:ie) = 0
767+
rhoSurf(:) = 0.0
760768
do i=is,ie
761769
dZ(i) = 0.5 * dZ_2d(i,1) ! Depth of center of surface layer
762-
if (dZ(i) >= hRef_MLD(i)) ks(i) = k
770+
if (dZ(i) >= hRef_MLD(i)) then
771+
call calculate_density(tv%T(i,j,1), tv%S(i,j,1), pRef_MLD(i), rhoSurf_k, tv%eqn_of_state)
772+
rhoSurf(i) = rhoSurf_k
773+
endif
763774
enddo
764775
do k=2,nz
765776
do i=is,ie
766-
if (ks(i)==0 .and. k<nz) then ! still havent found 10 m depth
767-
if (dZ(i) >= hRef_MLD(i)) then ! check whether the previous layer was closer
768-
errZ = dZ(i) - hRef_MLD(i)
769-
errZm1 = hRef_MLD(i) - dZm1(i)
770-
if (errZ <= errZm1) then; ks(i) = k; z_ref_diag(i,j)=dZ(i); ! set the k reference
771-
elseif (errZm1 < errZ) then; ks(i) = k-1; z_ref_diag(i,j)=dZm1(i); endif
772-
elseif (dZ(i) < hRef_MLD(i)) then ! go to the next layer
773-
dZm1(i) = dZ(i) ! Depth of center of layer K-1
774-
dZ(i) = dZ(i) + 0.5 * ( dZ_2d(i,k) + dZ_2d(i,k-1) ) ! Depth of center of layer K
775-
endif
776-
elseif (ks(i)==0 .and. k==nz) then
777-
ks(i)=1
777+
dZm1(i) = dZ(i) ! Depth of center of layer K-1
778+
dZ(i) = dZ(i) + 0.5 * ( dZ_2d(i,k) + dZ_2d(i,k-1) ) ! Depth of center of layer K
779+
dddpth = dZ(i) - dZm1(i)
780+
!if ((rhoSurf(i) == 0.) .and. (dddpth > thresh) .and. &
781+
if ((rhoSurf(i) == 0.) .and. &
782+
(dZm1(i) < hRef_MLD(i)) .and. (dZ(i) >= hRef_MLD(i))) then
783+
aFac = ( hRef_MLD(i) - dZm1(i) ) / dddpth
784+
z_ref_diag(i,j) = (dZ(i) * aFac + dZm1(i) * (1. - aFac))
785+
call calculate_density(tv%T(i,j,k) , tv%S(i,j,k) , pRef_MLD, rhoSurf_k, tv%eqn_of_state)
786+
call calculate_density(tv%T(i,j,k-1), tv%S(i,j,k-1), pRef_MLD, rhoSurf_km1, tv%eqn_of_state)
787+
rhoSurf(i) = (rhoSurf_k * aFac + rhoSurf_km1 * (1. - aFac))
788+
!elseif (dddpth<=thresh) then
789+
elseif ((rhoSurf(i) == 0.) .and. (k>=nz)) then
790+
call calculate_density(tv%T(i,j,1), tv%S(i,j,1), pRef_MLD(i), rhoSurf_k, tv%eqn_of_state)
791+
rhoSurf(i) = rhoSurf_k
778792
endif
779793
enddo
780794
enddo
781795
do i=is,ie
782-
call calculate_density(tv%T(i,j,ks(i)), tv%S(i,j,ks(i)), pRef_MLD(i), rhoSurf_ref, tv%eqn_of_state)
783-
rhoSurf(i) = rhoSurf_ref
796+
rhoSurf_2d(i,j)=rhoSurf(i)
784797
deltaRhoAtK(i) = 0.
785798
MLD(i,j) = 0.
786799
if (id_N2>0) then
@@ -790,6 +803,42 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
790803
N2_region_set(i) = (G%mask2dT(i,j)<0.5) ! Only need to work on ocean points.
791804
endif
792805
enddo
806+
do i=is,ie ; dZ(i) = 0.5 * dZ_2d(i,1) ; enddo ! Reset the depth of center to the surface layer
807+
!!!!!!!!!!!!!!
808+
!ks(is:ie) = 0
809+
!do i=is,ie
810+
! dZ(i) = 0.5 * dZ_2d(i,1) ! Depth of center of surface layer
811+
! if (dZ(i) >= hRef_MLD(i)) ks(i) = k
812+
!enddo
813+
!do k=2,nz
814+
! do i=is,ie
815+
! if (ks(i)==0 .and. k<nz) then ! still havent found 10 m depth
816+
! if (dZ(i) >= hRef_MLD(i)) then ! check whether the previous layer was closer
817+
! errZ = dZ(i) - hRef_MLD(i)
818+
! errZm1 = hRef_MLD(i) - dZm1(i)
819+
! if (errZ <= errZm1) then; ks(i) = k; z_ref_diag(i,j)=dZ(i); ! set the k reference
820+
! elseif (errZm1 < errZ) then; ks(i) = k-1; z_ref_diag(i,j)=dZm1(i); endif
821+
! elseif (dZ(i) < hRef_MLD(i)) then ! go to the next layer
822+
! dZm1(i) = dZ(i) ! Depth of center of layer K-1
823+
! dZ(i) = dZ(i) + 0.5 * ( dZ_2d(i,k) + dZ_2d(i,k-1) ) ! Depth of center of layer K
824+
! endif
825+
! elseif (ks(i)==0 .and. k==nz) then
826+
! ks(i)=1
827+
! endif
828+
! enddo
829+
!enddo
830+
!do i=is,ie
831+
! call calculate_density(tv%T(i,j,ks(i)), tv%S(i,j,ks(i)), pRef_MLD(i), rhoSurf_ref, tv%eqn_of_state)
832+
! rhoSurf(i) = rhoSurf_ref
833+
! deltaRhoAtK(i) = 0.
834+
! MLD(i,j) = 0.
835+
! if (id_N2>0) then
836+
! subMLN2(i,j) = 0.0
837+
! H_subML(i) = h(i,j,1) ; dH_N2(i) = 0.0 ; dZ_N2(i) = 0.0 ! may need to cahnge the 1 on this line?
838+
! T_subML(i) = 0.0 ; S_subML(i) = 0.0 ; T_deeper(i) = 0.0 ; S_deeper(i) = 0.0
839+
! N2_region_set(i) = (G%mask2dT(i,j)<0.5) ! Only need to work on ocean points.
840+
! endif
841+
!enddo
793842
!!!!!!!!!!!!
794843
!do i=is,ie
795844
! ks(i) = 0
@@ -825,9 +874,11 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
825874
! endif
826875
!enddo
827876
elseif (pRef_MLD(is) == 0.0) then
877+
rhoSurf(:) = 0.0
828878
do i=is,ie ; dZ(i) = 0.5 * dZ_2d(i,1) ; enddo ! Depth of center of surface layer
829879
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom)
830880
do i=is,ie
881+
rhoSurf_2d(i,j)=rhoSurf(i)
831882
deltaRhoAtK(i) = 0.
832883
MLD(i,j) = 0.
833884
if (id_N2>0) then
@@ -909,6 +960,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US,
909960
if (id_SQ > 0) call post_data(id_SQ, MLD2, diagPtr)
910961

911962
if ((id_ref_z > 0) .and. (pRef_MLD(is).ne.0.)) call post_data(id_ref_z, z_ref_diag , diagPtr)
963+
if (id_ref_rho > 0) call post_data(id_ref_rho, rhoSurf_2d , diagPtr)
912964

913965
end subroutine diagnoseMLDbyDensityDifference
914966

@@ -2051,5 +2103,4 @@ end subroutine diabatic_aux_end
20512103
!! salinities due to the application of the surface forcing. It may also calculate the implied
20522104
!! turbulent kinetic energy requirements for this forcing to be mixed over the model's finite
20532105
!! vertical resolution in the surface layers.
2054-
20552106
end module MOM_diabatic_aux

‎src/parameterizations/vertical/MOM_diabatic_driver.F90

+13-12
Original file line numberDiff line numberDiff line change
@@ -175,8 +175,8 @@ module MOM_diabatic_driver
175175
real :: dz_subML_N2 !< The distance over which to calculate a diagnostic of the
176176
!! average stratification at the base of the mixed layer [Z ~> m].
177177
real :: MLD_En_vals(3) !< Energy values for energy mixed layer diagnostics [R Z3 T-2 ~> J m-2]
178-
real :: ref_p_mld = 0.0 !< A referece pressure for density used in a difference mixed based
179-
!! MLD calculation [R L2 T-2 ~> Pa].
178+
real :: ref_h_mld = 0.0 !< A reference depthfor density used in a difference mixed based
179+
!! MLD calculation [Z ~> m].
180180

181181
!>@{ Diagnostic IDs
182182
integer :: id_ea = -1, id_eb = -1 ! used by layer diabatic
@@ -185,7 +185,7 @@ module MOM_diabatic_driver
185185
integer :: id_Tdif = -1, id_Sdif = -1, id_Tadv = -1, id_Sadv = -1
186186
! These are handles to diagnostics related to the mixed layer properties.
187187
integer :: id_MLD_003 = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_mlotstsq = -1
188-
integer :: id_MLD_003_zr = -1
188+
integer :: id_MLD_003_zr = -1, id_MLD_003_rr = -1
189189
integer :: id_MLD_EN1 = -1, id_MLD_EN2 = -1, id_MLD_EN3 = -1, id_subMLN2 = -1
190190

191191
! These are handles to diagnostics that are only available in non-ALE layered mode.
@@ -483,15 +483,14 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
483483
if (CS%id_MLD_003 > 0 .or. CS%id_subMLN2 > 0 .or. CS%id_mlotstsq > 0) then
484484
call diagnoseMLDbyDensityDifference(CS%id_MLD_003, h, tv, 0.03*US%kg_m3_to_R, G, GV, US, CS%diag, &
485485
id_N2subML=CS%id_subMLN2, id_MLDsq=CS%id_mlotstsq, dz_subML=CS%dz_subML_N2, &
486-
ref_P_MLD=CS%ref_p_mld, id_ref_z=CS%id_MLD_003_zr)
486+
ref_H_MLD=CS%ref_h_mld, id_ref_z=CS%id_MLD_003_zr, id_ref_rho=CS%id_MLD_003_rr)
487487
endif
488488
if (CS%id_MLD_0125 > 0) then
489489
call diagnoseMLDbyDensityDifference(CS%id_MLD_0125, h, tv, 0.125*US%kg_m3_to_R, G, GV, US, CS%diag, &
490-
ref_P_MLD=CS%ref_p_mld)
490+
ref_H_MLD=CS%ref_h_mld)
491491
endif
492492
if (CS%id_MLD_user > 0) then
493-
call diagnoseMLDbyDensityDifference(CS%id_MLD_user, h, tv, CS%MLDdensityDifference, G, GV, US, CS%diag, &
494-
ref_P_MLD=CS%ref_p_mld)
493+
call diagnoseMLDbyDensityDifference(CS%id_MLD_user, h, tv, CS%MLDdensityDifference, G, GV, US, CS%diag)
495494
endif
496495
if ((CS%id_MLD_EN1 > 0) .or. (CS%id_MLD_EN2 > 0) .or. (CS%id_MLD_EN3 > 0)) then
497496
call diagnoseMLDbyEnergy((/CS%id_MLD_EN1, CS%id_MLD_EN2, CS%id_MLD_EN3/),&
@@ -3225,6 +3224,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
32253224
cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t')
32263225
CS%id_MLD_003_zr = register_diag_field('ocean_model', 'MLD_003_refZ', diag%axesT1, Time, &
32273226
'Depth of reference density for MLD (delta rho = 0.03)', units='m', conversion=US%Z_to_m)
3227+
CS%id_MLD_003_rr = register_diag_field('ocean_model', 'MLD_003_refRho', diag%axesT1, Time, &
3228+
'Reference density for MLD (delta rho = 0.03)', units='kg/m3', conversion=US%R_to_kg_m3)
32283229
CS%id_mlotstsq = register_diag_field('ocean_model', 'mlotstsq', diag%axesT1, Time, &
32293230
long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', &
32303231
standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', &
@@ -3256,12 +3257,12 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
32563257
'Squared buoyancy frequency below mixed layer', units='s-2', conversion=US%s_to_T**2)
32573258
CS%id_MLD_user = register_diag_field('ocean_model', 'MLD_user', diag%axesT1, Time, &
32583259
'Mixed layer depth (used defined)', units='m', conversion=US%Z_to_m)
3259-
if (CS%id_MLD_003 > 0 .or. CS%id_MLD_0125 > 0 .or. CS%id_MLD_0125 > 0) then
3260-
call get_param(param_file, mdl, "PREF_FOR_MLD", CS%ref_p_mld, &
3261-
"Refernced pressure used to calculate the potential density used to find the mixed layer depth "//&
3260+
if (CS%id_MLD_003 > 0 .or. CS%id_MLD_0125 > 0 .or. CS%id_MLD_user > 0) then
3261+
call get_param(param_file, mdl, "HREF_FOR_MLD", CS%ref_h_mld, &
3262+
"Refernced depth used to calculate the potential density used to find the mixed layer depth "//&
32623263
"based on a density difference.", &
3263-
units='Pa', default=0.0, scale=US%Pa_to_RL2_T2)
3264-
!! MLD calculation [R L2 T-2 ~> Pa].
3264+
units='m', default=0.0, scale=US%m_to_Z)
3265+
!! MLD calculation [R L2 T-2 ~> Pa].
32653266
endif
32663267
endif
32673268
call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", CS%MLDdensityDifference, &

0 commit comments

Comments
 (0)