Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Added diagnostic for sub-mixed layer stratification. #21

Merged
merged 3 commits into from
May 22, 2014
Merged

Added diagnostic for sub-mixed layer stratification. #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
825b608
Added diagnostic for sub-mixed layer stratification.
MJHarrison-GFDL May 20, 2014
0981861
Fixed MLD diagnostic and refactored subML_N2
adcroft May 22, 2014
3d8913b
Added subML_N2 to two example diag_tables
adcroft May 22, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions examples/ocean_SIS/MOM6z_SIS_025/diag_table.MOM6
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@
"ocean_model", "KPP_OBLdepth", "KPP_OBLdepth_min", "ocean_month", "all", "min", "none",2
"ocean_model", "MLD_0125", "MLD_0125", "ocean_month", "all", "mean", "none",2
"ocean_model", "MLD_user", "MLD_user", "ocean_month", "all", "mean", "none",2
"ocean_model", "subML_N2", "subML_N2", "ocean_month", "all", "mean", "none",2
"ocean_model", "temp", "temp", "ocean_month", "all", "mean", "none",2
"ocean_model", "salt", "salt", "ocean_month", "all", "mean", "none",2
"ocean_model", "u", "u", "ocean_month", "all", "mean", "none",2
Expand Down
2 changes: 2 additions & 0 deletions examples/solo_ocean/global_ALE/common/diag_table
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,6 +169,8 @@
"ocean_model","KPP_Ksalt","KPP_Ksalt","visc","all",.true.,"none",2
"ocean_model","KPP_NLtransport_heat","KPP_NLtransport_heat","visc","all",.true.,"none",2
"ocean_model","KPP_NLtransport_salt","KPP_NLtransport_salt","visc","all",.true.,"none",2
"ocean_model","MLD_0125","MLD_0125","visc","all",.true.,"none",2
"ocean_model","subML_N2","subML_N2","visc","all",.true.,"none",2

#
# Kinetic Energy Balance Terms:
Expand Down
78 changes: 61 additions & 17 deletions src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ module MOM_diabatic_driver
integer :: id_ea = -1 , id_eb = -1, id_Kd_z = -1, id_Kd_interface = -1
integer :: id_Tdif_z = -1, id_Tadv_z = -1, id_Sdif_z = -1, id_Sadv_z = -1
integer :: id_Tdif = -1, id_Tadv = -1, id_Sdif = -1, id_Sadv = -1
integer :: id_createdH = -1, id_MLD_0125 = -1, id_MLD_user = -1
integer :: id_createdH = -1, id_MLD_0125 = -1, id_MLD_user = -1, id_subMLN2 = -1

type(entrain_diffusive_CS), pointer :: entrain_diffusive_CSp => NULL()
type(bulkmixedlayer_CS), pointer :: bulkmixedlayer_CSp => NULL()
Expand Down Expand Up @@ -1086,8 +1086,8 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS)
if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, ADp%du_dt_dia, CS%diag)
if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, ADp%dv_dt_dia, CS%diag)
if (CS%id_wd > 0) call post_data(CS%id_wd, CDp%diapyc_vel, CS%diag)
if (CS%id_MLD_0125 > 0) then
call diagnoseMLDbyDensityDifference(CS%id_MLD_0125, h, tv, 0.125, G, CS%diag)
if (CS%id_MLD_0125 > 0 .or. CS%id_subMLN2 > 0) then
call diagnoseMLDbyDensityDifference(CS%id_MLD_0125, h, tv, 0.125, G, CS%diag, CS%id_subMLN2)
endif
if (CS%id_MLD_user > 0) then
call diagnoseMLDbyDensityDifference(CS%id_MLD_user, h, tv, CS%MLDdensityDifference, G, CS%diag)
Expand Down Expand Up @@ -1660,6 +1660,8 @@ subroutine diabatic_driver_init(Time, G, param_file, useALEalgorithm, diag, &
if (CS%id_createdH>0) allocate(CS%createdH(isd:ied,jsd:jed))
CS%id_MLD_0125 = register_diag_field('ocean_model','MLD_0125',diag%axesT1,Time, &
'Mixed layer depth (delta rho = 0.125)', 'meter')
CS%id_subMLN2 = register_diag_field('ocean_model','subML_N2',diag%axesT1,Time, &
'Squared buoyancy frequency below Mixed layer ', 's-2')
CS%id_MLD_user = register_diag_field('ocean_model','MLD_user',diag%axesT1,Time, &
'Mixed layer depth (used defined)', 'meter')
call get_param(param_file, mod, "DIAG_MLD_DENSITY_DIFF", CS%MLDdensityDifference, &
Expand Down Expand Up @@ -1884,30 +1886,65 @@ subroutine find_uv_at_h(u, v, h, u_h, v_h, G, ea, eb)
end subroutine find_uv_at_h

!> Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface.
subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, diagPtr)
subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, diagPtr, id_N2subML)
integer, intent(in) :: id_MLD !< Handle (ID) of MLD diagnostic
real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h !< Layer thickness
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics type
real, intent(in) :: densityDiff !< Density difference to determine MLD (kg/m3)
type(ocean_grid_type), intent(in) :: G !< Grid type
type(diag_ctrl), pointer :: diagPtr !< Diagnostics structure
integer, optional, intent(in) :: id_N2subML !< Optional handle (ID) of subML stratification
! Local variables
real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK, dK, dKm1, pRef
real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK, dK, dKm1, pRef_MLD ! Used for MLD
real, dimension(SZI_(G)) :: rhoAtK, rho1, d1, pRef_N2 ! Used for N2
real, dimension(SZI_(G), SZJ_(G)) :: MLD ! Diagnosed mixed layer depth
integer :: i, j, is, ie, js, je, k, nz
real, dimension(SZI_(G), SZJ_(G)) :: subMLN2 ! Diagnosed stratification below ML
real, parameter :: dz_subML = 50. ! Depth below ML over which to diagnose stratification (m)
integer :: i, j, is, ie, js, je, k, nz, id_N2
real :: aFac, ddRho

id_N2 = -1
if (PRESENT(id_N2subML)) id_N2 = id_N2subML

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
pRef(:) = 0.
pRef_MLD(:) = 0. ; pRef_N2(:) = 0.
do j = js, je
dK(:) = 0.5 * h(:,j,1) ! Center of surface layer
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef, rhoSurf, is, ie-is+1, tv%eqn_of_state)
dK(:) = 0.5 * h(:,j,1) * G%H_to_m ! Depth of center of surface layer
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, is, ie-is+1, tv%eqn_of_state)
deltaRhoAtK(:) = 0.
MLD(:,j) = 0.
if (id_N2>0) then
subMLN2(:,j) = 0.
rho1(:) = 0.
d1(:) = 0.
pRef_N2(:) = G%g_Earth * G%Rho0 * h(:,j,1) * G%H_to_m ! Boussinesq approximation!!!! ?????
endif
do k = 2, nz
dKm1(:) = dK(:) ! Center of layer K-1
dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) ! Center of layer K
deltaRhoAtKm1(:) = deltaRhoAtK(:)
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef, deltaRhoAtK, is, ie-is+1, tv%eqn_of_state)
dKm1(:) = dK(:) ! Depth of center of layer K-1
dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) * G%H_to_m ! Depth of center of layer K

! Stratification, N2, immediately below the mixed layer, averaged over at least 50 m.
if (id_N2>0) then
pRef_N2(:) = pRef_N2(:) + G%g_Earth * G%Rho0 * h(:,j,k) * G%H_to_m ! Boussinesq approximation!!!! ?????
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_N2, rhoAtK, is, ie-is+1, tv%eqn_of_state)
do i = is, ie
if (MLD(i,j)>0. .and. subMLN2(i,j)==0.) then ! This block is below the mixed layer
if (d1(i)==0.) then ! Record the density, depth and pressure, immediately below the ML
rho1(i) = rhoAtK(i)
d1(i) = dK(i)
! Use pressure at the bottom of the upper layer used in calculating d/dz rho
pRef_N2(i) = pRef_N2(i) + G%g_Earth * G%Rho0 * h(i,j,k) * G%H_to_m ! Boussinesq approximation!!!! ?????
endif
if (d1(i)>0. .and. dK(i)-d1(i)>=dz_subML) then
subMLN2(i,j) = G%g_Earth/ G%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i))
endif
endif
enddo ! i-loop
endif ! id_N2>0

! Mixed-layer depth, using sigma-0 (surface reference pressure)
deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, is, ie-is+1, tv%eqn_of_state)
deltaRhoAtK(:) = deltaRhoAtK(:) - rhoSurf(:) ! Density difference between layer K and surface
do i = is, ie
ddRho = deltaRhoAtK(i) - deltaRhoAtKm1(i)
Expand All @@ -1916,13 +1953,20 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, diagPtr
aFac = ( densityDiff - deltaRhoAtKm1(i) ) / ddRho
MLD(i,j) = dK(i) * aFac + dKm1(i) * (1. - aFac)
endif
enddo
enddo
enddo ! i-loop

enddo ! k-loop
do i = is, ie
if ((MLD(i,j)==0.) .and. (deltaRhoAtK(i)<densityDiff)) MLD(i,j) = dK(i) ! Assume mixing to the bottom
! if (id_N2>0 .and. subMLN2(i,j)==0. .and. d1(i)>0. .and. dK(i)-d1(i)>0.) then
! ! Use what ever stratification we can, measured over what ever distance is available
! subMLN2(i,j) = G%g_Earth/ G%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i))
! endif
enddo
enddo
if (id_MLD > 0) call post_data(id_MLD, MLD * G%H_to_m, diagPtr) ! Convert to meters
enddo ! j-loop
if (id_MLD > 0) call post_data(id_MLD, MLD, diagPtr)
if (id_N2 > 0) call post_data(id_N2, subMLN2 , diagPtr)

end subroutine diagnoseMLDbyDensityDifference

subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv)
Expand Down