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