MOM6
MOM_opacity.F90
Go to the documentation of this file.
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
8 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
11 use mom_grid, only : ocean_grid_type
15 
16 implicit none ; private
17 
18 #include <MOM_memory.h>
19 
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 [R-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 [R 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 R2 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%L_to_Z**2*gv%g_Earth * gv%H_to_RZ) * gv%H_to_RZ
622  else
623  g_hconv2 = us%L_to_Z**2*gv%g_Earth * gv%H_to_RZ**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  do i=is,ie
845  netpen(i,1) = 0.
846  do n=1,max(nsw,1)
847  netpen(i,1) = netpen(i,1) + pen_sw_bnd(n,i) ! Surface interface
848  enddo
849  enddo
850 
851  ! Apply penetrating SW radiation to remaining parts of layers.
852  ! Excessively thin layers are not heated to avoid runaway temps.
853  do k=1,nz
854 
855  do i=is,ie
856  netpen(i,k+1) = 0.
857 
858  if (h(i,k) > 0.0) then
859  do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
860  ! SW_trans is the SW that is transmitted THROUGH the layer
861  opt_depth = h(i,k)*gv%H_to_m * optics%opacity_band(n,i,j,k)
862  exp_od = exp(-opt_depth)
863  sw_trans = exp_od
864 
865  ! Heating at a very small rate can be absorbed by a sufficiently thick layer or several
866  ! thin layers without further penetration.
867  if (optics%answers_2018) then
868  if (nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat*min(1.0, i_habs*h(i,k)) ) sw_trans = 0.0
869  elseif ((nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat) .and. (h(i,k) > h_min_heat)) then
870  if (nsw*pen_sw_bnd(n,i) <= min_sw_heat * (i_habs*(h(i,k) - h_min_heat))) then
871  sw_trans = 0.0
872  else
873  sw_trans = min(sw_trans, &
874  1.0 - (min_sw_heat*(i_habs*(h(i,k) - h_min_heat))) / (nsw*pen_sw_bnd(n,i)))
875  endif
876  endif
877 
878  pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
879  netpen(i,k+1) = netpen(i,k+1) + pen_sw_bnd(n,i)
880  endif ; enddo
881  endif ! h(i,k) > 0.0
882 
883  ! Add to the accumulated thickness above that could be heated.
884  ! Only layers greater than h_min_heat thick should get heated.
885  if (h(i,k) >= 2.0*h_min_heat) then
886  h_heat(i) = h_heat(i) + h(i,k)
887  elseif (h(i,k) > h_min_heat) then
888  h_heat(i) = h_heat(i) + (2.0*h(i,k) - 2.0*h_min_heat)
889  endif
890  enddo ! i loop
891  enddo ! k loop
892 
893  if (absorballsw) then
894 
895  ! If there is still shortwave radiation at this point, it could go into
896  ! the bottom (with a bottom mud model), or it could be redistributed back
897  ! through the water column.
898  do i=is,ie
899  pen_sw_rem(i) = pen_sw_bnd(1,i)
900  do n=2,nsw ; pen_sw_rem(i) = pen_sw_rem(i) + pen_sw_bnd(n,i) ; enddo
901  enddo
902  do i=is,ie ; if (pen_sw_rem(i) > 0.0) sw_remains = .true. ; enddo
903 
904  ih_limit = 1.0 / h_limit_fluxes
905  do i=is,ie ; if ((pen_sw_rem(i) > 0.0) .and. (h_heat(i) > 0.0)) then
906  if (h_heat(i)*ih_limit < 1.0) then
907  unabsorbed = 1.0 - h_heat(i)*ih_limit
908  else
909  unabsorbed = 0.0
910  endif
911  do n=1,nsw ; pen_sw_bnd(n,i) = unabsorbed * pen_sw_bnd(n,i) ; enddo
912  endif ; enddo
913 
914  endif ! absorbAllSW
915 
916 end subroutine sumswoverbands
917 
918 
919 
920 !> This routine initalizes the opacity module, including an optics_type.
921 subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
922  type(time_type), target, intent(in) :: time !< The current model time.
923  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
924  type(verticalgrid_type), intent(in) :: gv !< model vertical grid structure
925  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
926  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
927  !! parameters.
928  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
929  !! output.
930  type(opacity_cs), pointer :: cs !< A pointer that is set to point to the control
931  !! structure for this module.
932  type(optics_type), pointer :: optics !< An optics structure that has parameters
933  !! set and arrays allocated here.
934  ! Local variables
935  character(len=200) :: tmpstr
936  character(len=40) :: mdl = "MOM_opacity"
937  character(len=40) :: bandnum, shortname
938  character(len=200) :: longname
939  character(len=40) :: scheme_string
940  ! This include declares and sets the variable "version".
941 # include "version_variable.h"
942  real :: pensw_absorb_minthick ! A thickness that is used to absorb the remaining shortwave heat
943  ! flux when that flux drops below PEN_SW_FLUX_ABSORB [m].
944  real :: pensw_minthick_dflt ! The default for PenSW_absorb_minthick [m]
945  logical :: default_2018_answers
946  logical :: use_scheme
947  integer :: isd, ied, jsd, jed, nz, n
948  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
949 
950  if (associated(cs)) then
951  call mom_error(warning, "opacity_init called with an associated"// &
952  "associated control structure.")
953  return
954  else ; allocate(cs) ; endif
955 
956  cs%diag => diag
957 
958  ! Read all relevant parameters and write them to the model log.
959  call log_version(param_file, mdl, version, '')
960 
961 ! parameters for CHL_A routines
962  call get_param(param_file, mdl, "VAR_PEN_SW", cs%var_pen_sw, &
963  "If true, use one of the CHL_A schemes specified by "//&
964  "OPACITY_SCHEME to determine the e-folding depth of "//&
965  "incoming short wave radiation.", default=.false.)
966 
967  cs%opacity_scheme = no_scheme ; scheme_string = ''
968  if (cs%var_pen_sw) then
969  call get_param(param_file, mdl, "OPACITY_SCHEME", tmpstr, &
970  "This character string specifies how chlorophyll "//&
971  "concentrations are translated into opacities. Currently "//&
972  "valid options include:\n"//&
973  " \t\t MANIZZA_05 - Use Manizza et al., GRL, 2005. \n"//&
974  " \t\t MOREL_88 - Use Morel, JGR, 1988.", &
975  default=manizza_05_string)
976  if (len_trim(tmpstr) > 0) then
977  tmpstr = uppercase(tmpstr)
978  select case (tmpstr)
979  case (manizza_05_string)
980  cs%opacity_scheme = manizza_05 ; scheme_string = manizza_05_string
981  case (morel_88_string)
982  cs%opacity_scheme = morel_88 ; scheme_string = morel_88_string
983  case default
984  call mom_error(fatal, "opacity_init: #DEFINE OPACITY_SCHEME "//&
985  trim(tmpstr) // "in input file is invalid.")
986  end select
987  call mom_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5)
988  endif
989  if (cs%opacity_scheme == no_scheme) then
990  call mom_error(warning, "opacity_init: No scheme has successfully "//&
991  "been specified for the opacity. Using the default MANIZZA_05.")
992  cs%opacity_scheme = manizza_05 ; scheme_string = manizza_05_string
993  endif
994 
995  call get_param(param_file, mdl, "BLUE_FRAC_SW", cs%blue_frac, &
996  "The fraction of the penetrating shortwave radiation "//&
997  "that is in the blue band.", default=0.5, units="nondim")
998  else
999  call get_param(param_file, mdl, "EXP_OPACITY_SCHEME", tmpstr, &
1000  "This character string specifies which exponential "//&
1001  "opacity scheme to utilize. Currently "//&
1002  "valid options include:\n"//&
1003  " \t\t SINGLE_EXP - Single Exponent decay. \n"//&
1004  " \t\t DOUBLE_EXP - Double Exponent decay.", &
1005  default=single_exp_string)!New default for "else" above (non-Chl scheme)
1006  if (len_trim(tmpstr) > 0) then
1007  tmpstr = uppercase(tmpstr)
1008  select case (tmpstr)
1009  case (single_exp_string)
1010  cs%opacity_scheme = single_exp ; scheme_string = single_exp_string
1011  case (double_exp_string)
1012  cs%opacity_scheme = double_exp ; scheme_string = double_exp_string
1013  end select
1014  call mom_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5)
1015  endif
1016  call get_param(param_file, mdl, "PEN_SW_SCALE", cs%pen_sw_scale, &
1017  "The vertical absorption e-folding depth of the "//&
1018  "penetrating shortwave radiation.", units="m", default=0.0)
1019  !BGR/ Added for opacity_scheme==double_exp read in 2nd exp-decay and fraction
1020  if (cs%Opacity_scheme == double_exp ) then
1021  call get_param(param_file, mdl, "PEN_SW_SCALE_2ND", cs%pen_sw_scale_2nd, &
1022  "The (2nd) vertical absorption e-folding depth of the "//&
1023  "penetrating shortwave radiation "//&
1024  "(use if SW_EXP_MODE==double.)",&
1025  units="m", default=0.0)
1026  call get_param(param_file, mdl, "SW_1ST_EXP_RATIO", cs%sw_1st_exp_ratio, &
1027  "The fraction of 1st vertical absorption e-folding depth "//&
1028  "penetrating shortwave radiation if SW_EXP_MODE==double.",&
1029  units="m", default=0.0)
1030  elseif (cs%OPACITY_SCHEME == single_exp) then
1031  !/Else disable 2nd_exp scheme
1032  cs%pen_sw_scale_2nd = 0.0
1033  cs%sw_1st_exp_ratio = 1.0
1034  endif
1035  call get_param(param_file, mdl, "PEN_SW_FRAC", cs%pen_sw_frac, &
1036  "The fraction of the shortwave radiation that penetrates "//&
1037  "below the surface.", units="nondim", default=0.0)
1038 
1039  endif
1040  call get_param(param_file, mdl, "PEN_SW_NBANDS", optics%nbands, &
1041  "The number of bands of penetrating shortwave radiation.", &
1042  default=1)
1043  if (cs%Opacity_scheme == double_exp ) then
1044  if (optics%nbands /= 2) call mom_error(fatal, &
1045  "set_opacity: \Cannot use a double_exp opacity scheme with nbands!=2.")
1046  elseif (cs%Opacity_scheme == single_exp ) then
1047  if (optics%nbands /= 1) call mom_error(fatal, &
1048  "set_opacity: \Cannot use a single_exp opacity scheme with nbands!=1.")
1049  endif
1050 
1051  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1052  "This sets the default value for the various _2018_ANSWERS parameters.", &
1053  default=.true.)
1054  call get_param(param_file, mdl, "OPTICS_2018_ANSWERS", optics%answers_2018, &
1055  "If true, use the order of arithmetic and expressions that recover the "//&
1056  "answers from the end of 2018. Otherwise, use updated expressions for "//&
1057  "handling the absorption of small remaining shortwave fluxes.", &
1058  default=default_2018_answers)
1059 
1060  call get_param(param_file, mdl, "PEN_SW_FLUX_ABSORB", optics%PenSW_flux_absorb, &
1061  "A minimum remaining shortwave heating rate that will be simply absorbed in "//&
1062  "the next sufficiently thick layers for computational efficiency, instead of "//&
1063  "continuing to penetrate. The default, 2.5e-11 degC m s-1, is about 1e-4 W m-2 "//&
1064  "or 0.08 degC m century-1, but 0 is also a valid value.", &
1065  default=2.5e-11, units="degC m s-1", scale=gv%m_to_H*us%T_to_s)
1066 
1067  if (optics%answers_2018) then ; pensw_minthick_dflt = 0.001 ; else ; pensw_minthick_dflt = 1.0 ; endif
1068  call get_param(param_file, mdl, "PEN_SW_ABSORB_MINTHICK", pensw_absorb_minthick, &
1069  "A thickness that is used to absorb the remaining penetrating shortwave heat "//&
1070  "flux when it drops below PEN_SW_FLUX_ABSORB.", &
1071  default=pensw_minthick_dflt, units="m", scale=gv%m_to_H)
1072  optics%PenSW_absorb_Invlen = 1.0 / (pensw_absorb_minthick + gv%H_subroundoff)
1073 
1074  if (.not.associated(optics%min_wavelength_band)) &
1075  allocate(optics%min_wavelength_band(optics%nbands))
1076  if (.not.associated(optics%max_wavelength_band)) &
1077  allocate(optics%max_wavelength_band(optics%nbands))
1078 
1079  if (cs%opacity_scheme == manizza_05) then
1080  optics%min_wavelength_band(1) =0
1081  optics%max_wavelength_band(1) =550
1082  if (optics%nbands >= 2) then
1083  optics%min_wavelength_band(2)=550
1084  optics%max_wavelength_band(2)=700
1085  endif
1086  if (optics%nbands > 2) then
1087  do n=3,optics%nbands
1088  optics%min_wavelength_band(n) =700
1089  optics%max_wavelength_band(n) =2800
1090  enddo
1091  endif
1092  endif
1093 
1094  call get_param(param_file, mdl, "OPACITY_LAND_VALUE", cs%opacity_land_value, &
1095  "The value to use for opacity over land. The default is "//&
1096  "10 m-1 - a value for muddy water.", units="m-1", default=10.0)
1097 
1098  if (.not.associated(optics%opacity_band)) &
1099  allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz))
1100  if (.not.associated(optics%sw_pen_band)) &
1101  allocate(optics%sw_pen_band(optics%nbands,isd:ied,jsd:jed))
1102  allocate(cs%id_opacity(optics%nbands)) ; cs%id_opacity(:) = -1
1103 
1104  cs%id_sw_pen = register_diag_field('ocean_model', 'SW_pen', diag%axesT1, time, &
1105  'Penetrating shortwave radiation flux into ocean', 'W m-2')
1106  cs%id_sw_vis_pen = register_diag_field('ocean_model', 'SW_vis_pen', diag%axesT1, time, &
1107  'Visible penetrating shortwave radiation flux into ocean', 'W m-2')
1108  do n=1,optics%nbands
1109  write(bandnum,'(i3)') n
1110  shortname = 'opac_'//trim(adjustl(bandnum))
1111  longname = 'Opacity for shortwave radiation in band '//trim(adjustl(bandnum)) &
1112  // ', saved as L^-1 tanh(Opacity * L) for L = 10^-10 m'
1113  cs%id_opacity(n) = register_diag_field('ocean_model', shortname, diag%axesTL, time, &
1114  longname, 'm-1')
1115  enddo
1116 
1117 end subroutine opacity_init
1118 
1119 
1120 subroutine opacity_end(CS, optics)
1121  type(opacity_cs), pointer :: cs !< An opacity control structure that should be deallocated.
1122  type(optics_type), optional, pointer :: optics !< An optics type structure that should be deallocated.
1123 
1124  if (associated(cs%id_opacity)) deallocate(cs%id_opacity)
1125  if (associated(cs)) deallocate(cs)
1126 
1127  if (present(optics)) then ; if (associated(optics)) then
1128  if (associated(optics%opacity_band)) deallocate(optics%opacity_band)
1129  if (associated(optics%sw_pen_band)) deallocate(optics%sw_pen_band)
1130  endif ; endif
1131 
1132 end subroutine opacity_end
1133 
1134 !> \namespace mom_opacity
1135 !!
1136 !! opacity_from_chl:
1137 !! In this routine, the Morel (modified) or Manizza (modified)
1138 !! schemes use the "blue" band in the paramterizations to determine
1139 !! the e-folding depth of the incoming shortwave attenuation. The red
1140 !! portion is lumped into the net heating at the surface.
1141 !!
1142 !! Morel, A., 1988: Optical modeling of the upper ocean in relation
1143 !! to its biogenous matter content (case-i waters)., J. Geo. Res.,
1144 !! 93, 10,749-10,768.
1145 !!
1146 !! Manizza, M., C. LeQuere, A. J. Watson, and E. T. Buitenhuis, 2005:
1147 !! Bio-optical feedbacks among phytoplankton, upper ocean physics
1148 !! and sea-ice in a global model, Geophys. Res. Let., 32, L05603,
1149 !! doi:10.1029/2004GL020778.
1150 
1151 end module mom_opacity
mom_opacity::manizza_05_string
character *(10), parameter manizza_05_string
String to specify the opacity scheme.
Definition: MOM_opacity.F90:80
mom_opacity::sumswoverbands
subroutine, public sumswoverbands(G, GV, US, h, nsw, optics, j, dt, H_limit_fluxes, absorbAllSW, iPen_SW_bnd, netPen)
This subroutine calculates the total shortwave heat flux integrated over bands as a function of depth...
Definition: MOM_opacity.F90:783
mom_opacity::set_opacity
subroutine, public set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, G, GV, CS, chl_2d, chl_3d)
This sets the opacity of sea water based based on one of several different schemes.
Definition: MOM_opacity.F90:93
mom_opacity::absorbremainingsw
subroutine, public absorbremainingsw(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_limit_fluxes, adjustAbsorptionProfile, absorbAllSW, T, Pen_SW_bnd, eps, ksort, htot, Ttot, TKE, dSV_dT)
Apply shortwave heating below the boundary layer (when running with the bulk mixed layer inhereted fr...
Definition: MOM_opacity.F90:512
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Call this subroutine to determine whether the averaging is currently enabled. .true....
Definition: MOM_diag_mediator.F90:1850
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_opacity::opacity_init
subroutine, public opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
This routine initalizes the opacity module, including an optics_type.
Definition: MOM_opacity.F90:922
mom_opacity::op_diag_len
real, parameter op_diag_len
Lengthscale L used to remap opacity from op to 1/L * tanh(op * L)
Definition: MOM_opacity.F90:85
mom_opacity::morel_88_string
character *(10), parameter morel_88_string
String to specify the opacity scheme.
Definition: MOM_opacity.F90:81
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:78
mom_opacity::extract_optics_slice
subroutine, public extract_optics_slice(optics, j, G, GV, opacity, opacity_scale, penSW_top, penSW_scale)
This subroutine returns a 2-d slice at constant j of fields from an optics_type, with the potential f...
Definition: MOM_opacity.F90:447
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_string_functions
Handy functions for manipulating strings.
Definition: MOM_string_functions.F90:2
mom_opacity::manizza_05
integer, parameter manizza_05
Coded integers to specify the opacity scheme.
Definition: MOM_opacity.F90:77
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
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_string_functions::uppercase
character(len=len(input_string)) function, public uppercase(input_string)
Return a string in which all uppercase letters have been replaced by their lowercase counterparts.
Definition: MOM_string_functions.F90:42
mom_opacity::double_exp
integer, parameter double_exp
Coded integers to specify the opacity scheme.
Definition: MOM_opacity.F90:77
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_opacity::double_exp_string
character *(10), parameter double_exp_string
String to specify the opacity scheme.
Definition: MOM_opacity.F90:83
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_opacity::opacity_manizza
real function, public opacity_manizza(chl_data)
This sets the blue-wavelength opacity according to the scheme proposed by Manizza,...
Definition: MOM_opacity.F90:437
mom_opacity::sw_pen_frac_morel
real function sw_pen_frac_morel(chl_data)
This sets the penetrating shortwave fraction according to the scheme proposed by Morel and Antoine (1...
Definition: MOM_opacity.F90:417
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_opacity::single_exp_string
character *(10), parameter single_exp_string
String to specify the opacity scheme.
Definition: MOM_opacity.F90:82
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_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_opacity::no_scheme
integer, parameter no_scheme
Coded integers to specify the opacity scheme.
Definition: MOM_opacity.F90:77
mom_opacity::opacity_morel
real function, public opacity_morel(chl_data)
This sets the blue-wavelength opacity according to the scheme proposed by Morel and Antoine (1994).
Definition: MOM_opacity.F90:397
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::opacity_from_chl
subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, G, GV, CS, chl_2d, chl_3d)
This sets the "blue" band opacity based on chloophyll A concencentrations The red portion is lumped i...
Definition: MOM_opacity.F90:222
mom_opacity::opacity_end
subroutine, public opacity_end(CS, optics)
Definition: MOM_opacity.F90:1121
mom_opacity::optics_nbands
integer function, public optics_nbands(optics)
Return the number of bands of penetrating shortwave radiation.
Definition: MOM_opacity.F90:495
mom_opacity::extract_optics_fields
subroutine, public extract_optics_fields(optics, nbands)
Set arguments to fields from the optics type.
Definition: MOM_opacity.F90:485
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_opacity::morel_88
integer, parameter morel_88
Coded integers to specify the opacity scheme.
Definition: MOM_opacity.F90:77
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:239
mom_opacity::single_exp
integer, parameter single_exp
Coded integers to specify the opacity scheme.
Definition: MOM_opacity.F90:77