MOM_diabatic_aux.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 functions for some diabatic processes such as frazil, brine rejection,
6!! tendency due to surface flux divergence.
8
9use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
10use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
11use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
12use mom_diag_mediator, only : diag_ctrl, time_type
13use mom_eos, only : calculate_density, calculate_tfreeze, eos_domain
14use mom_eos, only : calculate_specific_vol_derivs, calculate_density_derivs
15use mom_error_handler, only : mom_error, fatal, warning, calltree_showquery
17use mom_file_parser, only : get_param, log_param, log_version, param_file_type
19use mom_grid, only : ocean_grid_type
20use mom_interface_heights, only : thickness_to_dz
21use mom_interpolate, only : init_external_field, time_interp_external, time_interp_external_init
22use mom_interpolate, only : external_field
23use mom_io, only : slasher
30
31implicit none ; private
32
33#include <MOM_memory.h>
34
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
45type, 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
100end type diabatic_aux_cs
101
102!>@{ CPU time clock IDs
103integer :: id_clock_uv_at_h, id_clock_frazil
104!>@}
105
106contains
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.
112subroutine 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 fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].
128 t_freeze, & ! The freezing potential temperature at the current salinity [C ~> degC].
129 ps ! Surface pressure [R L2 T-2 ~> Pa]
130 real, dimension(SZI_(G),SZK_(GV)) :: &
131 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
139 if (present(halo)) then
140 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
141 endif
142
143 call cpu_clock_begin(id_clock_frazil)
144
145 if (.not.cs%pressure_dependent_frazil) then
146 do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo
147 else
148 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 do j=js,je
153 ps(:) = 0.0
154 if (PRESENT(p_surf)) then ; do i=is,ie
155 ps(i) = p_surf(i,j)
156 enddo ; endif
157
158 do i=is,ie ; fraz_col(i) = 0.0 ; enddo
159
160 if (cs%pressure_dependent_frazil) then
161 do i=is,ie
162 pressure(i,1) = ps(i) + (0.5*h_to_rl2_t2)*h(i,j,1)
163 enddo
164 do k=2,nz ; do i=is,ie
165 pressure(i,k) = pressure(i,k-1) + (0.5*h_to_rl2_t2) * (h(i,j,k) + h(i,j,k-1))
166 enddo ; enddo
167 endif
168
169 if (cs%reclaim_frazil) then
170 t_fr_set = .false.
171 do i=is,ie ; if (tv%frazil(i,j) > 0.0) then
172 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 tv%eqn_of_state)
175 t_fr_set = .true.
176 endif
177
178 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 hc = (tv%C_p*gv%H_to_RZ) * h(i,j,1)
182 if (tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i)) <= 0.0) then
183 tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j) / hc
184 tv%frazil(i,j) = 0.0
185 else
186 tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i))
187 tv%T(i,j,1) = t_freeze(i)
188 endif
189 endif
190 endif ; enddo
191 endif
192
193 do k=nz,1,-1
194 t_fr_set = .false.
195 do i=is,ie
196 if ((g%mask2dT(i,j) > 0.0) .and. &
197 ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then
198 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 tv%eqn_of_state)
201 t_fr_set = .true.
202 endif
203
204 hc = (tv%C_p*gv%H_to_RZ) * h(i,j,k)
205 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 if (tv%T(i,j,k) < t_freeze(i)) then
208 fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
209 tv%T(i,j,k) = t_freeze(i)
210 endif
211 elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < t_freeze(i))) then
212 if (fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k)) < 0.0) then
213 tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc
214 fraz_col(i) = 0.0
215 else
216 fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
217 tv%T(i,j,k) = t_freeze(i)
218 endif
219 endif
220 endif
221 enddo
222 enddo
223 do i=is,ie
224 tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i)
225 enddo
226 enddo
227
228 tv%frazil_was_reset = .false.
229
230 call cpu_clock_end(id_clock_frazil)
231
232end subroutine make_frazil
233
234!> This subroutine applies double diffusion to T & S, assuming no diapycnal mass
235!! fluxes, using a simple tridiagonal solver.
236subroutine 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 b1_t, b1_s, & ! Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].
260 d1_t, d1_s ! Variables used by the tridiagonal solvers [nondim].
261 real, dimension(SZI_(G),SZK_(GV)) :: &
262 dz, & ! Height change across layers [Z ~> m]
263 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 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
278 h_neglect = gv%H_subroundoff
279 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 do j=js,je
283
284 ! Find the vertical distances across layers.
285 call thickness_to_dz(h, tv, dz, j, g, gv)
286
287 do i=is,ie
288 i_dz_int = 1.0 / (0.5 * (dz(i,1) + dz(i,2)) + dz_neglect)
289 mix_t(i,2) = (dt * kd_t(i,j,2)) * i_dz_int
290 mix_s(i,2) = (dt * kd_s(i,j,2)) * i_dz_int
291
292 h_tr = h(i,j,1) + h_neglect
293 b1_t(i) = 1.0 / (h_tr + mix_t(i,2))
294 b1_s(i) = 1.0 / (h_tr + mix_s(i,2))
295 d1_t(i) = h_tr * b1_t(i)
296 d1_s(i) = h_tr * b1_s(i)
297 t(i,j,1) = (b1_t(i)*h_tr)*t(i,j,1)
298 s(i,j,1) = (b1_s(i)*h_tr)*s(i,j,1)
299 enddo
300 do k=2,nz-1 ; do i=is,ie
301 ! Calculate the mixing across the interface below this layer.
302 i_dz_int = 1.0 / (0.5 * (dz(i,k) + dz(i,k+1)) + dz_neglect)
303 mix_t(i,k+1) = ((dt * kd_t(i,j,k+1))) * i_dz_int
304 mix_s(i,k+1) = ((dt * kd_s(i,j,k+1))) * i_dz_int
305
306 c1_t(i,k) = mix_t(i,k) * b1_t(i)
307 c1_s(i,k) = mix_s(i,k) * b1_s(i)
308
309 h_tr = h(i,j,k) + h_neglect
310 b_denom_t = h_tr + d1_t(i)*mix_t(i,k)
311 b_denom_s = h_tr + d1_s(i)*mix_s(i,k)
312 b1_t(i) = 1.0 / (b_denom_t + mix_t(i,k+1))
313 b1_s(i) = 1.0 / (b_denom_s + mix_s(i,k+1))
314 d1_t(i) = b_denom_t * b1_t(i)
315 d1_s(i) = b_denom_s * b1_s(i)
316
317 t(i,j,k) = b1_t(i) * (h_tr*t(i,j,k) + mix_t(i,k)*t(i,j,k-1))
318 s(i,j,k) = b1_s(i) * (h_tr*s(i,j,k) + mix_s(i,k)*s(i,j,k-1))
319 enddo ; enddo
320 do i=is,ie
321 c1_t(i,nz) = mix_t(i,nz) * b1_t(i)
322 c1_s(i,nz) = mix_s(i,nz) * b1_s(i)
323
324 h_tr = h(i,j,nz) + h_neglect
325 b1_t(i) = 1.0 / (h_tr + d1_t(i)*mix_t(i,nz))
326 b1_s(i) = 1.0 / (h_tr + d1_s(i)*mix_s(i,nz))
327
328 t(i,j,nz) = b1_t(i) * (h_tr*t(i,j,nz) + mix_t(i,nz)*t(i,j,nz-1))
329 s(i,j,nz) = b1_s(i) * (h_tr*s(i,j,nz) + mix_s(i,nz)*s(i,j,nz-1))
330 enddo
331 do k=nz-1,1,-1 ; do i=is,ie
332 t(i,j,k) = t(i,j,k) + c1_t(i,k+1)*t(i,j,k+1)
333 s(i,j,k) = s(i,j,k) + c1_s(i,k+1)*s(i,j,k+1)
334 enddo ; enddo
335 enddo
336end 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.
341subroutine 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 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 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 s_min = tv%min_salinity
362
363 salt_add_col(:,:) = 0.0
364
365 !$OMP parallel do default(shared) private(mc)
366 do j=js,je
367 do k=nz,1,-1 ; do i=is,ie
368 if ( (g%mask2dT(i,j) > 0.0) .and. &
369 ((tv%S(i,j,k) < s_min) .or. (salt_add_col(i,j) > 0.0)) ) then
370 mc = gv%H_to_RZ * h(i,j,k)
371 if (h(i,j,k) <= 10.0*gv%Angstrom_H) then
372 ! Very thin layers should not be adjusted by the salt flux
373 if (tv%S(i,j,k) < s_min) then
374 salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
375 tv%S(i,j,k) = s_min
376 endif
377 elseif (salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k)) <= 0.0) then
378 tv%S(i,j,k) = tv%S(i,j,k) - salt_add_col(i,j) / mc
379 salt_add_col(i,j) = 0.0
380 else
381 salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
382 tv%S(i,j,k) = s_min
383 endif
384 endif
385 enddo ; enddo
386 do i=is,ie
387 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
392end 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.
396subroutine 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 real :: b1(szib_(g)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
414 real :: d1(szib_(g)) ! A variable used by the tridiagonal solver [nondim].
415 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 do j=js,je
421 do i=is,ie
422 h_tr = hold(i,j,1) + gv%H_subroundoff
423 b1(i) = 1.0 / (h_tr + eb(i,j,1))
424 d1(i) = h_tr * b1(i)
425 t(i,j,1) = (b1(i)*h_tr)*t(i,j,1)
426 s(i,j,1) = (b1(i)*h_tr)*s(i,j,1)
427 enddo
428 do k=2,gv%ke ; do i=is,ie
429 c1(i,k) = eb(i,j,k-1) * b1(i)
430 h_tr = hold(i,j,k) + gv%H_subroundoff
431 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
432 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
433 d1(i) = b_denom_1 * b1(i)
434 t(i,j,k) = b1(i) * (h_tr*t(i,j,k) + ea(i,j,k)*t(i,j,k-1))
435 s(i,j,k) = b1(i) * (h_tr*s(i,j,k) + ea(i,j,k)*s(i,j,k-1))
436 enddo ; enddo
437 do k=gv%ke-1,1,-1 ; do i=is,ie
438 t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1)
439 s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1)
440 enddo ; enddo
441 enddo
442end 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.
446subroutine 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 real :: b1(szib_(g)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
462 real :: d1(szib_(g)) ! A variable used by the tridiagonal solver [nondim].
463 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 do j=js,je
469 do i=is,ie
470 h_tr = hold(i,j,1) + gv%H_subroundoff
471 b1(i) = 1.0 / (h_tr + ent(i,j,2))
472 d1(i) = h_tr * b1(i)
473 t(i,j,1) = (b1(i)*h_tr)*t(i,j,1)
474 s(i,j,1) = (b1(i)*h_tr)*s(i,j,1)
475 enddo
476 do k=2,gv%ke ; do i=is,ie
477 c1(i,k) = ent(i,j,k) * b1(i)
478 h_tr = hold(i,j,k) + gv%H_subroundoff
479 b_denom_1 = h_tr + d1(i)*ent(i,j,k)
480 b1(i) = 1.0 / (b_denom_1 + ent(i,j,k+1))
481 d1(i) = b_denom_1 * b1(i)
482 t(i,j,k) = b1(i) * (h_tr*t(i,j,k) + ent(i,j,k)*t(i,j,k-1))
483 s(i,j,k) = b1(i) * (h_tr*s(i,j,k) + ent(i,j,k)*s(i,j,k-1))
484 enddo ; enddo
485 do k=gv%ke-1,1,-1 ; do i=is,ie
486 t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1)
487 s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1)
488 enddo ; enddo
489 enddo
490end 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.
495subroutine 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 real :: b1(szi_(g)) ! A thickness used in the tridiagonal solver [H ~> m or kg m-2]
526 real :: c1(szi_(g),szk_(gv)) ! A variable used in the tridiagonal solver [nondim]
527 real :: d1(szi_(g)) ! The complement of c1 [nondim]
528 ! Fractional weights of the neighboring velocity points, ~1/2 in the open ocean.
529 real :: a_n(szi_(g)), a_s(szi_(g)) ! Fractional weights of the neighboring velocity points [nondim]
530 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
536 call cpu_clock_begin(id_clock_uv_at_h)
537 h_neglect = gv%H_subroundoff
538
539 mix_vertically = present(ea)
540 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 "in call to find_uv_at_h.")
543 zero_mixing = .false. ; if (present(zero_mix)) zero_mixing = zero_mix
544 if (zero_mixing) mix_vertically = .false.
545 !$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,zero_mixing,h, &
546 !$OMP h_neglect,ea,eb,u_h,u,v_h,v,nz) &
547 !$OMP private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1)
548 do j=js,je
549 do i=is,ie
550 sum_area = g%areaCu(i-1,j) + g%areaCu(i,j)
551 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.5*sum_area*G%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 idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
559 a_w(i) = g%areaCu(i-1,j) * idenom
560 a_e(i) = g%areaCu(i,j) * idenom
561 else
562 a_w(i) = 0.0 ; a_e(i) = 0.0
563 endif
564
565 sum_area = g%areaCv(i,j-1) + g%areaCv(i,j)
566 if (sum_area > 0.0) then
567 idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
568 a_s(i) = g%areaCv(i,j-1) * idenom
569 a_n(i) = g%areaCv(i,j) * idenom
570 else
571 a_s(i) = 0.0 ; a_n(i) = 0.0
572 endif
573 enddo
574
575 if (mix_vertically) then
576 do i=is,ie
577 b_denom_1 = h(i,j,1) + h_neglect
578 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
579 d1(i) = b_denom_1 * b1(i)
580 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 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 do k=2,nz ; do i=is,ie
584 c1(i,k) = eb(i,j,k-1) * b1(i)
585 b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
586 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
587 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 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 ea(i,j,k)*v_h(i,j,k-1))*b1(i)
592 enddo ; enddo
593 do k=nz-1,1,-1 ; do i=is,ie
594 u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
595 v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
596 enddo ; enddo
597 elseif (zero_mixing) then
598 do i=is,ie
599 b1(i) = 1.0 / (h(i,j,1) + h_neglect)
600 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 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 do k=2,nz ; do i=is,ie
604 b1(i) = 1.0 / (h(i,j,k) + h_neglect)
605 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 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 ; enddo
608 else
609 do k=1,nz ; do i=is,ie
610 u_h(i,j,k) = (a_e(i)*u(i,j,k)) + (a_w(i)*u(i-1,j,k))
611 v_h(i,j,k) = (a_n(i)*v(i,j,k)) + (a_s(i)*v(i,j-1,k))
612 enddo ; enddo
613 endif
614 enddo
615
616 call cpu_clock_end(id_clock_uv_at_h)
617end subroutine find_uv_at_h
618
619!> Estimate the optical properties of the water column and determine the penetrating shortwave
620!! radiation by band, extracting the relevant information from the fluxes type and storing it
621!! in the optics type for later application. This routine is effectively a wrapper for
622!! set_opacity with added error handling and diagnostics.
623subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity, tracer_flow_CSp)
624 type(optics_type), pointer :: optics !< An optics structure that has will contain
625 !! information about shortwave fluxes and absorption.
626 type(forcing), intent(inout) :: fluxes !< points to forcing fields
627 !! unused fields have NULL pointers
628 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
629 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
630 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
631 type(diabatic_aux_cs), pointer :: cs !< Control structure for diabatic_aux
632 type(opacity_cs) :: opacity !< The control structure for the opacity module.
633 type(tracer_flow_control_cs), pointer :: tracer_flow_csp !< A pointer to the control structure
634 !! organizing the tracer modules.
635
636 ! Local variables
637 real, dimension(SZI_(G),SZJ_(G)) :: chl_2d !< Vertically uniform chlorophyll-A concentrations [mg m-3]
638 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: chl_3d !< The chlorophyll-A concentrations of each layer [mg m-3]
639 character(len=128) :: mesg
640 integer :: i, j, is, ie, js, je
641 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
642
643 if (.not.associated(optics)) return
644
645 if (cs%var_pen_sw) then
646 if (cs%chl_from_file) then
647 ! Only the 2-d surface chlorophyll can be read in from a file. The
648 ! same value is assumed for all layers.
649 call time_interp_external(cs%sbc_chl, cs%Time, chl_2d, turns=g%HI%turns)
650 do j=js,je ; do i=is,ie
651 if ((g%mask2dT(i,j) > 0.0) .and. (chl_2d(i,j) < 0.0)) then
652 write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
653 & I0,", ",I0," lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
654 chl_2d(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
655 call mom_error(fatal, "MOM_diabatic_aux set_pen_shortwave: "//trim(mesg))
656 endif
657 enddo ; enddo
658
659 if (cs%id_chl > 0) call post_data(cs%id_chl, chl_2d, cs%diag)
660
661 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
662 fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity, chl_2d=chl_2d)
663 else
664 if (.not.associated(tracer_flow_csp)) call mom_error(fatal, &
665 "The tracer flow control structure must be associated when the model sets "//&
666 "the chlorophyll internally in set_pen_shortwave.")
667 call get_chl_from_model(chl_3d, g, gv, tracer_flow_csp)
668
669 if (cs%id_chl > 0) call post_data(cs%id_chl, chl_3d(:,:,1), cs%diag)
670
671 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
672 fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity, chl_3d=chl_3d)
673 endif
674 else
675 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
676 fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity)
677 endif
678
679end subroutine set_pen_shortwave
680
681!> Update the thickness, temperature, and salinity due to thermodynamic
682!! boundary forcing (contained in fluxes type) applied to h, tv%T and tv%S,
683!! and calculate the TKE implications of this heating.
684subroutine applyboundaryfluxesinout(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, &
685 aggregate_FW_forcing, evap_CFL_limit, &
686 minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
687 SkinBuoyFlux, MLD_h)
688 type(diabatic_aux_cs), pointer :: cs !< Control structure for diabatic_aux
689 type(ocean_grid_type), intent(in) :: g !< Grid structure
690 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
691 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
692 real, intent(in) :: dt !< Time-step over which forcing is applied [T ~> s]
693 type(forcing), intent(inout) :: fluxes !< Surface fluxes container
694 type(optics_type), pointer :: optics !< Optical properties container
695 integer, intent(in) :: nsw !< The number of frequency bands of penetrating
696 !! shortwave radiation
697 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
698 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
699 type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
700 !! available thermodynamic fields.
701 logical, intent(in) :: aggregate_fw_forcing !< If False, treat in/out fluxes separately.
702 real, intent(in) :: evap_cfl_limit !< The largest fraction of a layer that
703 !! can be evaporated in one time-step [nondim].
704 real, intent(in) :: minimum_forcing_depth !< The smallest depth over which
705 !! heat and freshwater fluxes is applied [H ~> m or kg m-2].
706 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
707 optional, intent(out) :: ctke !< Turbulent kinetic energy requirement to mix
708 !! forcing through each layer [R Z3 T-2 ~> J m-2]
709 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
710 optional, intent(out) :: dsv_dt !< Partial derivative of specific volume with
711 !! potential temperature [R-1 C-1 ~> m3 kg-1 degC-1].
712 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
713 optional, intent(out) :: dsv_ds !< Partial derivative of specific volume with
714 !! salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
715 real, dimension(SZI_(G),SZJ_(G)), &
716 optional, intent(out) :: skinbuoyflux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].
717 real, dimension(:,:), &
718 optional, pointer :: mld_h !< Mixed layer thickness for brine plumes [H ~> m or kg m-2]
719
720 ! Local variables
721 integer, parameter :: maxgroundings = 5
722 integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
723 real :: h_limit_fluxes ! Surface fluxes are scaled down fluxes when the total depth of the ocean
724 ! drops below this value [H ~> m or kg m-2]
725 real :: iforcingdepthscale ! The inverse of the layer thickness below which mass losses are
726 ! shifted to the next deeper layer [H ~> m or kg m-2]
727 real :: idt ! The inverse of the timestep [T-1 ~> s-1]
728 real :: dthickness ! The change in layer thickness [H ~> m or kg m-2]
729 real :: dtemp ! The integrated change in layer temperature [C H ~> degC m or degC kg m-2]
730 real :: dsalt ! The integrated change in layer salinity [S H ~> ppt m or ppt kg m-2]
731 real :: fractionofforcing ! THe fraction of the remaining forcing applied to a layer [nondim]
732 real :: hold ! The original thickness of a layer [H ~> m or kg m-2]
733 real :: ithickness ! The inverse of the new layer thickness [H-1 ~> m-1 or m2 kg-1]
734 real :: rivermixconst ! A constant used in implementing river mixing [R Z2 T-1 ~> Pa s].
735 real :: enthalpyconst ! A constant used to control the enthalpy calculation [nondim]
736 ! By default EnthalpyConst = 1.0. If fluxes%heat_content_evap
737 ! is associated enthalpy is provided via coupler and EnthalpyConst = 0.0.
738 real, dimension(SZI_(G)) :: &
739 d_pres, & ! pressure change across a layer [R L2 T-2 ~> Pa]
740 p_lay, & ! average pressure in a layer [R L2 T-2 ~> Pa]
741 pres, & ! pressure at an interface [R L2 T-2 ~> Pa]
742 netmassinout, & ! surface water fluxes [H ~> m or kg m-2] over time step
743 netmassin, & ! mass entering ocean surface [H ~> m or kg m-2] over a time step
744 netmassout, & ! mass leaving ocean surface [H ~> m or kg m-2] over a time step
745 netheat, & ! heat via surface fluxes excluding Pen_SW_bnd and netMassOut
746 ! [C H ~> degC m or degC kg m-2]
747 netsalt, & ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
748 ! [S H ~> ppt m or ppt kg m-2]
749 nonpensw, & ! non-downwelling SW, which is absorbed at ocean surface
750 ! [C H ~> degC m or degC kg m-2]
751 surfpressure, & ! Surface pressure (approximated as 0.0) [R L2 T-2 ~> Pa]
752 drhodt, & ! change in density per change in temperature [R C-1 ~> kg m-3 degC-1]
753 drhods, & ! change in density per change in salinity [R S-1 ~> kg m-3 ppt-1]
754 dspv_dt, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
755 dspv_ds, & ! Partial derivative of specific volume with to salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
756 netheat_rate, & ! netheat but for dt=1 [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
757 netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate)
758 ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
759 netmassinout_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
760 mixing_depth, & ! The mixing depth for brine plumes [H ~> m or kg m-2]
761 total_h ! Total thickness of the water column [H ~> m or kg m-2]
762 real, dimension(SZI_(G), SZK_(GV)) :: &
763 h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
764 ! dz, & ! Layer thicknesses in depth units [Z ~> m]
765 t2d, & ! A 2-d copy of the layer temperatures [C ~> degC]
766 pen_tke_2d, & ! The TKE required to homogenize the heating by shortwave radiation within
767 ! a layer [R Z3 T-2 ~> J m-2]
768 dsv_dt_2d ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
769 real, dimension(SZI_(G)) :: &
770 netpen_rate ! The surface penetrative shortwave heating rate summed over all bands
771 ! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
772 real, dimension(max(nsw,1),SZI_(G)) :: &
773 pen_sw_bnd, & ! The penetrative shortwave heating integrated over a timestep by band
774 ! [C H ~> degC m or degC kg m-2]
775 pen_sw_bnd_rate ! The penetrative shortwave heating rate by band
776 ! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
777 real, dimension(max(nsw,1),SZI_(G),SZK_(GV)) :: &
778 opacityband ! The opacity (inverse of the exponential absorption length) of each frequency
779 ! band of shortwave radiation in each layer [H-1 ~> m-1 or m2 kg-1]
780 real, dimension(maxGroundings) :: hgrounding ! Thickness added by each grounding event [H ~> m or kg m-2]
781 real :: temp_in ! The initial temperature of a layer [C ~> degC]
782 real :: salin_in ! The initial salinity of a layer [S ~> ppt]
783 real :: g_hconv2 ! A conversion factor for use in the TKE calculation
784 ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
785 real :: gorho ! g_Earth times a unit conversion factor divided by density
786 ! [Z T-2 R-1 ~> m4 s-2 kg-1]
787 real :: g_conv ! The gravitational acceleration times the conversion factors from non-Boussinesq
788 ! thickness units to mass per units area [R Z2 H-1 T-2 ~> kg m-2 s-2 or m s-2]
789 logical :: calculate_energetics ! If true, calculate the energy required to mix the newly added
790 ! water over the topmost grid cell, assuming that the fluxes of heat and salt
791 ! and rejected brine are initially applied in vanishingly thin layers at the
792 ! top of the layer before being mixed throughout the layer.
793 logical :: calculate_buoyancy ! If true, calculate the surface buoyancy flux.
794 real :: dk(szi_(g)) ! Depth of the layer center in thickness units [H ~> m or kg m-2]
795 real :: a_brine(szi_(g)) ! Constant [H-(n+1) ~> m-(n+1) or m(2n+2) kg-(n+1)].
796 real :: fraction_left_brine ! Fraction of the brine that has not been applied yet [nondim]
797 real :: plume_fraction ! Fraction of the brine that is applied to a layer [nondim]
798 real :: plume_flux ! Brine flux to move downwards [S H ~> ppt m or ppt kg m-2]
799 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
800 integer :: i, j, is, ie, js, je, k, nz, nb
801 character(len=45) :: mesg
802
803 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
804
805 idt = 1.0 / dt
806 plume_flux = 0.0
807
808 calculate_energetics = (present(ctke) .and. present(dsv_dt) .and. present(dsv_ds))
809 calculate_buoyancy = present(skinbuoyflux)
810 if (calculate_buoyancy) skinbuoyflux(:,:) = 0.0
811 if (present(ctke)) ctke(:,:,:) = 0.0
812 g_hconv2 = (gv%g_Earth_Z_T2 * gv%H_to_RZ) * gv%H_to_RZ
813 eosdom(:) = eos_domain(g%HI)
814
815 ! Only apply forcing if fluxes%sw is associated.
816 if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return
817
818 enthalpyconst = 1.0
819 if (associated(fluxes%heat_content_evap)) enthalpyconst = 0.0
820
821 if (calculate_buoyancy) then
822 surfpressure(:) = 0.0
823 gorho = gv%g_Earth_Z_T2 / gv%Rho0
824 endif
825
826 if (cs%do_brine_plume .and. .not.present(mld_h)) then
827 call mom_error(fatal, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
828 "Brine plume parameterization requires a mixed-layer depth argument,\n"//&
829 "currently coming from the energetic PBL scheme.")
830 endif
831 if (cs%do_brine_plume .and. .not.associated(mld_h)) then
832 call mom_error(fatal, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
833 "Brine plume parameterization requires an associated mixed-layer depth.")
834 endif
835 if (cs%do_brine_plume .and. .not. associated(fluxes%salt_left_behind)) then
836 call mom_error(fatal, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
837 "Brine plume parameterization requires DO_BRINE_PLUME\n"//&
838 "to be turned on in SIS2 as well as MOM6.")
839 endif
840
841 ! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total
842 ! depth of the ocean is vanishing. It does not (yet) handle a value of zero.
843 ! To accommodate vanishing upper layers, we need to allow for an instantaneous
844 ! distribution of forcing over some finite vertical extent. The bulk mixed layer
845 ! code handles this issue properly.
846 h_limit_fluxes = max(gv%Angstrom_H, gv%H_subroundoff)
847
848 ! diagnostic to see if need to create mass to avoid grounding
849 if (cs%id_createdH>0) cs%createdH(:,:) = 0.
850 numberofgroundings = 0
851
852 !$OMP parallel do default(none) shared(is,ie,js,je,nz,h,tv,nsw,G,GV,US,optics,fluxes, &
853 !$OMP H_limit_fluxes,numberOfGroundings,iGround,jGround,&
854 !$OMP nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing, &
855 !$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, &
856 !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,&
857 !$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, &
858 !$OMP EnthalpyConst,MLD_h) &
859 !$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &
860 !$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
861 !$OMP IforcingDepthScale,g_conv,dSpV_dT,dSpV_dS, &
862 !$OMP dThickness,dTemp,dSalt,hOld,Ithickness, &
863 !$OMP netMassIn,pres,d_pres,p_lay,dSV_dT_2d, &
864 !$OMP netmassinout_rate,netheat_rate,netsalt_rate, &
865 !$OMP drhodt,drhods,pen_sw_bnd_rate, &
866 !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, &
867 !$OMP mixing_depth,A_brine,fraction_left_brine, &
868 !$OMP plume_fraction,dK,total_h) &
869 !$OMP firstprivate(SurfPressure,plume_flux)
870 do j=js,je
871 ! Work in vertical slices for efficiency
872
873 ! Copy state into 2D-slice arrays
874 do k=1,nz ; do i=is,ie
875 h2d(i,k) = h(i,j,k)
876 t2d(i,k) = tv%T(i,j,k)
877 enddo ; enddo
878
879 if (calculate_energetics) then
880 ! The partial derivatives of specific volume with temperature and
881 ! salinity need to be precalculated to avoid having heating of
882 ! tiny layers give nonsensical values.
883 if (associated(tv%p_surf)) then
884 do i=is,ie ; pres(i) = tv%p_surf(i,j) ; enddo
885 else
886 do i=is,ie ; pres(i) = 0.0 ; enddo
887 endif
888 do k=1,nz
889 do i=is,ie
890 d_pres(i) = (gv%g_Earth * gv%H_to_RZ) * h2d(i,k)
891 p_lay(i) = pres(i) + 0.5*d_pres(i)
892 pres(i) = pres(i) + d_pres(i)
893 enddo
894 call calculate_specific_vol_derivs(t2d(:,k), tv%S(:,j,k), p_lay(:), &
895 dsv_dt(:,j,k), dsv_ds(:,j,k), tv%eqn_of_state, eosdom)
896 do i=is,ie ; dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; enddo
897 enddo
898 pen_tke_2d(:,:) = 0.0
899 endif
900
901 ! Nothing more is done on this j-slice if there is no buoyancy forcing.
902 if (.not.associated(fluxes%sw)) cycle
903
904 if (nsw>0) then
905 if (gv%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
906 call extract_optics_slice(optics, j, g, gv, opacity=opacityband, opacity_scale=gv%H_to_Z)
907 else
908 call extract_optics_slice(optics, j, g, gv, opacity=opacityband, opacity_scale=gv%H_to_RZ, &
909 spv_avg=tv%SpV_avg)
910 endif
911 endif
912
913 ! The surface forcing is contained in the fluxes type.
914 ! We aggregate the thermodynamic forcing for a time step into the following:
915 ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step
916 ! = lprec + fprec + vprec + evap + lrunoff + frunoff
917 ! note that lprec generally has sea ice melt/form included.
918 ! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.
919 ! netMassOut < 0 means mass leaves ocean.
920 ! netHeat = heat via surface fluxes [C H ~> degC m or degC kg m-2], excluding the part
921 ! contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
922 ! netSalt = surface salt fluxes [S H ~> ppt m or gSalt m-2]
923 ! Pen_SW_bnd = components to penetrative shortwave radiation split according to bands.
924 ! This field provides that portion of SW from atmosphere that in fact
925 ! enters to the ocean and participates in penetrative SW heating.
926 ! nonpenSW = non-downwelling SW flux, which is absorbed in ocean surface
927 ! (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.
928
929 !----------------------------------------------------------------------------------------
930 !BGR-June 26, 2017{
931 !Temporary action to preserve answers while fixing a bug.
932 ! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns
933 ! the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning
934 ! a second call to extractfluxes1d (causing the diagnostic net_heat to be incorrect).
935 ! Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,
936 ! which means it used the T/S fields after this routine. Therefore, the surface
937 ! buoyancy flux is computed here at the very end of this routine for legacy reasons.
938 ! A few specific notes follow:
939 ! 1) The old method did not included river/calving contributions to heat flux. This
940 ! is kept consistent here via commenting code in the present extractFluxes1d <_rate>
941 ! outputs, but we may reconsider this approach.
942 ! 2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead
943 ! of computing the integrated value (and dividing by dt). Hence the required
944 ! additional outputs from extractFluxes1d.
945 ! *** This is because: A*dt/dt =/= A due to round off.
946 ! 3) The old method computed buoyancy flux after this routine, meaning the returned
947 ! surface fluxes (from extractfluxes1d) must be recorded for use later in the code.
948 ! We could (and maybe should) move that loop up to before the surface fluxes are
949 ! applied, but this will change answers.
950 ! For all these reasons we compute additional values of <_rate> which are preserved
951 ! for the buoyancy flux calculation and reproduce the old answers.
952 ! In the future this needs more detailed investigation to make sure everything is
953 ! consistent and correct. These details should not significantly effect climate,
954 ! but do change answers.
955 !-----------------------------------------------------------------------------------------
956 if (calculate_buoyancy) then
957 call extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt, &
958 h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
959 h2d, t2d, netmassinout, netmassout, netheat, netsalt, &
960 pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw, &
961 net_heat_rate=netheat_rate, net_salt_rate=netsalt_rate, &
962 netmassinout_rate=netmassinout_rate, pen_sw_bnd_rate=pen_sw_bnd_rate)
963 else
964 call extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt, &
965 h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
966 h2d, t2d, netmassinout, netmassout, netheat, netsalt, &
967 pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw)
968 endif
969 ! ea is for passive tracers
970 do i=is,ie
971 ! ea(i,j,1) = netMassInOut(i)
972 if (aggregate_fw_forcing) then
973 netmassout(i) = netmassinout(i)
974 netmassin(i) = 0.
975 else
976 netmassin(i) = netmassinout(i) - netmassout(i)
977 endif
978 if (g%mask2dT(i,j) > 0.0) then
979 fluxes%netMassOut(i,j) = netmassout(i)
980 fluxes%netMassIn(i,j) = netmassin(i)
981 else
982 fluxes%netMassOut(i,j) = 0.0
983 fluxes%netMassIn(i,j) = 0.0
984 endif
985 enddo
986
987 ! Apply the surface boundary fluxes in three steps:
988 ! A/ update mass, temp, and salinity due to all terms except mass leaving
989 ! ocean (and corresponding outward heat content), and ignoring penetrative SW.
990 ! B/ update mass, salt, temp from mass leaving ocean.
991 ! C/ update temp due to penetrative SW
992 if (cs%do_brine_plume) then
993 ! Find the plume mixing depth.
994 do i=is,ie ; total_h(i) = 0.0 ; enddo
995 do k=1,nz ; do i=is,ie ; total_h(i) = total_h(i) + h(i,j,k) ; enddo ; enddo
996 do i=is,ie
997 mixing_depth(i) = min( max(mld_h(i,j) - minimum_forcing_depth, minimum_forcing_depth), &
998 max(total_h(i), gv%angstrom_h) ) + gv%H_subroundoff
999 a_brine(i) = (cs%brine_plume_n + 1) / (mixing_depth(i) ** (cs%brine_plume_n + 1))
1000 enddo
1001 endif
1002
1003 do i=is,ie
1004 if (g%mask2dT(i,j) > 0.) then
1005
1006 ! A/ Update mass, temp, and salinity due to incoming mass flux.
1007 do k=1,1
1008
1009 ! Change in state due to forcing
1010 dthickness = netmassin(i) ! Since we are adding mass, we can use all of it
1011 dtemp = 0.
1012 dsalt = 0.
1013
1014 ! Update the forcing by the part to be consumed within the present k-layer.
1015 ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
1016 netmassin(i) = netmassin(i) - dthickness
1017 ! This line accounts for the temperature of the mass exchange
1018 temp_in = t2d(i,k)
1019 salin_in = 0.0
1020 dtemp = dtemp + dthickness*temp_in*enthalpyconst
1021
1022 ! Diagnostics of heat content associated with mass fluxes
1023 if (.not. associated(fluxes%heat_content_evap)) then
1024 if (associated(fluxes%heat_content_massin)) &
1025 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1026 t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * tv%C_p * idt
1027 if (associated(fluxes%heat_content_massout)) &
1028 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1029 t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * tv%C_p * idt
1030 if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1031 t2d(i,k) * dthickness * gv%H_to_RZ
1032 endif
1033
1034 ! Determine the energetics of river mixing before updating the state.
1035 if (calculate_energetics .and. associated(fluxes%lrunoff) .and. cs%do_rivermix) then
1036 ! Here we add an additional source of TKE to the mixed layer where river
1037 ! is present to simulate unresolved estuaries. The TKE input, TKE_river in
1038 ! [Z3 T-3 ~> m3 s-3], is diagnosed as follows:
1039 ! TKE_river = 0.5*rivermix_depth*g*(1/rho)*drho_ds*
1040 ! River*(Samb - Sriver) = CS%mstar*U_star^3
1041 ! where River is in units of [Z T-1 ~> m s-1].
1042 ! Samb = Ambient salinity at the mouth of the estuary
1043 ! rivermix_depth = The prescribed depth over which to mix river inflow
1044 ! drho_ds = The derivative of density with salt at the ambient surface salinity.
1045 ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
1046 if (gv%Boussinesq) then
1047 rivermixconst = -0.5*(cs%rivermix_depth*dt) * gv%g_Earth_Z_T2 * gv%Rho0
1048 elseif (allocated(tv%SpV_avg)) then
1049 rivermixconst = -0.5*(cs%rivermix_depth*dt) * gv%g_Earth_Z_T2 / tv%SpV_avg(i,j,1)
1050 else
1051 rivermixconst = -0.5*(cs%rivermix_depth*dt) * gv%Rho0 * gv%g_Earth_Z_T2
1052 endif
1053 ctke(i,j,k) = ctke(i,j,k) + max(0.0, rivermixconst*dsv_ds(i,j,1) * &
1054 ((fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) + &
1055 (fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j))) * tv%S(i,j,1))
1056 endif
1057
1058 ! Update state
1059 hold = h2d(i,k) ! Keep original thickness in hand
1060 h2d(i,k) = h2d(i,k) + dthickness ! New thickness
1061 if (h2d(i,k) > 0.0) then
1062 if (calculate_energetics .and. (dthickness > 0.)) then
1063 ! Calculate the energy required to mix the newly added water over
1064 ! the topmost grid cell.
1065 ctke(i,j,k) = ctke(i,j,k) + 0.5*g_hconv2*(hold*dthickness) * &
1066 ((t2d(i,k) - temp_in) * dsv_dt(i,j,k) + (tv%S(i,j,k) - salin_in) * dsv_ds(i,j,k))
1067 endif
1068 ithickness = 1.0/h2d(i,k) ! Inverse new thickness
1069 ! The "if"s below avoid changing T/S by roundoff unnecessarily
1070 if (dthickness /= 0. .or. dtemp /= 0.) t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness
1071 if (dthickness /= 0. .or. dsalt /= 0.) tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness
1072
1073 endif
1074
1075 enddo ! k=1,1
1076
1077 ! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.
1078 fraction_left_brine = 1.0
1079 do k=1,nz
1080 ! Place forcing into this layer if this layer has nontrivial thickness.
1081 ! For layers thin relative to 1/IforcingDepthScale, then distribute
1082 ! forcing into deeper layers.
1083 iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth - netmassout(i) )
1084 ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
1085 fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)
1086
1087 ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
1088 ! limit the forcing applied to this cell, leaving the remaining forcing to
1089 ! be distributed downwards.
1090 if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k)) then
1091 fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
1092 endif
1093
1094 if (cs%do_brine_plume .and. associated(fluxes%salt_left_behind)) then
1095 if (fluxes%salt_left_behind(i,j) > 0 .and. fraction_left_brine > 0.0) then
1096 ! Place forcing into this layer by depth for brine plume parameterization.
1097 if (k == 1) then
1098 dk(i) = 0.5 * h(i,j,k) ! Depth of center of layer K
1099 plume_flux = - (1000.0*us%ppt_to_S * (cs%plume_strength * fluxes%salt_left_behind(i,j))) * gv%RZ_to_H
1100 plume_fraction = 1.0
1101 else
1102 dk(i) = dk(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) ! Depth of center of layer K
1103 plume_flux = 0.0
1104 endif
1105 if (dk(i) <= mixing_depth(i) .and. fraction_left_brine > 0.0) then
1106 plume_fraction = min(fraction_left_brine, (a_brine(i) * dk(i)**cs%brine_plume_n) * h(i,j,k))
1107 else
1108 iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth - netmassout(i) )
1109 ! plume_fraction = fraction_left_brine, unless h2d is less than IforcingDepthScale.
1110 plume_fraction = min(fraction_left_brine, h2d(i,k)*iforcingdepthscale)
1111 endif
1112 fraction_left_brine = fraction_left_brine - plume_fraction
1113 plume_flux = plume_flux + plume_fraction * (1000.0*us%ppt_to_S * (cs%plume_strength * &
1114 fluxes%salt_left_behind(i,j))) * gv%RZ_to_H
1115 else
1116 plume_flux = 0.0
1117 endif
1118 endif
1119
1120 ! Change in state due to forcing
1121
1122 dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
1123 dtemp = fractionofforcing*netheat(i)
1124 dsalt = max( fractionofforcing*netsalt(i), -cs%dSalt_frac_max * h2d(i,k) * tv%S(i,j,k))
1125
1126 ! Update the forcing by the part to be consumed within the present k-layer.
1127 ! If fractionOfForcing = 1, then new netMassOut vanishes.
1128 netmassout(i) = netmassout(i) - dthickness
1129 netheat(i) = netheat(i) - dtemp
1130 netsalt(i) = netsalt(i) - dsalt
1131
1132 ! This line accounts for the temperature of the mass exchange
1133 dtemp = dtemp + dthickness*t2d(i,k)*enthalpyconst
1134
1135 ! Diagnostics of heat content associated with mass fluxes
1136 if (.not. associated(fluxes%heat_content_evap)) then
1137 if (associated(fluxes%heat_content_massin)) &
1138 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
1139 t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * tv%C_p * idt
1140 if (associated(fluxes%heat_content_massout)) &
1141 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
1142 t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * tv%C_p * idt
1143 if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
1144 t2d(i,k) * dthickness * gv%H_to_RZ
1145 endif
1146
1147 ! Update state by the appropriate increment.
1148 hold = h2d(i,k) ! Keep original thickness in hand
1149 h2d(i,k) = h2d(i,k) + dthickness ! New thickness
1150
1151 if (h2d(i,k) > 0.) then
1152 if (calculate_energetics) then
1153 ! Calculate the energy required to mix the newly added water over the topmost grid
1154 ! cell, assuming that the fluxes of heat and salt and rejected brine are initially
1155 ! applied in vanishingly thin layers at the top of the layer before being mixed
1156 ! throughout the layer. Note that dThickness is always <= 0 here, and that
1157 ! negative cTKE is a deficit that will need to be filled later.
1158 ctke(i,j,k) = ctke(i,j,k) - (0.5*h2d(i,k)*g_hconv2) * &
1159 ((dtemp - dthickness*t2d(i,k)) * dsv_dt(i,j,k) + &
1160 (dsalt - dthickness*tv%S(i,j,k)) * dsv_ds(i,j,k))
1161 endif
1162 ithickness = 1.0/h2d(i,k) ! Inverse of new thickness
1163 t2d(i,k) = (hold*t2d(i,k) + dtemp)*ithickness
1164 tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt + plume_flux)*ithickness
1165 elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
1166 call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (h<0)')
1167 write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
1168 write(0,*) 'applyBoundaryFluxesInOut(): netT,netS,netH=', &
1169 us%C_to_degC*netheat(i), us%S_to_ppt*netsalt(i), netmassinout(i)
1170 write(0,*) 'applyBoundaryFluxesInOut(): dT,dS,dH=', &
1171 us%C_to_degC*dtemp, us%S_to_ppt*dsalt, dthickness
1172 write(0,*) 'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hold,h2d(i,k),k
1173 call mom_error(fatal, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
1174 "Complete mass loss in column!")
1175 endif
1176
1177 enddo ! k
1178
1179 ! Check if trying to apply fluxes over land points
1180 elseif ((abs(netheat(i)) + abs(netsalt(i)) + abs(netmassin(i)) + abs(netmassout(i))) > 0.) then
1181
1182 if (.not. cs%ignore_fluxes_over_land) then
1183 call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (land)')
1184 write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
1185 write(0,*) 'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
1186 us%C_to_degC*netheat(i), us%S_to_ppt*netsalt(i), netmassin(i), netmassout(i)
1187
1188 call mom_error(fatal, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
1189 "Mass loss over land?")
1190 endif
1191
1192 endif
1193
1194 ! If anything remains after the k-loop, then we have grounded out, which is a problem.
1195 if (netmassin(i)+netmassout(i) /= 0.0) then
1196!$OMP critical
1197 numberofgroundings = numberofgroundings +1
1198 if (numberofgroundings<=maxgroundings) then
1199 iground(numberofgroundings) = i ! Record i,j location of event for
1200 jground(numberofgroundings) = j ! warning message
1201 hgrounding(numberofgroundings) = netmassin(i)+netmassout(i)
1202 endif
1203!$OMP end critical
1204 if (cs%id_createdH>0) cs%createdH(i,j) = cs%createdH(i,j) - (netmassin(i)+netmassout(i))/dt
1205 endif
1206
1207 enddo ! i
1208
1209 ! Step C/ in the application of fluxes
1210 ! Heat by the convergence of penetrating SW.
1211 ! SW penetrative heating uses the updated thickness from above.
1212
1213 ! Save temperature before increment with SW heating
1214 ! and initialize CS%penSWflux_diag to zero.
1215 if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then
1216 do k=1,nz ; do i=is,ie
1217 cs%penSW_diag(i,j,k) = t2d(i,k)
1218 cs%penSWflux_diag(i,j,k) = 0.0
1219 enddo ; enddo
1220 k=nz+1 ; do i=is,ie
1221 cs%penSWflux_diag(i,j,k) = 0.0
1222 enddo
1223 endif
1224
1225 if (calculate_energetics) then
1226 call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
1227 .false., .true., t2d, pen_sw_bnd, tke=pen_tke_2d, dsv_dt=dsv_dt_2d)
1228 k = 1 ! For setting break-points.
1229 do k=1,nz ; do i=is,ie
1230 ctke(i,j,k) = ctke(i,j,k) + pen_tke_2d(i,k)
1231 enddo ; enddo
1232 else
1233 call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
1234 .false., .true., t2d, pen_sw_bnd)
1235 endif
1236
1237
1238 ! Step D/ copy updated thickness and temperature
1239 ! 2d slice now back into model state.
1240 do k=1,nz ; do i=is,ie
1241 h(i,j,k) = h2d(i,k)
1242 tv%T(i,j,k) = t2d(i,k)
1243 enddo ; enddo
1244
1245 ! Diagnose heating [Q R Z T-1 ~> W m-2] applied to a grid cell from SW penetration
1246 ! Also diagnose the penetrative SW heat flux at base of layer.
1247 if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then
1248
1249 ! convergence of SW into a layer
1250 do k=1,nz ; do i=is,ie
1251 ! Note that the units of penSW_diag change here, from [C ~> degC] to [Q R Z T-1 ~> W m-2].
1252 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
1253 enddo ; enddo
1254
1255 ! Perform a cumulative sum upwards from bottom to
1256 ! diagnose penetrative SW flux at base of tracer cell.
1257 ! CS%penSWflux_diag(i,j,k=1) is penetrative shortwave at top of ocean.
1258 ! CS%penSWflux_diag(i,j,k=kbot+1) is zero, since assume no SW penetrates rock.
1259 ! CS%penSWflux_diag = rsdo and CS%penSW_diag = rsdoabsorb
1260 ! rsdoabsorb(k) = rsdo(k) - rsdo(k+1), so that rsdo(k) = rsdo(k+1) + rsdoabsorb(k)
1261 if (cs%id_penSWflux_diag > 0) then
1262 do k=nz,1,-1 ; do i=is,ie
1263 cs%penSWflux_diag(i,j,k) = cs%penSW_diag(i,j,k) + cs%penSWflux_diag(i,j,k+1)
1264 enddo ; enddo
1265 endif
1266
1267 endif
1268
1269 ! Fill CS%nonpenSW_diag
1270 if (cs%id_nonpenSW_diag > 0) then
1271 do i=is,ie
1272 cs%nonpenSW_diag(i,j) = nonpensw(i) * idt * tv%C_p * gv%H_to_RZ
1273 enddo
1274 endif
1275
1276 ! BGR: Get buoyancy flux to return for ePBL
1277 ! We want the rate, so we use the rate values returned from extractfluxes1d.
1278 ! Note that the *dt values could be divided by dt here, but
1279 ! 1) Answers will change due to round-off
1280 ! 2) Be sure to save their values BEFORE fluxes are used.
1281 if (calculate_buoyancy) then
1282 netpen_rate(:) = 0.0
1283 ! Sum over bands and attenuate as a function of depth.
1284 ! netPen_rate is the netSW as a function of depth, but only the surface value is used here,
1285 ! in which case the values of dt, h, optics and H_limit_fluxes are irrelevant. Consider
1286 ! writing a shorter and simpler variant to handle this very limited case.
1287 ! Find the vertical distances across layers.
1288 ! call thickness_to_dz(h, tv, dz, j, G, GV)
1289 ! call sumSWoverBands(G, GV, US, h2d, dz, optics_nbands(optics), optics, j, dt, &
1290 ! H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
1291 do i=is,ie ; do nb=1,nsw ; netpen_rate(i) = netpen_rate(i) + pen_sw_bnd_rate(nb,i) ; enddo ; enddo
1292
1293 ! 1. Adjust netSalt to reflect dilution effect of FW flux
1294 ! 2. Add in the SW heating for purposes of calculating the net
1295 ! surface buoyancy flux affecting the top layer.
1296 ! 3. Convert to a buoyancy flux, excluding penetrating SW heating
1297 ! BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
1298 if (associated(tv%p_surf)) then ; do i=is,ie ; surfpressure(i) = tv%p_surf(i,j) ; enddo ; endif
1299
1300 if ((.not.gv%Boussinesq) .and. (.not.gv%semi_Boussinesq)) then
1301 g_conv = gv%g_Earth_Z_T2 * gv%H_to_RZ
1302
1303 ! Specific volume derivatives
1304 call calculate_specific_vol_derivs(t2d(:,1), tv%S(:,j,1), surfpressure, dspv_dt, dspv_ds, &
1305 tv%eqn_of_state, eos_domain(g%HI))
1306 do i=is,ie
1307 skinbuoyflux(i,j) = g_conv * &
1308 (dspv_ds(i) * ( netsalt_rate(i) - tv%S(i,j,1)*netmassinout_rate(i)) + &
1309 dspv_dt(i) * ( netheat_rate(i) + netpen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
1310 enddo
1311 else
1312 ! Density derivatives
1313 call calculate_density_derivs(t2d(:,1), tv%S(:,j,1), surfpressure, drhodt, drhods, &
1314 tv%eqn_of_state, eosdom)
1315 do i=is,ie
1316 skinbuoyflux(i,j) = - gorho * gv%H_to_Z * &
1317 (drhods(i) * ( netsalt_rate(i) - tv%S(i,j,1)*netmassinout_rate(i)) + &
1318 drhodt(i) * ( netheat_rate(i) + netpen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
1319 enddo
1320 endif
1321 endif
1322
1323 enddo ! j-loop finish
1324
1325 ! Post the diagnostics
1326 if (cs%id_createdH > 0) call post_data(cs%id_createdH , cs%createdH , cs%diag)
1327 if (cs%id_penSW_diag > 0) call post_data(cs%id_penSW_diag , cs%penSW_diag , cs%diag)
1328 if (cs%id_penSWflux_diag > 0) call post_data(cs%id_penSWflux_diag, cs%penSWflux_diag, cs%diag)
1329 if (cs%id_nonpenSW_diag > 0) call post_data(cs%id_nonpenSW_diag , cs%nonpenSW_diag , cs%diag)
1330
1331! The following check will be ignored if ignore_fluxes_over_land = true
1332 if ((numberofgroundings > 0) .and. .not.cs%ignore_fluxes_over_land) then
1333 do i = 1, min(numberofgroundings, maxgroundings)
1334 call forcing_singlepointprint(fluxes,g,iground(i),jground(i),'applyBoundaryFluxesInOut (grounding)')
1335 write(mesg(1:45),'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
1336 g%geoLatT( iground(i), jground(i)), hgrounding(i)*gv%H_to_m
1337 call mom_error(warning, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
1338 "Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
1339 enddo
1340
1341 if (numberofgroundings - maxgroundings > 0) then
1342 write(mesg, '(I0)') numberofgroundings - maxgroundings
1343 call mom_error(warning, "MOM_diabatic_aux:F90, applyBoundaryFluxesInOut(): "//&
1344 trim(mesg) // " groundings remaining")
1345 endif
1346 endif
1347
1348end subroutine applyboundaryfluxesinout
1349
1350!> This subroutine initializes the parameters and control structure of the diabatic_aux module.
1351subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm, use_ePBL)
1352 type(time_type), target, intent(in) :: time !< The current model time.
1353 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1354 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1355 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1356 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1357 type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output
1358 type(diabatic_aux_cs), pointer :: cs !< A pointer to the control structure for the
1359 !! diabatic_aux module, which is initialized here.
1360 logical, intent(in) :: usealealgorithm !< If true, use the ALE algorithm rather
1361 !! than layered mode.
1362 logical, intent(in) :: use_epbl !< If true, use the implicit energetics planetary
1363 !! boundary layer scheme to determine the diffusivity
1364 !! in the surface boundary layer.
1365
1366 ! This "include" declares and sets the variable "version".
1367# include "version_variable.h"
1368 character(len=40) :: mdl = "MOM_diabatic_aux" ! This module's name.
1369 character(len=200) :: inputdir ! The directory where NetCDF input files
1370 character(len=240) :: chl_filename ! A file from which chl_a concentrations are to be read.
1371 character(len=128) :: chl_file ! Data containing chl_a concentrations. Used
1372 ! when var_pen_sw is defined and reading from file.
1373 character(len=32) :: chl_varname ! Name of chl_a variable in chl_file.
1374 logical :: use_temperature ! True if thermodynamics are enabled.
1375 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
1376 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1377 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1378
1379 if (associated(cs)) then
1380 call mom_error(warning, "diabatic_aux_init called with an "// &
1381 "associated control structure.")
1382 return
1383 else
1384 allocate(cs)
1385 endif
1386
1387 cs%diag => diag
1388 cs%Time => time
1389
1390! Set default, read and log parameters
1391 call log_version(param_file, mdl, version, &
1392 "The following parameters are used for auxiliary diabatic processes.")
1393
1394 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
1395 "If true, temperature and salinity are used as state variables.", default=.true.)
1396
1397 call get_param(param_file, mdl, "RECLAIM_FRAZIL", cs%reclaim_frazil, &
1398 "If true, try to use any frazil heat deficit to cool any "//&
1399 "overlying layers down to the freezing point, thereby "//&
1400 "avoiding the creation of thin ice when the SST is above "//&
1401 "the freezing point.", default=.true., do_not_log=.not.use_temperature)
1402 call get_param(param_file, mdl, "SALT_EXTRACTION_LIMIT", cs%dSalt_frac_max, &
1403 "An upper limit on the fraction of the salt in a layer that can be lost to the "//&
1404 "net surface salt fluxes within a timestep.", &
1405 units="nondim", default=0.9999, do_not_log=.not.use_temperature)
1406 cs%dSalt_frac_max = max(min(cs%dSalt_frac_max, 1.0), 0.0)
1407 call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", cs%pressure_dependent_frazil, &
1408 "If true, use a pressure dependent freezing temperature "//&
1409 "when making frazil. The default is false, which will be "//&
1410 "faster but is inappropriate with ice-shelf cavities.", &
1411 default=.false., do_not_log=.not.use_temperature)
1412
1413 if (use_epbl) then
1414 call get_param(param_file, mdl, "IGNORE_FLUXES_OVER_LAND", cs%ignore_fluxes_over_land,&
1415 "If true, the model does not check if fluxes are being applied "//&
1416 "over land points. This is needed when the ocean is coupled "//&
1417 "with ice shelves and sea ice, since the sea ice mask needs to "//&
1418 "be different than the ocean mask to avoid sea ice formation "//&
1419 "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
1420 call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
1421 "If true, apply additional mixing wherever there is "//&
1422 "runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
1423 "if the ocean is that deep.", default=.false.)
1424 if (cs%do_rivermix) &
1425 call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
1426 "The depth to which rivers are mixed if DO_RIVERMIX is "//&
1427 "defined.", units="m", default=0.0, scale=us%m_to_Z)
1428 else
1429 cs%do_rivermix = .false. ; cs%rivermix_depth = 0.0 ; cs%ignore_fluxes_over_land = .false.
1430 endif
1431
1432 if (gv%nkml == 0) then
1433 call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
1434 "If true, use the fluxes%runoff_Hflx field to set the "//&
1435 "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
1436 default=.false., do_not_log=.not.use_temperature)
1437 call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
1438 "If true, use the fluxes%calving_Hflx field to set the "//&
1439 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
1440 default=.false., do_not_log=.not.use_temperature)
1441 else
1442 cs%use_river_heat_content = .false.
1443 cs%use_calving_heat_content = .false.
1444 endif
1445
1446 call get_param(param_file, mdl, "DO_BRINE_PLUME", cs%do_brine_plume, &
1447 "If true, use a brine plume parameterization from "//&
1448 "Nguyen et al., 2009.", default=.false.)
1449 call get_param(param_file, mdl, "BRINE_PLUME_EXPONENT", cs%brine_plume_n, &
1450 "If using the brine plume parameterization, set the integer exponent.", &
1451 default=5, do_not_log=.not.cs%do_brine_plume)
1452 call get_param(param_file, mdl, "BRINE_PLUME_FRACTION", cs%plume_strength, &
1453 "Fraction of the available brine to mix down using the brine plume parameterization.", &
1454 units="nondim", default=1.0, do_not_log=.not.cs%do_brine_plume)
1455
1456 if (usealealgorithm) then
1457 cs%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
1458 time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
1459 "m s-1", conversion=gv%H_to_m*us%s_to_T)
1460 if (cs%id_createdH>0) allocate(cs%createdH(isd:ied,jsd:jed))
1461
1462 ! diagnostic for heating of a grid cell from convergence of SW heat into the cell
1463 cs%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb', &
1464 diag%axesTL, time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
1465 'W m-2', conversion=us%QRZ_T_to_W_m2, &
1466 standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer', v_extensive=.true.)
1467
1468 ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)
1469 ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).
1470 cs%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo', &
1471 diag%axesTi, time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
1472 'W m-2', conversion=us%QRZ_T_to_W_m2, standard_name='downwelling_shortwave_flux_in_sea_water')
1473
1474 ! need both arrays for the SW diagnostics (one for flux, one for convergence)
1475 if (cs%id_penSW_diag>0 .or. cs%id_penSWflux_diag>0) then
1476 allocate(cs%penSW_diag(isd:ied,jsd:jed,nz), source=0.0)
1477 allocate(cs%penSWflux_diag(isd:ied,jsd:jed,nz+1), source=0.0)
1478 endif
1479
1480 ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)
1481 cs%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW', &
1482 diag%axesT1, time, &
1483 'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
1484 'W m-2', conversion=us%QRZ_T_to_W_m2, &
1485 standard_name='nondownwelling_shortwave_flux_in_sea_water')
1486 if (cs%id_nonpenSW_diag > 0) then
1487 allocate(cs%nonpenSW_diag(isd:ied,jsd:jed), source=0.0)
1488 endif
1489 endif
1490
1491 if (use_temperature) then
1492 call get_param(param_file, mdl, "VAR_PEN_SW", cs%var_pen_sw, &
1493 "If true, use one of the CHL_A schemes specified by "//&
1494 "OPACITY_SCHEME to determine the e-folding depth of "//&
1495 "incoming short wave radiation.", default=.false.)
1496 if (cs%var_pen_sw) then
1497
1498 call get_param(param_file, mdl, "CHL_FROM_FILE", cs%chl_from_file, &
1499 "If true, chl_a is read from a file.", default=.true.)
1500 if (cs%chl_from_file) then
1501 call time_interp_external_init()
1502
1503 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1504 call get_param(param_file, mdl, "CHL_FILE", chl_file, &
1505 "CHL_FILE is the file containing chl_a concentrations in "//&
1506 "the variable CHL_A. It is used when VAR_PEN_SW and "//&
1507 "CHL_FROM_FILE are true.", fail_if_missing=.true.)
1508 chl_filename = trim(slasher(inputdir))//trim(chl_file)
1509 call log_param(param_file, mdl, "INPUTDIR/CHL_FILE", chl_filename)
1510 call get_param(param_file, mdl, "CHL_VARNAME", chl_varname, &
1511 "Name of CHL_A variable in CHL_FILE.", default='CHL_A')
1512 if (modulo(g%Domain%turns, 4) /= 0) then
1513 cs%sbc_chl = init_external_field(chl_filename, trim(chl_varname), mom_domain=g%Domain%domain_in)
1514 else
1515 cs%sbc_chl = init_external_field(chl_filename, trim(chl_varname), mom_domain=g%Domain)
1516 endif
1517 endif
1518
1519 cs%id_chl = register_diag_field('ocean_model', 'Chl_opac', diag%axesT1, time, &
1520 'Surface chlorophyll A concentration used to find opacity', 'mg m-3')
1521 endif
1522 endif
1523
1524 id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=clock_routine)
1525 id_clock_frazil = cpu_clock_id('(Ocean frazil)', grain=clock_routine)
1526
1527end subroutine diabatic_aux_init
1528
1529!> This subroutine initializes the control structure and any related memory
1530!! for the diabatic_aux module.
1531subroutine diabatic_aux_end(CS)
1532 type(diabatic_aux_cs), pointer :: cs !< The control structure returned by a previous
1533 !! call to diabatic_aux_init; it is deallocated here.
1534
1535 if (.not.associated(cs)) return
1536
1537 if (cs%id_createdH >0) deallocate(cs%createdH)
1538 if (cs%id_penSW_diag >0) deallocate(cs%penSW_diag)
1539 if (cs%id_penSWflux_diag >0) deallocate(cs%penSWflux_diag)
1540 if (cs%id_nonpenSW_diag >0) deallocate(cs%nonpenSW_diag)
1541
1542 if (associated(cs)) deallocate(cs)
1543
1544end subroutine diabatic_aux_end
1545
1546!> \namespace mom_diabatic_aux
1547!!
1548!! This module contains subroutines that apply various diabatic processes. Usually these
1549!! subroutines are called from the MOM_diabatic module. All of these routines use appropriate
1550!! limiters or logic to work properly with arbitrary layer thicknesses (including massless layers)
1551!! and an arbitrarily large timestep.
1552!!
1553!! The subroutine make_frazil facilitates the formation of frazil ice when the ocean water
1554!! drops below the in situ freezing point by heating the water to the freezing point and
1555!! accumulating the required heat for exchange with the sea-ice module.
1556!!
1557!! The subroutine adjust_salt adds salt as necessary to keep the salinity above a
1558!! specified minimum value, and keeps track of the cumulative additions. If the minimum
1559!! salinity is the natural value of 0, this routine should never do anything.
1560!!
1561!! The subroutine differential_diffuse_T_S solves a pair of tridiagonal equations for
1562!! the diffusion of temperatures and salinities with differing diffusivities.
1563!!
1564!! The subroutine triDiagTS solves a tridiagonal equations for the evolution of temperatures
1565!! and salinities due to net entrainment by layers and a diffusion with the same diffusivity.
1566!!
1567!! The subroutine triDiagTS_Eulerian solves a tridiagonal equations for the evolution of
1568!! temperatures and salinities due to diffusion with the same diffusivity, but no net entrainment.
1569!!
1570!! The subroutine find_uv_at_h interpolates velocities to thickness points, optionally also
1571!! using tridiagonal equations to solve for the impacts of net entrainment or mixing of
1572!! momentum between layers.
1573!!
1574!! The subroutine set_pen_shortwave determines the optical properties of the water column and
1575!! the net shortwave fluxes, and stores them in the optics type, working via calls to set_opacity.
1576!!
1577!! The subroutine applyBoundaryFluxesInOut updates the layer thicknesses, temperatures and
1578!! salinities due to the application of the surface forcing. It may also calculate the implied
1579!! turbulent kinetic energy requirements for this forcing to be mixed over the model's finite
1580!! vertical resolution in the surface layers.
1581end module mom_diabatic_aux