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