diff --git a/src/diagnostics/MOM_diagnose_KdWork.F90 b/src/diagnostics/MOM_diagnose_KdWork.F90 index b981da9af1..34d5407509 100644 --- a/src/diagnostics/MOM_diagnose_KdWork.F90 +++ b/src/diagnostics/MOM_diagnose_KdWork.F90 @@ -32,19 +32,9 @@ module MOM_diagnose_kdwork !! routine use. type vbf_CS ! 3d varying Kd contributions + ! The following are all allocatable arrays that store copies of process driven Kd, so that + ! the process driven buoyancy flux and work can be derived at the end of the time step. real, pointer, dimension(:,:,:) :: & - Bflx_salt => NULL(), & !< Salinity contribution to buoyancy flux at interfaces - !! [H Z T-3 ~> m2 s-3 or W m-3] - Bflx_temp => NULL(), & !< Temperature contribution to buoyancy flux at interfaces - !! [H Z T-3 ~> m2 s-3 or W m-3] - Bflx_salt_dz => NULL(), & !< Salinity contribution to integral of buoyancy flux over layer - !! [H Z2 T-3 ~> m3 s-3 or W m-2] - Bflx_temp_dz => NULL(), & !< Temperature contribution to integral of buoyancy flux over layer - !! [H Z2 T-3 ~> m3 s-3 or W m-2] - ! The following are all allocatable arrays that store copies of process driven Kd, so that - ! the process driven buoyancy flux and work can be derived at the end of the time step. - Kd_salt => NULL(), & !< total diapycnal diffusivity of salt at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] - Kd_temp => NULL(), & !< total diapycnal diffusivity of heat at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_BBL => NULL(), & !< diapycnal diffusivity due to BBL at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_ePBL => NULL(), & !< diapycnal diffusivity due to ePBL at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_KS => NULL(), & !< diapycnal diffusivity due to Kappa Shear at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] @@ -52,7 +42,7 @@ module MOM_diagnose_kdwork Kd_ddiff_S => NULL(), &!< diapycnal diffusivity due to double diffusion of salt at interfaces !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_ddiff_T => NULL(), &!< diapycnal diffusivity due to double diffusion of heat at interfaces - !![H Z T-1 ~> m2 s-1 or kg m-1 s-1] + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_leak => NULL(), & !< diapycnal diffusivity due to Kd_leak at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_quad => NULL(), & !< diapycnal diffusivity due to Kd_quad at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_itidal => NULL(), & !< diapycnal diffusivity due to Kd_itidal at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] @@ -104,30 +94,51 @@ module MOM_diagnose_kdwork contains !> Loop over all implemented diffusivities to diagnose and output Kd Work/buoyancy fluxes -subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) +subroutine KdWork_Diagnostics(G, GV, US, diag, VBF, Kd_salt, Kd_temp, N2_Salt, N2_Temp, dz) 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 type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output type (vbf_CS), intent(inout) :: VBF !< Vertical buoyancy flux structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & - intent(in) :: N2_Salt !< Buoyancy frequency [T-2 ~> s-2] + intent(in) :: Kd_salt !< Total diapycnal diffusivity of salt and passive + !! tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & + intent(in) :: Kd_temp !< Total diapycnal diffusivity of heat at interfaces + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & + intent(in) :: N2_Salt !< Buoyancy frequency contribution from salinity [T-2 ~> s-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & - intent(in) :: N2_Temp !< Buoyancy frequency [T-2 ~> s-2] + intent(in) :: N2_Temp !< Buoyancy frequency contribution from temperature [T-2 ~> s-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: dz !< Grid spacing [Z ~> m] ! Work arrays for computing buoyancy flux integrals - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: work3d_i - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work3d_l - real, dimension(SZI_(G),SZJ_(G)) :: work2d, work2d_salt, work2d_temp - real :: work, work_salt, work_temp + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & + Bflx_salt, & ! Salinity contribution to buoyancy flux at interfaces [H Z T-3 ~> m2 s-3 or W m-3] + Bflx_temp, & ! Temperature contribution to buoyancy flux at interfaces [H Z T-3 ~> m2 s-3 or W m-3] + work3d_i ! A work array for sums of buoyancy fluxes at interfaces [H Z T-3 ~> m2 s-3 or W m-3] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: & + Bflx_salt_dz, & ! Salinity contribution to integral of buoyancy flux over layer + ! [H Z2 T-3 ~> m3 s-3 or W m-2] + Bflx_temp_dz, & ! Temperature contribution to integral of buoyancy flux over layer + ! [H Z2 T-3 ~> m3 s-3 or W m-2] + work3d_l ! A work array for sums of buoyancy fluxes at layers [H Z2 T-3 ~> m3 s-3 or W m-2] + real, dimension(SZI_(G),SZJ_(G)) :: & + work2d, & ! Vertically integrated buoyancy fluxes [H Z2 T-3 ~> m3 s-3 or W m-2] + work2d_salt, & ! Vertically integrated buoyancy fluxes due to salt fluxes [H Z2 T-3 ~> m3 s-3 or W m-2] + work2d_temp ! Vertically integrated buoyancy fluxes due to heat fluxes [H Z2 T-3 ~> m3 s-3 or W m-2] + real :: work ! Global integrated buoyancy fluxes [H Z2 L2 T-3 ~> m5 s-3 or W] + real :: work_salt ! Global integrated buoyancy fluxes due to salt fluxes [H Z2 L2 T-3 ~> m5 s-3 or W] + real :: work_temp ! Global integrated buoyancy fluxes due to heat fluxes [H Z2 L2 T-3 ~> m5 s-3 or W] + real :: work_unscale ! The combination of dimensional scaling factors for vertically integrated + ! buoyancy fluxes [H Z2 s3 m-3 T-3 or H Z2 m2 T-3 W-1 ~> 1] integer :: i, j, k, nz, isc, iec, jsc, jec isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec - nz = GV%ke + work_unscale = GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3 ! Compute total fluxes if (VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_salt_dz>0 .or. VBF%id_Bdif_temp_dz>0 .or. & @@ -136,35 +147,33 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Do Salt if (VBF%id_Bdif_salt_dz>0 .or. VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_salt>0 .or. VBF%id_Bdif>0 .or. & VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_salt_idz>0 .or. VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_salt_idV>0) & - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_salt, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_salt, Kd_salt, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) ! Do Temp if (VBF%id_Bdif_temp_dz>0 .or. VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_temp>0 .or. VBF%id_Bdif>0 .or. & VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_temp_idz>0 .or. VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_temp_idV>0) & - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_temp, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_temp, Kd_temp, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_temp_idz>0 .or. VBF%id_Bdif_idz>0) then work2d_temp(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d_temp(i,j) = work2d_temp(i,j) + VBF%Bflx_temp_dz(i,j,k) + work2d_temp(i,j) = work2d_temp(i,j) + Bflx_temp_dz(i,j,k) enddo ; enddo ; enddo endif if (VBF%id_Bdif_temp_idV>0 .or. VBF%id_Bdif_idV>0) then work_temp = 0.0 do k = 1,nz - work_temp = work_temp + global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, & - tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + work_temp = work_temp + global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) enddo endif if (VBF%id_Bdif_salt_idz>0 .or. VBF%id_Bdif_idz>0) then work2d_salt(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d_salt(i,j) = work2d_salt(i,j) + VBF%Bflx_salt_dz(i,j,k) + work2d_salt(i,j) = work2d_salt(i,j) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo endif if (VBF%id_Bdif_salt_idV>0 .or. VBF%id_Bdif_idV>0) then work_salt = 0.0 do k = 1,nz - work_salt = work_salt + global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, & - tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + work_salt = work_salt + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale) enddo endif work = work_temp + work_salt @@ -174,26 +183,26 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) elseif (VBF%id_Bdif>0 .or. VBF%id_Bdif_salt>0 .or. VBF%id_Bdif_temp>0) then ! Not doing vertical integrals ! Do Salt if (VBF%id_Bdif_salt>0 .or. VBF%id_Bdif>0) & - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_salt, VBF%Bflx_salt) + call diagnoseKdWork(G, GV, N2_salt, Kd_salt, Bflx_salt) if (VBF%id_Bdif_temp>0 .or. VBF%id_Bdif>0) & - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_temp, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_temp, Kd_temp, Bflx_temp) endif ! Post total fluxes - if (VBF%id_Bdif_salt>0) call post_data(VBF%id_Bdif_salt, VBF%Bflx_salt, diag) - if (VBF%id_Bdif_temp>0) call post_data(VBF%id_Bdif_temp, VBF%Bflx_temp, diag) + if (VBF%id_Bdif_salt>0) call post_data(VBF%id_Bdif_salt, Bflx_salt, diag) + if (VBF%id_Bdif_temp>0) call post_data(VBF%id_Bdif_temp, Bflx_temp, diag) if (VBF%id_Bdif>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif, work3d_i, diag) endif - if (VBF%id_Bdif_salt_dz>0) call post_data(VBF%id_Bdif_salt_dz, VBF%Bflx_salt_dz, diag) - if (VBF%id_Bdif_temp_dz>0) call post_data(VBF%id_Bdif_temp_dz, VBF%Bflx_temp_dz, diag) + if (VBF%id_Bdif_salt_dz>0) call post_data(VBF%id_Bdif_salt_dz, Bflx_salt_dz, diag) + if (VBF%id_Bdif_temp_dz>0) call post_data(VBF%id_Bdif_temp_dz, Bflx_temp_dz, diag) if (VBF%id_Bdif_dz>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz, work3d_l, diag) endif @@ -206,38 +215,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute ePBL fluxes if (VBF%id_Bdif_dz_ePBL>0.or.VBF%id_Bdif_idz_ePBL>0.or.VBF%id_Bdif_idV_ePBL>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ePBL, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ePBL, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ePBL, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ePBL, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_ePBL>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_ePBL>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_ePBL>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ePBL, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ePBL, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ePBL, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ePBL, Bflx_temp) endif ! Post ePBL fluxes if (VBF%id_Bdif_ePBL>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_ePBL, work3d_i, diag) endif if (VBF%id_Bdif_dz_ePBL>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_ePBL, work3d_l, diag) endif @@ -246,38 +254,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute BBL fluxes if (VBF%id_Bdif_dz_BBL>0.or.VBF%id_Bdif_idz_BBL>0.or.VBF%id_Bdif_idV_BBL>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_BBL, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_BBL, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_BBL, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_BBL, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_BBL>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_BBL>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_BBL>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_BBL, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_BBL, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_BBL, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_BBL, Bflx_temp) endif ! Post BBL fluxes if (VBF%id_Bdif_BBL>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_BBL, work3d_i, diag) endif if (VBF%id_Bdif_dz_BBL>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_BBL, work3d_l, diag) endif @@ -286,38 +293,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kappa Shear fluxes if (VBF%id_Bdif_dz_KS>0.or.VBF%id_Bdif_idz_KS>0.or.VBF%id_Bdif_idV_KS>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_KS, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_KS, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_KS, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_KS, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_KS>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_KS>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_KS>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_KS, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_KS, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_KS, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_KS, Bflx_temp) endif ! Post Kappa Shear fluxes if (VBF%id_Bdif_KS>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_KS, work3d_i, diag) endif if (VBF%id_Bdif_dz_KS>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_KS, work3d_l, diag) endif @@ -326,38 +332,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute bkgnd fluxes if (VBF%id_Bdif_dz_bkgnd>0.or.VBF%id_Bdif_idz_bkgnd>0.or.VBF%id_Bdif_idV_bkgnd>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_bkgnd, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_bkgnd, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_bkgnd, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_bkgnd, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_bkgnd>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_bkgnd>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_bkgnd>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_bkgnd, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_bkgnd, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_bkgnd, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_bkgnd, Bflx_temp) endif ! Post bkgnd fluxes if (VBF%id_Bdif_bkgnd>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_bkgnd, work3d_i, diag) endif if (VBF%id_Bdif_dz_bkgnd>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_bkgnd, work3d_l, diag) endif @@ -366,85 +371,82 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute double diffusion fluxes if (VBF%id_Bdif_dz_ddiff_temp>0.or.VBF%id_Bdif_idz_ddiff_temp>0.or.VBF%id_Bdif_idV_ddiff_temp>0) then - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ddiff_T, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ddiff_T, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_ddiff_temp>0) then work2d_temp(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d_temp(i,j) = work2d_temp(i,j) + VBF%Bflx_temp_dz(i,j,k) + work2d_temp(i,j) = work2d_temp(i,j) + Bflx_temp_dz(i,j,k) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_ddiff_temp>0) then work_temp = 0.0 do k = 1,nz - work_temp = work_temp + global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, & - tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + work_temp = work_temp + global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) enddo endif elseif (VBF%id_Bdif_ddiff_temp>0) then - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ddiff_T, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_ddiff_T, Bflx_temp) endif if (VBF%id_Bdif_dz_ddiff_salt>0.or.VBF%id_Bdif_idz_ddiff_salt>0.or.VBF%id_Bdif_idV_ddiff_salt>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ddiff_S, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ddiff_S, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) if (VBF%id_Bdif_idz_ddiff_salt>0) then work2d_salt(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d_salt(i,j) = work2d_salt(i,j) + VBF%Bflx_salt_dz(i,j,k) + work2d_salt(i,j) = work2d_salt(i,j) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_ddiff_salt>0) then work_salt = 0.0 do k = 1,nz - work_salt = work_salt + global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, & - tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + work_salt = work_salt + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale) enddo endif elseif (VBF%id_Bdif_ddiff_salt>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ddiff_S, VBF%Bflx_salt) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_ddiff_S, Bflx_salt) endif ! Post double diffusion fluxes - if (VBF%id_Bdif_ddiff_temp>0) call post_data(VBF%id_Bdif_ddiff_temp, VBF%Bflx_temp, diag) - if (VBF%id_Bdif_dz_ddiff_temp>0) call post_data(VBF%id_Bdif_dz_ddiff_temp, VBF%Bflx_temp_dz, diag) + if (VBF%id_Bdif_ddiff_temp>0) call post_data(VBF%id_Bdif_ddiff_temp, Bflx_temp, diag) + if (VBF%id_Bdif_dz_ddiff_temp>0) call post_data(VBF%id_Bdif_dz_ddiff_temp, Bflx_temp_dz, diag) if (VBF%id_Bdif_idz_ddiff_temp>0) call post_data(VBF%id_Bdif_idz_ddiff_temp, work2d_temp, diag) if (VBF%id_Bdif_idV_ddiff_temp>0) call post_data(VBF%id_Bdif_idV_ddiff_temp, work_temp, diag) - if (VBF%id_Bdif_ddiff_salt>0) call post_data(VBF%id_Bdif_ddiff_salt, VBF%Bflx_salt, diag) - if (VBF%id_Bdif_dz_ddiff_salt>0) call post_data(VBF%id_Bdif_dz_ddiff_salt, VBF%Bflx_salt_dz, diag) + if (VBF%id_Bdif_ddiff_salt>0) call post_data(VBF%id_Bdif_ddiff_salt, Bflx_salt, diag) + if (VBF%id_Bdif_dz_ddiff_salt>0) call post_data(VBF%id_Bdif_dz_ddiff_salt, Bflx_salt_dz, diag) if (VBF%id_Bdif_idz_ddiff_salt>0) call post_data(VBF%id_Bdif_idz_ddiff_salt, work2d_salt, diag) if (VBF%id_Bdif_idV_ddiff_salt>0) call post_data(VBF%id_Bdif_idV_ddiff_salt, work_salt, diag) ! Compute Kd_leak fluxes if (VBF%id_Bdif_dz_leak>0.or.VBF%id_Bdif_idz_leak>0.or.VBF%id_Bdif_idV_leak>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_leak, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_leak, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_leak, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_leak, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_leak>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_leak>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_leak>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_leak, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_leak, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_leak, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_leak, Bflx_temp) endif ! Post Kd_leak fluxes if (VBF%id_Bdif_leak>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_leak, work3d_i, diag) endif if (VBF%id_Bdif_dz_leak>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_leak, work3d_l, diag) endif @@ -453,38 +455,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kd_quad fluxes if (VBF%id_Bdif_dz_quad>0.or.VBF%id_Bdif_idz_quad>0.or.VBF%id_Bdif_idV_quad>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_quad, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_quad, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_quad, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_quad, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_quad>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_quad>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_quad>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_quad, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_quad, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_quad, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_quad, Bflx_temp) endif ! Post Kd_quad fluxes if (VBF%id_Bdif_quad>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_quad, work3d_i, diag) endif if (VBF%id_Bdif_dz_quad>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_quad, work3d_l, diag) endif @@ -493,38 +494,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kd_itidal fluxes if (VBF%id_Bdif_dz_itidal>0.or.VBF%id_Bdif_idz_itidal>0.or.VBF%id_Bdif_idV_itidal>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itidal, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itidal, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itidal, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itidal, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_itidal>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_itidal>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_itidal>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itidal, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itidal, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itidal, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itidal, Bflx_temp) endif ! Post Kd_itidal fluxes if (VBF%id_Bdif_itidal>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_itidal, work3d_i, diag) endif if (VBF%id_Bdif_dz_itidal>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k)+VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k)+Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_itidal, work3d_l, diag) endif @@ -533,38 +533,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kd_Froude fluxes if (VBF%id_Bdif_dz_Froude>0.or.VBF%id_Bdif_idz_Froude>0.or.VBF%id_Bdif_idV_Froude>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Froude, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Froude, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Froude, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Froude, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_Froude>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_Froude>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_Froude>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Froude, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Froude, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Froude, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Froude, Bflx_temp) endif ! Post Kd_Froude fluxes if (VBF%id_Bdif_Froude>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_Froude, work3d_i, diag) endif if (VBF%id_Bdif_dz_Froude>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_Froude, work3d_l, diag) endif @@ -573,38 +572,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kd_slope fluxes if (VBF%id_Bdif_dz_slope>0.or.VBF%id_Bdif_idz_slope>0.or.VBF%id_Bdif_idV_slope>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_slope, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_slope, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_slope, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_slope, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_slope>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_slope>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_slope>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_slope, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_slope, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_slope, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_slope, Bflx_temp) endif ! Post Kd_slope fluxes if (VBF%id_Bdif_slope>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_slope, work3d_i, diag) endif if (VBF%id_Bdif_dz_slope>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_slope, work3d_l, diag) endif @@ -613,38 +611,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kd_lowmode fluxes if (VBF%id_Bdif_dz_lowmode>0.or.VBF%id_Bdif_idz_lowmode>0.or.VBF%id_Bdif_idV_lowmode>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_lowmode, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_lowmode, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_lowmode, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_lowmode, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_lowmode>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_lowmode>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_lowmode>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_lowmode, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_lowmode, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_lowmode, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_lowmode, Bflx_temp) endif ! Post Kd_lowmode fluxes if (VBF%id_Bdif_lowmode>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_lowmode, work3d_i, diag) endif if (VBF%id_Bdif_dz_lowmode>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_lowmode, work3d_l, diag) endif @@ -653,38 +650,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kd_Niku fluxes if (VBF%id_Bdif_dz_Niku>0 .or. VBF%id_Bdif_idz_Niku>0 .or. VBF%id_Bdif_idV_Niku>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Niku, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Niku, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Niku, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Niku, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_Niku>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_Niku>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_Niku>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Niku, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Niku, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_Niku, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_Niku, Bflx_temp) endif ! Post Kd_Niku fluxes if (VBF%id_Bdif_Niku>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_lowmode, work3d_i, diag) endif if (VBF%id_Bdif_dz_Niku>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_Niku, work3d_l, diag) endif @@ -693,38 +689,37 @@ subroutine KdWork_Diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz) ! Compute Kd_itides fluxes if (VBF%id_Bdif_dz_itides>0 .or. VBF%id_Bdif_idz_itides>0 .or. VBF%id_Bdif_idV_itides>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itides, VBF%Bflx_salt, dz=dz, Bdif_flx_dz=VBF%Bflx_salt_dz) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itides, VBF%Bflx_temp, dz=dz, Bdif_flx_dz=VBF%Bflx_temp_dz) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itides, Bflx_salt, dz=dz, Bdif_flx_dz=Bflx_salt_dz) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itides, Bflx_temp, dz=dz, Bdif_flx_dz=Bflx_temp_dz) if (VBF%id_Bdif_idz_itides>0) then work2d(:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work2d(i,j) = work2d(i,j) + (VBF%Bflx_salt_dz(i,j,k) + VBF%Bflx_temp_dz(i,j,k)) + work2d(i,j) = work2d(i,j) + (Bflx_salt_dz(i,j,k) + Bflx_temp_dz(i,j,k)) enddo ; enddo ; enddo endif if (VBF%id_Bdif_idV_itides>0) then work = 0.0 do k = 1,nz - work = work + & - (global_area_integral(VBF%Bflx_temp_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) + & - global_area_integral(VBF%Bflx_salt_dz(:,:,k), G, tmp_scale=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3)) + work = work + (global_area_integral(Bflx_temp_dz(:,:,k), G, tmp_scale=work_unscale) + & + global_area_integral(Bflx_salt_dz(:,:,k), G, tmp_scale=work_unscale)) enddo endif elseif (VBF%id_Bdif_itides>0) then - call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itides, VBF%Bflx_salt) - call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itides, VBF%Bflx_temp) + call diagnoseKdWork(G, GV, N2_salt, VBF%Kd_itides, Bflx_salt) + call diagnoseKdWork(G, GV, N2_temp, VBF%Kd_itides, Bflx_temp) endif ! Post Kd_itides fluxes if (VBF%id_Bdif_itides>0) then work3d_i(:,:,:) = 0.0 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec - work3d_i(i,j,k) = VBF%Bflx_temp(i,j,k) + VBF%Bflx_salt(i,j,k) + work3d_i(i,j,k) = Bflx_temp(i,j,k) + Bflx_salt(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_itides, work3d_i, diag) endif if (VBF%id_Bdif_dz_itides>0) then work3d_l(:,:,:) = 0.0 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec - work3d_l(i,j,k) = VBF%Bflx_temp_dz(i,j,k) + VBF%Bflx_salt_dz(i,j,k) + work3d_l(i,j,k) = Bflx_temp_dz(i,j,k) + Bflx_salt_dz(i,j,k) enddo ; enddo ; enddo call post_data(VBF%id_Bdif_dz_itides, work3d_l, diag) endif @@ -776,22 +771,6 @@ subroutine Allocate_VBF_CS(G, GV, VBF) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke - if (VBF%do_bflx_salt) & - allocate(VBF%Bflx_salt(isd:ied,jsd:jed,nz+1), source=0.0) - if (VBF%do_bflx_salt_dz) & - allocate(VBF%Bflx_salt_dz(isd:ied,jsd:jed,nz), source=0.0) - if (VBF%do_bflx_temp) & - allocate(VBF%Bflx_temp(isd:ied,jsd:jed,nz+1), source=0.0) - if (VBF%do_bflx_temp_dz) & - allocate(VBF%Bflx_temp_dz(isd:ied,jsd:jed,nz), source=0.0) - - if (VBF%id_Bdif_salt_dz>0 .or. VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_salt>0 .or. VBF%id_Bdif>0 .or. & - VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_salt_idz>0 .or. VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_salt_idV>0) & - allocate(VBF%Kd_salt(isd:ied,jsd:jed,nz+1), source=0.0) - if (VBF%id_Bdif_temp_dz>0 .or. VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_temp>0 .or. VBF%id_Bdif>0 .or. & - VBF%id_Bdif_idz>0 .or. VBF%id_Bdif_temp_idz>0 .or. VBF%id_Bdif_idV>0 .or. VBF%id_Bdif_temp_idV>0) & - allocate(VBF%Kd_temp(isd:ied,jsd:jed,nz+1), source=0.0) - if (VBF%id_Bdif_BBL>0 .or. VBF%id_Bdif_dz_BBL>0 .or. VBF%id_Bdif_idz_BBL>0 .or. VBF%id_Bdif_idV_BBL>0) & allocate(VBF%Kd_BBL(isd:ied,jsd:jed,nz+1), source=0.0) if (VBF%id_Bdif_ePBL>0 .or. VBF%id_Bdif_dz_ePBL>0 .or. VBF%id_Bdif_idz_ePBL>0 .or. VBF%id_Bdif_idV_ePBL>0) & @@ -827,18 +806,6 @@ end subroutine Allocate_VBF_CS subroutine Deallocate_VBF_CS(VBF) type (vbf_CS), intent(inout) :: VBF !< Vertical buoyancy flux structure - if (associated(VBF%Bflx_salt)) & - deallocate(VBF%Bflx_salt) - if (associated(VBF%Bflx_temp)) & - deallocate(VBF%Bflx_temp) - if (associated(VBF%Bflx_salt_dz)) & - deallocate(VBF%Bflx_salt_dz) - if (associated(VBF%Bflx_temp_dz)) & - deallocate(VBF%Bflx_temp_dz) - if (associated(VBF%Kd_salt)) & - deallocate(VBF%Kd_salt) - if (associated(VBF%Kd_temp)) & - deallocate(VBF%Kd_temp) if (associated(VBF%Kd_BBL)) & deallocate(VBF%Kd_BBL) if (associated(VBF%Kd_ePBL)) & @@ -871,242 +838,259 @@ subroutine Deallocate_VBF_CS(VBF) end subroutine Deallocate_VBF_CS !> Handles all KdWork diagnostics and flags which calculations should be done. -subroutine KdWork_init(Time, G,GV,US,diag,VBF,Use_KdWork_diag) - type(time_type), target :: Time !< model time +subroutine KdWork_init(Time, G, GV, US, use_EOS, diag, VBF, Use_KdWork_diag) + type(time_type), target :: Time !< model time type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + logical, intent(in) :: use_EOS !< If true, density is calculated from temperature + !! and salinity using an equation of state. type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output type (vbf_CS), pointer, intent(inout) :: VBF !< Vertical buoyancy flux structure logical, intent(out) :: Use_KdWork_diag !< Flag if any output was turned on + real :: bflx_unscale ! The combination of dimensional scaling factors for buoyancy + ! fluxes [H Z s3 m-2 T-3 or H Z m3 T-3 W-1 ~> 1] + real :: work_unscale ! The combination of dimensional scaling factors for vertically integrated + ! buoyancy fluxes [H Z2 s3 m-3 T-3 or H Z2 m2 T-3 W-1 ~> 1] + allocate(VBF) + bflx_unscale = GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3 + work_unscale = GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3 + VBF%do_bflx_salt = .false. VBF%do_bflx_salt_dz = .false. VBF%do_bflx_temp = .false. VBF%do_bflx_temp_dz = .false. - VBF%id_Bdif = register_diag_field('ocean_model',"Bflx_dia_diff", diag%axesTi, & + if (.not.use_EOS) then + ! None of the diagnostics currently offered here make sense without temperature and salinity + ! as state variables and an equation of state. + Use_KdWork_diag = .false. + return + endif + + VBF%id_Bdif = register_diag_field('ocean_model', "Bflx_dia_diff", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz = register_diag_field('ocean_model',"Bflx_dia_diff_dz", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz = register_diag_field('ocean_model', "Bflx_dia_diff_dz", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz = register_diag_field('ocean_model',"Bflx_dia_diff_idz", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz = register_diag_field('ocean_model', "Bflx_dia_diff_idz", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV = register_scalar_field('ocean_model',"Bflx_dia_diff_idV", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV = register_scalar_field('ocean_model', "Bflx_dia_diff_idV", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_salt = register_diag_field('ocean_model',"Bflx_salt_dia_diff", diag%axesTi, & + VBF%id_Bdif_salt = register_diag_field('ocean_model', "Bflx_salt_dia_diff", diag%axesTi, & Time, "Salinity contribution to diffusive diapycnal buoyancy flux across interfaces", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_salt_dz = register_diag_field('ocean_model',"Bflx_salt_dia_diff_dz", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_salt_dz = register_diag_field('ocean_model', "Bflx_salt_dia_diff_dz", diag%axesTl, & Time, "Salinity contribution to layer integral of diffusive diapycnal buoyancy flux.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_salt_idz = register_diag_field('ocean_model',"Bflx_salt_dia_diff_idz", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_salt_idz = register_diag_field('ocean_model', "Bflx_salt_dia_diff_idz", diag%axesT1, & Time, "Salinity contribution to layer integrated diffusive diapycnal buoyancy flux.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_salt_idV = register_scalar_field('ocean_model',"Bflx_salt_dia_diff_idV", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_salt_idV = register_scalar_field('ocean_model', "Bflx_salt_dia_diff_idV", Time, diag, & "Salinity contribution to global integrated diffusive diapycnal buoyancy flux.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_temp = register_diag_field('ocean_model',"Bflx_temp_dia_diff", diag%axesTi, & + VBF%id_Bdif_temp = register_diag_field('ocean_model', "Bflx_temp_dia_diff", diag%axesTi, & Time, "Temperature contribution to diffusive diapycnal buoyancy flux across interfaces", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_temp_dz = register_diag_field('ocean_model',"Bflx_temp_dia_diff_dz", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_temp_dz = register_diag_field('ocean_model', "Bflx_temp_dia_diff_dz", diag%axesTl, & Time, "Temperature contribution to layer integral of diffusive diapycnal buoyancy flux.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_temp_idz = register_diag_field('ocean_model',"Bflx_temp_dia_diff_idz", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_temp_idz = register_diag_field('ocean_model', "Bflx_temp_dia_diff_idz", diag%axesT1, & Time, "Temperature contribution to layer integrated diffusive diapycnal buoyancy flux.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_temp_idV = register_scalar_field('ocean_model',"Bflx_temp_dia_diff_idV", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_temp_idV = register_scalar_field('ocean_model', "Bflx_temp_dia_diff_idV", Time, diag, & "Temperature contribution to global integrated diffusive diapycnal buoyancy flux.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_BBL", diag%axesTi, & + VBF%id_Bdif_BBL = register_diag_field('ocean_model', "Bflx_dia_diff_BBL", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to the BBL parameterization.", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_dz_BBL", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_BBL = register_diag_field('ocean_model', "Bflx_dia_diff_dz_BBL", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to the BBL parameterization.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_idz_BBL", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_BBL = register_diag_field('ocean_model', "Bflx_dia_diff_idz_BBL", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to the BBL parameterization.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_BBL = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_BBL", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_BBL = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_BBL", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to BBL.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_ePBL", diag%axesTi, & + VBF%id_Bdif_ePBL = register_diag_field('ocean_model', "Bflx_dia_diff_ePBL", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to ePBL", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ePBL", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_ePBL = register_diag_field('ocean_model', "Bflx_dia_diff_dz_ePBL", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to ePBL.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ePBL", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_ePBL = register_diag_field('ocean_model', "Bflx_dia_diff_idz_ePBL", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to ePBL.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_ePBL = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ePBL", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_ePBL = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_ePBL", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to ePBL.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_KS = register_diag_field('ocean_model',"Bflx_dia_diff_KS", diag%axesTi, & + VBF%id_Bdif_KS = register_diag_field('ocean_model', "Bflx_dia_diff_KS", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kappa Shear", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_KS = register_diag_field('ocean_model',"Bflx_dia_diff_dz_KS", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_KS = register_diag_field('ocean_model', "Bflx_dia_diff_dz_KS", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to Kappa Shear.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_KS = register_diag_field('ocean_model',"Bflx_dia_diff_idz_KS", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_KS = register_diag_field('ocean_model', "Bflx_dia_diff_idz_KS", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to Kappa Shear.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_KS = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_KS", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_KS = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_KS", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kappa Shear.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_bkgnd", diag%axesTi, & + VBF%id_Bdif_bkgnd = register_diag_field('ocean_model', "Bflx_dia_diff_bkgnd", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to bkgnd mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_dz_bkgnd", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_bkgnd = register_diag_field('ocean_model', "Bflx_dia_diff_dz_bkgnd", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_idz_bkgnd", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_bkgnd = register_diag_field('ocean_model', "Bflx_dia_diff_idz_bkgnd", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_bkgnd = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_bkgnd", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_bkgnd = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_bkgnd", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_bkgnd.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_ddiff_heat", diag%axesTi, & + VBF%id_Bdif_ddiff_temp = register_diag_field('ocean_model', "Bflx_dia_diff_ddiff_heat", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to double diffusion of heat", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ddiff_heat", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_ddiff_temp = register_diag_field('ocean_model', "Bflx_dia_diff_dz_ddiff_heat", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to double diffusion of heat.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ddiff_heat", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_ddiff_temp = register_diag_field('ocean_model', "Bflx_dia_diff_idz_ddiff_heat", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to double diffusion of heat.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_ddiff_temp = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ddiff_heat", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_ddiff_temp = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_ddiff_heat", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to double diffusion of heat.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_ddiff_salt", diag%axesTi, & + VBF%id_Bdif_ddiff_salt = register_diag_field('ocean_model', "Bflx_dia_diff_ddiff_salt", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to double diffusion of salt", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ddiff_salt", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_ddiff_salt = register_diag_field('ocean_model', "Bflx_dia_diff_dz_ddiff_salt", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to double diffusion of salt.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ddiff_salt", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_ddiff_salt = register_diag_field('ocean_model', "Bflx_dia_diff_idz_ddiff_salt", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to double diffusion of salt.", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_ddiff_salt = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ddiff_salt", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_ddiff_salt = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_ddiff_salt", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to double diffusion of salt.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_leak = register_diag_field('ocean_model',"Bflx_dia_diff_leak", diag%axesTi, & + VBF%id_Bdif_leak = register_diag_field('ocean_model', "Bflx_dia_diff_leak", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_leak mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_leak = register_diag_field('ocean_model',"Bflx_dia_diff_dz_leak", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_leak = register_diag_field('ocean_model', "Bflx_dia_diff_dz_leak", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_leak = register_diag_field('ocean_model',"Bflx_dia_diff_idz_leak", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_leak = register_diag_field('ocean_model', "Bflx_dia_diff_idz_leak", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_leak = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_leak", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_leak = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_leak", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_leak.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_quad = register_diag_field('ocean_model',"Bflx_dia_diff_quad", diag%axesTi, & + VBF%id_Bdif_quad = register_diag_field('ocean_model', "Bflx_dia_diff_quad", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_quad mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_quad = register_diag_field('ocean_model',"Bflx_dia_diff_dz_quad", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_quad = register_diag_field('ocean_model', "Bflx_dia_diff_dz_quad", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_quad = register_diag_field('ocean_model',"Bflx_dia_diff_idz_quad", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_quad = register_diag_field('ocean_model', "Bflx_dia_diff_idz_quad", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_quad = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_quad", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_quad = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_quad", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_quad.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_itidal", diag%axesTi, & + VBF%id_Bdif_itidal = register_diag_field('ocean_model', "Bflx_dia_diff_itidal", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_itidal mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_dz_itidal", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_itidal = register_diag_field('ocean_model', "Bflx_dia_diff_dz_itidal", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_idz_itidal", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_itidal = register_diag_field('ocean_model', "Bflx_dia_diff_idz_itidal", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_itidal = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_itidal", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_itidal = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_itidal", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_itidal.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_Froude", diag%axesTi, & + VBF%id_Bdif_Froude = register_diag_field('ocean_model', "Bflx_dia_diff_Froude", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_Froude mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_dz_Froude", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_Froude = register_diag_field('ocean_model', "Bflx_dia_diff_dz_Froude", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_idz_Froude", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_Froude = register_diag_field('ocean_model', "Bflx_dia_diff_idz_Froude", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_Froude = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_Froude", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_Froude = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_Froude", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_Froude.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_slope = register_diag_field('ocean_model',"Bflx_dia_diff_slope", diag%axesTi, & + VBF%id_Bdif_slope = register_diag_field('ocean_model', "Bflx_dia_diff_slope", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_slope mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_slope = register_diag_field('ocean_model',"Bflx_dia_diff_dz_slope", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_slope = register_diag_field('ocean_model', "Bflx_dia_diff_dz_slope", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_slope = register_diag_field('ocean_model',"Bflx_dia_diff_idz_slope", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_slope = register_diag_field('ocean_model', "Bflx_dia_diff_idz_slope", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_slope = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_slope", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_slope = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_slope", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_slope.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_lowmode", diag%axesTi, & + VBF%id_Bdif_lowmode = register_diag_field('ocean_model', "Bflx_dia_diff_lowmode", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_lowmode mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_dz_lowmode", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_lowmode = register_diag_field('ocean_model', "Bflx_dia_diff_dz_lowmode", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_idz_lowmode", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_lowmode = register_diag_field('ocean_model', "Bflx_dia_diff_idz_lowmode", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_lowmode = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_lowmode", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_lowmode = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_lowmode", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_lowmode.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_Niku", diag%axesTi, & + VBF%id_Bdif_Niku = register_diag_field('ocean_model', "Bflx_dia_diff_Niku", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_Niku mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_dz_Niku", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_Niku = register_diag_field('ocean_model', "Bflx_dia_diff_dz_Niku", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_idz_Niku", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_Niku = register_diag_field('ocean_model', "Bflx_dia_diff_idz_Niku", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_Niku = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_Niku", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_Niku = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_Niku", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_Niku.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) - VBF%id_Bdif_itides = register_diag_field('ocean_model',"Bflx_dia_diff_itides", diag%axesTi, & + VBF%id_Bdif_itides = register_diag_field('ocean_model', "Bflx_dia_diff_itides", diag%axesTi, & Time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_itides mixing", & - "W m-3", conversion=GV%H_to_kg_m2*US%Z_to_m*US%s_to_T**3) - VBF%id_Bdif_dz_itides = register_diag_field('ocean_model',"Bflx_dia_diff_dz_itides", diag%axesTl, & + units="W m-3", conversion=bflx_unscale) + VBF%id_Bdif_dz_itides = register_diag_field('ocean_model', "Bflx_dia_diff_dz_itides", diag%axesTl, & Time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idz_itides = register_diag_field('ocean_model',"Bflx_dia_diff_idz_itides", diag%axesT1, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idz_itides = register_diag_field('ocean_model', "Bflx_dia_diff_idz_itides", diag%axesT1, & Time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", & - "W m-2", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3) - VBF%id_Bdif_idV_itides = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_itides", Time, diag, & + units="W m-2", conversion=work_unscale) + VBF%id_Bdif_idV_itides = register_scalar_field('ocean_model', "Bflx_dia_diff_idV_itides", Time, diag, & "Global integrated diffusive diapycnal buoyancy flux due to Kd_itides.", & - units="W", conversion=GV%H_to_kg_m2*US%Z_to_m**2*US%s_to_T**3*US%L_to_m**2) + units="W", conversion=work_unscale*US%L_to_m**2) if (VBF%id_Bdif_dz>0 .or. VBF%id_Bdif_salt_dz>0 .or. VBF%id_Bdif_dz_BBL>0 .or. & VBF%id_Bdif_dz_ePBL>0 .or. VBF%id_Bdif_dz_KS>0 .or. VBF%id_Bdif_dz_bkgnd>0 .or. & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 55179734d4..f0d348f227 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -604,9 +604,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim] KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3] Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] - Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - N2_salt, & !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2] - N2_temp !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2] + Sdif_flx ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] real, dimension(SZI_(G),SZJ_(G)) :: & U_star, & ! The friction velocity [Z T-1 ~> m s-1]. @@ -616,15 +614,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim SkinBuoyFlux, & ! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL BBL_BuoyFlux ! 2d bottom buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL - real, dimension(SZI_(G)) :: & - p_i ,& ! Pressure at the interface [R L2 T-2 ~> Pa] - T_i, & ! Temperature at the interface [C ~> degC] - S_i, & ! Salinity at the interface [S ~> ppt] - drhodS, & ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1] - drhodT, & ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1] - dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] - dSpV_dS ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1] - logical, dimension(SZI_(G)) :: & in_boundary ! True if there are no massive layers below, where massive is defined as ! sufficiently thick that the no-flux boundary conditions have not restricted @@ -637,7 +626,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2] real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2] real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1] - real :: I_h ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1] real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is ! coupled to the bottom within a timestep [H ~> m or kg m-2] real :: Kd_add_here ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. @@ -645,12 +633,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim real :: Ent_int ! The diffusive entrainment rate at an interface [H ~> m or kg m-2] real :: Idt ! The inverse time step [T-1 ~> s-1] - real :: g_Rho0 ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2] - real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] - real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1] - logical :: nonBous ! True if not using the Boussinesq approximation - - integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state logical :: showCallTree ! If true, show the call tree integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -660,11 +642,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect h_neglect = GV%H_subroundoff - nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq) - g_Rho0 = GV%g_Earth_Z_T2 / GV%H_to_RZ - H_to_pres = GV%H_to_RZ * GV%g_Earth - alt_H_to_pres = H_to_pres * US%L_to_Z**2 * GV%Z_to_H - Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 showCallTree = callTree_showQuery() @@ -1084,9 +1061,6 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim ! Set diffusivities for VBF diagnostics if enabled if (CS%use_energetic_PBL .and. associated(CS%VBF%Kd_ePBL)) CS%VBF%Kd_ePBL(:,:,:) = Kd_ePBL(:,:,:) - if (associated(CS%VBF%Kd_temp)) CS%VBF%Kd_temp(:,:,:) = Kd_heat(:,:,:) - if (associated(CS%VBF%Kd_salt)) CS%VBF%Kd_salt(:,:,:) = Kd_salt(:,:,:) - ! Diagnose the diapycnal diffusivities and other related quantities. if (CS%id_Kd_int > 0) call post_data(CS%id_Kd_int, Kd_int, CS%diag) @@ -1123,65 +1097,8 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim endif if (CS%Use_KdWork_diag .or. CS%Use_N2_diag) then - N2_salt(:,:,:) = 0.0 - N2_temp(:,:,:) = 0.0 - !Compute N2 and don't mask negatives here - EOSdom(:) = EOS_domain(G%HI) - if (nonBous) then - !$OMP parallel do default(shared) - do j=js,je - if (associated(tv%p_surf)) then - do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo - else - do i=is,ie ; p_i(i) = 0.0 ; enddo - endif - do K=2,nz - do i=is,ie - p_i(i) = p_i(i) + H_to_pres * h(i,j,k-1) - enddo - T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k)) - S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k)) - call calculate_specific_vol_derivs(T_i, S_i, p_i, dSpV_dT, dSpV_dS, tv%eqn_of_state, EOSdom) - do i=is,ie - I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k))) - N2_salt(i,j,K) = (tv%S(i,j,k-1) - tv%S(i,j,k)) * (dSpv_dS(i) * (alt_H_to_pres * I_dzval)) - N2_temp(i,j,K) = (tv%T(i,j,k-1) - tv%T(i,j,k)) * (dSpV_dT(i) * (alt_H_to_pres * I_dzval)) - enddo - enddo - enddo - else - !$OMP parallel do default(shared) - do j=js,je - if (associated(tv%p_surf)) then - do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo - else - do i=is,ie ; p_i(i) = 0.0 ; enddo - endif - do K=2,nz - do i=is,ie - p_i(i) = p_i(i) + H_to_pres* h(i,j,k-1) - enddo - T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k)) - S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k)) - call calculate_density_derivs(T_i, S_i, p_i, dRhodT, dRhodS, tv%eqn_of_state, EOSdom) - do i=is,ie - I_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - N2_salt(i,j,K) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (dRhodS(i) * (g_rho0 * I_h)) - N2_temp(i,j,K) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (dRhodT(i) * (g_rho0 * I_h)) - enddo - enddo - enddo - endif - if (CS%id_N2_dd>0) call post_data(CS%id_N2_dd, N2_salt(:,:,:)+N2_temp(:,:,:), CS%diag) - if (CS%id_N2_salt_dd>0) call post_data(CS%id_N2_salt_dd, N2_salt, CS%diag) - if (CS%id_N2_temp_dd>0) call post_data(CS%id_N2_temp_dd, N2_temp, CS%diag) - - if (CS%Use_KdWork_diag) then - call KdWork_diagnostics(G,GV,US,CS%diag,CS%VBF,N2_salt,N2_temp,dz) - endif - - call deallocate_VBF_CS(CS%VBF) - + ! Diagnose contributions to stratification and the work done by diapycnal mixing. + call diagnose_strat_Kd_work(tv, h, dz, Kd_heat, Kd_salt, dt, G, GV, US, CS) endif ! mixing of passive tracers from massless boundary layers to interior @@ -1339,9 +1256,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, KPP_NLTscalar, & ! KPP non-local transport for scalars [nondim] KPP_buoy_flux, & ! KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3] Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] - Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] - N2_salt, & !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2] - N2_temp !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2] + Sdif_flx ! diffusive diapycnal salt flux across interfaces [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1] real, dimension(SZI_(G),SZJ_(G)) :: & U_star, & ! The friction velocity [Z T-1 ~> m s-1]. @@ -1356,14 +1271,6 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, ! sufficiently thick that the no-flux boundary conditions have not restricted ! the entrainment - usually sqrt(Kd*dt). - real, dimension(SZI_(G)) :: & - p_i ,& ! Pressure at the interface [R L2 T-2 ~> Pa] - T_i, & ! Temperature at the interface [C ~> degC] - S_i, & ! Salinity at the interface [S ~> ppt] - drhodS, & ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1] - drhodT, & ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1] - dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] - dSpV_dS ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1] real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2] @@ -1372,18 +1279,11 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, real :: dz_neglect2 ! dz_neglect^2 [Z2 ~> m2] real :: add_ent ! Entrainment that needs to be added when mixing tracers [H ~> m or kg m-2] real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1] - real :: I_h ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1] real :: Tr_ea_BBL ! The diffusive tracer thickness in the BBL that is ! coupled to the bottom within a timestep [H ~> m or kg m-2] real :: htot(SZIB_(G)) ! The summed thickness from the bottom [H ~> m or kg m-2]. real :: Kd_add_here ! An added diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. real :: Idt ! The inverse time step [T-1 ~> s-1] - real :: g_Rho0 ! G_Earth/Rho0 [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2] - real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] - real :: alt_H_to_pres! A conversion factor from thicknesses to pressure w/ alternative scaling [R Z T-2 ~> Pa m-1] - logical :: nonBous ! True if not using the Boussinesq approximation - - integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state logical :: showCallTree ! If true, show the call tree integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, Isq, Ieq, Jsq, Jeq, nz @@ -1394,11 +1294,6 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, dz_neglect = GV%dZ_subroundoff ; dz_neglect2 = dz_neglect*dz_neglect h_neglect = GV%H_subroundoff - nonBous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq) - g_Rho0 = GV%g_Earth_Z_T2 / GV%H_to_RZ - H_to_pres = GV%H_to_RZ * GV%g_Earth - alt_H_to_pres = H_to_pres * US%L_to_Z**2 * GV%Z_to_H - Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 ent_s(:,:,:) = 0.0 ; ent_t(:,:,:) = 0.0 @@ -1751,8 +1646,6 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, ! Set diffusivities for VBF diagnostics if enabled if (CS%use_energetic_PBL .and. associated(CS%VBF%Kd_ePBL)) CS%VBF%Kd_ePBL(:,:,:) = Kd_ePBL(:,:,:) - if (associated(CS%VBF%Kd_salt)) CS%VBF%Kd_salt(:,:,:) = Kd_salt(:,:,:) - if (associated(CS%VBF%Kd_temp)) CS%VBF%Kd_temp(:,:,:) = Kd_heat(:,:,:) ! Diagnose the diapycnal diffusivities and other related quantities. if (CS%id_Kd_heat > 0) call post_data(CS%id_Kd_heat, Kd_heat, CS%diag) @@ -1796,65 +1689,8 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, endif if (CS%Use_KdWork_diag .or. CS%Use_N2_diag) then - N2_salt(:,:,:) = 0.0 - N2_temp(:,:,:) = 0.0 - !Compute N2 and don't mask negatives here - EOSdom(:) = EOS_domain(G%HI) - if (nonBous) then - !$OMP parallel do default(shared) - do j=js,je - if (associated(tv%p_surf)) then - do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo - else - do i=is,ie ; p_i(i) = 0.0 ; enddo - endif - do K=2,nz - do i=is,ie - p_i(i) = p_i(i) + H_to_pres * h(i,j,k-1) - enddo - T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k)) - S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k)) - call calculate_specific_vol_derivs(T_i, S_i, p_i, dSpV_dT, dSpV_dS, tv%eqn_of_state, EOSdom) - do i=is,ie - I_dzval = 1.0 / (dz_neglect + 0.5*(dz(i,j,k-1) + dz(i,j,k))) - N2_salt(i,j,K) = (tv%S(i,j,k-1) - tv%S(i,j,k)) * (dSpv_dS(i) * (alt_H_to_pres * I_dzval)) - N2_temp(i,j,K) = (tv%T(i,j,k-1) - tv%T(i,j,k)) * (dSpV_dT(i) * (alt_H_to_pres * I_dzval)) - enddo - enddo - enddo - else - !$OMP parallel do default(shared) - do j=js,je - if (associated(tv%p_surf)) then - do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo - else - do i=is,ie ; p_i(i) = 0.0 ; enddo - endif - do K=2,nz - do i=is,ie - p_i(i) = p_i(i) + H_to_pres* h(i,j,k-1) - enddo - T_i = 0.5*(tv%T(:,j,k-1)+tv%T(:,j,k)) - S_i = 0.5*(tv%S(:,j,k-1)+tv%S(:,j,k)) - call calculate_density_derivs(T_i, S_i, p_i, dRhodT, dRhodS, tv%eqn_of_state, EOSdom) - do i=is,ie - I_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - N2_salt(i,j,K) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (dRhodS(i) * (g_rho0 * I_h)) - N2_temp(i,j,K) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (dRhodT(i) * (g_rho0 * I_h)) - enddo - enddo - enddo - endif - if (CS%id_N2_dd>0) call post_data(CS%id_N2_dd, N2_salt(:,:,:)+N2_temp(:,:,:), CS%diag) - if (CS%id_N2_salt_dd>0) call post_data(CS%id_N2_salt_dd, N2_salt, CS%diag) - if (CS%id_N2_temp_dd>0) call post_data(CS%id_N2_temp_dd, N2_temp, CS%diag) - - if (CS%Use_KdWork_diag) then - call KdWork_diagnostics(G,GV,US,CS%diag,CS%VBF,N2_salt,N2_temp,dz) - endif - - call deallocate_VBF_CS(CS%VBF) - + ! Diagnose contributions to stratification and the work done by diapycnal mixing. + call diagnose_strat_Kd_work(tv, h, dz, Kd_heat, Kd_salt, dt, G, GV, US, CS) endif ! mixing of passive tracers from massless boundary layers to interior @@ -2071,6 +1907,8 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep call enable_averages(dt, Time_end, CS%diag) + if (CS%Use_KdWork_diag) call Allocate_VBF_CS(G, GV, CS%VBF) + if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then halo = CS%halo_TS_diff !$OMP parallel do default(shared) @@ -2319,17 +2157,6 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e call cpu_clock_end(id_clock_differential_diff) if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G, GV, US) - - ! increment heat and salt diffusivity. - ! CS%useKPP==.true. already has extra_T and extra_S included - if (.not. CS%useKPP) then - !$OMP parallel do default(shared) - do K=2,nz ; do j=js,je ; do i=is,ie - Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_extra_T(i,j,K) - Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_extra_S(i,j,K) - enddo ; enddo ; enddo - endif - endif ! Calculate layer entrainments and detrainments from diffusivities and differences between @@ -2642,15 +2469,40 @@ subroutine layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_e enddo ; enddo ; enddo endif - ! mixing of passive tracers from massless boundary layers to interior - call cpu_clock_begin(id_clock_tracers) + if ((.not. CS%useKPP) .and. associated(tv%T)) then + if ((CS%id_Kd_heat > 0) .or. (CS%id_Kd_salt > 0) .or. CS%Use_KdWork_diag) then + ! In layered mode without KPP, Kd_heat and Kd_salt are only used for diagnostics. + ! CS%useKPP==.true. already has set Kd_heat and Kd_salt, including extra_T and extra_S. + if (CS%double_diffuse) then + !$OMP parallel do default(shared) + do K=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_heat(i,j,K) = Kd_int(i,j,K) + Kd_extra_T(i,j,K) + Kd_salt(i,j,K) = Kd_int(i,j,K) + Kd_extra_S(i,j,K) + enddo ; enddo ; enddo + else + !$OMP parallel do default(shared) + do K=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_salt(i,j,K) = Kd_int(i,j,K) + Kd_heat(i,j,K) = Kd_int(i,j,K) + enddo ; enddo ; enddo + endif + endif + endif ! Find the vertical distances across layers. - if (CS%mix_boundary_tracers .or. CS%double_diffuse) & + if (CS%mix_boundary_tracers .or. CS%double_diffuse .or. CS%Use_KdWork_diag) & call thickness_to_dz(h, tv, dz, G, GV, US) if (CS%double_diffuse) & call thickness_to_dz(hold, tv, dz_old, G, GV, US) + if (CS%Use_KdWork_diag .or. CS%Use_N2_diag) then + ! Diagnose contributions to stratification and the work done by diapycnal mixing. + call diagnose_strat_Kd_work(tv, h, dz, Kd_heat, Kd_salt, dt, G, GV, US, CS) + endif + + call cpu_clock_begin(id_clock_tracers) + + ! mixing of passive tracers from massless boundary layers to interior if (CS%mix_boundary_tracers) then Tr_ea_BBL = sqrt(dt * CS%Kd_BBL_tr) !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) @@ -3218,6 +3070,95 @@ subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, US, CS) end subroutine diagnose_frazil_tendency +!> This routine diagnoses contributions to stratification and the work done by diapyncal mixing. +subroutine diagnose_strat_Kd_work(tv, h, dz, Kd_heat, Kd_salt, dt, G, GV, US, CS) + type(ocean_grid_type), intent(in) :: G !< ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(thermo_var_ptrs), intent(in) :: tv !< points to updated thermodynamic fields + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< thickness [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: Kd_heat !< diapycnal diffusivity of heat + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: Kd_salt !< diapycnal diffusivity of salt and passive + !! tracers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real, intent(in) :: dt !< time step [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(diabatic_CS), pointer :: CS !< module control structure + + ! Local variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & + N2_salt, & !< Salinity contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2] + N2_temp, & !< Temperature contribution to squared buoyancy frequency at interfaces [T-2 ~> s-2] + N2_tot !< Squared buoyancy frequency at interfaces [T-2 ~> s-2] + real, dimension(SZI_(G)) :: & + p_i ,& ! Pressure at the interface [R L2 T-2 ~> Pa] + T_i, & ! Temperature at the interface [C ~> degC] + S_i, & ! Salinity at the interface [S ~> ppt] + drhodS, & ! Local change in density w.r.t. salinity using model EOS & state [R C-1 ~> kg m-3 ppt-1] + drhodT, & ! Local change in density w.r.t. temperature using model EOS & state [R C-1 ~> kg m-3 degC-1] + dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1] + dSpV_dS ! Partial derivative of specific volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1] + + real :: I_dzval ! The inverse of the thicknesses averaged to interfaces [Z-1 ~> m-1] + real :: I_h ! The inverse of the thicknesses averaged to interfaces [H-1 ~> m-1 or m2 kg-1] + real :: g_Rho ! G_Earth/Rho0 when Boussinesq, or just G_Earth times dimensional rescaling + ! factors in non-Boussinesq mode [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2] + real :: H_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected [H ~> m or kg m-2] + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state + integer :: i, j, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + h_neglect = GV%H_subroundoff + + g_Rho = GV%g_Earth_Z_T2 / GV%H_to_RZ + H_to_pres = GV%H_to_RZ * GV%g_Earth + + N2_salt(:,:,:) = 0.0 + N2_temp(:,:,:) = 0.0 + + ! Compute N2 and don't mask negatives here + EOSdom(:) = EOS_domain(G%HI) + + !$OMP parallel do default(shared) + do j=js,je + if (associated(tv%p_surf)) then + do i=is,ie ; p_i(i) = tv%p_surf(i,j) ; enddo + else + do i=is,ie ; p_i(i) = 0.0 ; enddo + endif + do K=2,nz + do i=is,ie + p_i(i) = p_i(i) + H_to_pres * h(i,j,k-1) + T_i(i) = 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k)) + S_i(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k)) + enddo + call calculate_density_derivs(T_i, S_i, p_i, dRhodT, dRhodS, tv%eqn_of_state, EOSdom) + do i=is,ie + I_h = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) + N2_salt(i,j,K) = -(tv%S(i,j,k-1) - tv%S(i,j,k)) * (dRhodS(i) * (g_rho * I_h)) + N2_temp(i,j,K) = -(tv%T(i,j,k-1) - tv%T(i,j,k)) * (dRhodT(i) * (g_rho * I_h)) + enddo + enddo + enddo + + if (CS%id_N2_dd>0) then + do K=1,nz+1 ; do j=js,je ; do i=is,ie + N2_tot(i,j,K) = N2_salt(i,j,K) + N2_temp(i,j,K) + enddo ; enddo ; enddo + call post_data(CS%id_N2_dd, N2_tot, CS%diag) + endif + if (CS%id_N2_salt_dd>0) call post_data(CS%id_N2_salt_dd, N2_salt, CS%diag) + if (CS%id_N2_temp_dd>0) call post_data(CS%id_N2_temp_dd, N2_temp, CS%diag) + + if (CS%Use_KdWork_diag) then + call KdWork_diagnostics(G, GV, US, CS%diag, CS%VBF, Kd_salt, Kd_heat, N2_salt, N2_temp, dz) + endif + + call deallocate_VBF_CS(CS%VBF) + +end subroutine diagnose_strat_Kd_work !> A simplified version of diabatic_driver_init that will allow !! tracer column functions to be called without allowing any @@ -3302,7 +3243,9 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! Local variables real :: Kd ! A diffusivity used in the default for other tracer diffusivities [Z2 T-1 ~> m2 s-1] - logical :: use_temperature, use_MARBL_tracers + logical :: use_temperature ! If true, temperature and salinity are used as state variables. + logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. + logical :: use_MARBL_tracers character(len=20) :: EN1, EN2, EN3 ! This "include" declares and sets the variable "version". @@ -3458,7 +3401,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di call get_param(param_file, mdl, "DO_BRINE_PLUME", do_brine_plume, & "If true, enables a brine plume parameterizations (not logged here)", & - do_not_log=.true.,default=.false.) + do_not_log=.true., default=.false.) if (do_brine_plume) then call get_param(param_file, mdl, "BRINE_PLUME_MLD_DEF", brine_plume_mld_def, & "A string that determines which mixed/mixing depth is used in setting "//& @@ -3617,11 +3560,13 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di endif endif - call KdWork_init(Time, G,GV,US,diag,CS%VBF,CS%Use_KdWork_diag) - if (CS%Use_KdWork_diag.and.(.not.useALEalgorithm)) & - call MOM_error(WARNING,"The KdWork diagnostics are not fully implemented for use in layer mode.") - if (CS%Use_KdWork_diag.and.(CS%use_legacy_diabatic)) & - call MOM_error(WARNING,"The KdWork diagnostics are only approximate with the legacy diabatic driver.") + call get_param(param_file, mdl, "USE_EOS", use_EOS, & + default=use_temperature, do_not_log=.true.) + call KdWork_init(Time, G, GV, US, use_EOS, diag, CS%VBF, CS%Use_KdWork_diag) + if (CS%Use_KdWork_diag .and. (.not.useALEalgorithm)) & + call MOM_error(WARNING, "The KdWork diagnostics are only approximate in layer mode.") + if (CS%Use_KdWork_diag .and. (CS%use_legacy_diabatic)) & + call MOM_error(WARNING, "The KdWork diagnostics are only approximate with the legacy diabatic driver.") call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", CS%MLDdensityDifference, & "The density difference used to determine a diagnostic mixed "//& @@ -3658,16 +3603,16 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) endif - CS%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, Time, & - 'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s, & - cmor_field_name='difvho', & - cmor_standard_name='ocean_vertical_heat_diffusivity', & - cmor_long_name='Ocean vertical heat diffusivity') - CS%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, Time, & - 'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s, & - cmor_field_name='difvso', & - cmor_standard_name='ocean_vertical_salt_diffusivity', & - cmor_long_name='Ocean vertical salt diffusivity') + if (use_temperature) then + CS%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, Time, & + 'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s, & + cmor_field_name='difvho', cmor_standard_name='ocean_vertical_heat_diffusivity', & + cmor_long_name='Ocean vertical heat diffusivity') + CS%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, Time, & + 'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1', conversion=GV%HZ_T_to_m2_s, & + cmor_field_name='difvso', cmor_standard_name='ocean_vertical_salt_diffusivity', & + cmor_long_name='Ocean vertical salt diffusivity') + endif ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise. ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive