From 825b608a663e6ddac280798aa3653b6904ac564b Mon Sep 17 00:00:00 2001 From: Matthew Harrison Date: Tue, 20 May 2014 11:17:31 -0400 Subject: [PATCH 1/3] Added diagnostic for sub-mixed layer stratification. + New diagnostic is "subML_N2" + only calculated if "MLD_0125" diagnostic is enabled + Based on stratification measured from ML base to + 50m below the mixed layer base. --- .../vertical/MOM_diabatic_driver.F90 | 47 ++++++++++++++++--- 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 3e8c507915..fd71817e28 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -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() @@ -1087,7 +1087,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, CS) 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) + 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) @@ -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, & @@ -1884,18 +1886,32 @@ 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, intent(in), optional :: id_N2subML !< Optional handle (ID) of subML stratification ! Local variables real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK, dK, dKm1, pRef real, dimension(SZI_(G), SZJ_(G)) :: MLD ! Diagnosed mixed layer depth + real, dimension(SZI_(G), SZJ_(G)) :: subMLN2 ! Diagnosed stratification below ML + + real, parameter :: dz_subML = 50. ! Depth below ML over which to diagnose stratification + integer :: i, j, is, ie, js, je, k, nz - real :: aFac, ddRho + real :: aFac, ddRho,dzML,dzMLm1,ddZ + integer :: idn2 + + + idn2 = -1 + if (PRESENT(id_N2subML)) then + idn2 = id_N2subML + subMLN2(:,:)=0. + endif + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke pRef(:) = 0. do j = js, je @@ -1903,6 +1919,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, diagPtr call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef, rhoSurf, is, ie-is+1, tv%eqn_of_state) deltaRhoAtK(:) = 0. MLD(:,j) = 0. + if (idn2>0) subMLN2(:,j) = 0. 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 @@ -1917,12 +1934,28 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, diagPtr MLD(i,j) = dK(i) * aFac + dKm1(i) * (1. - aFac) endif enddo - enddo - do i = is, ie - if ((MLD(i,j)==0.) .and. (deltaRhoAtK(i) 0)) then + do i = is, ie + if (MLD(i,j)>0.) then + dzMLm1 = (dKm1(i) - MLD(i,j)) * G%H_to_m + dzML = (dK(i)-MLD(i,j)) * G%H_to_m + if (dzML >= dz_subML .and. dzMLm1 < dz_subML) then + ddRho = (deltaRhoAtK(i) - deltaRhoAtKm1(i))/(dK(i)-dKm1(i))*(dz_subML-dzMLm1) + subMLN2(i,j) = G%g_Earth/ G%Rho0 * (deltaRhoAtKm1(i)-densityDiff+ddRho) / dz_subML + endif + endif + enddo + endif + + do i = is, ie + if ((MLD(i,j)==0.) .and. (deltaRhoAtK(i) 0) call post_data(id_MLD, MLD * G%H_to_m, diagPtr) ! Convert to meters + if (idn2 > 0) call post_data(idn2, subMLN2 , diagPtr) + end subroutine diagnoseMLDbyDensityDifference subroutine applyBoundaryFluxes(CS, G, dt, fluxes, optics, ea, h, tv) From 098186140a77ab14f6310d307ba583b92c6420e6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 22 May 2014 11:36:18 -0400 Subject: [PATCH 2/3] Fixed MLD diagnostic and refactored subML_N2 - Adding subML_N2 diagnostic inadvertently brok the MLD diagnostic. - Refactored subML_N2 diagnostic. - Added pressure dependence in N2 calculation (using pressure below ML). - No answer changes. --- .../vertical/MOM_diabatic_driver.F90 | 99 ++++++++++--------- 1 file changed, 55 insertions(+), 44 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index fd71817e28..7725b339af 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1086,7 +1086,7 @@ 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 + 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 @@ -1893,38 +1893,58 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, diagPtr 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, intent(in), optional :: id_N2subML !< Optional handle (ID) of subML stratification + 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 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 - real, parameter :: dz_subML = 50. ! Depth below ML over which to diagnose stratification - - integer :: i, j, is, ie, js, je, k, nz - real :: aFac, ddRho,dzML,dzMLm1,ddZ - integer :: idn2 - - - idn2 = -1 - if (PRESENT(id_N2subML)) then - idn2 = id_N2subML - subMLN2(:,:)=0. - endif + 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 (idn2>0) subMLN2(:,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) @@ -1933,28 +1953,19 @@ 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 - - if ((idn2 > 0)) then - do i = is, ie - if (MLD(i,j)>0.) then - dzMLm1 = (dKm1(i) - MLD(i,j)) * G%H_to_m - dzML = (dK(i)-MLD(i,j)) * G%H_to_m - if (dzML >= dz_subML .and. dzMLm1 < dz_subML) then - ddRho = (deltaRhoAtK(i) - deltaRhoAtKm1(i))/(dK(i)-dKm1(i))*(dz_subML-dzMLm1) - subMLN2(i,j) = G%g_Earth/ G%Rho0 * (deltaRhoAtKm1(i)-densityDiff+ddRho) / dz_subML - endif - endif - enddo - endif - - do i = is, ie - if ((MLD(i,j)==0.) .and. (deltaRhoAtK(i)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 - if (idn2 > 0) call post_data(idn2, subMLN2 , diagPtr) + 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 From 3d8913b70b3e45c01bae416878ef8b72774cffc3 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 22 May 2014 11:39:58 -0400 Subject: [PATCH 3/3] Added subML_N2 to two example diag_tables - Always good to exercise new code! --- examples/ocean_SIS/MOM6z_SIS_025/diag_table.MOM6 | 1 + examples/solo_ocean/global_ALE/common/diag_table | 2 ++ 2 files changed, 3 insertions(+) diff --git a/examples/ocean_SIS/MOM6z_SIS_025/diag_table.MOM6 b/examples/ocean_SIS/MOM6z_SIS_025/diag_table.MOM6 index f600997fef..6961908d25 100644 --- a/examples/ocean_SIS/MOM6z_SIS_025/diag_table.MOM6 +++ b/examples/ocean_SIS/MOM6z_SIS_025/diag_table.MOM6 @@ -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 diff --git a/examples/solo_ocean/global_ALE/common/diag_table b/examples/solo_ocean/global_ALE/common/diag_table index 891158ecfe..a565a7d69d 100644 --- a/examples/solo_ocean/global_ALE/common/diag_table +++ b/examples/solo_ocean/global_ALE/common/diag_table @@ -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: