← back to index

src/parameterizations/vertical/MOM_diabatic_aux.F90

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

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Provides functions for some diabatic processes such as frazil, brine rejection,
6!! tendency due to surface flux divergence.
7module MOM_diabatic_aux
8
9use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
10use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
11use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
12use MOM_diag_mediator, only : diag_ctrl, time_type
13use MOM_EOS, only : calculate_density, calculate_TFreeze, EOS_domain
14use MOM_EOS, only : calculate_specific_vol_derivs, calculate_density_derivs
15use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery
16use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
17use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
18use MOM_forcing_type, only : forcing, extractFluxes1d, forcing_SinglePointPrint
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
24use MOM_opacity, only : set_opacity, opacity_CS, extract_optics_slice, extract_optics_fields
25use MOM_opacity, only : optics_type, optics_nbands, absorbRemainingSW, sumSWoverBands
26use MOM_tracer_flow_control, only : get_chl_from_model, tracer_flow_control_CS
27use MOM_unit_scaling, only : unit_scale_type
28use MOM_variables, only : thermo_var_ptrs
29use MOM_verticalGrid, only : verticalGrid_type
30
31implicit none ; private
32
33#include <MOM_memory.h>
34
35public diabatic_aux_init, diabatic_aux_end
36public make_frazil, adjust_salt, differential_diffuse_T_S, triDiagTS, triDiagTS_Eulerian
37public find_uv_at_h, applyBoundaryFluxesInOut, set_pen_shortwave
38
39! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
40! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
41! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
42! vary with the Boussinesq approximation, the Boussinesq variant is given first.
43
44!> Control structure for diabatic_aux
45type, public :: diabatic_aux_CS ; private
46 logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff at the
47 !! river mouths to a depth of "rivermix_depth"
48 real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m].
49 real :: dSalt_frac_max !< An upper limit on the fraction of the salt in a layer that can be
50 !! lost to the net surface salt fluxes within a timestep [nondim]
51 logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to
52 !! to cool the topmost layer down to the freezing
53 !! point. The default is true.
54 logical :: pressure_dependent_frazil !< If true, use a pressure dependent
55 !! freezing temperature when making frazil. The
56 !! default is false, which will be faster but is
57 !! inappropriate with ice-shelf cavities.
58 logical :: ignore_fluxes_over_land !< If true, the model does not check
59 !! if fluxes are applied over land points. This
60 !! flag must be used when the ocean is coupled with
61 !! sea ice and ice shelves and use_ePBL = true.
62 logical :: use_river_heat_content !< If true, assumes that ice-ocean boundary
63 !! has provided a river heat content. Otherwise, runoff
64 !! is added with a temperature of the local SST.
65 logical :: use_calving_heat_content !< If true, assumes that ice-ocean boundary
66 !! has provided a calving heat content. Otherwise, calving
67 !! is added with a temperature of the local SST.
68 logical :: var_pen_sw !< If true, use one of the CHL_A schemes to determine the
69 !! e-folding depth of incoming shortwave radiation.
70 type(external_field) :: sbc_chl !< A handle used in time interpolation of
71 !! chlorophyll read from a file.
72 logical :: chl_from_file !< If true, chl_a is read from a file.
73 logical :: do_brine_plume !< If true, insert salt flux below the surface according to
74 !! a parameterization by \cite Nguyen2009.
75 integer :: brine_plume_n !< The exponent in the brine plume parameterization.
76 real :: plume_strength !< Fraction of the available brine to take to the bottom of the mixed
77 !! layer [nondim].
78
79 type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
80 type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output
81
82 ! Diagnostic handles
83 integer :: id_createdH = -1 !< Diagnostic ID of mass added to avoid grounding
84 integer :: id_brine_lay = -1 !< Diagnostic ID of which layer receives the brine
85 integer :: id_penSW_diag = -1 !< Diagnostic ID of Penetrative shortwave heating (flux convergence)
86 integer :: id_penSWflux_diag = -1 !< Diagnostic ID of Penetrative shortwave flux
87 integer :: id_nonpenSW_diag = -1 !< Diagnostic ID of Non-penetrative shortwave heating
88 integer :: id_Chl = -1 !< Diagnostic ID of chlorophyll-A handles for opacity
89
90 ! Optional diagnostic arrays
91 real, allocatable, dimension(:,:) :: createdH !< The amount of volume added in order to
92 !! avoid grounding [H T-1 ~> m s-1]
93 real, allocatable, dimension(:,:,:) :: penSW_diag !< Heating in a layer from convergence of
94 !! penetrative SW [Q R Z T-1 ~> W m-2]
95 real, allocatable, dimension(:,:,:) :: penSWflux_diag !< Penetrative SW flux at base of grid
96 !! layer [Q R Z T-1 ~> W m-2]
97 real, allocatable, dimension(:,:) :: nonpenSW_diag !< Non-downwelling SW radiation at ocean
98 !! surface [Q R Z T-1 ~> W m-2]
99
100end type diabatic_aux_CS
101
102!>@{ CPU time clock IDs
103integer :: id_clock_uv_at_h, id_clock_frazil
104!>@}
105
106contains
107
108!> Frazil formation keeps the temperature above the freezing point.
109!! This subroutine warms any water that is colder than the (currently
110!! surface) freezing point up to the freezing point and accumulates
111!! the required heat (in [Q R Z ~> J m-2]) in tv%frazil.
1120subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
113 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
114 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
115 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
116 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
117 type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any available
118 !! thermodynamic fields.
119 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
120 type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous
121 !! call to diabatic_aux_init.
122 real, dimension(SZI_(G),SZJ_(G)), &
123 optional, intent(in) :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa].
124 integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil
125 ! Local variables
126 real, dimension(SZI_(G)) :: &
12748 fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].
12848 T_freeze, & ! The freezing potential temperature at the current salinity [C ~> degC].
12948 ps ! Surface pressure [R L2 T-2 ~> Pa]
130 real, dimension(SZI_(G),SZK_(GV)) :: &
13148 pressure ! The pressure at the middle of each layer [R L2 T-2 ~> Pa].
132 real :: H_to_RL2_T2 ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
133 real :: hc ! A layer's heat capacity [Q R Z C-1 ~> J m-2 degC-1].
134 logical :: T_fr_set ! True if the freezing point has been calculated for a
135 ! row of points.
136 integer :: i, j, k, is, ie, js, je, nz
137
13824 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
13924 if (present(halo)) then
14012 is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo
141 endif
142
14324 call cpu_clock_begin(id_clock_frazil)
144
14524 if (.not.CS%pressure_dependent_frazil) then
1460 do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo
147 else
14824 H_to_RL2_T2 = GV%H_to_RZ * GV%g_Earth
149 endif
150 !$OMP parallel do default(shared) private(fraz_col,T_fr_set,T_freeze,hc,ps) &
151 !$OMP firstprivate(pressure) ! pressure might be set above, so should be firstprivate
1521488 do j=js,je
153188856 ps(:) = 0.0
1541464 if (PRESENT(p_surf)) then ; do i=is,ie
1550 ps(i) = p_surf(i,j)
156 enddo ; endif
157
158178632 do i=is,ie ; fraz_col(i) = 0.0 ; enddo
159
1601464 if (CS%pressure_dependent_frazil) then
161178632 do i=is,ie
162178632 pressure(i,1) = ps(i) + (0.5*H_to_RL2_T2)*h(i,j,1)
163 enddo
16413220232 do k=2,nz ; do i=is,ie
16513218768 pressure(i,k) = pressure(i,k-1) + (0.5*H_to_RL2_T2) * (h(i,j,k) + h(i,j,k-1))
166 enddo ; enddo
167 endif
168
1691464 if (CS%reclaim_frazil) then
1701464 T_fr_set = .false.
171178632 do i=is,ie ; if (tv%frazil(i,j) > 0.0) then
1720 if (.not.T_fr_set) then
173 call calculate_TFreeze(tv%S(i:ie,j,1), pressure(i:ie,1), T_freeze(i:ie), &
1740 tv%eqn_of_state)
1750 T_fr_set = .true.
176 endif
177
1780 if (tv%T(i,j,1) > T_freeze(i)) then
179 ! If frazil had previously been formed, but the surface temperature is now
180 ! above freezing, cool the surface layer with the frazil heat deficit.
1810 hc = (tv%C_p*GV%H_to_RZ) * h(i,j,1)
1820 if (tv%frazil(i,j) - hc * (tv%T(i,j,1) - T_freeze(i)) <= 0.0) then
1830 tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j) / hc
1840 tv%frazil(i,j) = 0.0
185 else
1860 tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - T_freeze(i))
1870 tv%T(i,j,1) = T_freeze(i)
188 endif
189 endif
190 endif ; enddo
191 endif
192
193111264 do k=nz,1,-1
194109800 T_fr_set = .false.
19513398864 do i=is,ie
19613287600 if ((G%mask2dT(i,j) > 0.0) .and. &
197109800 ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then
1980 if (.not.T_fr_set) then
199 call calculate_TFreeze(tv%S(i:ie,j,k), pressure(i:ie,k), T_freeze(i:ie), &
2000 tv%eqn_of_state)
2010 T_fr_set = .true.
202 endif
203
2040 hc = (tv%C_p*GV%H_to_RZ) * h(i,j,k)
2050 if (h(i,j,k) <= 10.0*(GV%Angstrom_H + GV%H_subroundoff)) then
206 ! Very thin layers should not be cooled by the frazil flux.
2070 if (tv%T(i,j,k) < T_freeze(i)) then
2080 fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k))
2090 tv%T(i,j,k) = T_freeze(i)
210 endif
2110 elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < T_freeze(i))) then
2120 if (fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k)) < 0.0) then
2130 tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc
2140 fraz_col(i) = 0.0
215 else
2160 fraz_col(i) = fraz_col(i) + hc * (T_freeze(i) - tv%T(i,j,k))
2170 tv%T(i,j,k) = T_freeze(i)
218 endif
219 endif
220 endif
221 enddo
222 enddo
223178656 do i=is,ie
224178632 tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i)
225 enddo
226 enddo
227
22824 tv%frazil_was_reset = .false.
229
23024 call cpu_clock_end(id_clock_frazil)
231
23224end subroutine make_frazil
233
234!> This subroutine applies double diffusion to T & S, assuming no diapycnal mass
235!! fluxes, using a simple tridiagonal solver.
2360subroutine differential_diffuse_T_S(h, T, S, Kd_T, Kd_S, tv, dt, G, GV)
237 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
238 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
239 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
240 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
241 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
242 intent(inout) :: T !< Potential temperature [C ~> degC].
243 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
244 intent(inout) :: S !< Salinity [PSU] or [gSalt/kg], generically [S ~> ppt].
245 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
246 intent(in) :: Kd_T !< The extra diffusivity of temperature due to
247 !! double diffusion relative to the diffusivity of
248 !! density [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
249 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
250 intent(in) :: Kd_S !< The extra diffusivity of salinity due to
251 !! double diffusion relative to the diffusivity of
252 !! density [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
253 type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
254 !! available thermodynamic fields.
255 real, intent(in) :: dt !< Time increment [T ~> s].
256
257 ! local variables
258 real, dimension(SZI_(G)) :: &
2590 b1_T, b1_S, & ! Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].
2600 d1_T, d1_S ! Variables used by the tridiagonal solvers [nondim].
261 real, dimension(SZI_(G),SZK_(GV)) :: &
2620 dz, & ! Height change across layers [Z ~> m]
2630 c1_T, c1_S ! Variables used by the tridiagonal solvers [H ~> m or kg m-2].
264 real, dimension(SZI_(G),SZK_(GV)+1) :: &
2650 mix_T, mix_S ! Mixing distances in both directions across each interface [H ~> m or kg m-2].
266 real :: h_tr ! h_tr is h at tracer points with a tiny thickness
267 ! added to ensure positive definiteness [H ~> m or kg m-2].
268 real :: h_neglect ! A thickness that is so small it is usually lost
269 ! in roundoff and can be neglected [H ~> m or kg m-2].
270 real :: dz_neglect ! A vertical distance that is so small it is usually lost
271 ! in roundoff and can be neglected [Z ~> m].
272 real :: I_dz_int ! The inverse of the height scale associated with an interface [Z-1 ~> m-1].
273 real :: b_denom_T ! The first term in the denominator for the expression for b1_T [H ~> m or kg m-2].
274 real :: b_denom_S ! The first term in the denominator for the expression for b1_S [H ~> m or kg m-2].
275 integer :: i, j, k, is, ie, js, je, nz
276
2770 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2780 h_neglect = GV%H_subroundoff
2790 dz_neglect = GV%dZ_subroundoff
280
281 !$OMP parallel do default(private) shared(is,ie,js,je,h,h_neglect,dt,Kd_T,Kd_S,G,GV,T,S,nz)
2820 do j=js,je
283
284 ! Find the vertical distances across layers.
2850 call thickness_to_dz(h, tv, dz, j, G, GV)
286
2870 do i=is,ie
2880 I_dz_int = 1.0 / (0.5 * (dz(i,1) + dz(i,2)) + dz_neglect)
2890 mix_T(i,2) = (dt * Kd_T(i,j,2)) * I_dz_int
2900 mix_S(i,2) = (dt * Kd_S(i,j,2)) * I_dz_int
291
2920 h_tr = h(i,j,1) + h_neglect
2930 b1_T(i) = 1.0 / (h_tr + mix_T(i,2))
2940 b1_S(i) = 1.0 / (h_tr + mix_S(i,2))
2950 d1_T(i) = h_tr * b1_T(i)
2960 d1_S(i) = h_tr * b1_S(i)
2970 T(i,j,1) = (b1_T(i)*h_tr)*T(i,j,1)
2980 S(i,j,1) = (b1_S(i)*h_tr)*S(i,j,1)
299 enddo
3000 do k=2,nz-1 ; do i=is,ie
301 ! Calculate the mixing across the interface below this layer.
3020 I_dz_int = 1.0 / (0.5 * (dz(i,k) + dz(i,k+1)) + dz_neglect)
3030 mix_T(i,K+1) = ((dt * Kd_T(i,j,K+1))) * I_dz_int
3040 mix_S(i,K+1) = ((dt * Kd_S(i,j,K+1))) * I_dz_int
305
3060 c1_T(i,k) = mix_T(i,K) * b1_T(i)
3070 c1_S(i,k) = mix_S(i,K) * b1_S(i)
308
3090 h_tr = h(i,j,k) + h_neglect
3100 b_denom_T = h_tr + d1_T(i)*mix_T(i,K)
3110 b_denom_S = h_tr + d1_S(i)*mix_S(i,K)
3120 b1_T(i) = 1.0 / (b_denom_T + mix_T(i,K+1))
3130 b1_S(i) = 1.0 / (b_denom_S + mix_S(i,K+1))
3140 d1_T(i) = b_denom_T * b1_T(i)
3150 d1_S(i) = b_denom_S * b1_S(i)
316
3170 T(i,j,k) = b1_T(i) * (h_tr*T(i,j,k) + mix_T(i,K)*T(i,j,k-1))
3180 S(i,j,k) = b1_S(i) * (h_tr*S(i,j,k) + mix_S(i,K)*S(i,j,k-1))
319 enddo ; enddo
3200 do i=is,ie
3210 c1_T(i,nz) = mix_T(i,nz) * b1_T(i)
3220 c1_S(i,nz) = mix_S(i,nz) * b1_S(i)
323
3240 h_tr = h(i,j,nz) + h_neglect
3250 b1_T(i) = 1.0 / (h_tr + d1_T(i)*mix_T(i,nz))
3260 b1_S(i) = 1.0 / (h_tr + d1_S(i)*mix_S(i,nz))
327
3280 T(i,j,nz) = b1_T(i) * (h_tr*T(i,j,nz) + mix_T(i,nz)*T(i,j,nz-1))
3290 S(i,j,nz) = b1_S(i) * (h_tr*S(i,j,nz) + mix_S(i,nz)*S(i,j,nz-1))
330 enddo
3310 do k=nz-1,1,-1 ; do i=is,ie
3320 T(i,j,k) = T(i,j,k) + c1_T(i,k+1)*T(i,j,k+1)
3330 S(i,j,k) = S(i,j,k) + c1_S(i,k+1)*S(i,j,k+1)
334 enddo ; enddo
335 enddo
3360end subroutine differential_diffuse_T_S
337
338!> This subroutine keeps salinity from falling below a small but positive threshold.
339!! This usually occurs when the ice model attempts to extract more salt then
340!! is actually available to it from the ocean.
3410subroutine adjust_salt(h, tv, G, GV, CS)
342 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
343 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
344 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
345 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
346 type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
347 !! available thermodynamic fields.
348 type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous
349 !! call to diabatic_aux_init.
350
351 ! local variables
3520 real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement [S R Z ~> gSalt m-2]
353 real :: S_min !< The minimum salinity [S ~> ppt].
354 real :: mc !< A layer's mass [R Z ~> kg m-2].
355 integer :: i, j, k, is, ie, js, je, nz
356
3570 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
358
359! call cpu_clock_begin(id_clock_adjust_salt)
360
3610 S_min = tv%min_salinity
362
3630 salt_add_col(:,:) = 0.0
364
365 !$OMP parallel do default(shared) private(mc)
3660 do j=js,je
3670 do k=nz,1,-1 ; do i=is,ie
3680 if ( (G%mask2dT(i,j) > 0.0) .and. &
3690 ((tv%S(i,j,k) < S_min) .or. (salt_add_col(i,j) > 0.0)) ) then
3700 mc = GV%H_to_RZ * h(i,j,k)
3710 if (h(i,j,k) <= 10.0*GV%Angstrom_H) then
372 ! Very thin layers should not be adjusted by the salt flux
3730 if (tv%S(i,j,k) < S_min) then
3740 salt_add_col(i,j) = salt_add_col(i,j) + mc * (S_min - tv%S(i,j,k))
3750 tv%S(i,j,k) = S_min
376 endif
3770 elseif (salt_add_col(i,j) + mc * (S_min - tv%S(i,j,k)) <= 0.0) then
3780 tv%S(i,j,k) = tv%S(i,j,k) - salt_add_col(i,j) / mc
3790 salt_add_col(i,j) = 0.0
380 else
3810 salt_add_col(i,j) = salt_add_col(i,j) + mc * (S_min - tv%S(i,j,k))
3820 tv%S(i,j,k) = S_min
383 endif
384 endif
385 enddo ; enddo
3860 do i=is,ie
3870 tv%salt_deficit(i,j) = tv%salt_deficit(i,j) + salt_add_col(i,j)
388 enddo
389 enddo
390! call cpu_clock_end(id_clock_adjust_salt)
391
3920end subroutine adjust_salt
393
394!> This is a simple tri-diagonal solver for T and S.
395!! "Simple" means it only uses arrays hold, ea and eb.
3960subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S)
397 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
398 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
399 integer, intent(in) :: is !< The start i-index to work on.
400 integer, intent(in) :: ie !< The end i-index to work on.
401 integer, intent(in) :: js !< The start j-index to work on.
402 integer, intent(in) :: je !< The end j-index to work on.
403 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: hold !< The layer thicknesses before entrainment,
404 !! [H ~> m or kg m-2].
405 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: ea !< The amount of fluid entrained from the layer
406 !! above within this time step [H ~> m or kg m-2]
407 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: eb !< The amount of fluid entrained from the layer
408 !! below within this time step [H ~> m or kg m-2]
409 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: T !< Layer potential temperatures [C ~> degC].
410 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: S !< Layer salinities [S ~> ppt].
411
412 ! Local variables
4130 real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
4140 real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim].
4150 real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim].
416 real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2].
417 integer :: i, j, k
418
419 !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
4200 do j=js,je
4210 do i=is,ie
4220 h_tr = hold(i,j,1) + GV%H_subroundoff
4230 b1(i) = 1.0 / (h_tr + eb(i,j,1))
4240 d1(i) = h_tr * b1(i)
4250 T(i,j,1) = (b1(i)*h_tr)*T(i,j,1)
4260 S(i,j,1) = (b1(i)*h_tr)*S(i,j,1)
427 enddo
4280 do k=2,GV%ke ; do i=is,ie
4290 c1(i,k) = eb(i,j,k-1) * b1(i)
4300 h_tr = hold(i,j,k) + GV%H_subroundoff
4310 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
4320 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
4330 d1(i) = b_denom_1 * b1(i)
4340 T(i,j,k) = b1(i) * (h_tr*T(i,j,k) + ea(i,j,k)*T(i,j,k-1))
4350 S(i,j,k) = b1(i) * (h_tr*S(i,j,k) + ea(i,j,k)*S(i,j,k-1))
436 enddo ; enddo
4370 do k=GV%ke-1,1,-1 ; do i=is,ie
4380 T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1)
4390 S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1)
440 enddo ; enddo
441 enddo
4420end subroutine triDiagTS
443
444!> This is a simple tri-diagonal solver for T and S, with mixing across interfaces but no net
445!! transfer of mass.
44612subroutine triDiagTS_Eulerian(G, GV, is, ie, js, je, hold, ent, T, S)
447 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
448 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
449 integer, intent(in) :: is !< The start i-index to work on.
450 integer, intent(in) :: ie !< The end i-index to work on.
451 integer, intent(in) :: js !< The start j-index to work on.
452 integer, intent(in) :: je !< The end j-index to work on.
453 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: hold !< The layer thicknesses before entrainment,
454 !! [H ~> m or kg m-2].
455 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: ent !< The amount of fluid mixed across an interface
456 !! within this time step [H ~> m or kg m-2]
457 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: T !< Layer potential temperatures [C ~> degC].
458 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: S !< Layer salinities [S ~> ppt].
459
460 ! Local variables
46124 real :: b1(SZIB_(G)) ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
46224 real :: d1(SZIB_(G)) ! A variable used by the tridiagonal solver [nondim].
46324 real :: c1(SZIB_(G),SZK_(GV)) ! A variable used by the tridiagonal solver [nondim].
464 real :: h_tr, b_denom_1 ! Two temporary thicknesses [H ~> m or kg m-2].
465 integer :: i, j, k
466
467 !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
468732 do j=js,je
46987120 do i=is,ie
47086400 h_tr = hold(i,j,1) + GV%H_subroundoff
47186400 b1(i) = 1.0 / (h_tr + ent(i,j,2))
47286400 d1(i) = h_tr * b1(i)
47386400 T(i,j,1) = (b1(i)*h_tr)*T(i,j,1)
47487120 S(i,j,1) = (b1(i)*h_tr)*S(i,j,1)
475 enddo
4766447600 do k=2,GV%ke ; do i=is,ie
4776393600 c1(i,k) = ent(i,j,K) * b1(i)
4786393600 h_tr = hold(i,j,k) + GV%H_subroundoff
4796393600 b_denom_1 = h_tr + d1(i)*ent(i,j,K)
4806393600 b1(i) = 1.0 / (b_denom_1 + ent(i,j,K+1))
4816393600 d1(i) = b_denom_1 * b1(i)
4826393600 T(i,j,k) = b1(i) * (h_tr*T(i,j,k) + ent(i,j,K)*T(i,j,k-1))
4836446880 S(i,j,k) = b1(i) * (h_tr*S(i,j,k) + ent(i,j,K)*S(i,j,k-1))
484 enddo ; enddo
4856447612 do k=GV%ke-1,1,-1 ; do i=is,ie
4866393600 T(i,j,k) = T(i,j,k) + c1(i,k+1)*T(i,j,k+1)
4876446880 S(i,j,k) = S(i,j,k) + c1(i,k+1)*S(i,j,k+1)
488 enddo ; enddo
489 enddo
49012end subroutine triDiagTS_Eulerian
491
492
493!> This subroutine calculates u_h and v_h (velocities at thickness
494!! points), optionally using the entrainment amounts passed in as arguments.
4950subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb, zero_mix)
496 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
497 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
498 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
499 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
500 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
501 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
502 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
503 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
504 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
505 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
506 intent(out) :: u_h !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].
507 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
508 intent(out) :: v_h !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].
509 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
510 optional, intent(in) :: ea !< The amount of fluid entrained from the layer
511 !! above within this time step [H ~> m or kg m-2].
512 !! Omitting ea is the same as setting it to 0.
513 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
514 optional, intent(in) :: eb !< The amount of fluid entrained from the layer
515 !! below within this time step [H ~> m or kg m-2].
516 !! Omitting eb is the same as setting it to 0.
517 logical, optional, intent(in) :: zero_mix !< If true, do the calculation of u_h and
518 !! v_h as though ea and eb were being supplied with
519 !! uniformly zero values.
520
521 ! Local variables
522 real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
523 real :: h_neglect ! A thickness that is so small it is usually lost
524 ! in roundoff and can be neglected [H ~> m or kg m-2].
52548 real :: b1(SZI_(G)) ! A thickness used in the tridiagonal solver [H ~> m or kg m-2]
52648 real :: c1(SZI_(G),SZK_(GV)) ! A variable used in the tridiagonal solver [nondim]
52748 real :: d1(SZI_(G)) ! The complement of c1 [nondim]
528 ! Fractional weights of the neighboring velocity points, ~1/2 in the open ocean.
52948 real :: a_n(SZI_(G)), a_s(SZI_(G)) ! Fractional weights of the neighboring velocity points [nondim]
53048 real :: a_e(SZI_(G)), a_w(SZI_(G)) ! Fractional weights of the neighboring velocity points [nondim]
531 real :: sum_area ! A sum of adjacent areas [L2 ~> m2]
532 real :: Idenom ! The inverse of the denominator in a weighted average [L-2 ~> m-2]
533 logical :: mix_vertically, zero_mixing
534 integer :: i, j, k, is, ie, js, je, nz
53524 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
53624 call cpu_clock_begin(id_clock_uv_at_h)
53724 h_neglect = GV%H_subroundoff
538
53924 mix_vertically = present(ea)
54024 if (present(ea) .neqv. present(eb)) call MOM_error(FATAL, &
541 "find_uv_at_h: Either both ea and eb or neither one must be present "// &
5420 "in call to find_uv_at_h.")
54324 zero_mixing = .false. ; if (present(zero_mix)) zero_mixing = zero_mix
54424 if (zero_mixing) mix_vertically = .false.
545 !$omp target enter data map(alloc: a_w,a_e,a_s,a_n,b1,d1,c1)
546 !$omp target teams loop private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1) &
547 !$omp map(to: ea, eb, h) map(from: u_h, v_h)
5481464 do j=js,je
5491440 do concurrent (i=is:ie)
550172800 sum_area = G%areaCu(I-1,j) + G%areaCu(I,j)
551172800 if (sum_area > 0.0) then
552 ! If this were a simple area weighted average, this would just be I_denom = 1.0 / sum_area.
553 ! The other factor of sqrt(0.5*sum_area*G%IareaT(i,j)) is 1 for open ocean points on a
554 ! Cartesian grid. This construct predates the initial commit of the MOM6 code, and was
555 ! present in the GOLD code before February, 2010. I do not recall why this was added, and
556 ! the GOLD CVS server that contained the relevant history and logs appears to have been
557 ! decommissioned.
558120336 Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
559120336 a_w(i) = G%areaCu(I-1,j) * Idenom
560120336 a_e(i) = G%areaCu(I,j) * Idenom
561 else
56252464 a_w(i) = 0.0 ; a_e(i) = 0.0
563 endif
564
565172800 sum_area = G%areaCv(i,J-1) + G%areaCv(i,J)
566347040 if (sum_area > 0.0) then
567120336 Idenom = sqrt(0.5*G%IareaT(i,j) / sum_area)
568120336 a_s(i) = G%areaCv(i,J-1) * Idenom
569120336 a_n(i) = G%areaCv(i,J) * Idenom
570 else
57152464 a_s(i) = 0.0 ; a_n(i) = 0.0
572 endif
573 enddo
574
5751464 if (mix_vertically) then
5760 do concurrent (i=is:ie)
5770 b_denom_1 = h(i,j,1) + h_neglect
5780 b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
5790 d1(i) = b_denom_1 * b1(i)
5800 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)))
5810 v_h(i,j,1) = (h(i,j,1)*b1(i)) * ((a_n(i)*v(i,J,1)) + (a_s(i)*v(i,J-1,1)))
582 enddo
5830 do k=2,nz ; do concurrent (i=is:ie)
5840 c1(i,k) = eb(i,j,k-1) * b1(i)
5850 b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
5860 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
5870 d1(i) = b_denom_1 * b1(i)
588 u_h(i,j,k) = (h(i,j,k) * ((a_e(i)*u(I,j,k)) + (a_w(i)*u(I-1,j,k))) + &
5890 ea(i,j,k)*u_h(i,j,k-1))*b1(i)
590 v_h(i,j,k) = (h(i,j,k) * ((a_n(i)*v(i,J,k)) + (a_s(i)*v(i,J-1,k))) + &
5910 ea(i,j,k)*v_h(i,j,k-1))*b1(i)
592 enddo ; enddo
5930 do k=nz-1,1,-1 ; do concurrent (i=is:ie)
5940 u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
5950 v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
596 enddo ; enddo
5971440 elseif (zero_mixing) then
598720 do concurrent (i=is:ie)
59986400 b1(i) = 1.0 / (h(i,j,1) + h_neglect)
60086400 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)))
60187120 v_h(i,j,1) = (h(i,j,1)*b1(i)) * ((a_n(i)*v(i,J,1)) + (a_s(i)*v(i,J-1,1)))
602 enddo
603720 do concurrent (k=2:nz, i=is:ie)
6046393600 b1(i) = 1.0 / (h(i,j,k) + h_neglect)
6056393600 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)
6066480720 v_h(i,j,k) = (h(i,j,k) * ((a_n(i)*v(i,J,k)) + (a_s(i)*v(i,J-1,k)))) * b1(i)
607 enddo
608 else
609720 do concurrent (k=1:nz, i=is:ie)
6106480000 u_h(i,j,k) = (a_e(i)*u(I,j,k)) + (a_w(i)*u(I-1,j,k))
6116567120 v_h(i,j,k) = (a_n(i)*v(i,J,k)) + (a_s(i)*v(i,J-1,k))
612 enddo
613 endif
614 enddo
615 !$omp target exit data map(release: a_w,a_e,a_s,a_n,b1,d1,c1)
616
61724 call cpu_clock_end(id_clock_uv_at_h)
61824end subroutine find_uv_at_h
619
620!> Estimate the optical properties of the water column and determine the penetrating shortwave
621!! radiation by band, extracting the relevant information from the fluxes type and storing it
622!! in the optics type for later application. This routine is effectively a wrapper for
623!! set_opacity with added error handling and diagnostics.
62412subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity, tracer_flow_CSp)
625 type(optics_type), pointer :: optics !< An optics structure that has will contain
626 !! information about shortwave fluxes and absorption.
627 type(forcing), intent(inout) :: fluxes !< points to forcing fields
628 !! unused fields have NULL pointers
629 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
630 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
631 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
632 type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
633 type(opacity_CS) :: opacity !< The control structure for the opacity module.
634 type(tracer_flow_control_CS), pointer :: tracer_flow_CSp !< A pointer to the control structure
635 !! organizing the tracer modules.
636
637 ! Local variables
63812 real, dimension(SZI_(G),SZJ_(G)) :: chl_2d !< Vertically uniform chlorophyll-A concentrations [mg m-3]
63912 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: chl_3d !< The chlorophyll-A concentrations of each layer [mg m-3]
640 character(len=128) :: mesg
641 integer :: i, j, is, ie, js, je
64212 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
643
64412 if (.not.associated(optics)) return
645
64612 if (CS%var_pen_sw) then
6470 if (CS%chl_from_file) then
648 ! Only the 2-d surface chlorophyll can be read in from a file. The
649 ! same value is assumed for all layers.
6500 call time_interp_external(CS%sbc_chl, CS%Time, chl_2d, turns=G%HI%turns)
6510 do j=js,je ; do i=is,ie
6520 if ((G%mask2dT(i,j) > 0.0) .and. (chl_2d(i,j) < 0.0)) then
653 write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
654 & I0,", ",I0," lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
6550 chl_2d(i,j), i, j, G%geoLonT(i,j), G%geoLatT(i,j)
6560 call MOM_error(FATAL, "MOM_diabatic_aux set_pen_shortwave: "//trim(mesg))
657 endif
658 enddo ; enddo
659
6600 if (CS%id_chl > 0) call post_data(CS%id_chl, chl_2d, CS%diag)
661
662 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
6630 fluxes%sw_nir_dir, fluxes%sw_nir_dif, G, GV, US, opacity, chl_2d=chl_2d)
664 else
6650 if (.not.associated(tracer_flow_CSp)) call MOM_error(FATAL, &
666 "The tracer flow control structure must be associated when the model sets "//&
6670 "the chlorophyll internally in set_pen_shortwave.")
6680 call get_chl_from_model(chl_3d, G, GV, tracer_flow_CSp)
669
6700 if (CS%id_chl > 0) call post_data(CS%id_chl, chl_3d(:,:,1), CS%diag)
671
672 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
6730 fluxes%sw_nir_dir, fluxes%sw_nir_dif, G, GV, US, opacity, chl_3d=chl_3d)
674 endif
675 else
676 call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
67712 fluxes%sw_nir_dir, fluxes%sw_nir_dif, G, GV, US, opacity)
678 endif
679
680end subroutine set_pen_shortwave
681
682!> Update the thickness, temperature, and salinity due to thermodynamic
683!! boundary forcing (contained in fluxes type) applied to h, tv%T and tv%S,
684!! and calculate the TKE implications of this heating.
68524subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, &
686 aggregate_FW_forcing, evap_CFL_limit, &
68736 minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
68812 SkinBuoyFlux, MLD_h)
689 type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
690 type(ocean_grid_type), intent(in) :: G !< Grid structure
691 type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
692 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
693 real, intent(in) :: dt !< Time-step over which forcing is applied [T ~> s]
694 type(forcing), intent(inout) :: fluxes !< Surface fluxes container
695 type(optics_type), pointer :: optics !< Optical properties container
696 integer, intent(in) :: nsw !< The number of frequency bands of penetrating
697 !! shortwave radiation
698 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
699 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
700 type(thermo_var_ptrs), intent(inout) :: tv !< Structure containing pointers to any
701 !! available thermodynamic fields.
702 logical, intent(in) :: aggregate_FW_forcing !< If False, treat in/out fluxes separately.
703 real, intent(in) :: evap_CFL_limit !< The largest fraction of a layer that
704 !! can be evaporated in one time-step [nondim].
705 real, intent(in) :: minimum_forcing_depth !< The smallest depth over which
706 !! heat and freshwater fluxes is applied [H ~> m or kg m-2].
707 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
708 optional, intent(out) :: cTKE !< Turbulent kinetic energy requirement to mix
709 !! forcing through each layer [R Z3 T-2 ~> J m-2]
710 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
711 optional, intent(out) :: dSV_dT !< Partial derivative of specific volume with
712 !! potential temperature [R-1 C-1 ~> m3 kg-1 degC-1].
713 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
714 optional, intent(out) :: dSV_dS !< Partial derivative of specific volume with
715 !! salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
716 real, dimension(SZI_(G),SZJ_(G)), &
717 optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].
718 real, dimension(:,:), &
719 optional, pointer :: MLD_h !< Mixed layer thickness for brine plumes [H ~> m or kg m-2]
720
721 ! Local variables
722 integer, parameter :: maxGroundings = 5
723 integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings)
724 real :: H_limit_fluxes ! Surface fluxes are scaled down fluxes when the total depth of the ocean
725 ! drops below this value [H ~> m or kg m-2]
726 real :: IforcingDepthScale ! The inverse of the layer thickness below which mass losses are
727 ! shifted to the next deeper layer [H ~> m or kg m-2]
728 real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
729 real :: dThickness ! The change in layer thickness [H ~> m or kg m-2]
730 real :: dTemp ! The integrated change in layer temperature [C H ~> degC m or degC kg m-2]
731 real :: dSalt ! The integrated change in layer salinity [S H ~> ppt m or ppt kg m-2]
732 real :: fractionOfForcing ! THe fraction of the remaining forcing applied to a layer [nondim]
733 real :: hOld ! The original thickness of a layer [H ~> m or kg m-2]
734 real :: Ithickness ! The inverse of the new layer thickness [H-1 ~> m-1 or m2 kg-1]
735 real :: RivermixConst ! A constant used in implementing river mixing [R Z2 T-1 ~> Pa s].
736 real :: EnthalpyConst ! A constant used to control the enthalpy calculation [nondim]
737 ! By default EnthalpyConst = 1.0. If fluxes%heat_content_evap
738 ! is associated enthalpy is provided via coupler and EnthalpyConst = 0.0.
739 real, dimension(SZI_(G)) :: &
74024 d_pres, & ! pressure change across a layer [R L2 T-2 ~> Pa]
74124 p_lay, & ! average pressure in a layer [R L2 T-2 ~> Pa]
74224 pres, & ! pressure at an interface [R L2 T-2 ~> Pa]
74324 netMassInOut, & ! surface water fluxes [H ~> m or kg m-2] over time step
74424 netMassIn, & ! mass entering ocean surface [H ~> m or kg m-2] over a time step
74524 netMassOut, & ! mass leaving ocean surface [H ~> m or kg m-2] over a time step
74624 netHeat, & ! heat via surface fluxes excluding Pen_SW_bnd and netMassOut
747 ! [C H ~> degC m or degC kg m-2]
74824 netSalt, & ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
749 ! [S H ~> ppt m or ppt kg m-2]
75024 nonpenSW, & ! non-downwelling SW, which is absorbed at ocean surface
751 ! [C H ~> degC m or degC kg m-2]
75224 SurfPressure, & ! Surface pressure (approximated as 0.0) [R L2 T-2 ~> Pa]
75324 dRhodT, & ! change in density per change in temperature [R C-1 ~> kg m-3 degC-1]
75424 dRhodS, & ! change in density per change in salinity [R S-1 ~> kg m-3 ppt-1]
75524 dSpV_dT, & ! Partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
75624 dSpV_dS, & ! Partial derivative of specific volume with to salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
75724 netheat_rate, & ! netheat but for dt=1 [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
75824 netsalt_rate, & ! netsalt but for dt=1 (e.g. returns a rate)
759 ! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
76024 netMassInOut_rate, & ! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
76124 mixing_depth, & ! The mixing depth for brine plumes [H ~> m or kg m-2]
76224 total_h ! Total thickness of the water column [H ~> m or kg m-2]
763 real, dimension(SZI_(G), SZK_(GV)) :: &
76424 h2d, & ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
765 ! dz, & ! Layer thicknesses in depth units [Z ~> m]
76624 T2d, & ! A 2-d copy of the layer temperatures [C ~> degC]
76724 pen_TKE_2d, & ! The TKE required to homogenize the heating by shortwave radiation within
768 ! a layer [R Z3 T-2 ~> J m-2]
76924 dSV_dT_2d ! The partial derivative of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
770 real, dimension(SZI_(G)) :: &
77124 netPen_rate ! The surface penetrative shortwave heating rate summed over all bands
772 ! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
773 real, dimension(max(nsw,1),SZI_(G)) :: &
77424 Pen_SW_bnd, & ! The penetrative shortwave heating integrated over a timestep by band
775 ! [C H ~> degC m or degC kg m-2]
77624 Pen_SW_bnd_rate ! The penetrative shortwave heating rate by band
777 ! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
778 real, dimension(max(nsw,1),SZI_(G),SZK_(GV)) :: &
77924 opacityBand ! The opacity (inverse of the exponential absorption length) of each frequency
780 ! band of shortwave radiation in each layer [H-1 ~> m-1 or m2 kg-1]
781 real, dimension(maxGroundings) :: hGrounding ! Thickness added by each grounding event [H ~> m or kg m-2]
782 real :: Temp_in ! The initial temperature of a layer [C ~> degC]
783 real :: Salin_in ! The initial salinity of a layer [S ~> ppt]
784 real :: g_Hconv2 ! A conversion factor for use in the TKE calculation
785 ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
786 real :: GoRho ! g_Earth times a unit conversion factor divided by density
787 ! [Z T-2 R-1 ~> m4 s-2 kg-1]
788 real :: g_conv ! The gravitational acceleration times the conversion factors from non-Boussinesq
789 ! thickness units to mass per units area [R Z2 H-1 T-2 ~> kg m-2 s-2 or m s-2]
790 logical :: calculate_energetics ! If true, calculate the energy required to mix the newly added
791 ! water over the topmost grid cell, assuming that the fluxes of heat and salt
792 ! and rejected brine are initially applied in vanishingly thin layers at the
793 ! top of the layer before being mixed throughout the layer.
794 logical :: calculate_buoyancy ! If true, calculate the surface buoyancy flux.
79524 real :: dK(SZI_(G)) ! Depth of the layer center in thickness units [H ~> m or kg m-2]
79612 real :: A_brine(SZI_(G)) ! Constant [H-(n+1) ~> m-(n+1) or m(2n+2) kg-(n+1)].
797 real :: fraction_left_brine ! Fraction of the brine that has not been applied yet [nondim]
798 real :: plume_fraction ! Fraction of the brine that is applied to a layer [nondim]
799 real :: plume_flux ! Brine flux to move downwards [S H ~> ppt m or ppt kg m-2]
800 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
801 integer :: i, j, is, ie, js, je, k, nz, nb
802 character(len=45) :: mesg
803
80412 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
805
80612 Idt = 1.0 / dt
80712 plume_flux = 0.0
808
80912 calculate_energetics = (present(cTKE) .and. present(dSV_dT) .and. present(dSV_dS))
81012 calculate_buoyancy = present(SkinBuoyFlux)
811105276 if (calculate_buoyancy) SkinBuoyFlux(:,:) = 0.0
8127895712 if (present(cTKE)) cTKE(:,:,:) = 0.0
81312 g_Hconv2 = (GV%g_Earth_Z_T2 * GV%H_to_RZ) * GV%H_to_RZ
81412 EOSdom(:) = EOS_domain(G%HI)
815
816 ! Only apply forcing if fluxes%sw is associated.
81712 if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return
818
81912 EnthalpyConst = 1.0
82012 if (associated(fluxes%heat_content_evap)) EnthalpyConst = 0.0
821
82212 if (calculate_buoyancy) then
8231548 SurfPressure(:) = 0.0
82412 GoRho = GV%g_Earth_Z_T2 / GV%Rho0
825 endif
826
82712 if (CS%do_brine_plume .and. .not.present(MLD_h)) then
828 call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
829 "Brine plume parameterization requires a mixed-layer depth argument,\n"//&
8300 "currently coming from the energetic PBL scheme.")
831 endif
83212 if (CS%do_brine_plume .and. .not.associated(MLD_h)) then
833 call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
8340 "Brine plume parameterization requires an associated mixed-layer depth.")
835 endif
83612 if (CS%do_brine_plume .and. .not. associated(fluxes%salt_left_behind)) then
837 call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
838 "Brine plume parameterization requires DO_BRINE_PLUME\n"//&
8390 "to be turned on in SIS2 as well as MOM6.")
840 endif
841
842 ! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total
843 ! depth of the ocean is vanishing. It does not (yet) handle a value of zero.
844 ! To accommodate vanishing upper layers, we need to allow for an instantaneous
845 ! distribution of forcing over some finite vertical extent. The bulk mixed layer
846 ! code handles this issue properly.
84712 H_limit_fluxes = max(GV%Angstrom_H, GV%H_subroundoff)
848
849 ! diagnostic to see if need to create mass to avoid grounding
85012 if (CS%id_createdH>0) CS%createdH(:,:) = 0.
85112 numberOfGroundings = 0
852
853 !$OMP parallel do default(none) shared(is,ie,js,je,nz,h,tv,nsw,G,GV,US,optics,fluxes, &
854 !$OMP H_limit_fluxes,numberOfGroundings,iGround,jGround,&
855 !$OMP nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing, &
856 !$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, &
857 !$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,&
858 !$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, &
859 !$OMP EnthalpyConst,MLD_h) &
860 !$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &
861 !$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
862 !$OMP IforcingDepthScale,g_conv,dSpV_dT,dSpV_dS, &
863 !$OMP dThickness,dTemp,dSalt,hOld,Ithickness, &
864 !$OMP netMassIn,pres,d_pres,p_lay,dSV_dT_2d, &
865 !$OMP netmassinout_rate,netheat_rate,netsalt_rate, &
866 !$OMP drhodt,drhods,pen_sw_bnd_rate, &
867 !$OMP pen_TKE_2d,Temp_in,Salin_in,RivermixConst, &
868 !$OMP mixing_depth,A_brine,fraction_left_brine, &
869 !$OMP plume_fraction,dK,total_h) &
870 !$OMP firstprivate(SurfPressure,plume_flux)
871732 do j=js,je
872 ! Work in vertical slices for efficiency
873
874 ! Copy state into 2D-slice arrays
8756534720 do k=1,nz ; do i=is,ie
8766480000 h2d(i,k) = h(i,j,k)
8776534000 T2d(i,k) = tv%T(i,j,k)
878 enddo ; enddo
879
880720 if (calculate_energetics) then
881 ! The partial derivatives of specific volume with temperature and
882 ! salinity need to be precalculated to avoid having heating of
883 ! tiny layers give nonsensical values.
884720 if (associated(tv%p_surf)) then
88587120 do i=is,ie ; pres(i) = tv%p_surf(i,j) ; enddo
886 else
8870 do i=is,ie ; pres(i) = 0.0 ; enddo
888 endif
88954720 do k=1,nz
8906534000 do i=is,ie
8916480000 d_pres(i) = (GV%g_Earth * GV%H_to_RZ) * h2d(i,k)
8926480000 p_lay(i) = pres(i) + 0.5*d_pres(i)
8936534000 pres(i) = pres(i) + d_pres(i)
894 enddo
895 call calculate_specific_vol_derivs(T2d(:,k), tv%S(:,j,k), p_lay(:), &
89654000 dSV_dT(:,j,k), dSV_dS(:,j,k), tv%eqn_of_state, EOSdom)
8976534720 do i=is,ie ; dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; enddo
898 enddo
8996966720 pen_TKE_2d(:,:) = 0.0
900 endif
901
902 ! Nothing more is done on this j-slice if there is no buoyancy forcing.
903720 if (.not.associated(fluxes%sw)) cycle
904
905720 if (nsw>0) then
906720 if (GV%Boussinesq .or. (.not.allocated(tv%SpV_avg))) then
907720 call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=GV%H_to_Z)
908 else
909 call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=GV%H_to_RZ, &
9100 SpV_avg=tv%SpV_avg)
911 endif
912 endif
913
914 ! The surface forcing is contained in the fluxes type.
915 ! We aggregate the thermodynamic forcing for a time step into the following:
916 ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step
917 ! = lprec + fprec + vprec + evap + lrunoff + frunoff
918 ! note that lprec generally has sea ice melt/form included.
919 ! netMassOut = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.
920 ! netMassOut < 0 means mass leaves ocean.
921 ! netHeat = heat via surface fluxes [C H ~> degC m or degC kg m-2], excluding the part
922 ! contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
923 ! netSalt = surface salt fluxes [S H ~> ppt m or gSalt m-2]
924 ! Pen_SW_bnd = components to penetrative shortwave radiation split according to bands.
925 ! This field provides that portion of SW from atmosphere that in fact
926 ! enters to the ocean and participates in penetrative SW heating.
927 ! nonpenSW = non-downwelling SW flux, which is absorbed in ocean surface
928 ! (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.
929
930 !----------------------------------------------------------------------------------------
931 !BGR-June 26, 2017{
932 !Temporary action to preserve answers while fixing a bug.
933 ! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns
934 ! the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning
935 ! a second call to extractfluxes1d (causing the diagnostic net_heat to be incorrect).
936 ! Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,
937 ! which means it used the T/S fields after this routine. Therefore, the surface
938 ! buoyancy flux is computed here at the very end of this routine for legacy reasons.
939 ! A few specific notes follow:
940 ! 1) The old method did not included river/calving contributions to heat flux. This
941 ! is kept consistent here via commenting code in the present extractFluxes1d <_rate>
942 ! outputs, but we may reconsider this approach.
943 ! 2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead
944 ! of computing the integrated value (and dividing by dt). Hence the required
945 ! additional outputs from extractFluxes1d.
946 ! *** This is because: A*dt/dt =/= A due to round off.
947 ! 3) The old method computed buoyancy flux after this routine, meaning the returned
948 ! surface fluxes (from extractfluxes1d) must be recorded for use later in the code.
949 ! We could (and maybe should) move that loop up to before the surface fluxes are
950 ! applied, but this will change answers.
951 ! For all these reasons we compute additional values of <_rate> which are preserved
952 ! for the buoyancy flux calculation and reproduce the old answers.
953 ! In the future this needs more detailed investigation to make sure everything is
954 ! consistent and correct. These details should not significantly effect climate,
955 ! but do change answers.
956 !-----------------------------------------------------------------------------------------
957720 if (calculate_buoyancy) then
958 call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
959 H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
960 h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, &
961 Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW, &
962 net_Heat_rate=netheat_rate, net_salt_rate=netsalt_rate, &
963720 netmassinout_rate=netmassinout_rate, pen_sw_bnd_rate=pen_sw_bnd_rate)
964 else
965 call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, &
966 H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, &
967 h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, &
9680 Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW)
969 endif
970 ! ea is for passive tracers
97187120 do i=is,ie
972 ! ea(i,j,1) = netMassInOut(i)
97386400 if (aggregate_FW_forcing) then
97486400 netMassOut(i) = netMassInOut(i)
97586400 netMassIn(i) = 0.
976 else
9770 netMassIn(i) = netMassInOut(i) - netMassOut(i)
978 endif
97987120 if (G%mask2dT(i,j) > 0.0) then
98060168 fluxes%netMassOut(i,j) = netMassOut(i)
98160168 fluxes%netMassIn(i,j) = netMassIn(i)
982 else
98326232 fluxes%netMassOut(i,j) = 0.0
98426232 fluxes%netMassIn(i,j) = 0.0
985 endif
986 enddo
987
988 ! Apply the surface boundary fluxes in three steps:
989 ! A/ update mass, temp, and salinity due to all terms except mass leaving
990 ! ocean (and corresponding outward heat content), and ignoring penetrative SW.
991 ! B/ update mass, salt, temp from mass leaving ocean.
992 ! C/ update temp due to penetrative SW
993720 if (CS%do_brine_plume) then
994 ! Find the plume mixing depth.
9950 do i=is,ie ; total_h(i) = 0.0 ; enddo
9960 do k=1,nz ; do i=is,ie ; total_h(i) = total_h(i) + h(i,j,k) ; enddo ; enddo
9970 do i=is,ie
998 mixing_depth(i) = min( max(MLD_h(i,j) - minimum_forcing_depth, minimum_forcing_depth), &
9990 max(total_h(i), GV%angstrom_h) ) + GV%H_subroundoff
10000 A_brine(i) = (CS%brine_plume_n + 1) / (mixing_depth(i) ** (CS%brine_plume_n + 1))
1001 enddo
1002 endif
1003
100487120 do i=is,ie
100586400 if (G%mask2dT(i,j) > 0.) then
1006
1007 ! A/ Update mass, temp, and salinity due to incoming mass flux.
1008120336 do k=1,1
1009
1010 ! Change in state due to forcing
101160168 dThickness = netMassIn(i) ! Since we are adding mass, we can use all of it
101260168 dTemp = 0.
101360168 dSalt = 0.
1014
1015 ! Update the forcing by the part to be consumed within the present k-layer.
1016 ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
101760168 netMassIn(i) = netMassIn(i) - dThickness
1018 ! This line accounts for the temperature of the mass exchange
101960168 Temp_in = T2d(i,k)
102060168 Salin_in = 0.0
102160168 dTemp = dTemp + dThickness*Temp_in*EnthalpyConst
1022
1023 ! Diagnostics of heat content associated with mass fluxes
102460168 if (.not. associated(fluxes%heat_content_evap)) then
102560168 if (associated(fluxes%heat_content_massin)) &
1026 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
102760168 T2d(i,k) * max(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
102860168 if (associated(fluxes%heat_content_massout)) &
1029 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
103060168 T2d(i,k) * min(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
103160168 if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
103260168 T2d(i,k) * dThickness * GV%H_to_RZ
1033 endif
1034
1035 ! Determine the energetics of river mixing before updating the state.
103660168 if (calculate_energetics .and. associated(fluxes%lrunoff) .and. CS%do_rivermix) then
1037 ! Here we add an additional source of TKE to the mixed layer where river
1038 ! is present to simulate unresolved estuaries. The TKE input, TKE_river in
1039 ! [Z3 T-3 ~> m3 s-3], is diagnosed as follows:
1040 ! TKE_river = 0.5*rivermix_depth*g*(1/rho)*drho_ds*
1041 ! River*(Samb - Sriver) = CS%mstar*U_star^3
1042 ! where River is in units of [Z T-1 ~> m s-1].
1043 ! Samb = Ambient salinity at the mouth of the estuary
1044 ! rivermix_depth = The prescribed depth over which to mix river inflow
1045 ! drho_ds = The derivative of density with salt at the ambient surface salinity.
1046 ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
10470 if (GV%Boussinesq) then
10480 RivermixConst = -0.5*(CS%rivermix_depth*dt) * GV%g_Earth_Z_T2 * GV%Rho0
10490 elseif (allocated(tv%SpV_avg)) then
10500 RivermixConst = -0.5*(CS%rivermix_depth*dt) * GV%g_Earth_Z_T2 / tv%SpV_avg(i,j,1)
1051 else
10520 RivermixConst = -0.5*(CS%rivermix_depth*dt) * GV%Rho0 * GV%g_Earth_Z_T2
1053 endif
1054 cTKE(i,j,k) = cTKE(i,j,k) + max(0.0, RivermixConst*dSV_dS(i,j,1) * &
1055 ((fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) + &
10560 (fluxes%lrunoff_glc(i,j) + fluxes%frunoff_glc(i,j))) * tv%S(i,j,1))
1057 endif
1058
1059 ! Update state
106060168 hOld = h2d(i,k) ! Keep original thickness in hand
106160168 h2d(i,k) = h2d(i,k) + dThickness ! New thickness
1062120336 if (h2d(i,k) > 0.0) then
106360168 if (calculate_energetics .and. (dThickness > 0.)) then
1064 ! Calculate the energy required to mix the newly added water over
1065 ! the topmost grid cell.
1066 cTKE(i,j,k) = cTKE(i,j,k) + 0.5*g_Hconv2*(hOld*dThickness) * &
10670 ((T2d(i,k) - Temp_in) * dSV_dT(i,j,k) + (tv%S(i,j,k) - Salin_in) * dSV_dS(i,j,k))
1068 endif
106960168 Ithickness = 1.0/h2d(i,k) ! Inverse new thickness
1070 ! The "if"s below avoid changing T/S by roundoff unnecessarily
107160168 if (dThickness /= 0. .or. dTemp /= 0.) T2d(i,k) = (hOld*T2d(i,k) + dTemp)*Ithickness
107260168 if (dThickness /= 0. .or. dSalt /= 0.) tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt)*Ithickness
1073
1074 endif
1075
1076 enddo ! k=1,1
1077
1078 ! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.
107960168 fraction_left_brine = 1.0
10804572768 do k=1,nz
1081 ! Place forcing into this layer if this layer has nontrivial thickness.
1082 ! For layers thin relative to 1/IforcingDepthScale, then distribute
1083 ! forcing into deeper layers.
10844512600 IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) )
1085 ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
10864512600 fractionOfForcing = min(1.0, h2d(i,k)*IforcingDepthScale)
1087
1088 ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
1089 ! limit the forcing applied to this cell, leaving the remaining forcing to
1090 ! be distributed downwards.
10914512600 if (-fractionOfForcing*netMassOut(i) > evap_CFL_limit*h2d(i,k)) then
10920 fractionOfForcing = -evap_CFL_limit*h2d(i,k)/netMassOut(i)
1093 endif
1094
10954512600 if (CS%do_brine_plume .and. associated(fluxes%salt_left_behind)) then
10960 if (fluxes%salt_left_behind(i,j) > 0 .and. fraction_left_brine > 0.0) then
1097 ! Place forcing into this layer by depth for brine plume parameterization.
10980 if (k == 1) then
10990 dK(i) = 0.5 * h(i,j,k) ! Depth of center of layer K
11000 plume_flux = - (1000.0*US%ppt_to_S * (CS%plume_strength * fluxes%salt_left_behind(i,j))) * GV%RZ_to_H
11010 plume_fraction = 1.0
1102 else
11030 dK(i) = dK(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) ! Depth of center of layer K
11040 plume_flux = 0.0
1105 endif
11060 if (dK(i) <= mixing_depth(i) .and. fraction_left_brine > 0.0) then
11070 plume_fraction = min(fraction_left_brine, (A_brine(i) * dK(i)**CS%brine_plume_n) * h(i,j,k))
1108 else
11090 IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) )
1110 ! plume_fraction = fraction_left_brine, unless h2d is less than IforcingDepthScale.
11110 plume_fraction = min(fraction_left_brine, h2d(i,k)*IforcingDepthScale)
1112 endif
11130 fraction_left_brine = fraction_left_brine - plume_fraction
1114 plume_flux = plume_flux + plume_fraction * (1000.0*US%ppt_to_S * (CS%plume_strength * &
11150 fluxes%salt_left_behind(i,j))) * GV%RZ_to_H
1116 else
11170 plume_flux = 0.0
1118 endif
1119 endif
1120
1121 ! Change in state due to forcing
1122
11234512600 dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) )
11244512600 dTemp = fractionOfForcing*netHeat(i)
11254512600 dSalt = max( fractionOfForcing*netSalt(i), -CS%dSalt_frac_max * h2d(i,k) * tv%S(i,j,k))
1126
1127 ! Update the forcing by the part to be consumed within the present k-layer.
1128 ! If fractionOfForcing = 1, then new netMassOut vanishes.
11294512600 netMassOut(i) = netMassOut(i) - dThickness
11304512600 netHeat(i) = netHeat(i) - dTemp
11314512600 netSalt(i) = netSalt(i) - dSalt
1132
1133 ! This line accounts for the temperature of the mass exchange
11344512600 dTemp = dTemp + dThickness*T2d(i,k)*EnthalpyConst
1135
1136 ! Diagnostics of heat content associated with mass fluxes
11374512600 if (.not. associated(fluxes%heat_content_evap)) then
11384512600 if (associated(fluxes%heat_content_massin)) &
1139 fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
11404512600 T2d(i,k) * max(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
11414512600 if (associated(fluxes%heat_content_massout)) &
1142 fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
11434512600 T2d(i,k) * min(0.,dThickness) * GV%H_to_RZ * tv%C_p * Idt
11444512600 if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
11454512600 T2d(i,k) * dThickness * GV%H_to_RZ
1146 endif
1147
1148 ! Update state by the appropriate increment.
11494512600 hOld = h2d(i,k) ! Keep original thickness in hand
11504512600 h2d(i,k) = h2d(i,k) + dThickness ! New thickness
1151
11524572768 if (h2d(i,k) > 0.) then
11534512600 if (calculate_energetics) then
1154 ! Calculate the energy required to mix the newly added water over the topmost grid
1155 ! cell, assuming that the fluxes of heat and salt and rejected brine are initially
1156 ! applied in vanishingly thin layers at the top of the layer before being mixed
1157 ! throughout the layer. Note that dThickness is always <= 0 here, and that
1158 ! negative cTKE is a deficit that will need to be filled later.
1159 cTKE(i,j,k) = cTKE(i,j,k) - (0.5*h2d(i,k)*g_Hconv2) * &
1160 ((dTemp - dthickness*T2d(i,k)) * dSV_dT(i,j,k) + &
11614512600 (dSalt - dthickness*tv%S(i,j,k)) * dSV_dS(i,j,k))
1162 endif
11634512600 Ithickness = 1.0/h2d(i,k) ! Inverse of new thickness
11644512600 T2d(i,k) = (hOld*T2d(i,k) + dTemp)*Ithickness
11654512600 tv%S(i,j,k) = (hOld*tv%S(i,j,k) + dSalt + plume_flux)*Ithickness
11660 elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
11670 call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (h<0)')
11680 write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j)
11690 write(0,*) 'applyBoundaryFluxesInOut(): netT,netS,netH=', &
11700 US%C_to_degC*netHeat(i), US%S_to_ppt*netSalt(i), netMassInOut(i)
11710 write(0,*) 'applyBoundaryFluxesInOut(): dT,dS,dH=', &
11720 US%C_to_degC*dTemp, US%S_to_ppt*dSalt, dThickness
11730 write(0,*) 'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hOld,h2d(i,k),k
1174 call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
11750 "Complete mass loss in column!")
1176 endif
1177
1178 enddo ! k
1179
1180 ! Check if trying to apply fluxes over land points
118126232 elseif ((abs(netHeat(i)) + abs(netSalt(i)) + abs(netMassIn(i)) + abs(netMassOut(i))) > 0.) then
1182
11830 if (.not. CS%ignore_fluxes_over_land) then
11840 call forcing_SinglePointPrint(fluxes,G,i,j,'applyBoundaryFluxesInOut (land)')
11850 write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j)
11860 write(0,*) 'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
11870 US%C_to_degC*netHeat(i), US%S_to_ppt*netSalt(i), netMassIn(i), netMassOut(i)
1188
1189 call MOM_error(FATAL, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
11900 "Mass loss over land?")
1191 endif
1192
1193 endif
1194
1195 ! If anything remains after the k-loop, then we have grounded out, which is a problem.
119687120 if (netMassIn(i)+netMassOut(i) /= 0.0) then
1197!$OMP critical
11980 numberOfGroundings = numberOfGroundings +1
11990 if (numberOfGroundings<=maxGroundings) then
12000 iGround(numberOfGroundings) = i ! Record i,j location of event for
12010 jGround(numberOfGroundings) = j ! warning message
12020 hGrounding(numberOfGroundings) = netMassIn(i)+netMassOut(i)
1203 endif
1204!$OMP end critical
12050 if (CS%id_createdH>0) CS%createdH(i,j) = CS%createdH(i,j) - (netMassIn(i)+netMassOut(i))/dt
1206 endif
1207
1208 enddo ! i
1209
1210 ! Step C/ in the application of fluxes
1211 ! Heat by the convergence of penetrating SW.
1212 ! SW penetrative heating uses the updated thickness from above.
1213
1214 ! Save temperature before increment with SW heating
1215 ! and initialize CS%penSWflux_diag to zero.
1216720 if (CS%id_penSW_diag > 0 .or. CS%id_penSWflux_diag > 0) then
12170 do k=1,nz ; do i=is,ie
12180 CS%penSW_diag(i,j,k) = T2d(i,k)
12190 CS%penSWflux_diag(i,j,k) = 0.0
1220 enddo ; enddo
12210 k=nz+1 ; do i=is,ie
12220 CS%penSWflux_diag(i,j,k) = 0.0
1223 enddo
1224 endif
1225
1226720 if (calculate_energetics) then
1227 call absorbRemainingSW(G, GV, US, h2d, opacityBand, nsw, optics, j, dt, H_limit_fluxes, &
1228720 .false., .true., T2d, Pen_SW_bnd, TKE=pen_TKE_2d, dSV_dT=dSV_dT_2d)
1229720 k = 1 ! For setting break-points.
12306534720 do k=1,nz ; do i=is,ie
12316534000 cTKE(i,j,k) = cTKE(i,j,k) + pen_TKE_2d(i,k)
1232 enddo ; enddo
1233 else
1234 call absorbRemainingSW(G, GV, US, h2d, opacityBand, nsw, optics, j, dt, H_limit_fluxes, &
12350 .false., .true., T2d, Pen_SW_bnd)
1236 endif
1237
1238
1239 ! Step D/ copy updated thickness and temperature
1240 ! 2d slice now back into model state.
12416534720 do k=1,nz ; do i=is,ie
12426480000 h(i,j,k) = h2d(i,k)
12436534000 tv%T(i,j,k) = T2d(i,k)
1244 enddo ; enddo
1245
1246 ! Diagnose heating [Q R Z T-1 ~> W m-2] applied to a grid cell from SW penetration
1247 ! Also diagnose the penetrative SW heat flux at base of layer.
1248720 if (CS%id_penSW_diag > 0 .or. CS%id_penSWflux_diag > 0) then
1249
1250 ! convergence of SW into a layer
12510 do k=1,nz ; do i=is,ie
1252 ! Note that the units of penSW_diag change here, from [C ~> degC] to [Q R Z T-1 ~> W m-2].
12530 CS%penSW_diag(i,j,k) = (T2d(i,k)-CS%penSW_diag(i,j,k))*h(i,j,k) * Idt * tv%C_p * GV%H_to_RZ
1254 enddo ; enddo
1255
1256 ! Perform a cumulative sum upwards from bottom to
1257 ! diagnose penetrative SW flux at base of tracer cell.
1258 ! CS%penSWflux_diag(i,j,k=1) is penetrative shortwave at top of ocean.
1259 ! CS%penSWflux_diag(i,j,k=kbot+1) is zero, since assume no SW penetrates rock.
1260 ! CS%penSWflux_diag = rsdo and CS%penSW_diag = rsdoabsorb
1261 ! rsdoabsorb(k) = rsdo(k) - rsdo(k+1), so that rsdo(k) = rsdo(k+1) + rsdoabsorb(k)
12620 if (CS%id_penSWflux_diag > 0) then
12630 do k=nz,1,-1 ; do i=is,ie
12640 CS%penSWflux_diag(i,j,k) = CS%penSW_diag(i,j,k) + CS%penSWflux_diag(i,j,k+1)
1265 enddo ; enddo
1266 endif
1267
1268 endif
1269
1270 ! Fill CS%nonpenSW_diag
1271720 if (CS%id_nonpenSW_diag > 0) then
12720 do i=is,ie
12730 CS%nonpenSW_diag(i,j) = nonpenSW(i) * Idt * tv%C_p * GV%H_to_RZ
1274 enddo
1275 endif
1276
1277 ! BGR: Get buoyancy flux to return for ePBL
1278 ! We want the rate, so we use the rate values returned from extractfluxes1d.
1279 ! Note that the *dt values could be divided by dt here, but
1280 ! 1) Answers will change due to round-off
1281 ! 2) Be sure to save their values BEFORE fluxes are used.
1282732 if (Calculate_Buoyancy) then
128392880 netPen_rate(:) = 0.0
1284 ! Sum over bands and attenuate as a function of depth.
1285 ! netPen_rate is the netSW as a function of depth, but only the surface value is used here,
1286 ! in which case the values of dt, h, optics and H_limit_fluxes are irrelevant. Consider
1287 ! writing a shorter and simpler variant to handle this very limited case.
1288 ! Find the vertical distances across layers.
1289 ! call thickness_to_dz(h, tv, dz, j, G, GV)
1290 ! call sumSWoverBands(G, GV, US, h2d, dz, optics_nbands(optics), optics, j, dt, &
1291 ! H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
1292173520 do i=is,ie ; do nb=1,nsw ; netPen_rate(i) = netPen_rate(i) + pen_SW_bnd_rate(nb,i) ; enddo ; enddo
1293
1294 ! 1. Adjust netSalt to reflect dilution effect of FW flux
1295 ! 2. Add in the SW heating for purposes of calculating the net
1296 ! surface buoyancy flux affecting the top layer.
1297 ! 3. Convert to a buoyancy flux, excluding penetrating SW heating
1298 ! BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
129987120 if (associated(tv%p_surf)) then ; do i=is,ie ; SurfPressure(i) = tv%p_surf(i,j) ; enddo ; endif
1300
1301720 if ((.not.GV%Boussinesq) .and. (.not.GV%semi_Boussinesq)) then
13020 g_conv = GV%g_Earth_Z_T2 * GV%H_to_RZ
1303
1304 ! Specific volume derivatives
1305 call calculate_specific_vol_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, dSpV_dT, dSpV_dS, &
13060 tv%eqn_of_state, EOS_domain(G%HI))
13070 do i=is,ie
1308 SkinBuoyFlux(i,j) = g_conv * &
1309 (dSpV_dS(i) * ( netSalt_rate(i) - tv%S(i,j,1)*netMassInOut_rate(i)) + &
13100 dSpV_dT(i) * ( netHeat_rate(i) + netPen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
1311 enddo
1312 else
1313 ! Density derivatives
1314 call calculate_density_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, dRhodT, dRhodS, &
1315720 tv%eqn_of_state, EOSdom)
131687120 do i=is,ie
1317 SkinBuoyFlux(i,j) = - GoRho * GV%H_to_Z * &
1318 (dRhodS(i) * ( netSalt_rate(i) - tv%S(i,j,1)*netMassInOut_rate(i)) + &
131987120 dRhodT(i) * ( netHeat_rate(i) + netPen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
1320 enddo
1321 endif
1322 endif
1323
1324 enddo ! j-loop finish
1325
1326 ! Post the diagnostics
132712 if (CS%id_createdH > 0) call post_data(CS%id_createdH , CS%createdH , CS%diag)
132812 if (CS%id_penSW_diag > 0) call post_data(CS%id_penSW_diag , CS%penSW_diag , CS%diag)
132912 if (CS%id_penSWflux_diag > 0) call post_data(CS%id_penSWflux_diag, CS%penSWflux_diag, CS%diag)
133012 if (CS%id_nonpenSW_diag > 0) call post_data(CS%id_nonpenSW_diag , CS%nonpenSW_diag , CS%diag)
1331
1332! The following check will be ignored if ignore_fluxes_over_land = true
133312 if ((numberOfGroundings > 0) .and. .not.CS%ignore_fluxes_over_land) then
13340 do i = 1, min(numberOfGroundings, maxGroundings)
13350 call forcing_SinglePointPrint(fluxes,G,iGround(i),jGround(i),'applyBoundaryFluxesInOut (grounding)')
13360 write(mesg(1:45),'(3es15.3)') G%geoLonT( iGround(i), jGround(i) ), &
13370 G%geoLatT( iGround(i), jGround(i)), hGrounding(i)*GV%H_to_m
1338 call MOM_error(WARNING, "MOM_diabatic_aux.F90, applyBoundaryFluxesInOut(): "//&
13390 "Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
1340 enddo
1341
13420 if (numberOfGroundings - maxGroundings > 0) then
13430 write(mesg, '(I0)') numberOfGroundings - maxGroundings
1344 call MOM_error(WARNING, "MOM_diabatic_aux:F90, applyBoundaryFluxesInOut(): "//&
13450 trim(mesg) // " groundings remaining")
1346 endif
1347 endif
1348
134948end subroutine applyBoundaryFluxesInOut
1350
1351!> This subroutine initializes the parameters and control structure of the diabatic_aux module.
13521subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm, use_ePBL)
1353 type(time_type), target, intent(in) :: Time !< The current model time.
1354 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1355 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1356 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1357 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1358 type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output
1359 type(diabatic_aux_CS), pointer :: CS !< A pointer to the control structure for the
1360 !! diabatic_aux module, which is initialized here.
1361 logical, intent(in) :: useALEalgorithm !< If true, use the ALE algorithm rather
1362 !! than layered mode.
1363 logical, intent(in) :: use_ePBL !< If true, use the implicit energetics planetary
1364 !! boundary layer scheme to determine the diffusivity
1365 !! in the surface boundary layer.
1366
1367 ! This "include" declares and sets the variable "version".
1368# include "version_variable.h"
1369 character(len=40) :: mdl = "MOM_diabatic_aux" ! This module's name.
1370 character(len=200) :: inputdir ! The directory where NetCDF input files
1371 character(len=240) :: chl_filename ! A file from which chl_a concentrations are to be read.
1372 character(len=128) :: chl_file ! Data containing chl_a concentrations. Used
1373 ! when var_pen_sw is defined and reading from file.
1374 character(len=32) :: chl_varname ! Name of chl_a variable in chl_file.
1375 logical :: use_temperature ! True if thermodynamics are enabled.
1376 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
13771 isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
13781 IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
1379
13801 if (associated(CS)) then
1381 call MOM_error(WARNING, "diabatic_aux_init called with an "// &
13820 "associated control structure.")
13830 return
1384 else
13851 allocate(CS)
1386 endif
1387
13881 CS%diag => diag
13891 CS%Time => Time
1390
1391! Set default, read and log parameters
1392 call log_version(param_file, mdl, version, &
13931 "The following parameters are used for auxiliary diabatic processes.")
1394
1395 call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
13961 "If true, temperature and salinity are used as state variables.", default=.true.)
1397
1398 call get_param(param_file, mdl, "RECLAIM_FRAZIL", CS%reclaim_frazil, &
1399 "If true, try to use any frazil heat deficit to cool any "//&
1400 "overlying layers down to the freezing point, thereby "//&
1401 "avoiding the creation of thin ice when the SST is above "//&
14021 "the freezing point.", default=.true., do_not_log=.not.use_temperature)
1403 call get_param(param_file, mdl, "SALT_EXTRACTION_LIMIT", CS%dSalt_frac_max, &
1404 "An upper limit on the fraction of the salt in a layer that can be lost to the "//&
1405 "net surface salt fluxes within a timestep.", &
14061 units="nondim", default=0.9999, do_not_log=.not.use_temperature)
14071 CS%dSalt_frac_max = max(min(CS%dSalt_frac_max, 1.0), 0.0)
1408 call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", CS%pressure_dependent_frazil, &
1409 "If true, use a pressure dependent freezing temperature "//&
1410 "when making frazil. The default is false, which will be "//&
1411 "faster but is inappropriate with ice-shelf cavities.", &
14121 default=.false., do_not_log=.not.use_temperature)
1413
14141 if (use_ePBL) then
1415 call get_param(param_file, mdl, "IGNORE_FLUXES_OVER_LAND", CS%ignore_fluxes_over_land,&
1416 "If true, the model does not check if fluxes are being applied "//&
1417 "over land points. This is needed when the ocean is coupled "//&
1418 "with ice shelves and sea ice, since the sea ice mask needs to "//&
1419 "be different than the ocean mask to avoid sea ice formation "//&
14201 "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
1421 call get_param(param_file, mdl, "DO_RIVERMIX", CS%do_rivermix, &
1422 "If true, apply additional mixing wherever there is "//&
1423 "runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
14241 "if the ocean is that deep.", default=.false.)
14251 if (CS%do_rivermix) &
1426 call get_param(param_file, mdl, "RIVERMIX_DEPTH", CS%rivermix_depth, &
1427 "The depth to which rivers are mixed if DO_RIVERMIX is "//&
14280 "defined.", units="m", default=0.0, scale=US%m_to_Z)
1429 else
14300 CS%do_rivermix = .false. ; CS%rivermix_depth = 0.0 ; CS%ignore_fluxes_over_land = .false.
1431 endif
1432
14331 if (GV%nkml == 0) then
1434 call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", CS%use_river_heat_content, &
1435 "If true, use the fluxes%runoff_Hflx field to set the "//&
1436 "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
14371 default=.false., do_not_log=.not.use_temperature)
1438 call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", CS%use_calving_heat_content, &
1439 "If true, use the fluxes%calving_Hflx field to set the "//&
1440 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
14411 default=.false., do_not_log=.not.use_temperature)
1442 else
14430 CS%use_river_heat_content = .false.
14440 CS%use_calving_heat_content = .false.
1445 endif
1446
1447 call get_param(param_file, mdl, "DO_BRINE_PLUME", CS%do_brine_plume, &
1448 "If true, use a brine plume parameterization from "//&
14491 "Nguyen et al., 2009.", default=.false.)
1450 call get_param(param_file, mdl, "BRINE_PLUME_EXPONENT", CS%brine_plume_n, &
1451 "If using the brine plume parameterization, set the integer exponent.", &
14521 default=5, do_not_log=.not.CS%do_brine_plume)
1453 call get_param(param_file, mdl, "BRINE_PLUME_FRACTION", CS%plume_strength, &
1454 "Fraction of the available brine to mix down using the brine plume parameterization.", &
14551 units="nondim", default=1.0, do_not_log=.not.CS%do_brine_plume)
1456
14571 if (useALEalgorithm) then
1458 CS%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
1459 Time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
14601 "m s-1", conversion=GV%H_to_m*US%s_to_T)
14611 if (CS%id_createdH>0) allocate(CS%createdH(isd:ied,jsd:jed))
1462
1463 ! diagnostic for heating of a grid cell from convergence of SW heat into the cell
1464 CS%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb', &
1465 diag%axesTL, Time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
1466 'W m-2', conversion=US%QRZ_T_to_W_m2, &
14671 standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer', v_extensive=.true.)
1468
1469 ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)
1470 ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).
1471 CS%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo', &
1472 diag%axesTi, Time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
14731 'W m-2', conversion=US%QRZ_T_to_W_m2, standard_name='downwelling_shortwave_flux_in_sea_water')
1474
1475 ! need both arrays for the SW diagnostics (one for flux, one for convergence)
14761 if (CS%id_penSW_diag>0 .or. CS%id_penSWflux_diag>0) then
14770 allocate(CS%penSW_diag(isd:ied,jsd:jed,nz), source=0.0)
14780 allocate(CS%penSWflux_diag(isd:ied,jsd:jed,nz+1), source=0.0)
1479 endif
1480
1481 ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)
1482 CS%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW', &
1483 diag%axesT1, Time, &
1484 'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
1485 'W m-2', conversion=US%QRZ_T_to_W_m2, &
14861 standard_name='nondownwelling_shortwave_flux_in_sea_water')
14871 if (CS%id_nonpenSW_diag > 0) then
14880 allocate(CS%nonpenSW_diag(isd:ied,jsd:jed), source=0.0)
1489 endif
1490 endif
1491
14921 if (use_temperature) then
1493 call get_param(param_file, mdl, "VAR_PEN_SW", CS%var_pen_sw, &
1494 "If true, use one of the CHL_A schemes specified by "//&
1495 "OPACITY_SCHEME to determine the e-folding depth of "//&
14961 "incoming short wave radiation.", default=.false.)
14971 if (CS%var_pen_sw) then
1498
1499 call get_param(param_file, mdl, "CHL_FROM_FILE", CS%chl_from_file, &
15000 "If true, chl_a is read from a file.", default=.true.)
15010 if (CS%chl_from_file) then
15020 call time_interp_external_init()
1503
15040 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1505 call get_param(param_file, mdl, "CHL_FILE", chl_file, &
1506 "CHL_FILE is the file containing chl_a concentrations in "//&
1507 "the variable CHL_A. It is used when VAR_PEN_SW and "//&
15080 "CHL_FROM_FILE are true.", fail_if_missing=.true.)
15090 chl_filename = trim(slasher(inputdir))//trim(chl_file)
15100 call log_param(param_file, mdl, "INPUTDIR/CHL_FILE", chl_filename)
1511 call get_param(param_file, mdl, "CHL_VARNAME", chl_varname, &
15120 "Name of CHL_A variable in CHL_FILE.", default='CHL_A')
15130 if (modulo(G%Domain%turns, 4) /= 0) then
15140 CS%sbc_chl = init_external_field(chl_filename, trim(chl_varname), MOM_domain=G%Domain%domain_in)
1515 else
15160 CS%sbc_chl = init_external_field(chl_filename, trim(chl_varname), MOM_domain=G%Domain)
1517 endif
1518 endif
1519
1520 CS%id_chl = register_diag_field('ocean_model', 'Chl_opac', diag%axesT1, Time, &
15210 'Surface chlorophyll A concentration used to find opacity', 'mg m-3')
1522 endif
1523 endif
1524
15251 id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=CLOCK_ROUTINE)
15261 id_clock_frazil = cpu_clock_id('(Ocean frazil)', grain=CLOCK_ROUTINE)
1527
1528end subroutine diabatic_aux_init
1529
1530!> This subroutine initializes the control structure and any related memory
1531!! for the diabatic_aux module.
15321subroutine diabatic_aux_end(CS)
1533 type(diabatic_aux_CS), pointer :: CS !< The control structure returned by a previous
1534 !! call to diabatic_aux_init; it is deallocated here.
1535
15361 if (.not.associated(CS)) return
1537
15381 if (CS%id_createdH >0) deallocate(CS%createdH)
15391 if (CS%id_penSW_diag >0) deallocate(CS%penSW_diag)
15401 if (CS%id_penSWflux_diag >0) deallocate(CS%penSWflux_diag)
15411 if (CS%id_nonpenSW_diag >0) deallocate(CS%nonpenSW_diag)
1542
15431 if (associated(CS)) deallocate(CS)
1544
1545end subroutine diabatic_aux_end
1546
1547!> \namespace mom_diabatic_aux
1548!!
1549!! This module contains subroutines that apply various diabatic processes. Usually these
1550!! subroutines are called from the MOM_diabatic module. All of these routines use appropriate
1551!! limiters or logic to work properly with arbitrary layer thicknesses (including massless layers)
1552!! and an arbitrarily large timestep.
1553!!
1554!! The subroutine make_frazil facilitates the formation of frazil ice when the ocean water
1555!! drops below the in situ freezing point by heating the water to the freezing point and
1556!! accumulating the required heat for exchange with the sea-ice module.
1557!!
1558!! The subroutine adjust_salt adds salt as necessary to keep the salinity above a
1559!! specified minimum value, and keeps track of the cumulative additions. If the minimum
1560!! salinity is the natural value of 0, this routine should never do anything.
1561!!
1562!! The subroutine differential_diffuse_T_S solves a pair of tridiagonal equations for
1563!! the diffusion of temperatures and salinities with differing diffusivities.
1564!!
1565!! The subroutine triDiagTS solves a tridiagonal equations for the evolution of temperatures
1566!! and salinities due to net entrainment by layers and a diffusion with the same diffusivity.
1567!!
1568!! The subroutine triDiagTS_Eulerian solves a tridiagonal equations for the evolution of
1569!! temperatures and salinities due to diffusion with the same diffusivity, but no net entrainment.
1570!!
1571!! The subroutine find_uv_at_h interpolates velocities to thickness points, optionally also
1572!! using tridiagonal equations to solve for the impacts of net entrainment or mixing of
1573!! momentum between layers.
1574!!
1575!! The subroutine set_pen_shortwave determines the optical properties of the water column and
1576!! the net shortwave fluxes, and stores them in the optics type, working via calls to set_opacity.
1577!!
1578!! The subroutine applyBoundaryFluxesInOut updates the layer thicknesses, temperatures and
1579!! salinities due to the application of the surface forcing. It may also calculate the implied
1580!! turbulent kinetic energy requirements for this forcing to be mixed over the model's finite
1581!! vertical resolution in the surface layers.
15820end module MOM_diabatic_aux