MOM_diagnose_KdWork.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Provides diagnostics of work due to a given diffusivity
7
8use mom_diag_mediator, only : diag_ctrl, time_type, post_data, register_diag_field
10use mom_error_handler, only : mom_error, fatal, warning
11use mom_grid, only : ocean_grid_type
16
17implicit none ; private
18
19#include <MOM_memory.h>
20
21public vbf_cs
23public allocate_vbf_cs
25public kdwork_init
26public kdwork_end
27
28!> This structure has memory for used in calculating diagnostics of diffusivity
29!! many of the diffusivity diagnostics are copies of other 3d arrays. It could
30!! be written more efficiently, but it is less intrusive to copy into this structure
31!! and do all calculations in this module. These diagnostics may be expensive for
32!! routine use.
33type vbf_cs
34 ! 3d varying Kd contributions
35 real, pointer, dimension(:,:,:) :: &
36 bflx_salt => null(), & !< Salinity contribution to buoyancy flux at interfaces
37 !! [H Z T-3 ~> m2 s-3 or W m-3]
38 bflx_temp => null(), & !< Temperature contribution to buoyancy flux at interfaces
39 !! [H Z T-3 ~> m2 s-3 or W m-3]
40 bflx_salt_dz => null(), & !< Salinity contribution to integral of buoyancy flux over layer
41 !! [H Z2 T-3 ~> m3 s-3 or W m-2]
42 bflx_temp_dz => null(), & !< Temperature contribution to integral of buoyancy flux over layer
43 !! [H Z2 T-3 ~> m3 s-3 or W m-2]
44 ! The following are all allocatable arrays that store copies of process driven Kd, so that
45 ! the process driven buoyancy flux and work can be derived at the end of the time step.
46 kd_salt => null(), & !< total diapycnal diffusivity of salt at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
47 kd_temp => null(), & !< total diapycnal diffusivity of heat at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
48 kd_bbl => null(), & !< diapycnal diffusivity due to BBL at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
49 kd_epbl => null(), & !< diapycnal diffusivity due to ePBL at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
50 kd_ks => null(), & !< diapycnal diffusivity due to Kappa Shear at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
51 kd_bkgnd => null(), & !< diapycnal diffusivity due to Kd_bkgnd at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
52 kd_ddiff_s => null(), &!< diapycnal diffusivity due to double diffusion of salt at interfaces
53 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
54 kd_ddiff_t => null(), &!< diapycnal diffusivity due to double diffusion of heat at interfaces
55 !![H Z T-1 ~> m2 s-1 or kg m-1 s-1]
56 kd_leak => null(), & !< diapycnal diffusivity due to Kd_leak at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
57 kd_quad => null(), & !< diapycnal diffusivity due to Kd_quad at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
58 kd_itidal => null(), & !< diapycnal diffusivity due to Kd_itidal at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
59 kd_froude => null(), & !< diapycnal diffusivity due to Kd_Froude at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
60 kd_slope => null(), & !< diapycnal diffusivity due to Kd_slope at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
61 kd_lowmode => null(), &!< diapycnal diffusivity due to Kd_lowmode at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
62 kd_niku => null(), & !< diapycnal diffusivity due to Kd_Niku at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
63 kd_itides => null() !< diapycnal diffusivity due to Kd_itides at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
64
65 ! Constant Kd contributions
66 real :: kd_add !< spatially uniform additional diapycnal diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
67 !! a diagnostic for this diffusivity is not yet included, but this makes it straightforward to add
68
69 !>@{ Diagnostic IDs
70 integer :: id_bdif = -1, id_bdif_salt = -1, id_bdif_temp = -1
71 integer :: id_bdif_dz = -1, id_bdif_salt_dz = -1, id_bdif_temp_dz = -1
72 integer :: id_bdif_idz = -1, id_bdif_salt_idz = -1, id_bdif_temp_idz = -1
73 integer :: id_bdif_idv = -1, id_bdif_salt_idv = -1, id_bdif_temp_idv = -1
74 integer :: id_bdif_epbl = -1, id_bdif_dz_epbl = -1, id_bdif_idz_epbl = -1, id_bdif_idv_epbl = -1
75 integer :: id_bdif_bbl = -1, id_bdif_dz_bbl = -1, id_bdif_idz_bbl = -1, id_bdif_idv_bbl = -1
76 integer :: id_bdif_ks = -1, id_bdif_dz_ks = -1, id_bdif_idz_ks = -1, id_bdif_idv_ks = -1
77 integer :: id_bdif_bkgnd = -1, id_bdif_dz_bkgnd = -1, id_bdif_idz_bkgnd = -1, id_bdif_idv_bkgnd = -1
78 integer :: id_bdif_ddiff_temp = -1, id_bdif_ddiff_salt = -1
79 integer :: id_bdif_dz_ddiff_temp = -1, id_bdif_dz_ddiff_salt = -1
80 integer :: id_bdif_idz_ddiff_temp = -1, id_bdif_idz_ddiff_salt = -1
81 integer :: id_bdif_idv_ddiff_temp = -1, id_bdif_idv_ddiff_salt = -1
82 integer :: id_bdif_leak = -1, id_bdif_dz_leak = -1, id_bdif_idz_leak = -1, id_bdif_idv_leak = -1
83 integer :: id_bdif_quad = -1, id_bdif_dz_quad = -1, id_bdif_idz_quad = -1, id_bdif_idv_quad = -1
84 integer :: id_bdif_itidal = -1, id_bdif_dz_itidal = -1, id_bdif_idz_itidal = -1, id_bdif_idv_itidal = -1
85 integer :: id_bdif_froude = -1, id_bdif_dz_froude = -1, id_bdif_idz_froude = -1, id_bdif_idv_froude = -1
86 integer :: id_bdif_slope = -1, id_bdif_dz_slope = -1, id_bdif_idz_slope = -1, id_bdif_idv_slope = -1
87 integer :: id_bdif_lowmode = -1, id_bdif_dz_lowmode = -1, id_bdif_idz_lowmode = -1, id_bdif_idv_lowmode = -1
88 integer :: id_bdif_niku = -1, id_bdif_dz_niku = -1, id_bdif_idz_niku = -1, id_bdif_idv_niku = -1
89 integer :: id_bdif_itides = -1, id_bdif_dz_itides = -1, id_bdif_idz_itides = -1, id_bdif_idv_itides = -1
90 !>@}
91
92 logical :: do_bflx_salt = .false. !< Logical flag to indicate if N2_salt should be computed
93 logical :: do_bflx_temp = .false. !< Logical flag to indicate if N2_temp should be computed
94 logical :: do_bflx_salt_dz = .false. !< Logical flag to indicate if N2_salt should be computed
95 logical :: do_bflx_temp_dz = .false. !< Logical flag to indicate if N2_temp should be computed
96
97end type vbf_cs
98
99! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
100! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
101! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
102! vary with the Boussinesq approximation, the Boussinesq variant is given first.
103
104contains
105
106!> Loop over all implemented diffusivities to diagnose and output Kd Work/buoyancy fluxes
107subroutine kdwork_diagnostics(G,GV,US,diag,VBF,N2_Salt,N2_Temp,dz)
108 type(ocean_grid_type), intent(in) :: g !< Grid type
109 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
110 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
111 type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
112 type (vbf_cs), intent(inout) :: vbf !< Vertical buoyancy flux structure
113 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
114 intent(in) :: n2_salt !< Buoyancy frequency [T-2 ~> s-2]
115 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
116 intent(in) :: n2_temp !< Buoyancy frequency [T-2 ~> s-2]
117 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
118 intent(in) :: dz !< Grid spacing [Z ~> m]
119
120 ! Work arrays for computing buoyancy flux integrals
121 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: work3d_i
122 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work3d_l
123 real, dimension(SZI_(G),SZJ_(G)) :: work2d, work2d_salt, work2d_temp
124 real :: work, work_salt, work_temp
125
126 integer :: i, j, k, nz, isc, iec, jsc, jec
127
128 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
129
130 nz = gv%ke
131
132 ! Compute total fluxes
133 if (vbf%id_Bdif_dz>0 .or. vbf%id_Bdif_salt_dz>0 .or. vbf%id_Bdif_temp_dz>0 .or. &
134 vbf%id_Bdif_idz>0 .or. vbf%id_Bdif_salt_idz>0 .or. vbf%id_Bdif_temp_idz>0 .or. &
135 vbf%id_Bdif_idV>0 .or. vbf%id_Bdif_salt_idV>0 .or. vbf%id_Bdif_temp_idV>0 ) then ! Doing vertical integrals
136 ! Do Salt
137 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. &
138 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) &
139 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_salt, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
140 ! Do Temp
141 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. &
142 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) &
143 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_temp, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
144 if (vbf%id_Bdif_temp_idz>0 .or. vbf%id_Bdif_idz>0) then
145 work2d_temp(:,:) = 0.0
146 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
147 work2d_temp(i,j) = work2d_temp(i,j) + vbf%Bflx_temp_dz(i,j,k)
148 enddo ; enddo ; enddo
149 endif
150 if (vbf%id_Bdif_temp_idV>0 .or. vbf%id_Bdif_idV>0) then
151 work_temp = 0.0
152 do k = 1,nz
153 work_temp = work_temp + global_area_integral(vbf%Bflx_temp_dz(:,:,k), g, &
154 tmp_scale=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
155 enddo
156 endif
157 if (vbf%id_Bdif_salt_idz>0 .or. vbf%id_Bdif_idz>0) then
158 work2d_salt(:,:) = 0.0
159 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
160 work2d_salt(i,j) = work2d_salt(i,j) + vbf%Bflx_salt_dz(i,j,k)
161 enddo ; enddo ; enddo
162 endif
163 if (vbf%id_Bdif_salt_idV>0 .or. vbf%id_Bdif_idV>0) then
164 work_salt = 0.0
165 do k = 1,nz
166 work_salt = work_salt + global_area_integral(vbf%Bflx_salt_dz(:,:,k), g, &
167 tmp_scale=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
168 enddo
169 endif
170 work = work_temp + work_salt
171 do j = jsc,jec ; do i = isc,iec
172 work2d(i,j) = work2d_temp(i,j) + work2d_salt(i,j)
173 enddo ; enddo
174 elseif (vbf%id_Bdif>0 .or. vbf%id_Bdif_salt>0 .or. vbf%id_Bdif_temp>0) then ! Not doing vertical integrals
175 ! Do Salt
176 if (vbf%id_Bdif_salt>0 .or. vbf%id_Bdif>0) &
177 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_salt, vbf%Bflx_salt)
178 if (vbf%id_Bdif_temp>0 .or. vbf%id_Bdif>0) &
179 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_temp, vbf%Bflx_temp)
180 endif
181 ! Post total fluxes
182 if (vbf%id_Bdif_salt>0) call post_data(vbf%id_Bdif_salt, vbf%Bflx_salt, diag)
183 if (vbf%id_Bdif_temp>0) call post_data(vbf%id_Bdif_temp, vbf%Bflx_temp, diag)
184 if (vbf%id_Bdif>0) then
185 work3d_i(:,:,:) = 0.0
186 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
187 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
188 enddo ; enddo ; enddo
189 call post_data(vbf%id_Bdif, work3d_i, diag)
190 endif
191 if (vbf%id_Bdif_salt_dz>0) call post_data(vbf%id_Bdif_salt_dz, vbf%Bflx_salt_dz, diag)
192 if (vbf%id_Bdif_temp_dz>0) call post_data(vbf%id_Bdif_temp_dz, vbf%Bflx_temp_dz, diag)
193 if (vbf%id_Bdif_dz>0) then
194 work3d_l(:,:,:) = 0.0
195 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
196 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
197 enddo ; enddo ; enddo
198 call post_data(vbf%id_Bdif_dz, work3d_l, diag)
199 endif
200 if (vbf%id_Bdif_salt_idz>0) call post_data(vbf%id_Bdif_salt_idz, work2d_salt, diag)
201 if (vbf%id_Bdif_temp_idz>0) call post_data(vbf%id_Bdif_temp_idz, work2d_temp, diag)
202 if (vbf%id_Bdif_idz>0) call post_data(vbf%id_Bdif_idz, work2d, diag)
203 if (vbf%id_Bdif_salt_idV>0) call post_data(vbf%id_Bdif_salt_idV, work_salt, diag)
204 if (vbf%id_Bdif_temp_idV>0) call post_data(vbf%id_Bdif_temp_idV, work_temp, diag)
205 if (vbf%id_Bdif_idV>0) call post_data(vbf%id_Bdif_idV, work, diag)
206
207 ! Compute ePBL fluxes
208 if (vbf%id_Bdif_dz_ePBL>0.or.vbf%id_Bdif_idz_ePBL>0.or.vbf%id_Bdif_idV_ePBL>0) then
209 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_ePBL, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
210 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_ePBL, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
211 if (vbf%id_Bdif_idz_ePBL>0) then
212 work2d(:,:) = 0.0
213 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
214 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
215 enddo ; enddo ; enddo
216 endif
217 if (vbf%id_Bdif_idV_ePBL>0) then
218 work = 0.0
219 do k = 1,nz
220 work = work + &
221 (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) + &
222 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))
223 enddo
224 endif
225 elseif (vbf%id_Bdif_ePBL>0) then
226 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_ePBL, vbf%Bflx_salt)
227 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_ePBL, vbf%Bflx_temp)
228 endif
229 ! Post ePBL fluxes
230 if (vbf%id_Bdif_ePBL>0) then
231 work3d_i(:,:,:) = 0.0
232 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
233 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
234 enddo ; enddo ; enddo
235 call post_data(vbf%id_Bdif_ePBL, work3d_i, diag)
236 endif
237 if (vbf%id_Bdif_dz_ePBL>0) then
238 work3d_l(:,:,:) = 0.0
239 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
240 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
241 enddo ; enddo ; enddo
242 call post_data(vbf%id_Bdif_dz_ePBL, work3d_l, diag)
243 endif
244 if (vbf%id_Bdif_idz_ePBL>0) call post_data(vbf%id_Bdif_idz_ePBL, work2d, diag)
245 if (vbf%id_Bdif_idV_ePBL>0) call post_data(vbf%id_Bdif_idV_ePBL, work, diag)
246
247 ! Compute BBL fluxes
248 if (vbf%id_Bdif_dz_BBL>0.or.vbf%id_Bdif_idz_BBL>0.or.vbf%id_Bdif_idV_BBL>0) then
249 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_BBL, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
250 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_BBL, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
251 if (vbf%id_Bdif_idz_BBL>0) then
252 work2d(:,:) = 0.0
253 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
254 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
255 enddo ; enddo ; enddo
256 endif
257 if (vbf%id_Bdif_idV_BBL>0) then
258 work = 0.0
259 do k = 1,nz
260 work = work + &
261 (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) + &
262 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))
263 enddo
264 endif
265 elseif (vbf%id_Bdif_BBL>0) then
266 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_BBL, vbf%Bflx_salt)
267 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_BBL, vbf%Bflx_temp)
268 endif
269 ! Post BBL fluxes
270 if (vbf%id_Bdif_BBL>0) then
271 work3d_i(:,:,:) = 0.0
272 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
273 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
274 enddo ; enddo ; enddo
275 call post_data(vbf%id_Bdif_BBL, work3d_i, diag)
276 endif
277 if (vbf%id_Bdif_dz_BBL>0) then
278 work3d_l(:,:,:) = 0.0
279 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
280 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
281 enddo ; enddo ; enddo
282 call post_data(vbf%id_Bdif_dz_BBL, work3d_l, diag)
283 endif
284 if (vbf%id_Bdif_idz_BBL>0) call post_data(vbf%id_Bdif_idz_BBL, work2d, diag)
285 if (vbf%id_Bdif_idV_BBL>0) call post_data(vbf%id_Bdif_idV_BBL, work, diag)
286
287 ! Compute Kappa Shear fluxes
288 if (vbf%id_Bdif_dz_KS>0.or.vbf%id_Bdif_idz_KS>0.or.vbf%id_Bdif_idV_KS>0) then
289 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_KS, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
290 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_KS, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
291 if (vbf%id_Bdif_idz_KS>0) then
292 work2d(:,:) = 0.0
293 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
294 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
295 enddo ; enddo ; enddo
296 endif
297 if (vbf%id_Bdif_idV_KS>0) then
298 work = 0.0
299 do k = 1,nz
300 work = work + &
301 (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) + &
302 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))
303 enddo
304 endif
305 elseif (vbf%id_Bdif_KS>0) then
306 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_KS, vbf%Bflx_salt)
307 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_KS, vbf%Bflx_temp)
308 endif
309 ! Post Kappa Shear fluxes
310 if (vbf%id_Bdif_KS>0) then
311 work3d_i(:,:,:) = 0.0
312 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
313 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
314 enddo ; enddo ; enddo
315 call post_data(vbf%id_Bdif_KS, work3d_i, diag)
316 endif
317 if (vbf%id_Bdif_dz_KS>0) then
318 work3d_l(:,:,:) = 0.0
319 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
320 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
321 enddo ; enddo ; enddo
322 call post_data(vbf%id_Bdif_dz_KS, work3d_l, diag)
323 endif
324 if (vbf%id_Bdif_idz_KS>0) call post_data(vbf%id_Bdif_idz_KS, work2d, diag)
325 if (vbf%id_Bdif_idV_KS>0) call post_data(vbf%id_Bdif_idV_KS, work, diag)
326
327 ! Compute bkgnd fluxes
328 if (vbf%id_Bdif_dz_bkgnd>0.or.vbf%id_Bdif_idz_bkgnd>0.or.vbf%id_Bdif_idV_bkgnd>0) then
329 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_bkgnd, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
330 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_bkgnd, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
331 if (vbf%id_Bdif_idz_bkgnd>0) then
332 work2d(:,:) = 0.0
333 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
334 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
335 enddo ; enddo ; enddo
336 endif
337 if (vbf%id_Bdif_idV_bkgnd>0) then
338 work = 0.0
339 do k = 1,nz
340 work = work + &
341 (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) + &
342 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))
343 enddo
344 endif
345 elseif (vbf%id_Bdif_bkgnd>0) then
346 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_bkgnd, vbf%Bflx_salt)
347 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_bkgnd, vbf%Bflx_temp)
348 endif
349 ! Post bkgnd fluxes
350 if (vbf%id_Bdif_bkgnd>0) then
351 work3d_i(:,:,:) = 0.0
352 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
353 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
354 enddo ; enddo ; enddo
355 call post_data(vbf%id_Bdif_bkgnd, work3d_i, diag)
356 endif
357 if (vbf%id_Bdif_dz_bkgnd>0) then
358 work3d_l(:,:,:) = 0.0
359 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
360 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
361 enddo ; enddo ; enddo
362 call post_data(vbf%id_Bdif_dz_bkgnd, work3d_l, diag)
363 endif
364 if (vbf%id_Bdif_idz_bkgnd>0) call post_data(vbf%id_Bdif_idz_bkgnd, work2d, diag)
365 if (vbf%id_Bdif_idV_bkgnd>0) call post_data(vbf%id_Bdif_idV_bkgnd, work, diag)
366
367 ! Compute double diffusion fluxes
368 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
369 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_ddiff_T, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
370 if (vbf%id_Bdif_idz_ddiff_temp>0) then
371 work2d_temp(:,:) = 0.0
372 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
373 work2d_temp(i,j) = work2d_temp(i,j) + vbf%Bflx_temp_dz(i,j,k)
374 enddo ; enddo ; enddo
375 endif
376 if (vbf%id_Bdif_idV_ddiff_temp>0) then
377 work_temp = 0.0
378 do k = 1,nz
379 work_temp = work_temp + global_area_integral(vbf%Bflx_temp_dz(:,:,k), g, &
380 tmp_scale=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
381 enddo
382 endif
383 elseif (vbf%id_Bdif_ddiff_temp>0) then
384 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_ddiff_T, vbf%Bflx_temp)
385 endif
386 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
387 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_ddiff_S, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
388 if (vbf%id_Bdif_idz_ddiff_salt>0) then
389 work2d_salt(:,:) = 0.0
390 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
391 work2d_salt(i,j) = work2d_salt(i,j) + vbf%Bflx_salt_dz(i,j,k)
392 enddo ; enddo ; enddo
393 endif
394 if (vbf%id_Bdif_idV_ddiff_salt>0) then
395 work_salt = 0.0
396 do k = 1,nz
397 work_salt = work_salt + global_area_integral(vbf%Bflx_salt_dz(:,:,k), g, &
398 tmp_scale=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
399 enddo
400 endif
401 elseif (vbf%id_Bdif_ddiff_salt>0) then
402 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_ddiff_S, vbf%Bflx_salt)
403 endif
404 ! Post double diffusion fluxes
405 if (vbf%id_Bdif_ddiff_temp>0) call post_data(vbf%id_Bdif_ddiff_temp, vbf%Bflx_temp, diag)
406 if (vbf%id_Bdif_dz_ddiff_temp>0) call post_data(vbf%id_Bdif_dz_ddiff_temp, vbf%Bflx_temp_dz, diag)
407 if (vbf%id_Bdif_idz_ddiff_temp>0) call post_data(vbf%id_Bdif_idz_ddiff_temp, work2d_temp, diag)
408 if (vbf%id_Bdif_idV_ddiff_temp>0) call post_data(vbf%id_Bdif_idV_ddiff_temp, work_temp, diag)
409 if (vbf%id_Bdif_ddiff_salt>0) call post_data(vbf%id_Bdif_ddiff_salt, vbf%Bflx_salt, diag)
410 if (vbf%id_Bdif_dz_ddiff_salt>0) call post_data(vbf%id_Bdif_dz_ddiff_salt, vbf%Bflx_salt_dz, diag)
411 if (vbf%id_Bdif_idz_ddiff_salt>0) call post_data(vbf%id_Bdif_idz_ddiff_salt, work2d_salt, diag)
412 if (vbf%id_Bdif_idV_ddiff_salt>0) call post_data(vbf%id_Bdif_idV_ddiff_salt, work_salt, diag)
413
414 ! Compute Kd_leak fluxes
415 if (vbf%id_Bdif_dz_leak>0.or.vbf%id_Bdif_idz_leak>0.or.vbf%id_Bdif_idV_leak>0) then
416 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_leak, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
417 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_leak, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
418 if (vbf%id_Bdif_idz_leak>0) then
419 work2d(:,:) = 0.0
420 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
421 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
422 enddo ; enddo ; enddo
423 endif
424 if (vbf%id_Bdif_idV_leak>0) then
425 work = 0.0
426 do k = 1,nz
427 work = work + &
428 (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) + &
429 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))
430 enddo
431 endif
432 elseif (vbf%id_Bdif_leak>0) then
433 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_leak, vbf%Bflx_salt)
434 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_leak, vbf%Bflx_temp)
435 endif
436 ! Post Kd_leak fluxes
437 if (vbf%id_Bdif_leak>0) then
438 work3d_i(:,:,:) = 0.0
439 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
440 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
441 enddo ; enddo ; enddo
442 call post_data(vbf%id_Bdif_leak, work3d_i, diag)
443 endif
444 if (vbf%id_Bdif_dz_leak>0) then
445 work3d_l(:,:,:) = 0.0
446 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
447 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
448 enddo ; enddo ; enddo
449 call post_data(vbf%id_Bdif_dz_leak, work3d_l, diag)
450 endif
451 if (vbf%id_Bdif_idz_leak>0) call post_data(vbf%id_Bdif_idz_leak, work2d, diag)
452 if (vbf%id_Bdif_idV_leak>0) call post_data(vbf%id_Bdif_idV_leak, work, diag)
453
454 ! Compute Kd_quad fluxes
455 if (vbf%id_Bdif_dz_quad>0.or.vbf%id_Bdif_idz_quad>0.or.vbf%id_Bdif_idV_quad>0) then
456 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_quad, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
457 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_quad, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
458 if (vbf%id_Bdif_idz_quad>0) then
459 work2d(:,:) = 0.0
460 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
461 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
462 enddo ; enddo ; enddo
463 endif
464 if (vbf%id_Bdif_idV_quad>0) then
465 work = 0.0
466 do k = 1,nz
467 work = work + &
468 (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) + &
469 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))
470 enddo
471 endif
472 elseif (vbf%id_Bdif_quad>0) then
473 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_quad, vbf%Bflx_salt)
474 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_quad, vbf%Bflx_temp)
475 endif
476 ! Post Kd_quad fluxes
477 if (vbf%id_Bdif_quad>0) then
478 work3d_i(:,:,:) = 0.0
479 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
480 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
481 enddo ; enddo ; enddo
482 call post_data(vbf%id_Bdif_quad, work3d_i, diag)
483 endif
484 if (vbf%id_Bdif_dz_quad>0) then
485 work3d_l(:,:,:) = 0.0
486 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
487 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
488 enddo ; enddo ; enddo
489 call post_data(vbf%id_Bdif_dz_quad, work3d_l, diag)
490 endif
491 if (vbf%id_Bdif_idz_quad>0) call post_data(vbf%id_Bdif_idz_quad, work2d, diag)
492 if (vbf%id_Bdif_idV_quad>0) call post_data(vbf%id_Bdif_idV_quad, work, diag)
493
494 ! Compute Kd_itidal fluxes
495 if (vbf%id_Bdif_dz_itidal>0.or.vbf%id_Bdif_idz_itidal>0.or.vbf%id_Bdif_idV_itidal>0) then
496 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_itidal, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
497 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_itidal, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
498 if (vbf%id_Bdif_idz_itidal>0) then
499 work2d(:,:) = 0.0
500 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
501 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
502 enddo ; enddo ; enddo
503 endif
504 if (vbf%id_Bdif_idV_itidal>0) then
505 work = 0.0
506 do k = 1,nz
507 work = work + &
508 (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) + &
509 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))
510 enddo
511 endif
512 elseif (vbf%id_Bdif_itidal>0) then
513 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_itidal, vbf%Bflx_salt)
514 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_itidal, vbf%Bflx_temp)
515 endif
516 ! Post Kd_itidal fluxes
517 if (vbf%id_Bdif_itidal>0) then
518 work3d_i(:,:,:) = 0.0
519 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
520 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
521 enddo ; enddo ; enddo
522 call post_data(vbf%id_Bdif_itidal, work3d_i, diag)
523 endif
524 if (vbf%id_Bdif_dz_itidal>0) then
525 work3d_l(:,:,:) = 0.0
526 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
527 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k)+vbf%Bflx_salt_dz(i,j,k)
528 enddo ; enddo ; enddo
529 call post_data(vbf%id_Bdif_dz_itidal, work3d_l, diag)
530 endif
531 if (vbf%id_Bdif_idz_itidal>0) call post_data(vbf%id_Bdif_idz_itidal, work2d, diag)
532 if (vbf%id_Bdif_idV_itidal>0) call post_data(vbf%id_Bdif_idV_itidal, work, diag)
533
534 ! Compute Kd_Froude fluxes
535 if (vbf%id_Bdif_dz_Froude>0.or.vbf%id_Bdif_idz_Froude>0.or.vbf%id_Bdif_idV_Froude>0) then
536 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_Froude, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
537 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_Froude, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
538 if (vbf%id_Bdif_idz_Froude>0) then
539 work2d(:,:) = 0.0
540 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
541 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
542 enddo ; enddo ; enddo
543 endif
544 if (vbf%id_Bdif_idV_Froude>0) then
545 work = 0.0
546 do k = 1,nz
547 work = work + &
548 (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) + &
549 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))
550 enddo
551 endif
552 elseif (vbf%id_Bdif_Froude>0) then
553 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_Froude, vbf%Bflx_salt)
554 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_Froude, vbf%Bflx_temp)
555 endif
556 ! Post Kd_Froude fluxes
557 if (vbf%id_Bdif_Froude>0) then
558 work3d_i(:,:,:) = 0.0
559 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
560 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
561 enddo ; enddo ; enddo
562 call post_data(vbf%id_Bdif_Froude, work3d_i, diag)
563 endif
564 if (vbf%id_Bdif_dz_Froude>0) then
565 work3d_l(:,:,:) = 0.0
566 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
567 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
568 enddo ; enddo ; enddo
569 call post_data(vbf%id_Bdif_dz_Froude, work3d_l, diag)
570 endif
571 if (vbf%id_Bdif_idz_Froude>0) call post_data(vbf%id_Bdif_idz_Froude, work2d, diag)
572 if (vbf%id_Bdif_idV_Froude>0) call post_data(vbf%id_Bdif_idV_Froude, work, diag)
573
574 ! Compute Kd_slope fluxes
575 if (vbf%id_Bdif_dz_slope>0.or.vbf%id_Bdif_idz_slope>0.or.vbf%id_Bdif_idV_slope>0) then
576 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_slope, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
577 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_slope, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
578 if (vbf%id_Bdif_idz_slope>0) then
579 work2d(:,:) = 0.0
580 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
581 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
582 enddo ; enddo ; enddo
583 endif
584 if (vbf%id_Bdif_idV_slope>0) then
585 work = 0.0
586 do k = 1,nz
587 work = work + &
588 (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) + &
589 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))
590 enddo
591 endif
592 elseif (vbf%id_Bdif_slope>0) then
593 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_slope, vbf%Bflx_salt)
594 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_slope, vbf%Bflx_temp)
595 endif
596 ! Post Kd_slope fluxes
597 if (vbf%id_Bdif_slope>0) then
598 work3d_i(:,:,:) = 0.0
599 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
600 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
601 enddo ; enddo ; enddo
602 call post_data(vbf%id_Bdif_slope, work3d_i, diag)
603 endif
604 if (vbf%id_Bdif_dz_slope>0) then
605 work3d_l(:,:,:) = 0.0
606 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
607 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
608 enddo ; enddo ; enddo
609 call post_data(vbf%id_Bdif_dz_slope, work3d_l, diag)
610 endif
611 if (vbf%id_Bdif_idz_slope>0) call post_data(vbf%id_Bdif_idz_slope, work2d, diag)
612 if (vbf%id_Bdif_idV_slope>0) call post_data(vbf%id_Bdif_idV_slope, work, diag)
613
614 ! Compute Kd_lowmode fluxes
615 if (vbf%id_Bdif_dz_lowmode>0.or.vbf%id_Bdif_idz_lowmode>0.or.vbf%id_Bdif_idV_lowmode>0) then
616 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_lowmode, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
617 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_lowmode, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
618 if (vbf%id_Bdif_idz_lowmode>0) then
619 work2d(:,:) = 0.0
620 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
621 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
622 enddo ; enddo ; enddo
623 endif
624 if (vbf%id_Bdif_idV_lowmode>0) then
625 work = 0.0
626 do k = 1,nz
627 work = work + &
628 (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) + &
629 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))
630 enddo
631 endif
632 elseif (vbf%id_Bdif_lowmode>0) then
633 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_lowmode, vbf%Bflx_salt)
634 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_lowmode, vbf%Bflx_temp)
635 endif
636 ! Post Kd_lowmode fluxes
637 if (vbf%id_Bdif_lowmode>0) then
638 work3d_i(:,:,:) = 0.0
639 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
640 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
641 enddo ; enddo ; enddo
642 call post_data(vbf%id_Bdif_lowmode, work3d_i, diag)
643 endif
644 if (vbf%id_Bdif_dz_lowmode>0) then
645 work3d_l(:,:,:) = 0.0
646 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
647 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
648 enddo ; enddo ; enddo
649 call post_data(vbf%id_Bdif_dz_lowmode, work3d_l, diag)
650 endif
651 if (vbf%id_Bdif_idz_lowmode>0) call post_data(vbf%id_Bdif_idz_lowmode, work2d, diag)
652 if (vbf%id_Bdif_idV_lowmode>0) call post_data(vbf%id_Bdif_idV_lowmode, work, diag)
653
654 ! Compute Kd_Niku fluxes
655 if (vbf%id_Bdif_dz_Niku>0 .or. vbf%id_Bdif_idz_Niku>0 .or. vbf%id_Bdif_idV_Niku>0) then
656 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_Niku, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
657 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_Niku, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
658 if (vbf%id_Bdif_idz_Niku>0) then
659 work2d(:,:) = 0.0
660 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
661 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
662 enddo ; enddo ; enddo
663 endif
664 if (vbf%id_Bdif_idV_Niku>0) then
665 work = 0.0
666 do k = 1,nz
667 work = work + &
668 (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) + &
669 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))
670 enddo
671 endif
672 elseif (vbf%id_Bdif_Niku>0) then
673 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_Niku, vbf%Bflx_salt)
674 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_Niku, vbf%Bflx_temp)
675 endif
676 ! Post Kd_Niku fluxes
677 if (vbf%id_Bdif_Niku>0) then
678 work3d_i(:,:,:) = 0.0
679 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
680 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
681 enddo ; enddo ; enddo
682 call post_data(vbf%id_Bdif_lowmode, work3d_i, diag)
683 endif
684 if (vbf%id_Bdif_dz_Niku>0) then
685 work3d_l(:,:,:) = 0.0
686 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
687 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
688 enddo ; enddo ; enddo
689 call post_data(vbf%id_Bdif_dz_Niku, work3d_l, diag)
690 endif
691 if (vbf%id_Bdif_idz_Niku>0) call post_data(vbf%id_Bdif_idz_Niku, work2d, diag)
692 if (vbf%id_Bdif_idV_Niku>0) call post_data(vbf%id_Bdif_idV_Niku, work, diag)
693
694 ! Compute Kd_itides fluxes
695 if (vbf%id_Bdif_dz_itides>0 .or. vbf%id_Bdif_idz_itides>0 .or. vbf%id_Bdif_idV_itides>0) then
696 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_itides, vbf%Bflx_salt, dz=dz, bdif_flx_dz=vbf%Bflx_salt_dz)
697 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_itides, vbf%Bflx_temp, dz=dz, bdif_flx_dz=vbf%Bflx_temp_dz)
698 if (vbf%id_Bdif_idz_itides>0) then
699 work2d(:,:) = 0.0
700 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
701 work2d(i,j) = work2d(i,j) + (vbf%Bflx_salt_dz(i,j,k) + vbf%Bflx_temp_dz(i,j,k))
702 enddo ; enddo ; enddo
703 endif
704 if (vbf%id_Bdif_idV_itides>0) then
705 work = 0.0
706 do k = 1,nz
707 work = work + &
708 (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) + &
709 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))
710 enddo
711 endif
712 elseif (vbf%id_Bdif_itides>0) then
713 call diagnosekdwork(g, gv, n2_salt, vbf%Kd_itides, vbf%Bflx_salt)
714 call diagnosekdwork(g, gv, n2_temp, vbf%Kd_itides, vbf%Bflx_temp)
715 endif
716 ! Post Kd_itides fluxes
717 if (vbf%id_Bdif_itides>0) then
718 work3d_i(:,:,:) = 0.0
719 do k = 1,nz+1 ; do j = jsc,jec ; do i = isc,iec
720 work3d_i(i,j,k) = vbf%Bflx_temp(i,j,k) + vbf%Bflx_salt(i,j,k)
721 enddo ; enddo ; enddo
722 call post_data(vbf%id_Bdif_itides, work3d_i, diag)
723 endif
724 if (vbf%id_Bdif_dz_itides>0) then
725 work3d_l(:,:,:) = 0.0
726 do k = 1,nz ; do j = jsc,jec ; do i = isc,iec
727 work3d_l(i,j,k) = vbf%Bflx_temp_dz(i,j,k) + vbf%Bflx_salt_dz(i,j,k)
728 enddo ; enddo ; enddo
729 call post_data(vbf%id_Bdif_dz_itides, work3d_l, diag)
730 endif
731 if (vbf%id_Bdif_idz_itides>0) call post_data(vbf%id_Bdif_idz_itides, work2d, diag)
732 if (vbf%id_Bdif_idV_itides>0) call post_data(vbf%id_Bdif_idV_itides, work, diag)
733
734end subroutine kdwork_diagnostics
735
736!> Diagnose the implied "work", or buoyancy forcing & its integral, due to a given diffusivity and column state.
737subroutine diagnosekdwork(G, GV, N2, Kd, Bdif_flx, dz, Bdif_flx_dz)
738 type(ocean_grid_type), intent(in) :: G !< Grid type
739 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
740 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
741 intent(in) :: N2 !< Buoyancy frequency [T-2 ~> s-2]
742 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
743 intent(in) :: Kd !< Diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
744 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
745 intent(out) :: Bdif_flx !< Buoyancy flux [H Z T-3 ~> m2 s-3 or W m-3]
746 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
747 intent(in), optional :: dz !< Grid spacing [Z ~> m]
748 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
749 intent(out), optional :: Bdif_flx_dz !< Buoyancy flux over layer [H Z2 T-3 ~> m3 s-3 or W m-2]
750
751 integer :: i, j, k
752
753 bdif_flx(:,:,1) = 0.0
754 bdif_flx(:,:,gv%ke+1) = 0.0
755 !$OMP parallel do default(shared)
756 do k=2,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
757 bdif_flx(i,j,k) = - n2(i,j,k) * kd(i,j,k)
758 enddo ; enddo ; enddo
759
760 if (present(bdif_flx_dz) .and. present(dz)) then
761 !$OMP parallel do default(shared)
762 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
763 bdif_flx_dz(i,j,k) = 0.5*(bdif_flx(i,j,k)+bdif_flx(i,j,k+1))*dz(i,j,k)
764 enddo ; enddo ; enddo
765 endif
766
767end subroutine diagnosekdwork
768
769!> Allocates arrays only when needed
770subroutine allocate_vbf_cs(G, GV, VBF)
771 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
772 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
773 type (vbf_cs), intent(inout) :: vbf !< Vertical buoyancy flux structure
774
775 integer :: isd, ied, jsd, jed, nz
776
777 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
778
779 if (vbf%do_bflx_salt) &
780 allocate(vbf%Bflx_salt(isd:ied,jsd:jed,nz+1), source=0.0)
781 if (vbf%do_bflx_salt_dz) &
782 allocate(vbf%Bflx_salt_dz(isd:ied,jsd:jed,nz), source=0.0)
783 if (vbf%do_bflx_temp) &
784 allocate(vbf%Bflx_temp(isd:ied,jsd:jed,nz+1), source=0.0)
785 if (vbf%do_bflx_temp_dz) &
786 allocate(vbf%Bflx_temp_dz(isd:ied,jsd:jed,nz), source=0.0)
787
788 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. &
789 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) &
790 allocate(vbf%Kd_salt(isd:ied,jsd:jed,nz+1), source=0.0)
791 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. &
792 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) &
793 allocate(vbf%Kd_temp(isd:ied,jsd:jed,nz+1), source=0.0)
794
795 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) &
796 allocate(vbf%Kd_BBL(isd:ied,jsd:jed,nz+1), source=0.0)
797 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) &
798 allocate(vbf%Kd_ePBL(isd:ied,jsd:jed,nz+1), source=0.0)
799 if (vbf%id_Bdif_KS>0 .or. vbf%id_Bdif_dz_KS>0 .or. vbf%id_Bdif_idz_KS>0 .or. vbf%id_Bdif_idV_KS>0) &
800 allocate(vbf%Kd_KS(isd:ied,jsd:jed,nz+1), source=0.0)
801 if (vbf%id_Bdif_bkgnd>0 .or. vbf%id_Bdif_dz_bkgnd>0 .or. vbf%id_Bdif_idz_bkgnd>0 .or. vbf%id_Bdif_idV_bkgnd>0) &
802 allocate(vbf%Kd_bkgnd(isd:ied,jsd:jed,nz+1), source=0.0)
803 if (vbf%id_Bdif_ddiff_temp>0 .or. vbf%id_Bdif_dz_ddiff_temp>0 .or. vbf%id_Bdif_idz_ddiff_temp>0 &
804 .or. vbf%id_Bdif_idV_ddiff_temp>0) allocate(vbf%Kd_ddiff_T(isd:ied,jsd:jed,nz+1), source=0.0)
805 if (vbf%id_Bdif_ddiff_salt>0 .or. vbf%id_Bdif_dz_ddiff_salt>0 .or. vbf%id_Bdif_idV_ddiff_salt>0 &
806 .or. vbf%id_Bdif_idV_ddiff_salt>0) allocate(vbf%Kd_ddiff_S(isd:ied,jsd:jed,nz+1), source=0.0)
807 if (vbf%id_Bdif_leak>0 .or. vbf%id_Bdif_dz_leak>0 .or. vbf%id_Bdif_idz_leak>0 .or. vbf%id_Bdif_idV_leak>0) &
808 allocate(vbf%Kd_leak(isd:ied,jsd:jed,nz+1), source=0.0)
809 if (vbf%id_Bdif_quad>0 .or. vbf%id_Bdif_dz_quad>0 .or. vbf%id_Bdif_idz_quad>0 .or. vbf%id_Bdif_idV_quad>0) &
810 allocate(vbf%Kd_quad(isd:ied,jsd:jed,nz+1), source=0.0)
811 if (vbf%id_Bdif_itidal>0 .or. vbf%id_Bdif_dz_itidal>0 .or. vbf%id_Bdif_idz_itidal>0 .or. vbf%id_Bdif_idV_itidal>0) &
812 allocate(vbf%Kd_itidal(isd:ied,jsd:jed,nz+1), source=0.0)
813 if (vbf%id_Bdif_Froude>0 .or. vbf%id_Bdif_dz_Froude>0 .or. vbf%id_Bdif_idz_Froude>0 .or. vbf%id_Bdif_idV_Froude>0) &
814 allocate(vbf%Kd_Froude(isd:ied,jsd:jed,nz+1), source=0.0)
815 if (vbf%id_Bdif_slope>0 .or. vbf%id_Bdif_dz_slope>0 .or. vbf%id_Bdif_idz_slope>0 .or. vbf%id_Bdif_idV_slope>0) &
816 allocate(vbf%Kd_slope(isd:ied,jsd:jed,nz+1), source=0.0)
817 if (vbf%id_Bdif_lowmode>0 .or. vbf%id_Bdif_dz_lowmode>0 .or. vbf%id_Bdif_idz_lowmode>0 .or. &
818 vbf%id_Bdif_idV_lowmode>0) allocate(vbf%Kd_lowmode(isd:ied,jsd:jed,nz+1), source=0.0)
819 if (vbf%id_Bdif_Niku>0 .or. vbf%id_Bdif_dz_Niku>0 .or. vbf%id_Bdif_idz_Niku>0 .or. vbf%id_Bdif_idV_Niku>0) &
820 allocate(vbf%Kd_Niku(isd:ied,jsd:jed,nz+1), source=0.0)
821 if (vbf%id_Bdif_itides>0 .or. vbf%id_Bdif_dz_itides>0 .or. vbf%id_Bdif_idz_itides>0 .or. vbf%id_Bdif_idV_itides>0) &
822 allocate(vbf%Kd_itides(isd:ied,jsd:jed,nz+1), source=0.0)
823
824end subroutine allocate_vbf_cs
825
826!> Deallocate any arrays that were allocated
827subroutine deallocate_vbf_cs(VBF)
828 type (vbf_cs), intent(inout) :: vbf !< Vertical buoyancy flux structure
829
830 if (associated(vbf%Bflx_salt)) &
831 deallocate(vbf%Bflx_salt)
832 if (associated(vbf%Bflx_temp)) &
833 deallocate(vbf%Bflx_temp)
834 if (associated(vbf%Bflx_salt_dz)) &
835 deallocate(vbf%Bflx_salt_dz)
836 if (associated(vbf%Bflx_temp_dz)) &
837 deallocate(vbf%Bflx_temp_dz)
838 if (associated(vbf%Kd_salt)) &
839 deallocate(vbf%Kd_salt)
840 if (associated(vbf%Kd_temp)) &
841 deallocate(vbf%Kd_temp)
842 if (associated(vbf%Kd_BBL)) &
843 deallocate(vbf%Kd_BBL)
844 if (associated(vbf%Kd_ePBL)) &
845 deallocate(vbf%Kd_ePBL)
846 if (associated(vbf%Kd_KS)) &
847 deallocate(vbf%Kd_KS)
848 if (associated(vbf%Kd_bkgnd)) &
849 deallocate(vbf%Kd_bkgnd)
850 if (associated(vbf%Kd_ddiff_T)) &
851 deallocate(vbf%Kd_ddiff_T)
852 if (associated(vbf%Kd_ddiff_S)) &
853 deallocate(vbf%Kd_ddiff_S)
854 if (associated(vbf%Kd_leak)) &
855 deallocate(vbf%Kd_leak)
856 if (associated(vbf%Kd_quad)) &
857 deallocate(vbf%Kd_quad)
858 if (associated(vbf%Kd_itidal)) &
859 deallocate(vbf%Kd_itidal)
860 if (associated(vbf%Kd_Froude)) &
861 deallocate(vbf%Kd_Froude)
862 if (associated(vbf%Kd_slope)) &
863 deallocate(vbf%Kd_slope)
864 if (associated(vbf%Kd_lowmode)) &
865 deallocate(vbf%Kd_lowmode)
866 if (associated(vbf%Kd_Niku)) &
867 deallocate(vbf%Kd_Niku)
868 if (associated(vbf%Kd_itides)) &
869 deallocate(vbf%Kd_itides)
870
871end subroutine deallocate_vbf_cs
872
873!> Handles all KdWork diagnostics and flags which calculations should be done.
874subroutine kdwork_init(Time, G,GV,US,diag,VBF,Use_KdWork_diag)
875 type(time_type), target :: time !< model time
876 type(ocean_grid_type), intent(in) :: g !< ocean grid structure
877 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
878 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
879 type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
880 type (vbf_cs), pointer, intent(inout) :: vbf !< Vertical buoyancy flux structure
881 logical, intent(out) :: use_kdwork_diag !< Flag if any output was turned on
882
883 allocate(vbf)
884
885 vbf%do_bflx_salt = .false.
886 vbf%do_bflx_salt_dz = .false.
887 vbf%do_bflx_temp = .false.
888 vbf%do_bflx_temp_dz = .false.
889
890 vbf%id_Bdif = register_diag_field('ocean_model',"Bflx_dia_diff", diag%axesTi, &
891 time, "Diffusive diapycnal buoyancy flux across interfaces", &
892 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
893 vbf%id_Bdif_dz = register_diag_field('ocean_model',"Bflx_dia_diff_dz", diag%axesTl, &
894 time, "Layerwise integral of diffusive diapycnal buoyancy flux.", &
895 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
896 vbf%id_Bdif_idz = register_diag_field('ocean_model',"Bflx_dia_diff_idz", diag%axesT1, &
897 time, "Layer integrated diffusive diapycnal buoyancy flux.", &
898 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
899 vbf%id_Bdif_idV = register_scalar_field('ocean_model',"Bflx_dia_diff_idV", time, diag, &
900 "Global integrated diffusive diapycnal buoyancy flux.", &
901 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
902
903 vbf%id_Bdif_salt = register_diag_field('ocean_model',"Bflx_salt_dia_diff", diag%axesTi, &
904 time, "Salinity contribution to diffusive diapycnal buoyancy flux across interfaces", &
905 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
906 vbf%id_Bdif_salt_dz = register_diag_field('ocean_model',"Bflx_salt_dia_diff_dz", diag%axesTl, &
907 time, "Salinity contribution to layer integral of diffusive diapycnal buoyancy flux.", &
908 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
909 vbf%id_Bdif_salt_idz = register_diag_field('ocean_model',"Bflx_salt_dia_diff_idz", diag%axesT1, &
910 time, "Salinity contribution to layer integrated diffusive diapycnal buoyancy flux.", &
911 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
912 vbf%id_Bdif_salt_idV = register_scalar_field('ocean_model',"Bflx_salt_dia_diff_idV", time, diag, &
913 "Salinity contribution to global integrated diffusive diapycnal buoyancy flux.", &
914 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
915
916 vbf%id_Bdif_temp = register_diag_field('ocean_model',"Bflx_temp_dia_diff", diag%axesTi, &
917 time, "Temperature contribution to diffusive diapycnal buoyancy flux across interfaces", &
918 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
919 vbf%id_Bdif_temp_dz = register_diag_field('ocean_model',"Bflx_temp_dia_diff_dz", diag%axesTl, &
920 time, "Temperature contribution to layer integral of diffusive diapycnal buoyancy flux.", &
921 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
922 vbf%id_Bdif_temp_idz = register_diag_field('ocean_model',"Bflx_temp_dia_diff_idz", diag%axesT1, &
923 time, "Temperature contribution to layer integrated diffusive diapycnal buoyancy flux.", &
924 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
925 vbf%id_Bdif_temp_idV = register_scalar_field('ocean_model',"Bflx_temp_dia_diff_idV", time, diag, &
926 "Temperature contribution to global integrated diffusive diapycnal buoyancy flux.", &
927 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
928
929 vbf%id_Bdif_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_BBL", diag%axesTi, &
930 time, "Diffusive diapycnal buoyancy flux across interfaces due to the BBL parameterization.", &
931 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
932 vbf%id_Bdif_dz_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_dz_BBL", diag%axesTl, &
933 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to the BBL parameterization.", &
934 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
935 vbf%id_Bdif_idz_BBL = register_diag_field('ocean_model',"Bflx_dia_diff_idz_BBL", diag%axesT1, &
936 time, "Layer integrated diffusive diapycnal buoyancy flux due to the BBL parameterization.", &
937 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
938 vbf%id_Bdif_idV_BBL = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_BBL", time, diag, &
939 "Global integrated diffusive diapycnal buoyancy flux due to BBL.", &
940 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
941
942 vbf%id_Bdif_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_ePBL", diag%axesTi, &
943 time, "Diffusive diapycnal buoyancy flux across interfaces due to ePBL", &
944 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
945 vbf%id_Bdif_dz_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ePBL", diag%axesTl, &
946 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to ePBL.", &
947 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
948 vbf%id_Bdif_idz_ePBL = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ePBL", diag%axesT1, &
949 time, "Layer integrated diffusive diapycnal buoyancy flux due to ePBL.", &
950 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
951 vbf%id_Bdif_idV_ePBL = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ePBL", time, diag, &
952 "Global integrated diffusive diapycnal buoyancy flux due to ePBL.", &
953 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
954
955 vbf%id_Bdif_KS = register_diag_field('ocean_model',"Bflx_dia_diff_KS", diag%axesTi, &
956 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kappa Shear", &
957 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
958 vbf%id_Bdif_dz_KS = register_diag_field('ocean_model',"Bflx_dia_diff_dz_KS", diag%axesTl, &
959 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to Kappa Shear.", &
960 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
961 vbf%id_Bdif_idz_KS = register_diag_field('ocean_model',"Bflx_dia_diff_idz_KS", diag%axesT1, &
962 time, "Layer integrated diffusive diapycnal buoyancy flux due to Kappa Shear.", &
963 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
964 vbf%id_Bdif_idV_KS = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_KS", time, diag, &
965 "Global integrated diffusive diapycnal buoyancy flux due to Kappa Shear.", &
966 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
967
968 vbf%id_Bdif_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_bkgnd", diag%axesTi, &
969 time, "Diffusive diapycnal buoyancy flux across interfaces due to bkgnd mixing", &
970 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
971 vbf%id_Bdif_dz_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_dz_bkgnd", diag%axesTl, &
972 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
973 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
974 vbf%id_Bdif_idz_bkgnd = register_diag_field('ocean_model',"Bflx_dia_diff_idz_bkgnd", diag%axesT1, &
975 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
976 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
977 vbf%id_Bdif_idV_bkgnd = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_bkgnd", time, diag, &
978 "Global integrated diffusive diapycnal buoyancy flux due to Kd_bkgnd.", &
979 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
980
981 vbf%id_Bdif_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_ddiff_heat", diag%axesTi, &
982 time, "Diffusive diapycnal buoyancy flux across interfaces due to double diffusion of heat", &
983 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
984 vbf%id_Bdif_dz_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ddiff_heat", diag%axesTl, &
985 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to double diffusion of heat.", &
986 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
987 vbf%id_Bdif_idz_ddiff_temp = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ddiff_heat", diag%axesT1, &
988 time, "Layer integrated diffusive diapycnal buoyancy flux due to double diffusion of heat.", &
989 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
990 vbf%id_Bdif_idV_ddiff_temp = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ddiff_heat", time, diag, &
991 "Global integrated diffusive diapycnal buoyancy flux due to double diffusion of heat.", &
992 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
993
994 vbf%id_Bdif_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_ddiff_salt", diag%axesTi, &
995 time, "Diffusive diapycnal buoyancy flux across interfaces due to double diffusion of salt", &
996 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
997 vbf%id_Bdif_dz_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_dz_ddiff_salt", diag%axesTl, &
998 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to double diffusion of salt.", &
999 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1000 vbf%id_Bdif_idz_ddiff_salt = register_diag_field('ocean_model',"Bflx_dia_diff_idz_ddiff_salt", diag%axesT1, &
1001 time, "Layer integrated diffusive diapycnal buoyancy flux due to double diffusion of salt.", &
1002 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1003 vbf%id_Bdif_idV_ddiff_salt = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_ddiff_salt", time, diag, &
1004 "Global integrated diffusive diapycnal buoyancy flux due to double diffusion of salt.", &
1005 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1006
1007 vbf%id_Bdif_leak = register_diag_field('ocean_model',"Bflx_dia_diff_leak", diag%axesTi, &
1008 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_leak mixing", &
1009 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1010 vbf%id_Bdif_dz_leak = register_diag_field('ocean_model',"Bflx_dia_diff_dz_leak", diag%axesTl, &
1011 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1012 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1013 vbf%id_Bdif_idz_leak = register_diag_field('ocean_model',"Bflx_dia_diff_idz_leak", diag%axesT1, &
1014 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1015 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1016 vbf%id_Bdif_idV_leak = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_leak", time, diag, &
1017 "Global integrated diffusive diapycnal buoyancy flux due to Kd_leak.", &
1018 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1019
1020 vbf%id_Bdif_quad = register_diag_field('ocean_model',"Bflx_dia_diff_quad", diag%axesTi, &
1021 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_quad mixing", &
1022 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1023 vbf%id_Bdif_dz_quad = register_diag_field('ocean_model',"Bflx_dia_diff_dz_quad", diag%axesTl, &
1024 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1025 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1026 vbf%id_Bdif_idz_quad = register_diag_field('ocean_model',"Bflx_dia_diff_idz_quad", diag%axesT1, &
1027 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1028 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1029 vbf%id_Bdif_idV_quad = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_quad", time, diag, &
1030 "Global integrated diffusive diapycnal buoyancy flux due to Kd_quad.", &
1031 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1032
1033 vbf%id_Bdif_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_itidal", diag%axesTi, &
1034 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_itidal mixing", &
1035 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1036 vbf%id_Bdif_dz_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_dz_itidal", diag%axesTl, &
1037 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1038 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1039 vbf%id_Bdif_idz_itidal = register_diag_field('ocean_model',"Bflx_dia_diff_idz_itidal", diag%axesT1, &
1040 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1041 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1042 vbf%id_Bdif_idV_itidal = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_itidal", time, diag, &
1043 "Global integrated diffusive diapycnal buoyancy flux due to Kd_itidal.", &
1044 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1045
1046 vbf%id_Bdif_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_Froude", diag%axesTi, &
1047 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_Froude mixing", &
1048 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1049 vbf%id_Bdif_dz_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_dz_Froude", diag%axesTl, &
1050 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1051 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1052 vbf%id_Bdif_idz_Froude = register_diag_field('ocean_model',"Bflx_dia_diff_idz_Froude", diag%axesT1, &
1053 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1054 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1055 vbf%id_Bdif_idV_Froude = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_Froude", time, diag, &
1056 "Global integrated diffusive diapycnal buoyancy flux due to Kd_Froude.", &
1057 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1058
1059 vbf%id_Bdif_slope = register_diag_field('ocean_model',"Bflx_dia_diff_slope", diag%axesTi, &
1060 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_slope mixing", &
1061 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1062 vbf%id_Bdif_dz_slope = register_diag_field('ocean_model',"Bflx_dia_diff_dz_slope", diag%axesTl, &
1063 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1064 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1065 vbf%id_Bdif_idz_slope = register_diag_field('ocean_model',"Bflx_dia_diff_idz_slope", diag%axesT1, &
1066 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1067 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1068 vbf%id_Bdif_idV_slope = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_slope", time, diag, &
1069 "Global integrated diffusive diapycnal buoyancy flux due to Kd_slope.", &
1070 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1071
1072 vbf%id_Bdif_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_lowmode", diag%axesTi, &
1073 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_lowmode mixing", &
1074 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1075 vbf%id_Bdif_dz_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_dz_lowmode", diag%axesTl, &
1076 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1077 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1078 vbf%id_Bdif_idz_lowmode = register_diag_field('ocean_model',"Bflx_dia_diff_idz_lowmode", diag%axesT1, &
1079 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1080 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1081 vbf%id_Bdif_idV_lowmode = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_lowmode", time, diag, &
1082 "Global integrated diffusive diapycnal buoyancy flux due to Kd_lowmode.", &
1083 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1084
1085 vbf%id_Bdif_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_Niku", diag%axesTi, &
1086 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_Niku mixing", &
1087 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1088 vbf%id_Bdif_dz_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_dz_Niku", diag%axesTl, &
1089 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1090 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1091 vbf%id_Bdif_idz_Niku = register_diag_field('ocean_model',"Bflx_dia_diff_idz_Niku", diag%axesT1, &
1092 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1093 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1094 vbf%id_Bdif_idV_Niku = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_Niku", time, diag, &
1095 "Global integrated diffusive diapycnal buoyancy flux due to Kd_Niku.", &
1096 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1097
1098 vbf%id_Bdif_itides = register_diag_field('ocean_model',"Bflx_dia_diff_itides", diag%axesTi, &
1099 time, "Diffusive diapycnal buoyancy flux across interfaces due to Kd_itides mixing", &
1100 "W m-3", conversion=gv%H_to_kg_m2*us%Z_to_m*us%s_to_T**3)
1101 vbf%id_Bdif_dz_itides = register_diag_field('ocean_model',"Bflx_dia_diff_dz_itides", diag%axesTl, &
1102 time, "Layerwise integral of diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1103 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1104 vbf%id_Bdif_idz_itides = register_diag_field('ocean_model',"Bflx_dia_diff_idz_itides", diag%axesT1, &
1105 time, "Layer integrated diffusive diapycnal buoyancy flux due to bkgnd mixing", &
1106 "W m-2", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3)
1107 vbf%id_Bdif_idV_itides = register_scalar_field('ocean_model',"Bflx_dia_diff_idV_itides", time, diag, &
1108 "Global integrated diffusive diapycnal buoyancy flux due to Kd_itides.", &
1109 units="W", conversion=gv%H_to_kg_m2*us%Z_to_m**2*us%s_to_T**3*us%L_to_m**2)
1110
1111 if (vbf%id_Bdif_dz>0 .or. vbf%id_Bdif_salt_dz>0 .or. vbf%id_Bdif_dz_BBL>0 .or. &
1112 vbf%id_Bdif_dz_ePBL>0 .or. vbf%id_Bdif_dz_KS>0 .or. vbf%id_Bdif_dz_bkgnd>0 .or. &
1113 vbf%id_Bdif_dz_ddiff_salt>0 .or. vbf%id_Bdif_dz_leak>0 .or. vbf%id_Bdif_dz_quad>0 .or. &
1114 vbf%id_Bdif_dz_itidal>0 .or. vbf%id_Bdif_dz_Froude>0 .or. vbf%id_Bdif_dz_slope>0 .or. &
1115 vbf%id_Bdif_dz_lowmode>0 .or. vbf%id_Bdif_dz_Niku>0 .or. vbf%id_Bdif_dz_itides>0 .or. &
1116 vbf%id_Bdif_idV>0 .or. vbf%id_Bdif_salt_idV>0 .or. vbf%id_Bdif_idV_BBL>0 .or. &
1117 vbf%id_Bdif_idV_ePBL>0 .or. vbf%id_Bdif_idV_KS>0 .or. vbf%id_Bdif_idV_bkgnd>0 .or. &
1118 vbf%id_Bdif_idV_ddiff_salt>0 .or. vbf%id_Bdif_idV_leak>0 .or. vbf%id_Bdif_idV_quad>0 .or. &
1119 vbf%id_Bdif_idV_itidal>0 .or. vbf%id_Bdif_idV_Froude>0 .or. vbf%id_Bdif_idV_slope>0 .or. &
1120 vbf%id_Bdif_idV_lowmode>0 .or. vbf%id_Bdif_idV_Niku>0 .or. vbf%id_Bdif_idV_itides>0 .or. &
1121 vbf%id_Bdif_idz>0 .or. vbf%id_Bdif_salt_idz>0 .or. vbf%id_Bdif_idz_BBL>0 .or. &
1122 vbf%id_Bdif_idz_ePBL>0 .or. vbf%id_Bdif_idz_KS>0 .or. vbf%id_Bdif_idz_bkgnd>0 .or. &
1123 vbf%id_Bdif_idz_ddiff_salt>0 .or. vbf%id_Bdif_idz_leak>0 .or. vbf%id_Bdif_idz_quad>0 .or. &
1124 vbf%id_Bdif_idz_itidal>0 .or. vbf%id_Bdif_idz_Froude>0 .or. vbf%id_Bdif_idz_slope>0 .or. &
1125 vbf%id_Bdif_idz_lowmode>0 .or. vbf%id_Bdif_idz_Niku>0 .or. vbf%id_Bdif_idz_itides>0 ) then
1126 vbf%do_bflx_salt_dz = .true.
1127 endif
1128 if (vbf%id_Bdif_dz>0 .or. vbf%id_Bdif_temp_dz>0 .or. vbf%id_Bdif_dz_BBL>0 .or. &
1129 vbf%id_Bdif_dz_ePBL>0 .or. vbf%id_Bdif_dz_KS>0 .or. vbf%id_Bdif_dz_bkgnd>0 .or. &
1130 vbf%id_Bdif_dz_ddiff_temp>0 .or. vbf%id_Bdif_dz_leak>0 .or. vbf%id_Bdif_dz_quad>0 .or. &
1131 vbf%id_Bdif_dz_itidal>0 .or. vbf%id_Bdif_dz_Froude>0 .or. vbf%id_Bdif_dz_slope>0 .or. &
1132 vbf%id_Bdif_dz_lowmode>0 .or. vbf%id_Bdif_dz_Niku>0 .or. vbf%id_Bdif_dz_itides>0 .or. &
1133 vbf%id_Bdif_idV>0 .or. vbf%id_Bdif_temp_idV>0 .or. vbf%id_Bdif_idV_BBL>0 .or. &
1134 vbf%id_Bdif_idV_ePBL>0 .or. vbf%id_Bdif_idV_KS>0 .or. vbf%id_Bdif_idV_bkgnd>0 .or. &
1135 vbf%id_Bdif_idV_ddiff_temp>0 .or. vbf%id_Bdif_idV_leak>0 .or. vbf%id_Bdif_idV_quad>0 .or. &
1136 vbf%id_Bdif_idV_itidal>0 .or. vbf%id_Bdif_idV_Froude>0 .or. vbf%id_Bdif_idV_slope>0 .or. &
1137 vbf%id_Bdif_idV_lowmode>0 .or. vbf%id_Bdif_idV_Niku>0 .or. vbf%id_Bdif_idV_itides>0 .or. &
1138 vbf%id_Bdif_idz>0 .or. vbf%id_Bdif_temp_idz>0 .or. vbf%id_Bdif_idz_BBL>0 .or. &
1139 vbf%id_Bdif_idz_ePBL>0 .or. vbf%id_Bdif_idz_KS>0 .or. vbf%id_Bdif_idz_bkgnd>0 .or. &
1140 vbf%id_Bdif_idz_ddiff_temp>0 .or. vbf%id_Bdif_idz_leak>0 .or. vbf%id_Bdif_idz_quad>0 .or. &
1141 vbf%id_Bdif_idz_itidal>0 .or. vbf%id_Bdif_idz_Froude>0 .or. vbf%id_Bdif_idz_slope>0 .or. &
1142 vbf%id_Bdif_idz_lowmode>0 .or. vbf%id_Bdif_idz_Niku>0 .or. vbf%id_Bdif_idz_itides>0 ) then
1143 vbf%do_bflx_temp_dz = .true.
1144 endif
1145 if (vbf%id_Bdif>0 .or. vbf%id_Bdif_salt>0 .or. vbf%id_Bdif_BBL>0 .or. &
1146 vbf%id_Bdif_ePBL>0 .or. vbf%id_Bdif_KS>0 .or. vbf%id_Bdif_bkgnd>0 .or. &
1147 vbf%id_Bdif_ddiff_salt>0 .or. vbf%id_Bdif_leak>0 .or. vbf%id_Bdif_quad>0 .or. &
1148 vbf%id_Bdif_itidal>0 .or. vbf%id_Bdif_Froude>0 .or. vbf%id_Bdif_slope>0 .or. &
1149 vbf%id_Bdif_lowmode>0 .or. vbf%id_Bdif_Niku>0 .or. vbf%id_Bdif_itides>0 .or. &
1150 vbf%do_bflx_salt_dz) then
1151 vbf%do_bflx_salt = .true.
1152 endif
1153 if (vbf%id_Bdif>0 .or. vbf%id_Bdif_temp>0 .or. vbf%id_Bdif_BBL>0 .or. &
1154 vbf%id_Bdif_ePBL>0 .or. vbf%id_Bdif_KS>0 .or. vbf%id_Bdif_bkgnd>0 .or. &
1155 vbf%id_Bdif_ddiff_temp>0 .or. vbf%id_Bdif_leak>0 .or. vbf%id_Bdif_quad>0 .or. &
1156 vbf%id_Bdif_itidal>0 .or. vbf%id_Bdif_Froude>0 .or. vbf%id_Bdif_slope>0 .or. &
1157 vbf%id_Bdif_lowmode>0 .or. vbf%id_Bdif_Niku>0 .or. vbf%id_Bdif_itides>0 .or. &
1158 vbf%do_bflx_temp_dz) then
1159 vbf%do_bflx_temp = .true.
1160 endif
1161
1162 use_kdwork_diag = (vbf%do_bflx_salt .or. vbf%do_bflx_temp .or. vbf%do_bflx_salt_dz .or. vbf%do_bflx_temp_dz)
1163
1164end subroutine kdwork_init
1165
1166!> Deallocates control structrue
1167subroutine kdwork_end(VBF)
1168 type (vbf_cs), pointer, intent(inout) :: vbf !< Vertical buoyancy flux structure
1169
1170 if (associated(vbf)) deallocate(vbf)
1171
1172end subroutine kdwork_end
1173
1174!> \namespace mom_diagnose_kdwork
1175!!
1176!! The subroutine diagnoseKdWork diagnoses the energetics associated with various vertical diffusivities
1177!! inside MOM6 diabatic routines.
1178!!
1179
1180end module mom_diagnose_kdwork