MOM6
MOM_opacity.F90
1 !> Routines used to calculate the opacity of the ocean.
2 module mom_opacity
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data
7 use mom_diag_mediator, only : query_averaging_enabled, register_diag_field
8 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
9 use mom_file_parser, only : get_param, log_param, log_version, param_file_type
10 use mom_string_functions, only : uppercase
11 use mom_grid, only : ocean_grid_type
12 use mom_unit_scaling, only : unit_scale_type
13 use mom_variables, only : thermo_var_ptrs
14 use mom_verticalgrid, only : verticalgrid_type
15 
16 implicit none ; private
17 
18 #include <MOM_memory.h>
19 
20 public set_opacity, opacity_init, opacity_end, opacity_manizza, opacity_morel
21 public extract_optics_slice, extract_optics_fields, optics_nbands
22 public absorbremainingsw, sumswoverbands
23 
24 !> This type is used to store information about ocean optical properties
25 type, public :: optics_type
26  integer :: nbands !< The number of penetrating bands of SW radiation
27 
28  real, pointer, dimension(:,:,:,:) :: opacity_band => null() !< SW optical depth per unit thickness [m-1]
29  !! The number of radiation bands is most rapidly varying (first) index.
30 
31  real, pointer, dimension(:,:,:) :: sw_pen_band => null() !< shortwave radiation [W m-2] at the surface
32  !! in each of the nbands bands that penetrates beyond the surface.
33  !! The most rapidly varying dimension is the band.
34 
35  real, pointer, dimension(:) :: &
36  min_wavelength_band => null(), & !< The minimum wavelength in each band of penetrating shortwave radiation [nm]
37  max_wavelength_band => null() !< The maximum wavelength in each band of penetrating shortwave radiation [nm]
38 
39  real :: pensw_flux_absorb !< A heat flux that is small enough to be completely absorbed in the next
40  !! sufficiently thick layer [H degC T-1 ~> degC m s-1 or degC kg m-2 s-1].
41  real :: pensw_absorb_invlen !< The inverse of the thickness that is used to absorb the remaining
42  !! shortwave heat flux when it drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2].
43  logical :: answers_2018 !< If true, use the order of arithmetic and expressions that recover the
44  !! answers from the end of 2018. Otherwise, use updated and more robust
45  !! forms of the same expressions.
46 
47 end type optics_type
48 
49 !> The control structure with paramters for the MOM_opacity module
50 type, public :: opacity_cs ; private
51  logical :: var_pen_sw !< If true, use one of the CHL_A schemes (specified by OPACITY_SCHEME) to
52  !! determine the e-folding depth of incoming shortwave radiation.
53  integer :: opacity_scheme !< An integer indicating which scheme should be used to translate
54  !! water properties into the opacity (i.e., the e-folding depth) and
55  !! (perhaps) the number of bands of penetrating shortwave radiation to use.
56  real :: pen_sw_scale !< The vertical absorption e-folding depth of the
57  !! penetrating shortwave radiation [m].
58  real :: pen_sw_scale_2nd !< The vertical absorption e-folding depth of the
59  !! (2nd) penetrating shortwave radiation [m].
60  real :: sw_1st_exp_ratio !< Ratio for 1st exp decay in Two Exp decay opacity
61  real :: pen_sw_frac !< The fraction of shortwave radiation that is
62  !! penetrating with a constant e-folding approach.
63  real :: blue_frac !< The fraction of the penetrating shortwave
64  !! radiation that is in the blue band [nondim].
65  real :: opacity_land_value !< The value to use for opacity over land [m-1].
66  !! The default is 10 m-1 - a value for muddy water.
67  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
68  !! regulate the timing of diagnostic output.
69 
70  !>@{ Diagnostic IDs
71  integer :: id_sw_pen = -1, id_sw_vis_pen = -1
72  integer, pointer :: id_opacity(:) => null()
73  !!@}
74 end type opacity_cs
75 
76 !>@{ Coded integers to specify the opacity scheme
77 integer, parameter :: no_scheme = 0, manizza_05 = 1, morel_88 = 2, single_exp = 3, double_exp = 4
78 !!@}
79 
80 character*(10), parameter :: manizza_05_string = "MANIZZA_05" !< String to specify the opacity scheme
81 character*(10), parameter :: morel_88_string = "MOREL_88" !< String to specify the opacity scheme
82 character*(10), parameter :: single_exp_string = "SINGLE_EXP" !< String to specify the opacity scheme
83 character*(10), parameter :: double_exp_string = "DOUBLE_EXP" !< String to specify the opacity scheme
84 
85 real, parameter :: op_diag_len = 1e-10 !< Lengthscale L used to remap opacity
86  !! from op to 1/L * tanh(op * L)
87 
88 contains
89 
90 !> This sets the opacity of sea water based based on one of several different schemes.
91 subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
92  G, GV, CS, chl_2d, chl_3d)
93  type(optics_type), intent(inout) :: optics !< An optics structure that has values
94  !! set based on the opacities.
95  real, dimension(:,:), pointer :: sw_total !< Total shortwave flux into the ocean [W m-2]
96  real, dimension(:,:), pointer :: sw_vis_dir !< Visible, direct shortwave into the ocean [W m-2]
97  real, dimension(:,:), pointer :: sw_vis_dif !< Visible, diffuse shortwave into the ocean [W m-2]
98  real, dimension(:,:), pointer :: sw_nir_dir !< Near-IR, direct shortwave into the ocean [W m-2]
99  real, dimension(:,:), pointer :: sw_nir_dif !< Near-IR, diffuse shortwave into the ocean [W m-2]
100  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
101  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
102  type(opacity_cs), pointer :: cs !< The control structure earlier set up by
103  !! opacity_init.
104  real, dimension(SZI_(G),SZJ_(G)), &
105  optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentractions[mg m-3]
106  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
107  optional, intent(in) :: chl_3d !< The chlorophyll-A concentractions of each layer [mg m-3]
108 
109  ! Local variables
110  integer :: i, j, k, n, is, ie, js, je, nz
111  real :: inv_sw_pen_scale ! The inverse of the e-folding scale [m-1].
112  real :: inv_nbands ! The inverse of the number of bands of penetrating
113  ! shortwave radiation.
114  logical :: call_for_surface ! if horizontal slice is the surface layer
115  real :: tmp(szi_(g),szj_(g),szk_(gv)) ! A 3-d temporary array.
116  real :: chl(szi_(g),szj_(g),szk_(gv)) ! The concentration of chlorophyll-A [mg m-3].
117  real :: pen_sw_tot(szi_(g),szj_(g)) ! The penetrating shortwave radiation
118  ! summed across all bands [W m-2].
119  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
120 
121  if (.not. associated(cs)) call mom_error(fatal, "set_opacity: "// &
122  "Module must be initialized via opacity_init before it is used.")
123 
124  if (present(chl_2d) .or. present(chl_3d)) then
125  ! The optical properties are based on cholophyll concentrations.
126  call opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
127  g, gv, cs, chl_2d, chl_3d)
128  else ! Use sw e-folding scale set by MOM_input
129  if (optics%nbands <= 1) then ; inv_nbands = 1.0
130  else ; inv_nbands = 1.0 / real(optics%nbands) ; endif
131 
132  ! Make sure there is no division by 0.
133  inv_sw_pen_scale = 1.0 / max(cs%pen_sw_scale, 0.1*gv%Angstrom_m, &
134  gv%H_to_m*gv%H_subroundoff)
135  if ( cs%Opacity_scheme == double_exp ) then
136  !$OMP parallel do default(shared)
137  do k=1,nz ; do j=js,je ; do i=is,ie
138  optics%opacity_band(1,i,j,k) = inv_sw_pen_scale
139  optics%opacity_band(2,i,j,k) = 1.0 / max(cs%pen_sw_scale_2nd, &
140  0.1*gv%Angstrom_m,gv%H_to_m*gv%H_subroundoff)
141  enddo ; enddo ; enddo
142  if (.not.associated(sw_total) .or. (cs%pen_SW_scale <= 0.0)) then
143  !$OMP parallel do default(shared)
144  do j=js,je ; do i=is,ie ; do n=1,optics%nbands
145  optics%sw_pen_band(n,i,j) = 0.0
146  enddo ; enddo ; enddo
147  else
148  !$OMP parallel do default(shared)
149  do j=js,je ; do i=is,ie
150  optics%sw_pen_band(1,i,j) = (cs%SW_1st_EXP_RATIO) * sw_total(i,j)
151  optics%sw_pen_band(2,i,j) = (1.-cs%SW_1st_EXP_RATIO) * sw_total(i,j)
152  enddo ; enddo
153  endif
154  else
155  do k=1,nz ; do j=js,je ; do i=is,ie ; do n=1,optics%nbands
156  optics%opacity_band(n,i,j,k) = inv_sw_pen_scale
157  enddo ; enddo ; enddo ; enddo
158  if (.not.associated(sw_total) .or. (cs%pen_SW_scale <= 0.0)) then
159  !$OMP parallel do default(shared)
160  do j=js,je ; do i=is,ie ; do n=1,optics%nbands
161  optics%sw_pen_band(n,i,j) = 0.0
162  enddo ; enddo ; enddo
163  else
164  !$OMP parallel do default(shared)
165  do j=js,je ; do i=is,ie ; do n=1,optics%nbands
166  optics%sw_pen_band(n,i,j) = cs%pen_SW_frac * inv_nbands * sw_total(i,j)
167  enddo ; enddo ; enddo
168  endif
169  endif
170  endif
171 
172  if (query_averaging_enabled(cs%diag)) then
173  if (cs%id_sw_pen > 0) then
174  !$OMP parallel do default(shared)
175  do j=js,je ; do i=is,ie
176  pen_sw_tot(i,j) = 0.0
177  do n=1,optics%nbands
178  pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
179  enddo
180  enddo ; enddo
181  call post_data(cs%id_sw_pen, pen_sw_tot, cs%diag)
182  endif
183  if (cs%id_sw_vis_pen > 0) then
184  if (cs%opacity_scheme == manizza_05) then
185  !$OMP parallel do default(shared)
186  do j=js,je ; do i=is,ie
187  pen_sw_tot(i,j) = 0.0
188  do n=1,min(optics%nbands,2)
189  pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
190  enddo
191  enddo ; enddo
192  else
193  !$OMP parallel do default(shared)
194  do j=js,je ; do i=is,ie
195  pen_sw_tot(i,j) = 0.0
196  do n=1,optics%nbands
197  pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
198  enddo
199  enddo ; enddo
200  endif
201  call post_data(cs%id_sw_vis_pen, pen_sw_tot, cs%diag)
202  endif
203  do n=1,optics%nbands ; if (cs%id_opacity(n) > 0) then
204  !$OMP parallel do default(shared)
205  do k=1,nz ; do j=js,je ; do i=is,ie
206  ! Remap opacity (op) to 1/L * tanh(op * L) where L is one Angstrom.
207  ! This gives a nearly identical value when op << 1/L but allows one to
208  ! store the values when opacity is divergent (i.e. opaque).
209  tmp(i,j,k) = tanh(op_diag_len * optics%opacity_band(n,i,j,k)) / op_diag_len
210  enddo ; enddo ; enddo
211  call post_data(cs%id_opacity(n), tmp, cs%diag)
212  endif ; enddo
213  endif
214 
215 end subroutine set_opacity
216 
217 
218 !> This sets the "blue" band opacity based on chloophyll A concencentrations
219 !! The red portion is lumped into the net heating at the surface.
220 subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
221  G, GV, CS, chl_2d, chl_3d)
222  type(optics_type), intent(inout) :: optics !< An optics structure that has values
223  !! set based on the opacities.
224  real, dimension(:,:), pointer :: sw_total !< Total shortwave flux into the ocean [W m-2]
225  real, dimension(:,:), pointer :: sw_vis_dir !< Visible, direct shortwave into the ocean [W m-2]
226  real, dimension(:,:), pointer :: sw_vis_dif !< Visible, diffuse shortwave into the ocean [W m-2]
227  real, dimension(:,:), pointer :: sw_nir_dir !< Near-IR, direct shortwave into the ocean [W m-2]
228  real, dimension(:,:), pointer :: sw_nir_dif !< Near-IR, diffuse shortwave into the ocean [W m-2]
229  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
230  type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
231  type(opacity_cs), pointer :: CS !< The control structure.
232  real, dimension(SZI_(G),SZJ_(G)), &
233  optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentractions [mg m-3]
234  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
235  optional, intent(in) :: chl_3d !< A 3-d field of chlorophyll-A concentractions [mg m-3]
236 
237  real :: chl_data(SZI_(G),SZJ_(G)) ! The chlorophyll A concentrations in a layer [mg m-3].
238  real :: Inv_nbands ! The inverse of the number of bands of penetrating
239  ! shortwave radiation.
240  real :: Inv_nbands_nir ! The inverse of the number of bands of penetrating
241  ! near-infrafed radiation.
242  real :: SW_pen_tot ! The sum across the bands of the penetrating
243  ! shortwave radiation [W m-2].
244  real :: SW_vis_tot ! The sum across the visible bands of shortwave
245  ! radiation [W m-2].
246  real :: SW_nir_tot ! The sum across the near infrared bands of shortwave
247  ! radiation [W m-2].
248  type(time_type) :: day
249  character(len=128) :: mesg
250  integer :: i, j, k, n, is, ie, js, je, nz, nbands
251  logical :: multiband_vis_input, multiband_nir_input
252 
253  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
254 
255 ! In this model, the Morel (modified) and Manizza (modified) schemes
256 ! use the "blue" band in the parameterizations to determine the e-folding
257 ! depth of the incoming shortwave attenuation. The red portion is lumped
258 ! into the net heating at the surface.
259 !
260 ! Morel, A., Optical modeling of the upper ocean in relation to its biogenous
261 ! matter content (case-i waters).,J. Geo. Res., {93}, 10,749--10,768, 1988.
262 !
263 ! Manizza, M., C.~L. Quere, A.~Watson, and E.~T. Buitenhuis, Bio-optical
264 ! feedbacks among phytoplankton, upper ocean physics and sea-ice in a
265 ! global model, Geophys. Res. Let., , L05,603, 2005.
266 
267  nbands = optics%nbands
268 
269  if (nbands <= 1) then ; inv_nbands = 1.0
270  else ; inv_nbands = 1.0 / real(nbands) ; endif
271 
272  if (nbands <= 2) then ; inv_nbands_nir = 0.0
273  else ; inv_nbands_nir = 1.0 / real(nbands - 2.0) ; endif
274 
275  multiband_vis_input = (associated(sw_vis_dir) .and. &
276  associated(sw_vis_dif))
277  multiband_nir_input = (associated(sw_nir_dir) .and. &
278  associated(sw_nir_dif))
279 
280  chl_data(:,:) = 0.0
281  if (present(chl_3d)) then
282  do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,1) ; enddo ; enddo
283  do k=1,nz; do j=js,je ; do i=is,ie
284  if ((g%mask2dT(i,j) > 0.5) .and. (chl_3d(i,j,k) < 0.0)) then
285  write(mesg,'(" Negative chl_3d of ",(1pe12.4)," found at i,j,k = ", &
286  & 3(1x,i3), " lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
287  chl_3d(i,j,k), i, j, k, g%geoLonT(i,j), g%geoLatT(i,j)
288  call mom_error(fatal, "MOM_opacity opacity_from_chl: "//trim(mesg))
289  endif
290  enddo ; enddo ; enddo
291  elseif (present(chl_2d)) then
292  do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_2d(i,j) ; enddo ; enddo
293  do j=js,je ; do i=is,ie
294  if ((g%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then
295  write(mesg,'(" Negative chl_2d of ",(1pe12.4)," at i,j = ", &
296  & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
297  chl_data(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
298  call mom_error(fatal, "MOM_opacity opacity_from_chl: "//trim(mesg))
299  endif
300  enddo ; enddo
301  else
302  call mom_error(fatal, "Either chl_2d or chl_3d must be preesnt in a call to opacity_form_chl.")
303  endif
304 
305  select case (cs%opacity_scheme)
306  case (manizza_05)
307  !$OMP parallel do default(shared) private(SW_vis_tot,SW_nir_tot)
308  do j=js,je ; do i=is,ie
309  sw_vis_tot = 0.0 ; sw_nir_tot = 0.0
310  if (g%mask2dT(i,j) > 0.5) then
311  if (multiband_vis_input) then
312  sw_vis_tot = sw_vis_dir(i,j) + sw_vis_dif(i,j)
313  else ! Follow Manizza 05 in assuming that 42% of SW is visible.
314  sw_vis_tot = 0.42 * sw_total(i,j)
315  endif
316  if (multiband_nir_input) then
317  sw_nir_tot = sw_nir_dir(i,j) + sw_nir_dif(i,j)
318  else
319  sw_nir_tot = sw_total(i,j) - sw_vis_tot
320  endif
321  endif
322 
323  ! Band 1 is Manizza blue.
324  optics%sw_pen_band(1,i,j) = cs%blue_frac*sw_vis_tot
325  ! Band 2 (if used) is Manizza red.
326  if (nbands > 1) &
327  optics%sw_pen_band(2,i,j) = (1.0-cs%blue_frac)*sw_vis_tot
328  ! All remaining bands are NIR, for lack of something better to do.
329  do n=3,nbands
330  optics%sw_pen_band(n,i,j) = inv_nbands_nir * sw_nir_tot
331  enddo
332  enddo ; enddo
333  case (morel_88)
334  !$OMP parallel do default(shared) private(SW_pen_tot)
335  do j=js,je ; do i=is,ie
336  sw_pen_tot = 0.0
337  if (g%mask2dT(i,j) > 0.5) then ; if (multiband_vis_input) then
338  sw_pen_tot = sw_pen_frac_morel(chl_data(i,j)) * &
339  (sw_vis_dir(i,j) + sw_vis_dif(i,j))
340  else
341  sw_pen_tot = sw_pen_frac_morel(chl_data(i,j)) * &
342  0.5*sw_total(i,j)
343  endif ; endif
344 
345  do n=1,nbands
346  optics%sw_pen_band(n,i,j) = inv_nbands*sw_pen_tot
347  enddo
348  enddo ; enddo
349  case default
350  call mom_error(fatal, "opacity_from_chl: CS%opacity_scheme is not valid.")
351  end select
352 
353  !$OMP parallel do default(shared) firstprivate(chl_data)
354  do k=1,nz
355  if (present(chl_3d)) then
356  do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,k) ; enddo ; enddo
357  endif
358 
359  select case (cs%opacity_scheme)
360  case (manizza_05)
361  do j=js,je ; do i=is,ie
362  if (g%mask2dT(i,j) <= 0.5) then
363  do n=1,optics%nbands
364  optics%opacity_band(n,i,j,k) = cs%opacity_land_value
365  enddo
366  else
367  ! Band 1 is Manizza blue.
368  optics%opacity_band(1,i,j,k) = 0.0232 + 0.074*chl_data(i,j)**0.674
369  if (nbands >= 2) & ! Band 2 is Manizza red.
370  optics%opacity_band(2,i,j,k) = 0.225 + 0.037*chl_data(i,j)**0.629
371  ! All remaining bands are NIR, for lack of something better to do.
372  do n=3,nbands ; optics%opacity_band(n,i,j,k) = 2.86 ; enddo
373  endif
374  enddo ; enddo
375  case (morel_88)
376  do j=js,je ; do i=is,ie
377  optics%opacity_band(1,i,j,k) = cs%opacity_land_value
378  if (g%mask2dT(i,j) > 0.5) &
379  optics%opacity_band(1,i,j,k) = opacity_morel(chl_data(i,j))
380 
381  do n=2,optics%nbands
382  optics%opacity_band(n,i,j,k) = optics%opacity_band(1,i,j,k)
383  enddo
384  enddo ; enddo
385 
386  case default
387  call mom_error(fatal, "opacity_from_chl: CS%opacity_scheme is not valid.")
388  end select
389  enddo
390 
391 
392 end subroutine opacity_from_chl
393 
394 !> This sets the blue-wavelength opacity according to the scheme proposed by
395 !! Morel and Antoine (1994).
396 function opacity_morel(chl_data)
397  real, intent(in) :: chl_data !< The chlorophyll-A concentration in mg m-3.
398  real :: opacity_morel !< The returned opacity [m-1]
399 
400  ! The following are coefficients for the optical model taken from Morel and
401  ! Antoine (1994). These coeficients represent a non uniform distribution of
402  ! chlorophyll-a through the water column. Other approaches may be more
403  ! appropriate when using an interactive ecosystem model that predicts
404  ! three-dimensional chl-a values.
405  real, dimension(6), parameter :: &
406  z2_coef=(/7.925, -6.644, 3.662, -1.815, -0.218, 0.502/)
407  real :: chl, chl2 ! The log10 of chl_data (in mg m-3), and Chl^2.
408 
409  chl = log10(min(max(chl_data,0.02),60.0)) ; chl2 = chl*chl
410  opacity_morel = 1.0 / ( (z2_coef(1) + z2_coef(2)*chl) + chl2 * &
411  ((z2_coef(3) + chl*z2_coef(4)) + chl2*(z2_coef(5) + chl*z2_coef(6))) )
412 end function
413 
414 !> This sets the penetrating shortwave fraction according to the scheme proposed by
415 !! Morel and Antoine (1994).
416 function sw_pen_frac_morel(chl_data)
417  real, intent(in) :: chl_data !< The chlorophyll-A concentration in mg m-3.
418  real :: sw_pen_frac_morel !< The returned penetrating shortwave fraction [nondim]
419 
420  ! The following are coefficients for the optical model taken from Morel and
421  ! Antoine (1994). These coeficients represent a non uniform distribution of
422  ! chlorophyll-a through the water column. Other approaches may be more
423  ! appropriate when using an interactive ecosystem model that predicts
424  ! three-dimensional chl-a values.
425  real :: chl, chl2 ! The log10 of chl_data in mg m-3, and Chl^2.
426  real, dimension(6), parameter :: &
427  v1_coef=(/0.321, 0.008, 0.132, 0.038, -0.017, -0.007/)
428 
429  chl = log10(min(max(chl_data,0.02),60.0)) ; chl2 = chl*chl
430  sw_pen_frac_morel = 1.0 - ( (v1_coef(1) + v1_coef(2)*chl) + chl2 * &
431  ((v1_coef(3) + chl*v1_coef(4)) + chl2*(v1_coef(5) + chl*v1_coef(6))) )
432 end function sw_pen_frac_morel
433 
434 !> This sets the blue-wavelength opacity according to the scheme proposed by
435 !! Manizza, M. et al, 2005.
436 function opacity_manizza(chl_data)
437  real, intent(in) :: chl_data !< The chlorophyll-A concentration in mg m-3.
438  real :: opacity_manizza !< The returned opacity [m-1]
439 ! This sets the blue-wavelength opacity according to the scheme proposed by Manizza, M. et al, 2005.
440 
441  opacity_manizza = 0.0232 + 0.074*chl_data**0.674
442 end function
443 
444 !> This subroutine returns a 2-d slice at constant j of fields from an optics_type, with the potential
445 !! for rescaling these fields.
446 subroutine extract_optics_slice(optics, j, G, GV, opacity, opacity_scale, penSW_top, penSW_scale)
447  type(optics_type), intent(in) :: optics !< An optics structure that has values of opacities
448  !! and shortwave fluxes.
449  integer, intent(in) :: j !< j-index to extract
450  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
451  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
452  real, dimension(max(optics%nbands,1),SZI_(G),SZK_(GV)), &
453  optional, intent(out) :: opacity !< The opacity in each band, i-point, and layer
454  real, optional, intent(in) :: opacity_scale !< A factor by which to rescale the opacity.
455  real, dimension(max(optics%nbands,1),SZI_(G)), &
456  optional, intent(out) :: pensw_top !< The shortwave radiation [W m-2] at the surface
457  !! in each of the nbands bands that penetrates
458  !! beyond the surface skin layer.
459  real, optional, intent(in) :: pensw_scale !< A factor by which to rescale the shortwave flux.
460 
461  ! Local variables
462  real :: scale_opacity, scale_pensw ! Rescaling factors
463  integer :: i, is, ie, k, nz, n
464  is = g%isc ; ie = g%iec ; nz = g%ke
465 
466  scale_opacity = 1.0 ; if (present(opacity_scale)) scale_opacity = opacity_scale
467  scale_pensw = 1.0 ; if (present(pensw_scale)) scale_pensw = pensw_scale
468 
469  if (present(opacity)) then ; do k=1,nz ; do i=is,ie
470  do n=1,optics%nbands
471  opacity(n,i,k) = scale_opacity * optics%opacity_band(n,i,j,k)
472  enddo
473  enddo ; enddo ; endif
474 
475  if (present(pensw_top)) then ; do k=1,nz ; do i=is,ie
476  do n=1,optics%nbands
477  pensw_top(n,i) = scale_pensw * optics%SW_pen_band(n,i,j)
478  enddo
479  enddo ; enddo ; endif
480 
481 end subroutine extract_optics_slice
482 
483 !> Set arguments to fields from the optics type.
484 subroutine extract_optics_fields(optics, nbands)
485  type(optics_type), intent(in) :: optics !< An optics structure that has values of opacities
486  !! and shortwave fluxes.
487  integer, optional, intent(out) :: nbands !< The number of penetrating bands of SW radiation
488 
489  if (present(nbands)) nbands = optics%nbands
490 
491 end subroutine extract_optics_fields
492 
493 !> Return the number of bands of penetrating shortwave radiation.
494 function optics_nbands(optics)
495  type(optics_type), intent(in) :: optics !< An optics structure that has values of opacities
496  !! and shortwave fluxes.
497  integer :: optics_nbands !< The number of penetrating bands of SW radiation
498 
499  optics_nbands = optics%nbands
500 end function optics_nbands
501 
502 !> Apply shortwave heating below the boundary layer (when running with the bulk mixed layer inhereted
503 !! from GOLD) or throughout the water column.
504 !!
505 !! In addition, it causes all of the remaining SW radiation to be absorbed, provided that the total
506 !! water column thickness is greater than H_limit_fluxes.
507 !! For thinner water columns, the heating is scaled down proportionately, the assumption being that the
508 !! remaining heating (which is left in Pen_SW) should go into an (absent for now) ocean bottom sediment layer.
509 subroutine absorbremainingsw(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_limit_fluxes, &
510  adjustAbsorptionProfile, absorbAllSW, T, Pen_SW_bnd, &
511  eps, ksort, htot, Ttot, TKE, dSV_dT)
512 
513  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
514  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
515  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
516  integer, intent(in) :: nsw !< Number of bands of penetrating
517  !! shortwave radiation.
518  real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
519  real, dimension(max(1,nsw),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< Opacity in each band of penetrating
520  !! shortwave radiation [H-1 ~> m-1 or m2 kg-1].
521  !! The indicies are band, i, k.
522  type(optics_type), intent(in) :: optics !< An optics structure that has values of
523  !! opacities and shortwave fluxes.
524  integer, intent(in) :: j !< j-index to work on.
525  real, intent(in) :: dt !< Time step [T ~> s].
526  real, intent(in) :: h_limit_fluxes !< If the total ocean depth is
527  !! less than this, they are scaled away
528  !! to avoid numerical instabilities
529  !! [H ~> m or kg m-2]. This would
530  !! not be necessary if a finite heat
531  !! capacity mud-layer were added.
532  logical, intent(in) :: adjustabsorptionprofile !< If true, apply
533  !! heating above the layers in which it
534  !! should have occurred to get the
535  !! correct mean depth (and potential
536  !! energy change) of the shortwave that
537  !! should be absorbed by each layer.
538  logical, intent(in) :: absorballsw !< If true, apply heating above the
539  !! layers in which it should have occurred
540  !! to get the correct mean depth (and
541  !! potential energy change) of the
542  !! shortwave that should be absorbed by
543  !! each layer.
544  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: t !< Layer potential/conservative
545  !! temperatures [degC]
546  real, dimension(max(1,nsw),SZI_(G)), intent(inout) :: pen_sw_bnd !< Penetrating shortwave heating in
547  !! each band that hits the bottom and will
548  !! will be redistributed through the water
549  !! column [degC H ~> degC m or degC kg m-2],
550  !! size nsw x SZI_(G).
551  real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: eps !< Small thickness that must remain in
552  !! each layer, and which will not be
553  !! subject to heating [H ~> m or kg m-2]
554  integer, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: ksort !< Density-sorted k-indicies.
555  real, dimension(SZI_(G)), optional, intent(in) :: htot !< Total mixed layer thickness [H ~> m or kg m-2].
556  real, dimension(SZI_(G)), optional, intent(inout) :: ttot !< Depth integrated mixed layer
557  !! temperature [degC H ~> degC m or degC kg m-2]
558  real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: dsv_dt !< The partial derivative of specific
559  !! volume with temperature [m3 kg-1 degC-1].
560  real, dimension(SZI_(G),SZK_(GV)), optional, intent(inout) :: tke !< The TKE sink from mixing the heating
561  !! throughout a layer [kg m-3 Z3 T-2 ~> J m-2].
562 
563  ! Local variables
564  real, dimension(SZI_(G),SZK_(GV)) :: &
565  t_chg_above ! A temperature change that will be applied to all the thick
566  ! layers above a given layer [degC]. This is only nonzero if
567  ! adjustAbsorptionProfile is true, in which case the net
568  ! change in the temperature of a layer is the sum of the
569  ! direct heating of that layer plus T_chg_above from all of
570  ! the layers below, plus any contribution from absorbing
571  ! radiation that hits the bottom.
572  real, dimension(SZI_(G)) :: &
573  h_heat, & ! The thickness of the water column that will be heated by
574  ! any remaining shortwave radiation [H ~> m or kg m-2].
575  t_chg, & ! The temperature change of thick layers due to the remaining
576  ! shortwave radiation and contributions from T_chg_above [degC].
577  pen_sw_rem ! The sum across all wavelength bands of the penetrating shortwave
578  ! heating that hits the bottom and will be redistributed through
579  ! the water column [degC H ~> degC m or degC kg m-2]
580  real :: sw_trans ! fraction of shortwave radiation that is not
581  ! absorbed in a layer [nondim]
582  real :: unabsorbed ! fraction of the shortwave radiation that
583  ! is not absorbed because the layers are too thin
584  real :: ih_limit ! inverse of the total depth at which the
585  ! surface fluxes start to be limited [H-1 ~> m-1 or m2 kg-1]
586  real :: h_min_heat ! minimum thickness layer that should get heated [H ~> m or kg m-2]
587  real :: opt_depth ! optical depth of a layer [nondim]
588  real :: exp_od ! exp(-opt_depth) [nondim]
589  real :: heat_bnd ! heating due to absorption in the current
590  ! layer by the current band, including any piece that
591  ! is moved upward [degC H ~> degC m or degC kg m-2]
592  real :: swa ! fraction of the absorbed shortwave that is
593  ! moved to layers above with adjustAbsorptionProfile [nondim]
594  real :: coswa_frac ! The fraction of SWa that is actually moved upward.
595  real :: min_sw_heat ! A minimum remaining shortwave heating within a timestep that will be simply
596  ! absorbed in the next layer for computational efficiency, instead of
597  ! continuing to penetrate [degC H ~> degC m or degC kg m-2].
598  real :: i_habs ! The inverse of the absorption length for a minimal flux [H-1 ~> m-1 or m2 kg-1]
599  real :: epsilon ! A small thickness that must remain in each
600  ! layer, and which will not be subject to heating [H ~> m or kg m-2]
601  real :: g_hconv2 ! A conversion factor for use in the TKE calculation
602  ! in units of [Z3 kg2 m-6 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
603  logical :: sw_remains ! If true, some column has shortwave radiation that
604  ! was not entirely absorbed.
605  logical :: tke_calc ! If true, calculate the implications to the
606  ! TKE budget of the shortwave heating.
607  real :: c1_6, c1_60
608  integer :: is, ie, nz, i, k, ks, n
609  sw_remains = .false.
610 
611  min_sw_heat = optics%PenSW_flux_absorb * dt
612  i_habs = optics%PenSW_absorb_Invlen
613 
614  h_min_heat = 2.0*gv%Angstrom_H + gv%H_subroundoff
615  is = g%isc ; ie = g%iec ; nz = g%ke
616  c1_6 = 1.0 / 6.0 ; c1_60 = 1.0 / 60.0
617 
618  tke_calc = (present(tke) .and. present(dsv_dt))
619 
620  if (optics%answers_2018) then
621  g_hconv2 = (us%m_to_Z**2 * us%L_to_Z**2*gv%g_Earth * gv%H_to_kg_m2) * gv%H_to_kg_m2
622  else
623  g_hconv2 = us%m_to_Z**2 * us%L_to_Z**2*gv%g_Earth * gv%H_to_kg_m2**2
624  endif
625 
626  h_heat(:) = 0.0
627  if (present(htot)) then ; do i=is,ie ; h_heat(i) = htot(i) ; enddo ; endif
628 
629  ! Apply penetrating SW radiation to remaining parts of layers.
630  ! Excessively thin layers are not heated to avoid runaway temps.
631  do ks=1,nz ; do i=is,ie
632  k = ks
633  if (present(ksort)) then
634  if (ksort(i,ks) <= 0) cycle
635  k = ksort(i,ks)
636  endif
637  epsilon = 0.0 ; if (present(eps)) epsilon = eps(i,k)
638 
639  t_chg_above(i,k) = 0.0
640 
641  if (h(i,k) > 1.5*epsilon) then
642  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
643  ! SW_trans is the SW that is transmitted THROUGH the layer
644  opt_depth = h(i,k) * opacity_band(n,i,k)
645  exp_od = exp(-opt_depth)
646  sw_trans = exp_od
647 
648  ! Heating at a very small rate can be absorbed by a sufficiently thick layer or several
649  ! thin layers without further penetration.
650  if (optics%answers_2018) then
651  if (nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat*min(1.0, i_habs*h(i,k)) ) sw_trans = 0.0
652  elseif ((nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat) .and. (h(i,k) > h_min_heat)) then
653  if (nsw*pen_sw_bnd(n,i) <= min_sw_heat * (i_habs*(h(i,k) - h_min_heat))) then
654  sw_trans = 0.0
655  else
656  sw_trans = min(sw_trans, &
657  1.0 - (min_sw_heat*(i_habs*(h(i,k) - h_min_heat))) / (nsw*pen_sw_bnd(n,i)))
658  endif
659  endif
660 
661  heat_bnd = pen_sw_bnd(n,i) * (1.0 - sw_trans)
662  if (adjustabsorptionprofile .and. (h_heat(i) > 0.0)) then
663  ! In this case, a fraction of the heating is applied to the
664  ! overlying water so that the mean pressure at which the shortwave
665  ! heating occurs is exactly what it would have been with a careful
666  ! pressure-weighted averaging of the exponential heating profile,
667  ! hence there should be no TKE budget requirements due to this
668  ! layer. Very clever, but this is also limited so that the
669  ! water above is not heated at a faster rate than the layer
670  ! actually being heated, i.e., SWA <= h_heat / (h_heat + h(i,k))
671  ! and takes the energetics of the rest of the heating into account.
672  ! (-RWH, ~7 years later.)
673  if (opt_depth > 1e-5) then
674  swa = ((opt_depth + (opt_depth + 2.0)*exp_od) - 2.0) / &
675  ((opt_depth + opacity_band(n,i,k) * h_heat(i)) * &
676  (1.0 - exp_od))
677  else
678  ! Use Taylor series expansion of the expression above for a
679  ! more accurate form with very small layer optical depths.
680  swa = h(i,k) * (opt_depth * (1.0 - opt_depth)) / &
681  ((h_heat(i) + h(i,k)) * (6.0 - 3.0*opt_depth))
682  endif
683  coswa_frac = 0.0
684  if (swa*(h_heat(i) + h(i,k)) > h_heat(i)) then
685  coswa_frac = (swa*(h_heat(i) + h(i,k)) - h_heat(i) ) / &
686  (swa*(h_heat(i) + h(i,k)))
687  swa = h_heat(i) / (h_heat(i) + h(i,k))
688  endif
689 
690  t_chg_above(i,k) = t_chg_above(i,k) + (swa * heat_bnd) / h_heat(i)
691  t(i,k) = t(i,k) + ((1.0 - swa) * heat_bnd) / h(i,k)
692  else
693  coswa_frac = 1.0
694  t(i,k) = t(i,k) + pen_sw_bnd(n,i) * (1.0 - sw_trans) / h(i,k)
695  endif
696 
697  if (tke_calc) then
698  if (opt_depth > 1e-2) then
699  tke(i,k) = tke(i,k) - coswa_frac*heat_bnd*dsv_dt(i,k)* &
700  (0.5*h(i,k)*g_hconv2) * &
701  (opt_depth*(1.0+exp_od) - 2.0*(1.0-exp_od)) / (opt_depth*(1.0-exp_od))
702  else
703  ! Use Taylor series-derived approximation to the above expression
704  ! that is well behaved and more accurate when opt_depth is small.
705  tke(i,k) = tke(i,k) - coswa_frac*heat_bnd*dsv_dt(i,k)* &
706  (0.5*h(i,k)*g_hconv2) * &
707  (c1_6*opt_depth * (1.0 - c1_60*opt_depth**2))
708  endif
709  endif
710 
711  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
712  endif ; enddo
713  endif
714 
715  ! Add to the accumulated thickness above that could be heated.
716  ! Only layers greater than h_min_heat thick should get heated.
717  if (h(i,k) >= 2.0*h_min_heat) then
718  h_heat(i) = h_heat(i) + h(i,k)
719  elseif (h(i,k) > h_min_heat) then
720  h_heat(i) = h_heat(i) + (2.0*h(i,k) - 2.0*h_min_heat)
721  endif
722  enddo ; enddo ! i & k loops
723 
724 
725 ! if (.not.absorbAllSW .and. .not.adjustAbsorptionProfile) return
726 
727  ! Unless modified, there is no temperature change due to fluxes from the bottom.
728  do i=is,ie ; t_chg(i) = 0.0 ; enddo
729 
730  if (absorballsw) then
731  ! If there is still shortwave radiation at this point, it could go into
732  ! the bottom (with a bottom mud model), or it could be redistributed back
733  ! through the water column.
734  do i=is,ie
735  pen_sw_rem(i) = pen_sw_bnd(1,i)
736  do n=2,nsw ; pen_sw_rem(i) = pen_sw_rem(i) + pen_sw_bnd(n,i) ; enddo
737  enddo
738  do i=is,ie ; if (pen_sw_rem(i) > 0.0) sw_remains = .true. ; enddo
739 
740  ih_limit = 1.0 / h_limit_fluxes
741  do i=is,ie ; if ((pen_sw_rem(i) > 0.0) .and. (h_heat(i) > 0.0)) then
742  if (h_heat(i)*ih_limit >= 1.0) then
743  t_chg(i) = pen_sw_rem(i) / h_heat(i) ; unabsorbed = 0.0
744  else
745  t_chg(i) = pen_sw_rem(i) * ih_limit
746  unabsorbed = 1.0 - h_heat(i)*ih_limit
747  endif
748  do n=1,nsw ; pen_sw_bnd(n,i) = unabsorbed * pen_sw_bnd(n,i) ; enddo
749  endif ; enddo
750  endif ! absorbAllSW
751 
752  if (absorballsw .or. adjustabsorptionprofile) then
753  do ks=nz,1,-1 ; do i=is,ie
754  k = ks
755  if (present(ksort)) then
756  if (ksort(i,ks) <= 0) cycle
757  k = ksort(i,ks)
758  endif
759 
760  if (t_chg(i) > 0.0) then
761  ! Only layers greater than h_min_heat thick should get heated.
762  if (h(i,k) >= 2.0*h_min_heat) then ; t(i,k) = t(i,k) + t_chg(i)
763  elseif (h(i,k) > h_min_heat) then
764  t(i,k) = t(i,k) + t_chg(i) * (2.0 - 2.0*h_min_heat/h(i,k))
765  endif
766  endif
767  ! Increase the heating for layers above.
768  t_chg(i) = t_chg(i) + t_chg_above(i,k)
769  enddo ; enddo
770  if (present(htot) .and. present(ttot)) then
771  do i=is,ie ; ttot(i) = ttot(i) + t_chg(i) * htot(i) ; enddo
772  endif
773  endif ! absorbAllSW .or. adjustAbsorptionProfile
774 
775 end subroutine absorbremainingsw
776 
777 
778 !> This subroutine calculates the total shortwave heat flux integrated over
779 !! bands as a function of depth. This routine is only called for computing
780 !! buoyancy fluxes for use in KPP. This routine does not updat e the state.
781 subroutine sumswoverbands(G, GV, US, h, nsw, optics, j, dt, &
782  H_limit_fluxes, absorbAllSW, iPen_SW_bnd, netPen)
783  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
784  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
785  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
786  real, dimension(SZI_(G),SZK_(GV)), &
787  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
788  integer, intent(in) :: nsw !< The number of bands of penetrating shortwave
789  !! radiation, perhaps from optics_nbands(optics),
790  type(optics_type), intent(in) :: optics !< An optics structure that has values
791  !! set based on the opacities.
792  integer, intent(in) :: j !< j-index to work on.
793  real, intent(in) :: dt !< Time step [T ~> s].
794  real, intent(in) :: h_limit_fluxes !< the total depth at which the
795  !! surface fluxes start to be limited to avoid
796  !! excessive heating of a thin ocean [H ~> m or kg m-2]
797  logical, intent(in) :: absorballsw !< If true, ensure that all shortwave
798  !! radiation is absorbed in the ocean water column.
799  real, dimension(max(nsw,1),SZI_(G)), intent(in) :: ipen_sw_bnd !< The incident penetrating shortwave
800  !! heating in each band that hits the bottom and
801  !! will be redistributed through the water column
802  !! [degC H ~> degC m or degC kg m-2]; size nsw x SZI_(G).
803  real, dimension(SZI_(G),SZK_(GV)+1), &
804  intent(inout) :: netpen !< Net penetrating shortwave heat flux at each
805  !! interface, summed across all bands
806  !! [degC H ~> degC m or degC kg m-2].
807  ! Local variables
808  real :: h_heat(szi_(g)) ! thickness of the water column that receives
809  ! remaining shortwave radiation [H ~> m or kg m-2].
810  real :: pen_sw_rem(szi_(g)) ! sum across all wavelength bands of the
811  ! penetrating shortwave heating that hits the bottom
812  ! and will be redistributed through the water column
813  ! [degC H ~> degC m or degC kg m-2]
814 
815  real, dimension(max(nsw,1),SZI_(G)) :: pen_sw_bnd ! The remaining penetrating shortwave radiation
816  ! in each band, initially iPen_SW_bnd [degC H ~> degC m or degC kg m-2]
817  real :: sw_trans ! fraction of shortwave radiation not
818  ! absorbed in a layer [nondim]
819  real :: unabsorbed ! fraction of the shortwave radiation
820  ! not absorbed because the layers are too thin.
821  real :: ih_limit ! inverse of the total depth at which the
822  ! surface fluxes start to be limited [H-1 ~> m-1 or m2 kg-1]
823  real :: min_sw_heat ! A minimum remaining shortwave heating within a timestep that will be simply
824  ! absorbed in the next layer for computational efficiency, instead of
825  ! continuing to penetrate [degC H ~> degC m or degC kg m-2].
826  real :: i_habs ! The inverse of the absorption length for a minimal flux [H-1 ~> m-1 or m2 kg-1]
827  real :: h_min_heat ! minimum thickness layer that should get heated [H ~> m or kg m-2]
828  real :: opt_depth ! optical depth of a layer [nondim]
829  real :: exp_od ! exp(-opt_depth) [nondim]
830  logical :: sw_remains ! If true, some column has shortwave radiation that
831  ! was not entirely absorbed.
832 
833  integer :: is, ie, nz, i, k, ks, n
834  sw_remains = .false.
835 
836  min_sw_heat = optics%PenSW_flux_absorb*dt ! Default of 2.5e-11*US%T_to_s*GV%m_to_H
837  i_habs = 1e3*gv%H_to_m ! optics%PenSW_absorb_Invlen
838 
839  h_min_heat = 2.0*gv%Angstrom_H + gv%H_subroundoff
840  is = g%isc ; ie = g%iec ; nz = g%ke
841 
842  pen_sw_bnd(:,:) = ipen_sw_bnd(:,:)
843  do i=is,ie ; h_heat(i) = 0.0 ; enddo
844  netpen(:,1) = sum( pen_sw_bnd(:,:), dim=1 ) ! Surface interface
845 
846  ! Apply penetrating SW radiation to remaining parts of layers.
847  ! Excessively thin layers are not heated to avoid runaway temps.
848  do k=1,nz
849 
850  do i=is,ie
851  netpen(i,k+1) = 0.
852 
853  if (h(i,k) > 0.0) then
854  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
855  ! SW_trans is the SW that is transmitted THROUGH the layer
856  opt_depth = h(i,k)*gv%H_to_m * optics%opacity_band(n,i,j,k)
857  exp_od = exp(-opt_depth)
858  sw_trans = exp_od
859 
860  ! Heating at a very small rate can be absorbed by a sufficiently thick layer or several
861  ! thin layers without further penetration.
862  if (optics%answers_2018) then
863  if (nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat*min(1.0, i_habs*h(i,k)) ) sw_trans = 0.0
864  elseif ((nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat) .and. (h(i,k) > h_min_heat)) then
865  if (nsw*pen_sw_bnd(n,i) <= min_sw_heat * (i_habs*(h(i,k) - h_min_heat))) then
866  sw_trans = 0.0
867  else
868  sw_trans = min(sw_trans, &
869  1.0 - (min_sw_heat*(i_habs*(h(i,k) - h_min_heat))) / (nsw*pen_sw_bnd(n,i)))
870  endif
871  endif
872 
873  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
874  netpen(i,k+1) = netpen(i,k+1) + pen_sw_bnd(n,i)
875  endif ; enddo
876  endif ! h(i,k) > 0.0
877 
878  ! Add to the accumulated thickness above that could be heated.
879  ! Only layers greater than h_min_heat thick should get heated.
880  if (h(i,k) >= 2.0*h_min_heat) then
881  h_heat(i) = h_heat(i) + h(i,k)
882  elseif (h(i,k) > h_min_heat) then
883  h_heat(i) = h_heat(i) + (2.0*h(i,k) - 2.0*h_min_heat)
884  endif
885  enddo ! i loop
886  enddo ! k loop
887 
888  if (absorballsw) then
889 
890  ! If there is still shortwave radiation at this point, it could go into
891  ! the bottom (with a bottom mud model), or it could be redistributed back
892  ! through the water column.
893  do i=is,ie
894  pen_sw_rem(i) = pen_sw_bnd(1,i)
895  do n=2,nsw ; pen_sw_rem(i) = pen_sw_rem(i) + pen_sw_bnd(n,i) ; enddo
896  enddo
897  do i=is,ie ; if (pen_sw_rem(i) > 0.0) sw_remains = .true. ; enddo
898 
899  ih_limit = 1.0 / h_limit_fluxes
900  do i=is,ie ; if ((pen_sw_rem(i) > 0.0) .and. (h_heat(i) > 0.0)) then
901  if (h_heat(i)*ih_limit < 1.0) then
902  unabsorbed = 1.0 - h_heat(i)*ih_limit
903  else
904  unabsorbed = 0.0
905  endif
906  do n=1,nsw ; pen_sw_bnd(n,i) = unabsorbed * pen_sw_bnd(n,i) ; enddo
907  endif ; enddo
908 
909  endif ! absorbAllSW
910 
911 end subroutine sumswoverbands
912 
913 
914 
915 !> This routine initalizes the opacity module, including an optics_type.
916 subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
917  type(time_type), target, intent(in) :: time !< The current model time.
918  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
919  type(verticalgrid_type), intent(in) :: gv !< model vertical grid structure
920  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
921  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
922  !! parameters.
923  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
924  !! output.
925  type(opacity_cs), pointer :: cs !< A pointer that is set to point to the control
926  !! structure for this module.
927  type(optics_type), pointer :: optics !< An optics structure that has parameters
928  !! set and arrays allocated here.
929  ! Local variables
930  character(len=200) :: tmpstr
931  character(len=40) :: mdl = "MOM_opacity"
932  character(len=40) :: bandnum, shortname
933  character(len=200) :: longname
934  character(len=40) :: scheme_string
935  ! This include declares and sets the variable "version".
936 # include "version_variable.h"
937  real :: pensw_absorb_minthick ! A thickness that is used to absorb the remaining shortwave heat
938  ! flux when that flux drops below PEN_SW_FLUX_ABSORB [m].
939  real :: pensw_minthick_dflt ! The default for PenSW_absorb_minthick [m]
940  logical :: default_2018_answers
941  logical :: use_scheme
942  integer :: isd, ied, jsd, jed, nz, n
943  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
944 
945  if (associated(cs)) then
946  call mom_error(warning, "opacity_init called with an associated"// &
947  "associated control structure.")
948  return
949  else ; allocate(cs) ; endif
950 
951  cs%diag => diag
952 
953  ! Read all relevant parameters and write them to the model log.
954  call log_version(param_file, mdl, version, '')
955 
956 ! parameters for CHL_A routines
957  call get_param(param_file, mdl, "VAR_PEN_SW", cs%var_pen_sw, &
958  "If true, use one of the CHL_A schemes specified by "//&
959  "OPACITY_SCHEME to determine the e-folding depth of "//&
960  "incoming short wave radiation.", default=.false.)
961 
962  cs%opacity_scheme = no_scheme ; scheme_string = ''
963  if (cs%var_pen_sw) then
964  call get_param(param_file, mdl, "OPACITY_SCHEME", tmpstr, &
965  "This character string specifies how chlorophyll "//&
966  "concentrations are translated into opacities. Currently "//&
967  "valid options include:\n"//&
968  " \t\t MANIZZA_05 - Use Manizza et al., GRL, 2005. \n"//&
969  " \t\t MOREL_88 - Use Morel, JGR, 1988.", &
970  default=manizza_05_string)
971  if (len_trim(tmpstr) > 0) then
972  tmpstr = uppercase(tmpstr)
973  select case (tmpstr)
974  case (manizza_05_string)
975  cs%opacity_scheme = manizza_05 ; scheme_string = manizza_05_string
976  case (morel_88_string)
977  cs%opacity_scheme = morel_88 ; scheme_string = morel_88_string
978  case default
979  call mom_error(fatal, "opacity_init: #DEFINE OPACITY_SCHEME "//&
980  trim(tmpstr) // "in input file is invalid.")
981  end select
982  call mom_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5)
983  endif
984  if (cs%opacity_scheme == no_scheme) then
985  call mom_error(warning, "opacity_init: No scheme has successfully "//&
986  "been specified for the opacity. Using the default MANIZZA_05.")
987  cs%opacity_scheme = manizza_05 ; scheme_string = manizza_05_string
988  endif
989 
990  call get_param(param_file, mdl, "BLUE_FRAC_SW", cs%blue_frac, &
991  "The fraction of the penetrating shortwave radiation "//&
992  "that is in the blue band.", default=0.5, units="nondim")
993  else
994  call get_param(param_file, mdl, "EXP_OPACITY_SCHEME", tmpstr, &
995  "This character string specifies which exponential "//&
996  "opacity scheme to utilize. Currently "//&
997  "valid options include:\n"//&
998  " \t\t SINGLE_EXP - Single Exponent decay. \n"//&
999  " \t\t DOUBLE_EXP - Double Exponent decay.", &
1000  default=single_exp_string)!New default for "else" above (non-Chl scheme)
1001  if (len_trim(tmpstr) > 0) then
1002  tmpstr = uppercase(tmpstr)
1003  select case (tmpstr)
1004  case (single_exp_string)
1005  cs%opacity_scheme = single_exp ; scheme_string = single_exp_string
1006  case (double_exp_string)
1007  cs%opacity_scheme = double_exp ; scheme_string = double_exp_string
1008  end select
1009  call mom_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5)
1010  endif
1011  call get_param(param_file, mdl, "PEN_SW_SCALE", cs%pen_sw_scale, &
1012  "The vertical absorption e-folding depth of the "//&
1013  "penetrating shortwave radiation.", units="m", default=0.0)
1014  !BGR/ Added for opacity_scheme==double_exp read in 2nd exp-decay and fraction
1015  if (cs%Opacity_scheme == double_exp ) then
1016  call get_param(param_file, mdl, "PEN_SW_SCALE_2ND", cs%pen_sw_scale_2nd, &
1017  "The (2nd) vertical absorption e-folding depth of the "//&
1018  "penetrating shortwave radiation "//&
1019  "(use if SW_EXP_MODE==double.)",&
1020  units="m", default=0.0)
1021  call get_param(param_file, mdl, "SW_1ST_EXP_RATIO", cs%sw_1st_exp_ratio, &
1022  "The fraction of 1st vertical absorption e-folding depth "//&
1023  "penetrating shortwave radiation if SW_EXP_MODE==double.",&
1024  units="m", default=0.0)
1025  elseif (cs%OPACITY_SCHEME == single_exp) then
1026  !/Else disable 2nd_exp scheme
1027  cs%pen_sw_scale_2nd = 0.0
1028  cs%sw_1st_exp_ratio = 1.0
1029  endif
1030  call get_param(param_file, mdl, "PEN_SW_FRAC", cs%pen_sw_frac, &
1031  "The fraction of the shortwave radiation that penetrates "//&
1032  "below the surface.", units="nondim", default=0.0)
1033 
1034  endif
1035  call get_param(param_file, mdl, "PEN_SW_NBANDS", optics%nbands, &
1036  "The number of bands of penetrating shortwave radiation.", &
1037  default=1)
1038  if (cs%Opacity_scheme == double_exp ) then
1039  if (optics%nbands /= 2) call mom_error(fatal, &
1040  "set_opacity: \Cannot use a double_exp opacity scheme with nbands!=2.")
1041  elseif (cs%Opacity_scheme == single_exp ) then
1042  if (optics%nbands /= 1) call mom_error(fatal, &
1043  "set_opacity: \Cannot use a single_exp opacity scheme with nbands!=1.")
1044  endif
1045 
1046  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1047  "This sets the default value for the various _2018_ANSWERS parameters.", &
1048  default=.true.)
1049  call get_param(param_file, mdl, "OPTICS_2018_ANSWERS", optics%answers_2018, &
1050  "If true, use the order of arithmetic and expressions that recover the "//&
1051  "answers from the end of 2018. Otherwise, use updated expressions for "//&
1052  "handling the absorption of small remaining shortwave fluxes.", &
1053  default=default_2018_answers)
1054 
1055  call get_param(param_file, mdl, "PEN_SW_FLUX_ABSORB", optics%PenSW_flux_absorb, &
1056  "A minimum remaining shortwave heating rate that will be simply absorbed in "//&
1057  "the next sufficiently thick layers for computational efficiency, instead of "//&
1058  "continuing to penetrate. The default, 2.5e-11 degC m s-1, is about 1e-4 W m-2 "//&
1059  "or 0.08 degC m century-1, but 0 is also a valid value.", &
1060  default=2.5e-11, units="degC m s-1", scale=gv%m_to_H*us%T_to_s)
1061 
1062  if (optics%answers_2018) then ; pensw_minthick_dflt = 0.001 ; else ; pensw_minthick_dflt = 1.0 ; endif
1063  call get_param(param_file, mdl, "PEN_SW_ABSORB_MINTHICK", pensw_absorb_minthick, &
1064  "A thickness that is used to absorb the remaining penetrating shortwave heat "//&
1065  "flux when it drops below PEN_SW_FLUX_ABSORB.", &
1066  default=pensw_minthick_dflt, units="m", scale=gv%m_to_H)
1067  optics%PenSW_absorb_Invlen = 1.0 / (pensw_absorb_minthick + gv%H_subroundoff)
1068 
1069  if (.not.associated(optics%min_wavelength_band)) &
1070  allocate(optics%min_wavelength_band(optics%nbands))
1071  if (.not.associated(optics%max_wavelength_band)) &
1072  allocate(optics%max_wavelength_band(optics%nbands))
1073 
1074  if (cs%opacity_scheme == manizza_05) then
1075  optics%min_wavelength_band(1) =0
1076  optics%max_wavelength_band(1) =550
1077  if (optics%nbands >= 2) then
1078  optics%min_wavelength_band(2)=550
1079  optics%max_wavelength_band(2)=700
1080  endif
1081  if (optics%nbands > 2) then
1082  do n=3,optics%nbands
1083  optics%min_wavelength_band(n) =700
1084  optics%max_wavelength_band(n) =2800
1085  enddo
1086  endif
1087  endif
1088 
1089  call get_param(param_file, mdl, "OPACITY_LAND_VALUE", cs%opacity_land_value, &
1090  "The value to use for opacity over land. The default is "//&
1091  "10 m-1 - a value for muddy water.", units="m-1", default=10.0)
1092 
1093  if (.not.associated(optics%opacity_band)) &
1094  allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz))
1095  if (.not.associated(optics%sw_pen_band)) &
1096  allocate(optics%sw_pen_band(optics%nbands,isd:ied,jsd:jed))
1097  allocate(cs%id_opacity(optics%nbands)) ; cs%id_opacity(:) = -1
1098 
1099  cs%id_sw_pen = register_diag_field('ocean_model', 'SW_pen', diag%axesT1, time, &
1100  'Penetrating shortwave radiation flux into ocean', 'W m-2')
1101  cs%id_sw_vis_pen = register_diag_field('ocean_model', 'SW_vis_pen', diag%axesT1, time, &
1102  'Visible penetrating shortwave radiation flux into ocean', 'W m-2')
1103  do n=1,optics%nbands
1104  write(bandnum,'(i3)') n
1105  shortname = 'opac_'//trim(adjustl(bandnum))
1106  longname = 'Opacity for shortwave radiation in band '//trim(adjustl(bandnum)) &
1107  // ', saved as L^-1 tanh(Opacity * L) for L = 10^-10 m'
1108  cs%id_opacity(n) = register_diag_field('ocean_model', shortname, diag%axesTL, time, &
1109  longname, 'm-1')
1110  enddo
1111 
1112 end subroutine opacity_init
1113 
1114 
1115 subroutine opacity_end(CS, optics)
1116  type(opacity_cs), pointer :: cs !< An opacity control structure that should be deallocated.
1117  type(optics_type), optional, pointer :: optics !< An optics type structure that should be deallocated.
1118 
1119  if (associated(cs%id_opacity)) deallocate(cs%id_opacity)
1120  if (associated(cs)) deallocate(cs)
1121 
1122  if (present(optics)) then ; if (associated(optics)) then
1123  if (associated(optics%opacity_band)) deallocate(optics%opacity_band)
1124  if (associated(optics%sw_pen_band)) deallocate(optics%sw_pen_band)
1125  endif ; endif
1126 
1127 end subroutine opacity_end
1128 
1129 !> \namespace mom_opacity
1130 !!
1131 !! opacity_from_chl:
1132 !! In this routine, the Morel (modified) or Manizza (modified)
1133 !! schemes use the "blue" band in the paramterizations to determine
1134 !! the e-folding depth of the incoming shortwave attenuation. The red
1135 !! portion is lumped into the net heating at the surface.
1136 !!
1137 !! Morel, A., 1988: Optical modeling of the upper ocean in relation
1138 !! to its biogenous matter content (case-i waters)., J. Geo. Res.,
1139 !! 93, 10,749-10,768.
1140 !!
1141 !! Manizza, M., C. LeQuere, A. J. Watson, and E. T. Buitenhuis, 2005:
1142 !! Bio-optical feedbacks among phytoplankton, upper ocean physics
1143 !! and sea-ice in a global model, Geophys. Res. Let., 32, L05603,
1144 !! doi:10.1029/2004GL020778.
1145 
1146 end module mom_opacity
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_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.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_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
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_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
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_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_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