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