src/parameterizations/vertical/MOM_diabatic_aux.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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 functions for some diabatic processes such as frazil, brine rejection,
6		!! tendency due to surface flux divergence.
7		module MOM_diabatic_aux
8
9		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
10		use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
11		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
12		use MOM_diag_mediator, only : diag_ctrl, time_type
13		use MOM_EOS, only : calculate_density, calculate_TFreeze, EOS_domain
14		use MOM_EOS, only : calculate_specific_vol_derivs, calculate_density_derivs
15		use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery
16		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
17		use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
18		use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint
19		use MOM_grid, only : ocean_grid_type
20		use MOM_interface_heights, only : thickness_to_dz
21		use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init
22		use MOM_interpolate, only : external_field
23		use MOM_io, only : slasher
24		use MOM_opacity, only : set_opacity, opacity_CS, extract_optics_slice, extract_optics_fields
25		use MOM_opacity, only : optics_type, optics_nbands, absorbRemainingSW, sumSWoverBands
26		use MOM_tracer_flow_control, only : get_chl_from_model, tracer_flow_control_CS
27		use MOM_unit_scaling, only : unit_scale_type
28		use MOM_variables, only : thermo_var_ptrs
29		use MOM_verticalGrid, only : verticalGrid_type
30
31		implicit none ; private
32
33		#include <MOM_memory.h>
34
35		public diabatic_aux_init, diabatic_aux_end
36		public make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS, triDiagTS_Eulerian
37		public find_uv_at_h, applyBoundaryFluxesInOut, set_pen_shortwave
38
39		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
40		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
41		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
42		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
43
44		!> Control structure for diabatic_aux
45		type, public :: diabatic_aux_CS ; private
46		logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff at the
47		!! river mouths to a depth of "rivermix_depth"
48		real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m].
49		real :: dSalt_frac_max !< An upper limit on the fraction of the salt in a layer that can be
50		!! lost to the net surface salt fluxes within a timestep [nondim]
51		logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to
52		!! to cool the topmost layer down to the freezing
53		!! point. The default is true.
54		logical :: pressure_dependent_frazil !< If true, use a pressure dependent
55		!! freezing temperature when making frazil. The
56		!! default is false, which will be faster but is
57		!! inappropriate with ice-shelf cavities.
58		logical :: ignore_fluxes_over_land !< If true, the model does not check
59		!! if fluxes are applied over land points. This
60		!! flag must be used when the ocean is coupled with
61		!! sea ice and ice shelves and use_ePBL = true.
62		logical :: use_river_heat_content !< If true, assumes that ice-ocean boundary
63		!! has provided a river heat content. Otherwise, runoff
64		!! is added with a temperature of the local SST.
65		logical :: use_calving_heat_content !< If true, assumes that ice-ocean boundary
66		!! has provided a calving heat content. Otherwise, calving
67		!! is added with a temperature of the local SST.
68		logical :: var_pen_sw !< If true, use one of the CHL_A schemes to determine the
69		!! e-folding depth of incoming shortwave radiation.
70		type(external_field) :: sbc_chl !< A handle used in time interpolation of
71		!! chlorophyll read from a file.
72		logical :: chl_from_file !< If true, chl_a is read from a file.
73		logical :: do_brine_plume !< If true, insert salt flux below the surface according to
74		!! a parameterization by \cite Nguyen2009.
75		integer :: brine_plume_n !< The exponent in the brine plume parameterization.
76		real :: plume_strength !< Fraction of the available brine to take to the bottom of the mixed
77		!! layer [nondim].
78
79		type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
80		type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output
81
82		! Diagnostic handles
83		integer :: id_createdH = -1 !< Diagnostic ID of mass added to avoid grounding
84		integer :: id_brine_lay = -1 !< Diagnostic ID of which layer receives the brine
85		integer :: id_penSW_diag = -1 !< Diagnostic ID of Penetrative shortwave heating (flux convergence)
86		integer :: id_penSWflux_diag = -1 !< Diagnostic ID of Penetrative shortwave flux
87		integer :: id_nonpenSW_diag = -1 !< Diagnostic ID of Non-penetrative shortwave heating
88		integer :: id_Chl = -1 !< Diagnostic ID of chlorophyll-A handles for opacity
89
90		! Optional diagnostic arrays
91		real, allocatable, dimension(:,:) :: createdH !< The amount of volume added in order to
92		!! avoid grounding [H T-1 ~> m s-1]
93		real, allocatable, dimension(:,:,:) :: penSW_diag !< Heating in a layer from convergence of
94		!! penetrative SW [Q R Z T-1 ~> W m-2]
95		real, allocatable, dimension(:,:,:) :: penSWflux_diag !< Penetrative SW flux at base of grid
96		!! layer [Q R Z T-1 ~> W m-2]
97		real, allocatable, dimension(:,:) :: nonpenSW_diag !< Non-downwelling SW radiation at ocean
98		!! surface [Q R Z T-1 ~> W m-2]
99
100		end type diabatic_aux_CS
101
102		!>@{ CPU time clock IDs
103		integer :: id_clock_uv_at_h, id_clock_frazil
104		!>@}
105
106		contains
107
108		!> Frazil formation keeps the temperature above the freezing point.
109		!! This subroutine warms any water that is colder than the (currently
110		!! surface) freezing point up to the freezing point and accumulates
111		!! the required heat (in [Q R Z ~> J m-2]) in tv%frazil.
112	0	subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
113		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
114		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
115		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
116		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
117		type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any available
118		!! thermodynamic fields.
119		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
120		type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous
121		!! call to diabatic_aux_init.
122		real, dimension(SZI_(G),SZJ_(G)), &
123		optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa].
124		integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil
125		! Local variables
126		real, dimension(SZI_(G)) :: &
127	4	fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].
128	4	T_freeze, & ! The freezing potential temperature at the current salinity [C ~> degC].
129	4	ps ! Surface pressure [R L2 T-2 ~> Pa]
130		real, dimension(SZI_(G),SZK_(GV)) :: &
131	4	pressure ! The pressure at the middle of each layer [R L2 T-2 ~> Pa].
132		real :: H_to_RL2_T2 ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
133		real :: hc ! A layer's heat capacity [Q R Z C-1 ~> J m-2 degC-1].
134		logical :: T_fr_set ! True if the freezing point has been calculated for a
135		! row of points.
136		integer :: i, j, k, is, ie, js, je, nz
137
138	2	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
139	2	if (present(halo)) then
140	1	is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo
141		endif
142
143	2	call cpu_clock_begin(id_clock_frazil)
144
145	2	if (.not.CS%pressure_dependent_frazil) then
146	15886	do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo
147		else
148	0	H_to_RL2_T2 = GV%H_to_RZ * GV%g_Earth
149		endif
150		!$OMP parallel do default(shared) private(fraz_col,T_fr_set,T_freeze,hc,ps) &
151		!$OMP firstprivate(pressure) ! pressure might be set above, so should be firstprivate
152	362	do j=js,je
153	132840	ps(:) = 0.0
154	360	if (PRESENT(p_surf)) then ; do i=is,ie
155	0	ps(i) = p_surf(i,j)
156		enddo ; endif
157
158	129960	do i=is,ie ; fraz_col(i) = 0.0 ; enddo
159
160	360	if (CS%pressure_dependent_frazil) then
161	0	do i=is,ie
162	0	pressure(i,1) = ps(i) + (0.5H_to_RL2_T2)h(i,j,1)
163		enddo
164	0	do k=2,nz ; do i=is,ie
165	0	pressure(i,k) = pressure(i,k-1) + (0.5H_to_RL2_T2) (h(i,j,k) + h(i,j,k-1))
166		enddo ; enddo
167		endif
168
169	360	if (CS%reclaim_frazil) then
170	0	T_fr_set = .false.
171	0	do i=is,ie ; if (tv%frazil(i,j) > 0.0) then
172	0	if (.not.T_fr_set) then
173		call calculate_TFreeze(tv%S(i:ie,j,1), pressure(i:ie,1), T_freeze(i:ie), &
174	0	tv%eqn_of_state)
175	0	T_fr_set = .true.
176		endif
177
178	0	if (tv%T(i,j,1) > T_freeze(i)) then
179		! If frazil had previously been formed, but the surface temperature is now
180		! above freezing, cool the surface layer with the frazil heat deficit.
181	0	hc = (tv%C_pGV%H_to_RZ) h(i,j,1)
182	0	if (tv%frazil(i,j) - hc * (tv%T(i,j,1) - T_freeze(i)) <= 0.0) then
183	0	tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j) / hc
184	0	tv%frazil(i,j) = 0.0
185		else
186	0	tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - T_freeze(i))
187	0	tv%T(i,j,1) = T_freeze(i)
188		endif
189		endif
190		endif ; enddo
191		endif
192
193	8280	do k=nz,1,-1
194	7920	T_fr_set = .false.
195	2859480	do i=is,ie
196	2851200	if ((G%mask2dT(i,j) > 0.0) .and. &
197	7920	((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then
198	0	if (.not.T_fr_set) then
199		call calculate_TFreeze(tv%S(i:ie,j,k), pressure(i:ie,k), T_freeze(i:ie), &
200	0	tv%eqn_of_state)
201	0	T_fr_set = .true.
202		endif
203
204	0	hc = (tv%C_pGV%H_to_RZ) h(i,j,k)
205	0	if (h(i,j,k) <= 10.0*(GV%Angstrom_H + GV%H_subroundoff)) then
206		! Very thin layers should not be cooled by the frazil flux.
207	0	if (tv%T(i,j,k) < T_freeze(i)) then
208	0	fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k))
209	0	tv%T(i,j,k) = T_freeze(i)
210		endif
211	0	elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < T_freeze(i))) then
212	0	if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) < 0.0) then
213	0	tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc
214	0	fraz_col(i) = 0.0
215		else
216	0	fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k))
217	0	tv%T(i,j,k) = T_freeze(i)
218		endif
219		endif
220		endif
221		enddo
222		enddo
223	129962	do i=is,ie
224	129960	tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i)
225		enddo
226		enddo
227
228	2	tv%frazil_was_reset = .false.
229
230	2	call cpu_clock_end(id_clock_frazil)
231
232	2	end subroutine make_frazil
233
234		!> This subroutine applies double diffusion to T & S, assuming no diapycnal mass
235		!! fluxes, using a simple tridiagonal solver.
236	0	subroutine differential_diffuse_T_S(h, T, S, Kd_T, Kd_S, tv, dt, G, GV)
237		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
238		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
239		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
240		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
241		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
242		intent(inout) :: T !< Potential temperature [C ~> degC].
243		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
244		intent(inout) :: S !< Salinity [PSU] or [gSalt/kg], generically [S ~> ppt].
245		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
246		intent(in) :: Kd_T !< The extra diffusivity of temperature due to
247		!! double diffusion relative to the diffusivity of
248		!! density [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
249		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
250		intent(in) :: Kd_S !< The extra diffusivity of salinity due to
251		!! double diffusion relative to the diffusivity of
252		!! density [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
253		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
254		!! available thermodynamic fields.
255		real, intent(in) :: dt !< Time increment [T ~> s].
256
257		! local variables
258		real, dimension(SZI_(G)) :: &
259	0	b1_T, b1_S, & ! Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].
260	0	d1_T, d1_S ! Variables used by the tridiagonal solvers [nondim].
261		real, dimension(SZI_(G),SZK_(GV)) :: &
262	0	dz, & ! Height change across layers [Z ~> m]
263	0	c1_T, c1_S ! Variables used by the tridiagonal solvers [H ~> m or kg m-2].
264		real, dimension(SZI_(G),SZK_(GV)+1) :: &
265	0	mix_T, mix_S ! Mixing distances in both directions across each interface [H ~> m or kg m-2].
266		real :: h_tr ! h_tr is h at tracer points with a tiny thickness
267		! added to ensure positive definiteness [H ~> m or kg m-2].
268		real :: h_neglect ! A thickness that is so small it is usually lost
269		! in roundoff and can be neglected [H ~> m or kg m-2].
270		real :: dz_neglect ! A vertical distance that is so small it is usually lost
271		! in roundoff and can be neglected [Z ~> m].
272		real :: I_dz_int ! The inverse of the height scale associated with an interface [Z-1 ~> m-1].
273		real :: b_denom_T ! The first term in the denominator for the expression for b1_T [H ~> m or kg m-2].
274		real :: b_denom_S ! The first term in the denominator for the expression for b1_S [H ~> m or kg m-2].
275		integer :: i, j, k, is, ie, js, je, nz
276
277	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
278	0	h_neglect = GV%H_subroundoff
279	0	dz_neglect = GV%dZ_subroundoff
280
281		!$OMP parallel do default(private) shared(is,ie,js,je,h,h_neglect,dt,Kd_T,Kd_S,G,GV,T,S,nz)
282	0	do j=js,je
283
284		! Find the vertical distances across layers.
285	0	call thickness_to_dz(h, tv, dz, j, G, GV)
286
287	0	do i=is,ie
288	0	I_dz_int = 1.0 / (0.5 * (dz(i,1) + dz(i,2)) + dz_neglect)
289	0	mix_T(i,2) = (dt * Kd_T(i,j,2)) * I_dz_int
290	0	mix_S(i,2) = (dt * Kd_S(i,j,2)) * I_dz_int
291
292	0	h_tr = h(i,j,1) + h_neglect
293	0	b1_T(i) = 1.0 / (h_tr + mix_T(i,2))
294	0	b1_S(i) = 1.0 / (h_tr + mix_S(i,2))
295	0	d1_T(i) = h_tr * b1_T(i)
296	0	d1_S(i) = h_tr * b1_S(i)
297	0	T(i,j,1) = (b1_T(i)h_tr)T(i,j,1)
298	0	S(i,j,1) = (b1_S(i)h_tr)S(i,j,1)
299		enddo
300	0	do k=2,nz-1 ; do i=is,ie
301		! Calculate the mixing across the interface below this layer.
302	0	I_dz_int = 1.0 / (0.5 * (dz(i,k) + dz(i,k+1)) + dz_neglect)
303	0	mix_T(i,K+1) = ((dt * Kd_T(i,j,K+1))) * I_dz_int
304	0	mix_S(i,K+1) = ((dt * Kd_S(i,j,K+1))) * I_dz_int
305
306	0	c1_T(i,k) = mix_T(i,K) * b1_T(i)
307	0	c1_S(i,k) = mix_S(i,K) * b1_S(i)
308
309	0	h_tr = h(i,j,k) + h_neglect
310	0	b_denom_T = h_tr + d1_T(i)*mix_T(i,K)
311	0	b_denom_S = h_tr + d1_S(i)*mix_S(i,K)
312	0	b1_T(i) = 1.0 / (b_denom_T + mix_T(i,K+1))
313	0	b1_S(i) = 1.0 / (b_denom_S + mix_S(i,K+1))
314	0	d1_T(i) = b_denom_T * b1_T(i)
315	0	d1_S(i) = b_denom_S * b1_S(i)
316
317	0	T(i,j,k) = b1_T(i) * (h_trT(i,j,k) + mix_T(i,K)T(i,j,k-1))
318	0	S(i,j,k) = b1_S(i) * (h_trS(i,j,k) + mix_S(i,K)S(i,j,k-1))
319		enddo ; enddo
320	0	do i=is,ie
321	0	c1_T(i,nz) = mix_T(i,nz) * b1_T(i)
322	0	c1_S(i,nz) = mix_S(i,nz) * b1_S(i)
323
324	0	h_tr = h(i,j,nz) + h_neglect
325	0	b1_T(i) = 1.0 / (h_tr + d1_T(i)*mix_T(i,nz))
326	0	b1_S(i) = 1.0 / (h_tr + d1_S(i)*mix_S(i,nz))
327
328	0	T(i,j,nz) = b1_T(i) * (h_trT(i,j,nz) + mix_T(i,nz)T(i,j,nz-1))
329	0	S(i,j,nz) = b1_S(i) * (h_trS(i,j,nz) + mix_S(i,nz)S(i,j,nz-1))
330		enddo
331	0	do k=nz-1,1,-1 ; do i=is,ie
332	0	T(i,j,k) = T(i,j,k) + c1_T(i,k+1)*T(i,j,k+1)
333	0	S(i,j,k) = S(i,j,k) + c1_S(i,k+1)*S(i,j,k+1)
334		enddo ; enddo
335		enddo
336	0	end subroutine differential_diffuse_T_S
337
338		!> This subroutine keeps salinity from falling below a small but positive threshold.
339		!! This usually occurs when the ice model attempts to extract more salt then
340		!! is actually available to it from the ocean.
341	0	subroutine adjust_salt(h, tv, G, GV, CS)
342		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
343		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
344		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
345		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
346		type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
347		!! available thermodynamic fields.
348		type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous
349		!! call to diabatic_aux_init.
350
351		! local variables
352	0	real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement [S R Z ~> gSalt m-2]
353		real :: S_min !< The minimum salinity [S ~> ppt].
354		real :: mc !< A layer's mass [R Z ~> kg m-2].
355		integer :: i, j, k, is, ie, js, je, nz
356
357	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
358
359		! call cpu_clock_begin(id_clock_adjust_salt)
360
361	0	S_min = tv%min_salinity
362
363	0	salt_add_col(:,:) = 0.0
364
365		!$OMP parallel do default(shared) private(mc)
366	0	do j=js,je
367	0	do k=nz,1,-1 ; do i=is,ie
368	0	if ( (G%mask2dT(i,j) > 0.0) .and. &
369	0	((tv%S(i,j,k) < S_min) .or. (salt_add_col(i,j) > 0.0)) ) then
370	0	mc = GV%H_to_RZ * h(i,j,k)
371	0	if (h(i,j,k) <= 10.0*GV%Angstrom_H) then
372		! Very thin layers should not be adjusted by the salt flux
373	0	if (tv%S(i,j,k) < S_min) then
374	0	salt_add_col(i,j) = salt_add_col(i,j) + mc * (S_min - tv%S(i,j,k))
375	0	tv%S(i,j,k) = S_min
376		endif
377	0	elseif (salt_add_col(i,j) + mc * (S_min - tv%S(i,j,k)) <= 0.0) then
378	0	tv%S(i,j,k) = tv%S(i,j,k) - salt_add_col(i,j) / mc
379	0	salt_add_col(i,j) = 0.0
380		else
381	0	salt_add_col(i,j) = salt_add_col(i,j) + mc * (S_min - tv%S(i,j,k))
382	0	tv%S(i,j,k) = S_min
383		endif
384		endif
385		enddo ; enddo
386	0	do i=is,ie
387	0	tv%salt_deficit(i,j) = tv%salt_deficit(i,j) + salt_add_col(i,j)
388		enddo
389		enddo
390		! call cpu_clock_end(id_clock_adjust_salt)
391
392	0	end subroutine adjust_salt
393
394		!> This is a simple tri-diagonal solver for T and S.
395		!! "Simple" means it only uses arrays hold, ea and eb.
396	0	subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S)
397		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
398		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
399		integer, intent(in) :: is !< The start i-index to work on.
400		integer, intent(in) :: ie !< The end i-index to work on.
401		integer, intent(in) :: js !< The start j-index to work on.
402		integer, intent(in) :: je !< The end j-index to work on.
403		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: hold !< The layer thicknesses before entrainment,
404		!! [H ~> m or kg m-2].
405		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: ea !< The amount of fluid entrained from the layer
406		!! above within this time step [H ~> m or kg m-2]
407		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: eb !< The amount of fluid entrained from the layer
408		!! below within this time step [H ~> m or kg m-2]
409		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: T !< Layer potential temperatures [C ~> degC].
410		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: S !< Layer salinities [S ~> ppt].
411
412		! Local variables
413	0	real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
414	0	real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim].
415	0	real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim].
416		real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2].
417		integer :: i, j, k
418
419		!$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
420	0	do j=js,je
421	0	do i=is,ie
422	0	h_tr = hold(i,j,1) + GV%H_subroundoff
423	0	b1(i) = 1.0 / (h_tr + eb(i,j,1))
424	0	d1(i) = h_tr * b1(i)
425	0	T(i,j,1) = (b1(i)h_tr)T(i,j,1)
426	0	S(i,j,1) = (b1(i)h_tr)S(i,j,1)
427		enddo
428	0	do k=2,GV%ke ; do i=is,ie
429	0	c1(i,k) = eb(i,j,k-1) * b1(i)
430	0	h_tr = hold(i,j,k) + GV%H_subroundoff
431	0	b_denom_1 = h_tr + d1(i)*ea(i,j,k)
432	0	b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
433	0	d1(i) = b_denom_1 * b1(i)
434	0	T(i,j,k) = b1(i) * (h_trT(i,j,k) + ea(i,j,k)T(i,j,k-1))
435	0	S(i,j,k) = b1(i) * (h_trS(i,j,k) + ea(i,j,k)S(i,j,k-1))
436		enddo ; enddo
437	0	do k=GV%ke-1,1,-1 ; do i=is,ie
438	0	T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1)
439	0	S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1)
440		enddo ; enddo
441		enddo
442	0	end subroutine triDiagTS
443
444		!> This is a simple tri-diagonal solver for T and S, with mixing across interfaces but no net
445		!! transfer of mass.
446	0	subroutine triDiagTS_Eulerian(G, GV, is, ie, js, je, hold, ent, T, S)
447		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
448		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
449		integer, intent(in) :: is !< The start i-index to work on.
450		integer, intent(in) :: ie !< The end i-index to work on.
451		integer, intent(in) :: js !< The start j-index to work on.
452		integer, intent(in) :: je !< The end j-index to work on.
453		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: hold !< The layer thicknesses before entrainment,
454		!! [H ~> m or kg m-2].
455		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: ent !< The amount of fluid mixed across an interface
456		!! within this time step [H ~> m or kg m-2]
457		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: T !< Layer potential temperatures [C ~> degC].
458		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: S !< Layer salinities [S ~> ppt].
459
460		! Local variables
461	0	real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
462	0	real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim].
463	0	real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim].
464		real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2].
465		integer :: i, j, k
466
467		!$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
468	0	do j=js,je
469	0	do i=is,ie
470	0	h_tr = hold(i,j,1) + GV%H_subroundoff
471	0	b1(i) = 1.0 / (h_tr + ent(i,j,2))
472	0	d1(i) = h_tr * b1(i)
473	0	T(i,j,1) = (b1(i)h_tr)T(i,j,1)
474	0	S(i,j,1) = (b1(i)h_tr)S(i,j,1)
475		enddo
476	0	do k=2,GV%ke ; do i=is,ie
477	0	c1(i,k) = ent(i,j,K) * b1(i)
478	0	h_tr = hold(i,j,k) + GV%H_subroundoff
479	0	b_denom_1 = h_tr + d1(i)*ent(i,j,K)
480	0	b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1))
481	0	d1(i) = b_denom_1 * b1(i)
482	0	T(i,j,k) = b1(i) * (h_trT(i,j,k) + ent(i,j,K)T(i,j,k-1))
483	0	S(i,j,k) = b1(i) * (h_trS(i,j,k) + ent(i,j,K)S(i,j,k-1))
484		enddo ; enddo
485	0	do k=GV%ke-1,1,-1 ; do i=is,ie
486	0	T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1)
487	0	S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1)
488		enddo ; enddo
489		enddo
490	0	end subroutine triDiagTS_Eulerian
491
492
493		!> This subroutine calculates u_h and v_h (velocities at thickness
494		!! points), optionally using the entrainment amounts passed in as arguments.
495	2	subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix)
496		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
497		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
498		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
499		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
500		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
501		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
502		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
503		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
504		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
505		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
506		intent(out) :: u_h !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].
507		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
508		intent(out) :: v_h !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].
509		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
510		optional, intent(in) :: ea !< The amount of fluid entrained from the layer
511		!! above within this time step [H ~> m or kg m-2].
512		!! Omitting ea is the same as setting it to 0.
513		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
514		optional, intent(in) :: eb !< The amount of fluid entrained from the layer
515		!! below within this time step [H ~> m or kg m-2].
516		!! Omitting eb is the same as setting it to 0.
517		logical, optional, intent(in) :: zero_mix !< If true, do the calculation of u_h and
518		!! v_h as though ea and eb were being supplied with
519		!! uniformly zero values.
520
521		! Local variables
522		real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
523		real :: h_neglect ! A thickness that is so small it is usually lost
524		! in roundoff and can be neglected [H ~> m or kg m-2].
525	4	real :: b1(SZI_(G)) ! A thickness used in the tridiagonal solver [H ~> m or kg m-2]
526	4	real :: c1(SZI_(G),SZK_(GV)) ! A variable used in the tridiagonal solver [nondim]
527	4	real :: d1(SZI_(G)) ! The complement of c1 [nondim]
528		! Fractional weights of the neighboring velocity points, ~1/2 in the open ocean.
529	4	real :: a_n(SZI_(G)), a_s(SZI_(G)) ! Fractional weights of the neighboring velocity points [nondim]
530	4	real :: a_e(SZI_(G)), a_w(SZI_(G)) ! Fractional weights of the neighboring velocity points [nondim]
531		real :: sum_area ! A sum of adjacent areas [L2 ~> m2]
532		real :: Idenom ! The inverse of the denominator in a weighted average [L-2 ~> m-2]
533		logical :: mix_vertically, zero_mixing
534		integer :: i, j, k, is, ie, js, je, nz
535	2	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
536	2	call cpu_clock_begin(id_clock_uv_at_h)
537	2	h_neglect = GV%H_subroundoff
538
539	2	mix_vertically = present(ea)
540	2	if (present(ea) .neqv. present(eb)) call MOM_error(FATAL, &
541		"find_uv_at_h: Either both ea and eb or neither one must be present "// &
542	0	"in call to find_uv_at_h.")
543	2	zero_mixing = .false. ; if (present(zero_mix)) zero_mixing = zero_mix
544	2	if (zero_mixing) mix_vertically = .false.
545		!$omp target enter data map(alloc: a_w,a_e,a_s,a_n,b1,d1,c1)
546		!$omp target teams loop private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1) &
547		!$omp map(to: ea, eb, h) map(from: u_h, v_h)
548	362	do j=js,je
549	360	do concurrent (i=is:ie)
550	129600	sum_area = G%areaCu(I-1,j) + G%areaCu(I,j)
551	129600	if (sum_area > 0.0) then
552		! If this were a simple area weighted average, this would just be I_denom = 1.0 / sum_area.
553		! The other factor of sqrt(0.5sum_areaG%IareaT(i,j)) is 1 for open ocean points on a
554		! Cartesian grid. This construct predates the initial commit of the MOM6 code, and was
555		! present in the GOLD code before February, 2010. I do not recall why this was added, and
556		! the GOLD CVS server that contained the relevant history and logs appears to have been
557		! decommissioned.
558	92500	Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
559	92500	a_w(i) = G%areaCu(I-1,j) * Idenom
560	92500	a_e(i) = G%areaCu(I,j) * Idenom
561		else
562	37100	a_w(i) = 0.0 ; a_e(i) = 0.0
563		endif
564
565	129600	sum_area = G%areaCv(i,J-1) + G%areaCv(i,J)
566	259560	if (sum_area > 0.0) then
567	92500	Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
568	92500	a_s(i) = G%areaCv(i,J-1) * Idenom
569	92500	a_n(i) = G%areaCv(i,J) * Idenom
570		else
571	37100	a_s(i) = 0.0 ; a_n(i) = 0.0
572		endif
573		enddo
574
575	362	if (mix_vertically) then
576	180	do concurrent (i=is:ie)
577	64800	b_denom_1 = h(i,j,1) + h_neglect
578	64800	b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
579	64800	d1(i) = b_denom_1 * b1(i)
580	64800	u_h(i,j,1) = (h(i,j,1)b1(i)) ((a_e(i)u(I,j,1)) + (a_w(i)u(I-1,j,1)))
581	64980	v_h(i,j,1) = (h(i,j,1)b1(i)) ((a_n(i)v(i,J,1)) + (a_s(i)v(i,J-1,1)))
582		enddo
583	3960	do k=2,nz ; do concurrent (i=is:ie)
584	1360800	c1(i,k) = eb(i,j,k-1) * b1(i)
585	1360800	b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
586	1360800	b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
587	1360800	d1(i) = b_denom_1 * b1(i)
588		u_h(i,j,k) = (h(i,j,k) * ((a_e(i)u(I,j,k)) + (a_w(i)u(I-1,j,k))) + &
589	1360800	ea(i,j,k)u_h(i,j,k-1))b1(i)
590		v_h(i,j,k) = (h(i,j,k) * ((a_n(i)v(i,J,k)) + (a_s(i)v(i,J-1,k))) + &
591	1364580	ea(i,j,k)v_h(i,j,k-1))b1(i)
592		enddo ; enddo
593	3960	do k=nz-1,1,-1 ; do concurrent (i=is:ie)
594	1360800	u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
595	1364580	v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
596		enddo ; enddo
597	180	elseif (zero_mixing) then
598	0	do concurrent (i=is:ie)
599	0	b1(i) = 1.0 / (h(i,j,1) + h_neglect)
600	0	u_h(i,j,1) = (h(i,j,1)b1(i)) ((a_e(i)u(I,j,1)) + (a_w(i)u(I-1,j,1)))
601	0	v_h(i,j,1) = (h(i,j,1)b1(i)) ((a_n(i)v(i,J,1)) + (a_s(i)v(i,J-1,1)))
602		enddo
603	0	do concurrent (k=2:nz, i=is:ie)
604	0	b1(i) = 1.0 / (h(i,j,k) + h_neglect)
605	0	u_h(i,j,k) = (h(i,j,k) * ((a_e(i)u(I,j,k)) + (a_w(i)u(I-1,j,k)))) * b1(i)
606	0	v_h(i,j,k) = (h(i,j,k) * ((a_n(i)v(i,J,k)) + (a_s(i)v(i,J-1,k)))) * b1(i)
607		enddo
608		else
609	180	do concurrent (k=1:nz, i=is:ie)
610	1425600	u_h(i,j,k) = (a_e(i)u(I,j,k)) + (a_w(i)u(I-1,j,k))
611	1490580	v_h(i,j,k) = (a_n(i)v(i,J,k)) + (a_s(i)v(i,J-1,k))
612		enddo
613		endif
614		enddo
615		!$omp target exit data map(release: a_w,a_e,a_s,a_n,b1,d1,c1)
616
617	2	call cpu_clock_end(id_clock_uv_at_h)
618	3	end subroutine find_uv_at_h
619
620		!> Estimate the optical properties of the water column and determine the penetrating shortwave
621		!! radiation by band, extracting the relevant information from the fluxes type and storing it
622		!! in the optics type for later application. This routine is effectively a wrapper for
623		!! set_opacity with added error handling and diagnostics.
624	1	subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity, tracer_flow_CSp)
625		type(optics_type), pointer :: optics !< An optics structure that has will contain
626		!! information about shortwave fluxes and absorption.
627		type(forcing), intent(inout) :: fluxes !< points to forcing fields
628		!! unused fields have NULL pointers
629		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
630		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
631		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
632		type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
633		type(opacity_CS) :: opacity !< The control structure for the opacity module.
634		type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< A pointer to the control structure
635		!! organizing the tracer modules.
636
637		! Local variables
638	1	real, dimension(SZI_(G),SZJ_(G)) :: chl_2d !< Vertically uniform chlorophyll-A concentrations [mg m-3]
639	1	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: chl_3d !< The chlorophyll-A concentrations of each layer [mg m-3]
640		character(len=128) :: mesg
641		integer :: i, j, is, ie, js, je
642	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
643
644	1	if (.not.associated(optics)) return
645
646	1	if (CS%var_pen_sw) then
647	0	if (CS%chl_from_file) then
648		! Only the 2-d surface chlorophyll can be read in from a file. The
649		! same value is assumed for all layers.
650	0	call time_interp_external(CS%sbc_chl, CS%Time, chl_2d, turns=G%HI%turns)
651	0	do j=js,je ; do i=is,ie
652	0	if ((G%mask2dT(i,j) > 0.0) .and. (chl_2d(i,j) < 0.0)) then
653		write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
654		& I0,", ",I0," lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
655	0	chl_2d(i,j), i, j, G%geoLonT(i,j), G%geoLatT(i,j)
656	0	call MOM_error(FATAL, "MOM_diabatic_aux set_pen_shortwave: "//trim(mesg))
657		endif
658		enddo ; enddo
659
660	0	if (CS%id_chl > 0) call post_data(CS%id_chl, chl_2d, CS%diag)
661
662		call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
663	0	fluxes%sw_nir_dir, fluxes%sw_nir_dif, G, GV, US, opacity, chl_2d=chl_2d)
664		else
665	0	if (.not.associated(tracer_flow_CSp)) call MOM_error(FATAL, &
666		"The tracer flow control structure must be associated when the model sets "//&
667	0	"the chlorophyll internally in set_pen_shortwave.")
668	0	call get_chl_from_model(chl_3d, G, GV, tracer_flow_CSp)
669
670	0	if (CS%id_chl > 0) call post_data(CS%id_chl, chl_3d(:,:,1), CS%diag)
671
672		call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
673	0	fluxes%sw_nir_dir, fluxes%sw_nir_dif, G, GV, US, opacity, chl_3d=chl_3d)
674		endif
675		else
676		call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
677	1	fluxes%sw_nir_dir, fluxes%sw_nir_dif, G, GV, US, opacity)
678		endif
679
680		end subroutine set_pen_shortwave
681
682		!> Update the thickness, temperature, and salinity due to thermodynamic
683		!! boundary forcing (contained in fluxes type) applied to h, tv%T and tv%S,
684		!! and calculate the TKE implications of this heating.
685	0	subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, &
686		aggregate_FW_forcing, evap_CFL_limit, &
687	0	minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
688	0	SkinBuoyFlux, MLD_h)
689		type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
690		type(ocean_grid_type), intent(in) :: G !< Grid structure
691		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
692		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
693		real, intent(in) :: dt !< Time-step over which forcing is applied [T ~> s]
694		type(forcing), intent(inout) :: fluxes !< Surface fluxes container
695		type(optics_type), pointer :: optics !< Optical properties container
696		integer, intent(in) :: nsw !< The number of frequency bands of penetrating
697		!! shortwave radiation
698		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
699		intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
700		type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
701		!! available thermodynamic fields.
702		logical, intent(in) :: aggregate_FW_forcing !< If False, treat in/out fluxes separately.
703		real, intent(in) :: evap_CFL_limit !< The largest fraction of a layer that
704		!! can be evaporated in one time-step [nondim].
705		real, intent(in) :: minimum_forcing_depth !< The smallest depth over which
706		!! heat and freshwater fluxes is applied [H ~> m or kg m-2].
707		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
708		optional, intent(out) :: cTKE !< Turbulent kinetic energy requirement to mix
709		!! forcing through each layer [R Z3 T-2 ~> J m-2]
710		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
711		optional, intent(out) :: dSV_dT !< Partial derivative of specific volume with
712		!! potential temperature [R-1 C-1 ~> m3 kg-1 degC-1].
713		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
714		optional, intent(out) :: dSV_dS !< Partial derivative of specific volume with
715		!! salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
716		real, dimension(SZI_(G),SZJ_(G)), &
717		optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].
718		real, dimension(:,:), &
719		optional, pointer :: MLD_h !< Mixed layer thickness for brine plumes [H ~> m or kg m-2]
720
721		! Local variables
722		integer, parameter :: maxGroundings = 5
723		integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings)
724		real :: H_limit_fluxes ! Surface fluxes are scaled down fluxes when the total depth of the ocean
725		! drops below this value [H ~> m or kg m-2]
726		real :: IforcingDepthScale ! The inverse of the layer thickness below which mass losses are
727		! shifted to the next deeper layer [H ~> m or kg m-2]
728		real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
729		real :: dThickness ! The change in layer thickness [H ~> m or kg m-2]
730		real :: dTemp ! The integrated change in layer temperature [C H ~> degC m or degC kg m-2]
731		real :: dSalt ! The integrated change in layer salinity [S H ~> ppt m or ppt kg m-2]
732		real :: fractionOfForcing ! THe fraction of the remaining forcing applied to a layer [nondim]
733		real :: hOld ! The original thickness of a layer [H ~> m or kg m-2]
734		real :: Ithickness ! The inverse of the new layer thickness [H-1 ~> m-1 or m2 kg-1]
735		real :: RivermixConst ! A constant used in implementing river mixing [R Z2 T-1 ~> Pa s].
736		real :: EnthalpyConst ! A constant used to control the enthalpy calculation [nondim]
737		! By default EnthalpyConst = 1.0. If fluxes%heat_content_evap
738		! is associated enthalpy is provided via coupler and EnthalpyConst = 0.0.
739		real, dimension(SZI_(G)) :: &
740	0	d_pres, & ! pressure change across a layer [R L2 T-2 ~> Pa]
741	0	p_lay, & ! average pressure in a layer [R L2 T-2 ~> Pa]
742	0	pres, & ! pressure at an interface [R L2 T-2 ~> Pa]
743	0	netMassInOut, & ! surface water fluxes [H ~> m or kg m-2] over time step
744	0	netMassIn, & ! mass entering ocean surface [H ~> m or kg m-2] over a time step
745	0	netMassOut, & ! mass leaving ocean surface [H ~> m or kg m-2] over a time step
746	0	netHeat, & ! heat via surface fluxes excluding Pen_SW_bnd and netMassOut
747		! [C H ~> degC m or degC kg m-2]
748	0	netSalt, & ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
749		! [S H ~> ppt m or ppt kg m-2]
750	0	nonpenSW, & ! non-downwelling SW, which is absorbed at ocean surface
751		! [C H ~> degC m or degC kg m-2]
752	0	SurfPressure, & ! Surface pressure (approximated as 0.0) [R L2 T-2 ~> Pa]
753	0	dRhodT, & ! change in density per change in temperature [R C-1 ~> kg m-3 degC-1]
754	0	dRhodS, & ! change in density per change in salinity [R S-1 ~> kg m-3 ppt-1]
755	0	dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
756	0	dSpV_dS, & ! Partial derivative of specific volume with to salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
757	0	netheat_rate, & ! netheat but for dt=1 [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
758	0	netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate)
759		! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
760	0	netMassInOut_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
761	0	mixing_depth, & ! The mixing depth for brine plumes [H ~> m or kg m-2]
762	0	total_h ! Total thickness of the water column [H ~> m or kg m-2]
763		real, dimension(SZI_(G), SZK_(GV)) :: &
764	0	h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
765		! dz, & ! Layer thicknesses in depth units [Z ~> m]
766	0	T2d, & ! A 2-d copy of the layer temperatures [C ~> degC]
767	0	pen_TKE_2d, & ! The TKE required to homogenize the heating by shortwave radiation within
768		! a layer [R Z3 T-2 ~> J m-2]
769	0	dSV_dT_2d ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
770		real, dimension(SZI_(G)) :: &
771	0	netPen_rate ! The surface penetrative shortwave heating rate summed over all bands
772		! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
773		real, dimension(max(nsw,1),SZI_(G)) :: &
774	0	Pen_SW_bnd, & ! The penetrative shortwave heating integrated over a timestep by band
775		! [C H ~> degC m or degC kg m-2]
776	0	Pen_SW_bnd_rate ! The penetrative shortwave heating rate by band
777		! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
778		real, dimension(max(nsw,1),SZI_(G),SZK_(GV)) :: &
779	0	opacityBand ! The opacity (inverse of the exponential absorption length) of each frequency
780		! band of shortwave radiation in each layer [H-1 ~> m-1 or m2 kg-1]
781		real, dimension(maxGroundings) :: hGrounding ! Thickness added by each grounding event [H ~> m or kg m-2]
782		real :: Temp_in ! The initial temperature of a layer [C ~> degC]
783		real :: Salin_in ! The initial salinity of a layer [S ~> ppt]
784		real :: g_Hconv2 ! A conversion factor for use in the TKE calculation
785		! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
786		real :: GoRho ! g_Earth times a unit conversion factor divided by density
787		! [Z T-2 R-1 ~> m4 s-2 kg-1]
788		real :: g_conv ! The gravitational acceleration times the conversion factors from non-Boussinesq
789		! thickness units to mass per units area [R Z2 H-1 T-2 ~> kg m-2 s-2 or m s-2]
790		logical :: calculate_energetics ! If true, calculate the energy required to mix the newly added
791		! water over the topmost grid cell, assuming that the fluxes of heat and salt
792		! and rejected brine are initially applied in vanishingly thin layers at the
793		! top of the layer before being mixed throughout the layer.
794		logical :: calculate_buoyancy ! If true, calculate the surface buoyancy flux.
795	0	real :: dK(SZI_(G)) ! Depth of the layer center in thickness units [H ~> m or kg m-2]
796	0	real :: A_brine(SZI_(G)) ! Constant [H-(n+1) ~> m-(n+1) or m(2n+2) kg-(n+1)].
797		real :: fraction_left_brine ! Fraction of the brine that has not been applied yet [nondim]
798		real :: plume_fraction ! Fraction of the brine that is applied to a layer [nondim]
799		real :: plume_flux ! Brine flux to move downwards [S H ~> ppt m or ppt kg m-2]
800		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
801		integer :: i, j, is, ie, js, je, k, nz, nb
802		character(len=45) :: mesg
803
804	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
805
806	0	Idt = 1.0 / dt
807	0	plume_flux = 0.0
808
809	0	calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS))
810	0	calculate_buoyancy = present(SkinBuoyFlux)
811	0	if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0
812	0	if (present(cTKE)) cTKE(:,:,:) = 0.0
813	0	g_Hconv2 = (GV%g_Earth_Z_T2 * GV%H_to_RZ) * GV%H_to_RZ
814	0	EOSdom(:) = EOS_domain(G%HI)
815
816		! Only apply forcing if fluxes%sw is associated.
817	0	if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return
818
819	0	EnthalpyConst = 1.0
820	0	if (associated(fluxes%heat_content_evap)) EnthalpyConst = 0.0
821
822	0	if (calculate_buoyancy) then
823	0	SurfPressure(:) = 0.0
824	0	GoRho = GV%g_Earth_Z_T2 / GV%Rho0
825		endif
826
827	0	if (CS%do_brine_plume .and. .not.present(MLD_h)) then
828		call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
829		"Brine plume parameterization requires a mixed-layer depth argument,\n"//&
830	0	"currently coming from the energetic PBL scheme.")
831		endif
832	0	if (CS%do_brine_plume .and. .not.associated(MLD_h)) then
833		call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
834	0	"Brine plume parameterization requires an associated mixed-layer depth.")
835		endif
836	0	if (CS%do_brine_plume .and. .not. associated(fluxes%salt_left_behind)) then
837		call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
838		"Brine plume parameterization requires DO_BRINE_PLUME\n"//&
839	0	"to be turned on in SIS2 as well as MOM6.")
840		endif
841
842		! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total
843		! depth of the ocean is vanishing. It does not (yet) handle a value of zero.
844		! To accommodate vanishing upper layers, we need to allow for an instantaneous
845		! distribution of forcing over some finite vertical extent. The bulk mixed layer
846		! code handles this issue properly.
847	0	H_limit_fluxes = max(GV%Angstrom_H, GV%H_subroundoff)
848
849		! diagnostic to see if need to create mass to avoid grounding
850	0	if (CS%id_createdH>0) CS%createdH(:,:) = 0.
851	0	numberOfGroundings = 0
852
853		!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,tv,nsw,G,GV,US,optics,fluxes, &
854		!$OMP H_limit_fluxes,numberOfGroundings,iGround,jGround,&
855		!$OMP nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing, &
856		!$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, &
857		!$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,&
858		!$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, &
859		!$OMP EnthalpyConst,MLD_h) &
860		!$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &
861		!$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
862		!$OMP IforcingDepthScale,g_conv,dSpV_dT,dSpV_dS, &
863		!$OMP dThickness,dTemp,dSalt,hOld,Ithickness, &
864		!$OMP netMassIn,pres,d_pres,p_lay,dSV_dT_2d, &
865		!$OMP netmassinout_rate,netheat_rate,netsalt_rate, &
866		!$OMP drhodt,drhods,pen_sw_bnd_rate, &
867		!$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, &
868		!$OMP mixing_depth,A_brine,fraction_left_brine, &
869		!$OMP plume_fraction,dK,total_h) &
870		!$OMP firstprivate(SurfPressure,plume_flux)
871	0	do j=js,je
872		! Work in vertical slices for efficiency
873
874		! Copy state into 2D-slice arrays
875	0	do k=1,nz ; do i=is,ie
876	0	h2d(i,k) = h(i,j,k)
877	0	T2d(i,k) = tv%T(i,j,k)
878		enddo ; enddo
879
880	0	if (calculate_energetics) then
881		! The partial derivatives of specific volume with temperature and
882		! salinity need to be precalculated to avoid having heating of
883		! tiny layers give nonsensical values.
884	0	if (associated(tv%p_surf)) then
885	0	do i=is,ie ; pres(i) = tv%p_surf(i,j) ; enddo
886		else
887	0	do i=is,ie ; pres(i) = 0.0 ; enddo
888		endif
889	0	do k=1,nz
890	0	do i=is,ie
891	0	d_pres(i) = (GV%g_Earth * GV%H_to_RZ) * h2d(i,k)
892	0	p_lay(i) = pres(i) + 0.5*d_pres(i)
893	0	pres(i) = pres(i) + d_pres(i)
894		enddo
895		call calculate_specific_vol_derivs(T2d(:,k), tv%S(:,j,k), p_lay(:), &
896	0	dSV_dT(:,j,k), dSV_dS(:,j,k), tv%eqn_of_state, EOSdom)
897	0	do i=is,ie ; dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; enddo
898		enddo
899	0	pen_TKE_2d(:,:) = 0.0
900		endif
901
902		! Nothing more is done on this j-slice if there is no buoyancy forcing.
903	0	if (.not.associated(fluxes%sw)) cycle
904
905	0	if (nsw>0) then
906	0	if (GV%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
907	0	call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=GV%H_to_Z)
908		else
909		call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=GV%H_to_RZ, &
910	0	SpV_avg=tv%SpV_avg)
911		endif
912		endif
913
914		! The surface forcing is contained in the fluxes type.
915		! We aggregate the thermodynamic forcing for a time step into the following:
916		! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step
917		! = lprec + fprec + vprec + evap + lrunoff + frunoff
918		! note that lprec generally has sea ice melt/form included.
919		! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.
920		! netMassOut < 0 means mass leaves ocean.
921		! netHeat = heat via surface fluxes [C H ~> degC m or degC kg m-2], excluding the part
922		! contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
923		! netSalt = surface salt fluxes [S H ~> ppt m or gSalt m-2]
924		! Pen_SW_bnd = components to penetrative shortwave radiation split according to bands.
925		! This field provides that portion of SW from atmosphere that in fact
926		! enters to the ocean and participates in penetrative SW heating.
927		! nonpenSW = non-downwelling SW flux, which is absorbed in ocean surface
928		! (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.
929
930		!----------------------------------------------------------------------------------------
931		!BGR-June 26, 2017{
932		!Temporary action to preserve answers while fixing a bug.
933		! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns
934		! the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning
935		! a second call to extractfluxes1d (causing the diagnostic net_heat to be incorrect).
936		! Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,
937		! which means it used the T/S fields after this routine. Therefore, the surface
938		! buoyancy flux is computed here at the very end of this routine for legacy reasons.
939		! A few specific notes follow:
940		! 1) The old method did not included river/calving contributions to heat flux. This
941		! is kept consistent here via commenting code in the present extractFluxes1d <_rate>
942		! outputs, but we may reconsider this approach.
943		! 2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead
944		! of computing the integrated value (and dividing by dt). Hence the required
945		! additional outputs from extractFluxes1d.
946		! *** This is because: A*dt/dt =/= A due to round off.
947		! 3) The old method computed buoyancy flux after this routine, meaning the returned
948		! surface fluxes (from extractfluxes1d) must be recorded for use later in the code.
949		! We could (and maybe should) move that loop up to before the surface fluxes are
950		! applied, but this will change answers.
951		! For all these reasons we compute additional values of <_rate> which are preserved
952		! for the buoyancy flux calculation and reproduce the old answers.
953		! In the future this needs more detailed investigation to make sure everything is
954		! consistent and correct. These details should not significantly effect climate,
955		! but do change answers.
956		!-----------------------------------------------------------------------------------------
957	0	if (calculate_buoyancy) then
958		call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
959		H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
960		h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, &
961		Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW, &
962		net_Heat_rate=netheat_rate, net_salt_rate=netsalt_rate, &
963	0	netmassinout_rate=netmassinout_rate, pen_sw_bnd_rate=pen_sw_bnd_rate)
964		else
965		call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
966		H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
967		h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, &
968	0	Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW)
969		endif
970		! ea is for passive tracers
971	0	do i=is,ie
972		! ea(i,j,1) = netMassInOut(i)
973	0	if (aggregate_FW_forcing) then
974	0	netMassOut(i) = netMassInOut(i)
975	0	netMassIn(i) = 0.
976		else
977	0	netMassIn(i) = netMassInOut(i) - netMassOut(i)
978		endif
979	0	if (G%mask2dT(i,j) > 0.0) then
980	0	fluxes%netMassOut(i,j) = netMassOut(i)
981	0	fluxes%netMassIn(i,j) = netMassIn(i)
982		else
983	0	fluxes%netMassOut(i,j) = 0.0
984	0	fluxes%netMassIn(i,j) = 0.0
985		endif
986		enddo
987
988		! Apply the surface boundary fluxes in three steps:
989		! A/ update mass, temp, and salinity due to all terms except mass leaving
990		! ocean (and corresponding outward heat content), and ignoring penetrative SW.
991		! B/ update mass, salt, temp from mass leaving ocean.
992		! C/ update temp due to penetrative SW
993	0	if (CS%do_brine_plume) then
994		! Find the plume mixing depth.
995	0	do i=is,ie ; total_h(i) = 0.0 ; enddo
996	0	do k=1,nz ; do i=is,ie ; total_h(i) = total_h(i) + h(i,j,k) ; enddo ; enddo
997	0	do i=is,ie
998		mixing_depth(i) = min( max(MLD_h(i,j) - minimum_forcing_depth, minimum_forcing_depth), &
999	0	max(total_h(i), GV%angstrom_h) ) + GV%H_subroundoff
1000	0	A_brine(i) = (CS%brine_plume_n + 1) / (mixing_depth(i) ** (CS%brine_plume_n + 1))
1001		enddo
1002		endif
1003
1004	0	do i=is,ie
1005	0	if (G%mask2dT(i,j) > 0.) then
1006
1007		! A/ Update mass, temp, and salinity due to incoming mass flux.
1008	0	do k=1,1
1009
1010		! Change in state due to forcing
1011	0	dThickness = netMassIn(i) ! Since we are adding mass, we can use all of it
1012	0	dTemp = 0.
1013	0	dSalt = 0.
1014
1015		! Update the forcing by the part to be consumed within the present k-layer.
1016		! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
1017	0	netMassIn(i) = netMassIn(i) - dThickness
1018		! This line accounts for the temperature of the mass exchange
1019	0	Temp_in = T2d(i,k)
1020	0	Salin_in = 0.0
1021	0	dTemp = dTemp + dThicknessTemp_inEnthalpyConst
1022
1023		! Diagnostics of heat content associated with mass fluxes
1024	0	if (.not. associated(fluxes%heat_content_evap)) then
1025	0	if (associated(fluxes%heat_content_massin)) &
1026		fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1027	0	T2d(i,k) * max(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
1028	0	if (associated(fluxes%heat_content_massout)) &
1029		fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1030	0	T2d(i,k) * min(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
1031	0	if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1032	0	T2d(i,k) * dThickness * GV%H_to_RZ
1033		endif
1034
1035		! Determine the energetics of river mixing before updating the state.
1036	0	if (calculate_energetics .and. associated(fluxes%lrunoff) .and. CS%do_rivermix) then
1037		! Here we add an additional source of TKE to the mixed layer where river
1038		! is present to simulate unresolved estuaries. The TKE input, TKE_river in
1039		! [Z3 T-3 ~> m3 s-3], is diagnosed as follows:
1040		! TKE_river = 0.5rivermix_depthg(1/rho)drho_ds*
1041		! River(Samb - Sriver) = CS%mstarU_star^3
1042		! where River is in units of [Z T-1 ~> m s-1].
1043		! Samb = Ambient salinity at the mouth of the estuary
1044		! rivermix_depth = The prescribed depth over which to mix river inflow
1045		! drho_ds = The derivative of density with salt at the ambient surface salinity.
1046		! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
1047	0	if (GV%Boussinesq) then
1048	0	RivermixConst = -0.5(CS%rivermix_depthdt) * GV%g_Earth_Z_T2 * GV%Rho0
1049	0	elseif (allocated(tv%SpV_avg)) then
1050	0	RivermixConst = -0.5(CS%rivermix_depthdt) * GV%g_Earth_Z_T2 / tv%SpV_avg(i,j,1)
1051		else
1052	0	RivermixConst = -0.5(CS%rivermix_depthdt) * GV%Rho0 * GV%g_Earth_Z_T2
1053		endif
1054		cTKE(i,j,k) = cTKE(i,j,k) + max(0.0, RivermixConstdSV_dS(i,j,1) &
1055		((fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) + &
1056	0	(fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j))) * tv%S(i,j,1))
1057		endif
1058
1059		! Update state
1060	0	hOld = h2d(i,k) ! Keep original thickness in hand
1061	0	h2d(i,k) = h2d(i,k) + dThickness ! New thickness
1062	0	if (h2d(i,k) > 0.0) then
1063	0	if (calculate_energetics .and. (dThickness > 0.)) then
1064		! Calculate the energy required to mix the newly added water over
1065		! the topmost grid cell.
1066		cTKE(i,j,k) = cTKE(i,j,k) + 0.5g_Hconv2(hOlddThickness) &
1067	0	((T2d(i,k) - Temp_in) * dSV_dT(i,j,k) + (tv%S(i,j,k) - Salin_in) * dSV_dS(i,j,k))
1068		endif
1069	0	Ithickness = 1.0/h2d(i,k) ! Inverse new thickness
1070		! The "if"s below avoid changing T/S by roundoff unnecessarily
1071	0	if (dThickness /= 0. .or. dTemp /= 0.) T2d(i,k) = (hOldT2d(i,k) + dTemp)Ithickness
1072	0	if (dThickness /= 0. .or. dSalt /= 0.) tv%S(i,j,k) = (hOldtv%S(i,j,k) + dSalt)Ithickness
1073
1074		endif
1075
1076		enddo ! k=1,1
1077
1078		! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.
1079	0	fraction_left_brine = 1.0
1080	0	do k=1,nz
1081		! Place forcing into this layer if this layer has nontrivial thickness.
1082		! For layers thin relative to 1/IforcingDepthScale, then distribute
1083		! forcing into deeper layers.
1084	0	IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) )
1085		! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
1086	0	fractionOfForcing = min(1.0, h2d(i,k)*IforcingDepthScale)
1087
1088		! In the case with (-1)netMassOutfractionOfForcing greater than cfl*h, we
1089		! limit the forcing applied to this cell, leaving the remaining forcing to
1090		! be distributed downwards.
1091	0	if (-fractionOfForcingnetMassOut(i) > evap_CFL_limith2d(i,k)) then
1092	0	fractionOfForcing = -evap_CFL_limit*h2d(i,k)/netMassOut(i)
1093		endif
1094
1095	0	if (CS%do_brine_plume .and. associated(fluxes%salt_left_behind)) then
1096	0	if (fluxes%salt_left_behind(i,j) > 0 .and. fraction_left_brine > 0.0) then
1097		! Place forcing into this layer by depth for brine plume parameterization.
1098	0	if (k == 1) then
1099	0	dK(i) = 0.5 * h(i,j,k) ! Depth of center of layer K
1100	0	plume_flux = - (1000.0US%ppt_to_S (CS%plume_strength * fluxes%salt_left_behind(i,j))) * GV%RZ_to_H
1101	0	plume_fraction = 1.0
1102		else
1103	0	dK(i) = dK(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) ! Depth of center of layer K
1104	0	plume_flux = 0.0
1105		endif
1106	0	if (dK(i) <= mixing_depth(i) .and. fraction_left_brine > 0.0) then
1107	0	plume_fraction = min(fraction_left_brine, (A_brine(i) * dK(i)*CS%brine_plume_n) h(i,j,k))
1108		else
1109	0	IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) )
1110		! plume_fraction = fraction_left_brine, unless h2d is less than IforcingDepthScale.
1111	0	plume_fraction = min(fraction_left_brine, h2d(i,k)*IforcingDepthScale)
1112		endif
1113	0	fraction_left_brine = fraction_left_brine - plume_fraction
1114		plume_flux = plume_flux + plume_fraction * (1000.0US%ppt_to_S (CS%plume_strength * &
1115	0	fluxes%salt_left_behind(i,j))) * GV%RZ_to_H
1116		else
1117	0	plume_flux = 0.0
1118		endif
1119		endif
1120
1121		! Change in state due to forcing
1122
1123	0	dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) )
1124	0	dTemp = fractionOfForcing*netHeat(i)
1125	0	dSalt = max( fractionOfForcingnetSalt(i), -CS%dSalt_frac_max h2d(i,k) * tv%S(i,j,k))
1126
1127		! Update the forcing by the part to be consumed within the present k-layer.
1128		! If fractionOfForcing = 1, then new netMassOut vanishes.
1129	0	netMassOut(i) = netMassOut(i) - dThickness
1130	0	netHeat(i) = netHeat(i) - dTemp
1131	0	netSalt(i) = netSalt(i) - dSalt
1132
1133		! This line accounts for the temperature of the mass exchange
1134	0	dTemp = dTemp + dThicknessT2d(i,k)EnthalpyConst
1135
1136		! Diagnostics of heat content associated with mass fluxes
1137	0	if (.not. associated(fluxes%heat_content_evap)) then
1138	0	if (associated(fluxes%heat_content_massin)) &
1139		fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1140	0	T2d(i,k) * max(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
1141	0	if (associated(fluxes%heat_content_massout)) &
1142		fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1143	0	T2d(i,k) * min(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
1144	0	if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1145	0	T2d(i,k) * dThickness * GV%H_to_RZ
1146		endif
1147
1148		! Update state by the appropriate increment.
1149	0	hOld = h2d(i,k) ! Keep original thickness in hand
1150	0	h2d(i,k) = h2d(i,k) + dThickness ! New thickness
1151
1152	0	if (h2d(i,k) > 0.) then
1153	0	if (calculate_energetics) then
1154		! Calculate the energy required to mix the newly added water over the topmost grid
1155		! cell, assuming that the fluxes of heat and salt and rejected brine are initially
1156		! applied in vanishingly thin layers at the top of the layer before being mixed
1157		! throughout the layer. Note that dThickness is always <= 0 here, and that
1158		! negative cTKE is a deficit that will need to be filled later.
1159		cTKE(i,j,k) = cTKE(i,j,k) - (0.5h2d(i,k)g_Hconv2) * &
1160		((dTemp - dthicknessT2d(i,k)) dSV_dT(i,j,k) + &
1161	0	(dSalt - dthicknesstv%S(i,j,k)) dSV_dS(i,j,k))
1162		endif
1163	0	Ithickness = 1.0/h2d(i,k) ! Inverse of new thickness
1164	0	T2d(i,k) = (hOldT2d(i,k) + dTemp)Ithickness
1165	0	tv%S(i,j,k) = (hOldtv%S(i,j,k) + dSalt + plume_flux)Ithickness
1166	0	elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
1167	0	call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (h<0)')
1168	0	write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j)
1169	0	write(0,*) 'applyBoundaryFluxesInOut(): netT,netS,netH=', &
1170	0	US%C_to_degCnetHeat(i), US%S_to_pptnetSalt(i), netMassInOut(i)
1171	0	write(0,*) 'applyBoundaryFluxesInOut(): dT,dS,dH=', &
1172	0	US%C_to_degCdTemp, US%S_to_pptdSalt, dThickness
1173	0	write(0,*) 'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hOld,h2d(i,k),k
1174		call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
1175	0	"Complete mass loss in column!")
1176		endif
1177
1178		enddo ! k
1179
1180		! Check if trying to apply fluxes over land points
1181	0	elseif ((abs(netHeat(i)) + abs(netSalt(i)) + abs(netMassIn(i)) + abs(netMassOut(i))) > 0.) then
1182
1183	0	if (.not. CS%ignore_fluxes_over_land) then
1184	0	call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (land)')
1185	0	write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j)
1186	0	write(0,*) 'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
1187	0	US%C_to_degCnetHeat(i), US%S_to_pptnetSalt(i), netMassIn(i), netMassOut(i)
1188
1189		call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
1190	0	"Mass loss over land?")
1191		endif
1192
1193		endif
1194
1195		! If anything remains after the k-loop, then we have grounded out, which is a problem.
1196	0	if (netMassIn(i)+netMassOut(i) /= 0.0) then
1197		!$OMP critical
1198	0	numberOfGroundings = numberOfGroundings +1
1199	0	if (numberOfGroundings<=maxGroundings) then
1200	0	iGround(numberOfGroundings) = i ! Record i,j location of event for
1201	0	jGround(numberOfGroundings) = j ! warning message
1202	0	hGrounding(numberOfGroundings) = netMassIn(i)+netMassOut(i)
1203		endif
1204		!$OMP end critical
1205	0	if (CS%id_createdH>0) CS%createdH(i,j) = CS%createdH(i,j) - (netMassIn(i)+netMassOut(i))/dt
1206		endif
1207
1208		enddo ! i
1209
1210		! Step C/ in the application of fluxes
1211		! Heat by the convergence of penetrating SW.
1212		! SW penetrative heating uses the updated thickness from above.
1213
1214		! Save temperature before increment with SW heating
1215		! and initialize CS%penSWflux_diag to zero.
1216	0	if (CS%id_penSW_diag > 0 .or. CS%id_penSWflux_diag > 0) then
1217	0	do k=1,nz ; do i=is,ie
1218	0	CS%penSW_diag(i,j,k) = T2d(i,k)
1219	0	CS%penSWflux_diag(i,j,k) = 0.0
1220		enddo ; enddo
1221	0	k=nz+1 ; do i=is,ie
1222	0	CS%penSWflux_diag(i,j,k) = 0.0
1223		enddo
1224		endif
1225
1226	0	if (calculate_energetics) then
1227		call absorbRemainingSW(G, GV, US, h2d, opacityBand, nsw, optics, j, dt, H_limit_fluxes, &
1228	0	.false., .true., T2d, Pen_SW_bnd, TKE=pen_TKE_2d, dSV_dT=dSV_dT_2d)
1229	0	k = 1 ! For setting break-points.
1230	0	do k=1,nz ; do i=is,ie
1231	0	cTKE(i,j,k) = cTKE(i,j,k) + pen_TKE_2d(i,k)
1232		enddo ; enddo
1233		else
1234		call absorbRemainingSW(G, GV, US, h2d, opacityBand, nsw, optics, j, dt, H_limit_fluxes, &
1235	0	.false., .true., T2d, Pen_SW_bnd)
1236		endif
1237
1238
1239		! Step D/ copy updated thickness and temperature
1240		! 2d slice now back into model state.
1241	0	do k=1,nz ; do i=is,ie
1242	0	h(i,j,k) = h2d(i,k)
1243	0	tv%T(i,j,k) = T2d(i,k)
1244		enddo ; enddo
1245
1246		! Diagnose heating [Q R Z T-1 ~> W m-2] applied to a grid cell from SW penetration
1247		! Also diagnose the penetrative SW heat flux at base of layer.
1248	0	if (CS%id_penSW_diag > 0 .or. CS%id_penSWflux_diag > 0) then
1249
1250		! convergence of SW into a layer
1251	0	do k=1,nz ; do i=is,ie
1252		! Note that the units of penSW_diag change here, from [C ~> degC] to [Q R Z T-1 ~> W m-2].
1253	0	CS%penSW_diag(i,j,k) = (T2d(i,k)-CS%penSW_diag(i,j,k))h(i,j,k) Idt * tv%C_p * GV%H_to_RZ
1254		enddo ; enddo
1255
1256		! Perform a cumulative sum upwards from bottom to
1257		! diagnose penetrative SW flux at base of tracer cell.
1258		! CS%penSWflux_diag(i,j,k=1) is penetrative shortwave at top of ocean.
1259		! CS%penSWflux_diag(i,j,k=kbot+1) is zero, since assume no SW penetrates rock.
1260		! CS%penSWflux_diag = rsdo and CS%penSW_diag = rsdoabsorb
1261		! rsdoabsorb(k) = rsdo(k) - rsdo(k+1), so that rsdo(k) = rsdo(k+1) + rsdoabsorb(k)
1262	0	if (CS%id_penSWflux_diag > 0) then
1263	0	do k=nz,1,-1 ; do i=is,ie
1264	0	CS%penSWflux_diag(i,j,k) = CS%penSW_diag(i,j,k) + CS%penSWflux_diag(i,j,k+1)
1265		enddo ; enddo
1266		endif
1267
1268		endif
1269
1270		! Fill CS%nonpenSW_diag
1271	0	if (CS%id_nonpenSW_diag > 0) then
1272	0	do i=is,ie
1273	0	CS%nonpenSW_diag(i,j) = nonpenSW(i) * Idt * tv%C_p * GV%H_to_RZ
1274		enddo
1275		endif
1276
1277		! BGR: Get buoyancy flux to return for ePBL
1278		! We want the rate, so we use the rate values returned from extractfluxes1d.
1279		! Note that the *dt values could be divided by dt here, but
1280		! 1) Answers will change due to round-off
1281		! 2) Be sure to save their values BEFORE fluxes are used.
1282	0	if (Calculate_Buoyancy) then
1283	0	netPen_rate(:) = 0.0
1284		! Sum over bands and attenuate as a function of depth.
1285		! netPen_rate is the netSW as a function of depth, but only the surface value is used here,
1286		! in which case the values of dt, h, optics and H_limit_fluxes are irrelevant. Consider
1287		! writing a shorter and simpler variant to handle this very limited case.
1288		! Find the vertical distances across layers.
1289		! call thickness_to_dz(h, tv, dz, j, G, GV)
1290		! call sumSWoverBands(G, GV, US, h2d, dz, optics_nbands(optics), optics, j, dt, &
1291		! H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
1292	0	do i=is,ie ; do nb=1,nsw ; netPen_rate(i) = netPen_rate(i) + pen_SW_bnd_rate(nb,i) ; enddo ; enddo
1293
1294		! 1. Adjust netSalt to reflect dilution effect of FW flux
1295		! 2. Add in the SW heating for purposes of calculating the net
1296		! surface buoyancy flux affecting the top layer.
1297		! 3. Convert to a buoyancy flux, excluding penetrating SW heating
1298		! BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
1299	0	if (associated(tv%p_surf)) then ; do i=is,ie ; SurfPressure(i) = tv%p_surf(i,j) ; enddo ; endif
1300
1301	0	if ((.not.GV%Boussinesq) .and. (.not.GV%semi_Boussinesq)) then
1302	0	g_conv = GV%g_Earth_Z_T2 * GV%H_to_RZ
1303
1304		! Specific volume derivatives
1305		call calculate_specific_vol_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, dSpV_dT, dSpV_dS, &
1306	0	tv%eqn_of_state, EOS_domain(G%HI))
1307	0	do i=is,ie
1308		SkinBuoyFlux(i,j) = g_conv * &
1309		(dSpV_dS(i) * ( netSalt_rate(i) - tv%S(i,j,1)*netMassInOut_rate(i)) + &
1310	0	dSpV_dT(i) * ( netHeat_rate(i) + netPen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
1311		enddo
1312		else
1313		! Density derivatives
1314		call calculate_density_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, dRhodT, dRhodS, &
1315	0	tv%eqn_of_state, EOSdom)
1316	0	do i=is,ie
1317		SkinBuoyFlux(i,j) = - GoRho * GV%H_to_Z * &
1318		(dRhodS(i) * ( netSalt_rate(i) - tv%S(i,j,1)*netMassInOut_rate(i)) + &
1319	0	dRhodT(i) * ( netHeat_rate(i) + netPen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
1320		enddo
1321		endif
1322		endif
1323
1324		enddo ! j-loop finish
1325
1326		! Post the diagnostics
1327	0	if (CS%id_createdH > 0) call post_data(CS%id_createdH , CS%createdH , CS%diag)
1328	0	if (CS%id_penSW_diag > 0) call post_data(CS%id_penSW_diag , CS%penSW_diag , CS%diag)
1329	0	if (CS%id_penSWflux_diag > 0) call post_data(CS%id_penSWflux_diag, CS%penSWflux_diag, CS%diag)
1330	0	if (CS%id_nonpenSW_diag > 0) call post_data(CS%id_nonpenSW_diag , CS%nonpenSW_diag , CS%diag)
1331
1332		! The following check will be ignored if ignore_fluxes_over_land = true
1333	0	if ((numberOfGroundings > 0) .and. .not.CS%ignore_fluxes_over_land) then
1334	0	do i = 1, min(numberOfGroundings, maxGroundings)
1335	0	call forcing_SinglePointPrint(fluxes,G,iGround(i),jGround(i),'applyBoundaryFluxesInOut (grounding)')
1336	0	write(mesg(1:45),'(3es15.3)') G%geoLonT( iGround(i), jGround(i) ), &
1337	0	G%geoLatT( iGround(i), jGround(i)), hGrounding(i)*GV%H_to_m
1338		call MOM_error(WARNING, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
1339	0	"Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
1340		enddo
1341
1342	0	if (numberOfGroundings - maxGroundings > 0) then
1343	0	write(mesg, '(I0)') numberOfGroundings - maxGroundings
1344		call MOM_error(WARNING, "MOM_diabatic_aux:F90, applyBoundaryFluxesInOut(): "//&
1345	0	trim(mesg) // " groundings remaining")
1346		endif
1347		endif
1348
1349	0	end subroutine applyBoundaryFluxesInOut
1350
1351		!> This subroutine initializes the parameters and control structure of the diabatic_aux module.
1352	1	subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm, use_ePBL)
1353		type(time_type), target, intent(in) :: Time !< The current model time.
1354		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1355		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1356		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1357		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1358		type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output
1359		type(diabatic_aux_CS), pointer :: CS !< A pointer to the control structure for the
1360		!! diabatic_aux module, which is initialized here.
1361		logical, intent(in) :: useALEalgorithm !< If true, use the ALE algorithm rather
1362		!! than layered mode.
1363		logical, intent(in) :: use_ePBL !< If true, use the implicit energetics planetary
1364		!! boundary layer scheme to determine the diffusivity
1365		!! in the surface boundary layer.
1366
1367		! This "include" declares and sets the variable "version".
1368		# include "version_variable.h"
1369		character(len=40) :: mdl = "MOM_diabatic_aux" ! This module's name.
1370		character(len=200) :: inputdir ! The directory where NetCDF input files
1371		character(len=240) :: chl_filename ! A file from which chl_a concentrations are to be read.
1372		character(len=128) :: chl_file ! Data containing chl_a concentrations. Used
1373		! when var_pen_sw is defined and reading from file.
1374		character(len=32) :: chl_varname ! Name of chl_a variable in chl_file.
1375		logical :: use_temperature ! True if thermodynamics are enabled.
1376		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
1377	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
1378	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
1379
1380	1	if (associated(CS)) then
1381		call MOM_error(WARNING, "diabatic_aux_init called with an "// &
1382	0	"associated control structure.")
1383	0	return
1384		else
1385	1	allocate(CS)
1386		endif
1387
1388	1	CS%diag => diag
1389	1	CS%Time => Time
1390
1391		! Set default, read and log parameters
1392		call log_version(param_file, mdl, version, &
1393	1	"The following parameters are used for auxiliary diabatic processes.")
1394
1395		call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
1396	1	"If true, temperature and salinity are used as state variables.", default=.true.)
1397
1398		call get_param(param_file, mdl, "RECLAIM_FRAZIL", CS%reclaim_frazil, &
1399		"If true, try to use any frazil heat deficit to cool any "//&
1400		"overlying layers down to the freezing point, thereby "//&
1401		"avoiding the creation of thin ice when the SST is above "//&
1402	1	"the freezing point.", default=.true., do_not_log=.not.use_temperature)
1403		call get_param(param_file, mdl, "SALT_EXTRACTION_LIMIT", CS%dSalt_frac_max, &
1404		"An upper limit on the fraction of the salt in a layer that can be lost to the "//&
1405		"net surface salt fluxes within a timestep.", &
1406	1	units="nondim", default=0.9999, do_not_log=.not.use_temperature)
1407	1	CS%dSalt_frac_max = max(min(CS%dSalt_frac_max, 1.0), 0.0)
1408		call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", CS%pressure_dependent_frazil, &
1409		"If true, use a pressure dependent freezing temperature "//&
1410		"when making frazil. The default is false, which will be "//&
1411		"faster but is inappropriate with ice-shelf cavities.", &
1412	1	default=.false., do_not_log=.not.use_temperature)
1413
1414	1	if (use_ePBL) then
1415		call get_param(param_file, mdl, "IGNORE_FLUXES_OVER_LAND", CS%ignore_fluxes_over_land,&
1416		"If true, the model does not check if fluxes are being applied "//&
1417		"over land points. This is needed when the ocean is coupled "//&
1418		"with ice shelves and sea ice, since the sea ice mask needs to "//&
1419		"be different than the ocean mask to avoid sea ice formation "//&
1420	0	"under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
1421		call get_param(param_file, mdl, "DO_RIVERMIX", CS%do_rivermix, &
1422		"If true, apply additional mixing wherever there is "//&
1423		"runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
1424	0	"if the ocean is that deep.", default=.false.)
1425	0	if (CS%do_rivermix) &
1426		call get_param(param_file, mdl, "RIVERMIX_DEPTH", CS%rivermix_depth, &
1427		"The depth to which rivers are mixed if DO_RIVERMIX is "//&
1428	0	"defined.", units="m", default=0.0, scale=US%m_to_Z)
1429		else
1430	1	CS%do_rivermix = .false. ; CS%rivermix_depth = 0.0 ; CS%ignore_fluxes_over_land = .false.
1431		endif
1432
1433	1	if (GV%nkml == 0) then
1434		call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", CS%use_river_heat_content, &
1435		"If true, use the fluxes%runoff_Hflx field to set the "//&
1436		"heat carried by runoff, instead of using SSTCPliq_runoff.", &
1437	0	default=.false., do_not_log=.not.use_temperature)
1438		call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", CS%use_calving_heat_content, &
1439		"If true, use the fluxes%calving_Hflx field to set the "//&
1440		"heat carried by runoff, instead of using SSTCPfroz_runoff.", &
1441	0	default=.false., do_not_log=.not.use_temperature)
1442		else
1443	1	CS%use_river_heat_content = .false.
1444	1	CS%use_calving_heat_content = .false.
1445		endif
1446
1447		call get_param(param_file, mdl, "DO_BRINE_PLUME", CS%do_brine_plume, &
1448		"If true, use a brine plume parameterization from "//&
1449	1	"Nguyen et al., 2009.", default=.false.)
1450		call get_param(param_file, mdl, "BRINE_PLUME_EXPONENT", CS%brine_plume_n, &
1451		"If using the brine plume parameterization, set the integer exponent.", &
1452	1	default=5, do_not_log=.not.CS%do_brine_plume)
1453		call get_param(param_file, mdl, "BRINE_PLUME_FRACTION", CS%plume_strength, &
1454		"Fraction of the available brine to mix down using the brine plume parameterization.", &
1455	1	units="nondim", default=1.0, do_not_log=.not.CS%do_brine_plume)
1456
1457	1	if (useALEalgorithm) then
1458		CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
1459		Time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
1460	0	"m s-1", conversion=GV%H_to_m*US%s_to_T)
1461	0	if (CS%id_createdH>0) allocate(CS%createdH(isd:ied,jsd:jed))
1462
1463		! diagnostic for heating of a grid cell from convergence of SW heat into the cell
1464		CS%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb', &
1465		diag%axesTL, Time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
1466		'W m-2', conversion=US%QRZ_T_to_W_m2, &
1467	0	standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer', v_extensive=.true.)
1468
1469		! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)
1470		! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).
1471		CS%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo', &
1472		diag%axesTi, Time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
1473	0	'W m-2', conversion=US%QRZ_T_to_W_m2, standard_name='downwelling_shortwave_flux_in_sea_water')
1474
1475		! need both arrays for the SW diagnostics (one for flux, one for convergence)
1476	0	if (CS%id_penSW_diag>0 .or. CS%id_penSWflux_diag>0) then
1477	0	allocate(CS%penSW_diag(isd:ied,jsd:jed,nz), source=0.0)
1478	0	allocate(CS%penSWflux_diag(isd:ied,jsd:jed,nz+1), source=0.0)
1479		endif
1480
1481		! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)
1482		CS%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW', &
1483		diag%axesT1, Time, &
1484		'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
1485		'W m-2', conversion=US%QRZ_T_to_W_m2, &
1486	0	standard_name='nondownwelling_shortwave_flux_in_sea_water')
1487	0	if (CS%id_nonpenSW_diag > 0) then
1488	0	allocate(CS%nonpenSW_diag(isd:ied,jsd:jed), source=0.0)
1489		endif
1490		endif
1491
1492	1	if (use_temperature) then
1493		call get_param(param_file, mdl, "VAR_PEN_SW", CS%var_pen_sw, &
1494		"If true, use one of the CHL_A schemes specified by "//&
1495		"OPACITY_SCHEME to determine the e-folding depth of "//&
1496	1	"incoming short wave radiation.", default=.false.)
1497	1	if (CS%var_pen_sw) then
1498
1499		call get_param(param_file, mdl, "CHL_FROM_FILE", CS%chl_from_file, &
1500	0	"If true, chl_a is read from a file.", default=.true.)
1501	0	if (CS%chl_from_file) then
1502	0	call time_interp_external_init()
1503
1504	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1505		call get_param(param_file, mdl, "CHL_FILE", chl_file, &
1506		"CHL_FILE is the file containing chl_a concentrations in "//&
1507		"the variable CHL_A. It is used when VAR_PEN_SW and "//&
1508	0	"CHL_FROM_FILE are true.", fail_if_missing=.true.)
1509	0	chl_filename = trim(slasher(inputdir))//trim(chl_file)
1510	0	call log_param(param_file, mdl, "INPUTDIR/CHL_FILE", chl_filename)
1511		call get_param(param_file, mdl, "CHL_VARNAME", chl_varname, &
1512	0	"Name of CHL_A variable in CHL_FILE.", default='CHL_A')
1513	0	if (modulo(G%Domain%turns, 4) /= 0) then
1514	0	CS%sbc_chl = init_external_field(chl_filename, trim(chl_varname), MOM_domain=G%Domain%domain_in)
1515		else
1516	0	CS%sbc_chl = init_external_field(chl_filename, trim(chl_varname), MOM_domain=G%Domain)
1517		endif
1518		endif
1519
1520		CS%id_chl = register_diag_field('ocean_model', 'Chl_opac', diag%axesT1, Time, &
1521	0	'Surface chlorophyll A concentration used to find opacity', 'mg m-3')
1522		endif
1523		endif
1524
1525	1	id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=CLOCK_ROUTINE)
1526	1	id_clock_frazil = cpu_clock_id('(Ocean frazil)', grain=CLOCK_ROUTINE)
1527
1528		end subroutine diabatic_aux_init
1529
1530		!> This subroutine initializes the control structure and any related memory
1531		!! for the diabatic_aux module.
1532	1	subroutine diabatic_aux_end(CS)
1533		type(diabatic_aux_CS), pointer :: CS !< The control structure returned by a previous
1534		!! call to diabatic_aux_init; it is deallocated here.
1535
1536	1	if (.not.associated(CS)) return
1537
1538	1	if (CS%id_createdH >0) deallocate(CS%createdH)
1539	1	if (CS%id_penSW_diag >0) deallocate(CS%penSW_diag)
1540	1	if (CS%id_penSWflux_diag >0) deallocate(CS%penSWflux_diag)
1541	1	if (CS%id_nonpenSW_diag >0) deallocate(CS%nonpenSW_diag)
1542
1543	1	if (associated(CS)) deallocate(CS)
1544
1545		end subroutine diabatic_aux_end
1546
1547		!> \namespace mom_diabatic_aux
1548		!!
1549		!! This module contains subroutines that apply various diabatic processes. Usually these
1550		!! subroutines are called from the MOM_diabatic module. All of these routines use appropriate
1551		!! limiters or logic to work properly with arbitrary layer thicknesses (including massless layers)
1552		!! and an arbitrarily large timestep.
1553		!!
1554		!! The subroutine make_frazil facilitates the formation of frazil ice when the ocean water
1555		!! drops below the in situ freezing point by heating the water to the freezing point and
1556		!! accumulating the required heat for exchange with the sea-ice module.
1557		!!
1558		!! The subroutine adjust_salt adds salt as necessary to keep the salinity above a
1559		!! specified minimum value, and keeps track of the cumulative additions. If the minimum
1560		!! salinity is the natural value of 0, this routine should never do anything.
1561		!!
1562		!! The subroutine differential_diffuse_T_S solves a pair of tridiagonal equations for
1563		!! the diffusion of temperatures and salinities with differing diffusivities.
1564		!!
1565		!! The subroutine triDiagTS solves a tridiagonal equations for the evolution of temperatures
1566		!! and salinities due to net entrainment by layers and a diffusion with the same diffusivity.
1567		!!
1568		!! The subroutine triDiagTS_Eulerian solves a tridiagonal equations for the evolution of
1569		!! temperatures and salinities due to diffusion with the same diffusivity, but no net entrainment.
1570		!!
1571		!! The subroutine find_uv_at_h interpolates velocities to thickness points, optionally also
1572		!! using tridiagonal equations to solve for the impacts of net entrainment or mixing of
1573		!! momentum between layers.
1574		!!
1575		!! The subroutine set_pen_shortwave determines the optical properties of the water column and
1576		!! the net shortwave fluxes, and stores them in the optics type, working via calls to set_opacity.
1577		!!
1578		!! The subroutine applyBoundaryFluxesInOut updates the layer thicknesses, temperatures and
1579		!! salinities due to the application of the surface forcing. It may also calculate the implied
1580		!! turbulent kinetic energy requirements for this forcing to be mixed over the model's finite
1581		!! vertical resolution in the surface layers.
1582	0	end module MOM_diabatic_aux