MOM_opacity.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Routines used to calculate the opacity of the ocean.
6module mom_opacity
7
8use mom_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data
10use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
11use mom_file_parser, only : get_param, log_param, log_version, param_file_type
13use mom_grid, only : ocean_grid_type
17
18implicit none ; private
19
20#include <MOM_memory.h>
21
25
26!> This type is used to store information about ocean optical properties
27type, public :: optics_type
28 integer :: nbands !< The number of penetrating bands of SW radiation
29
30 real, allocatable :: opacity_band(:,:,:,:) !< SW optical depth per unit thickness [Z-1 ~> m-1]
31 !! The number of radiation bands is most rapidly varying (first) index.
32
33 real, allocatable :: sw_pen_band(:,:,:) !< shortwave radiation [Q R Z T-1 ~> W m-2]
34 !! at the surface in each of the nbands bands that penetrates beyond the surface.
35 !! The most rapidly varying dimension is the band.
36
37 real, allocatable :: min_wavelength_band(:)
38 !< The minimum wavelength in each band of penetrating shortwave radiation [nm]
39 real, allocatable :: max_wavelength_band(:)
40 !< The maximum wavelength in each band of penetrating shortwave radiation [nm]
41
42 real :: pensw_flux_absorb !< A heat flux that is small enough to be completely absorbed in the next
43 !! sufficiently thick layer [C H T-1 ~> degC m s-1 or degC kg m-2 s-1].
44 real :: pensw_absorb_invlen !< The inverse of the thickness that is used to absorb the remaining
45 !! shortwave heat flux when it drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2].
46
47 !! Lookup tables for Ohlmann solar penetration scheme
48 !! These would naturally exist as private module variables but that is prohibited in MOM6
49 real :: dlog10chl !< Chl increment within lookup table [log10 of Chl in mg m-3]
50 real :: chl_min !< Lower bound of Chl in lookup table [mg m-3]
51 real :: log10chl_min !< Lower bound of Chl in lookup table [log10 of Chl in mg m-3]
52 real :: log10chl_max !< Upper bound of Chl in lookup table [log10 of Chl in mg m-3]
53 real, allocatable, dimension(:) :: a1_lut,& !< Coefficient for band 1 [nondim]
54 & a2_lut,& !< Coefficient for band 2 [nondim]
55 & b1_lut,& !< Exponential decay scale for band 1 [Z-1 ~> m-1]
56 & b2_lut !< Exponential decay scale for band 2 [Z-1 ~> m-1]
57
58 integer :: answer_date !< The vintage of the order of arithmetic and expressions in the optics
59 !! calculations. Values below 20190101 recover the answers from the
60 !! end of 2018, while higher values use updated and more robust
61 !! forms of the same expressions.
62
63end type optics_type
64
65!> The control structure with parameters for the MOM_opacity module
66type, public :: opacity_cs ; private
67 logical :: var_pen_sw !< If true, use one of the CHL_A schemes (specified by OPACITY_SCHEME) to
68 !! determine the e-folding depth of incoming shortwave radiation.
69 integer :: opacity_scheme !< An integer indicating which scheme should be used to translate
70 !! water properties into the opacity (i.e., the e-folding depth) and
71 !! (perhaps) the number of bands of penetrating shortwave radiation to use.
72 real :: pen_sw_scale !< The vertical absorption e-folding depth of the
73 !! penetrating shortwave radiation [Z ~> m].
74 real :: pen_sw_scale_2nd !< The vertical absorption e-folding depth of the
75 !! (2nd) penetrating shortwave radiation [Z ~> m].
76 real :: sw_1st_exp_ratio !< Ratio for 1st exp decay in Two Exp decay opacity [nondim]
77 real :: pen_sw_frac !< The fraction of shortwave radiation that is
78 !! penetrating with a constant e-folding approach [nondim]
79 real :: blue_frac !< The fraction of the penetrating shortwave
80 !! radiation that is in the blue band [nondim].
81 real :: opacity_land_value !< The value to use for opacity over land [Z-1 ~> m-1].
82 !! The default is 10 m-1 - a value for muddy water.
83 real, allocatable, dimension(:,:) &
84 :: opacity_coef !< Groups of coefficients, in [Z-1 ~> m-1] or [Z ~> m] depending on the
85 !! scheme, in expressions for opacity, with the second index being the
86 !! wavelength band. For example, when OPACITY_SCHEME = MANIZZA_05,
87 !! these are coef_1 and coef_2 in the
88 !! expression opacity = coef_1 + coef_2 * chl**pow.
89 real, allocatable, dimension(:) &
90 :: sw_pen_frac_coef !< Coefficients in the expression for the penetrating shortwave
91 !! fracetion [nondim]
92 real, allocatable, dimension(:) &
93 :: chl_power !< Powers of chlorophyll [nondim] for each band for expressions for
94 !! opacity of the form opacity = coef_1 + coef_2 * chl**pow.
95 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
96 !! regulate the timing of diagnostic output.
97 integer :: chl_dep_bands !< The number of bands that depend on the Chlorophyll concentrations.
98 logical :: warning_issued !< A flag that is used to avoid repetitive warnings.
99
100 !>@{ Diagnostic IDs
101 integer :: id_sw_pen = -1, id_sw_vis_pen = -1
102 integer, allocatable :: id_opacity(:)
103 !>@}
104end type opacity_cs
105
106!>@{ Coded integers to specify the opacity scheme
107integer, parameter :: no_scheme = 0, manizza_05 = 1, morel_88 = 2, single_exp = 3, double_exp = 4,&
108 & ohlmann_03 = 5
109!>@}
110
111character*(10), parameter :: manizza_05_string = "MANIZZA_05" !< String to specify the opacity scheme
112character*(10), parameter :: morel_88_string = "MOREL_88" !< String to specify the opacity scheme
113character*(10), parameter :: ohlmann_03_string = "OHLMANN_03" !< String to specify the opacity scheme
114character*(10), parameter :: single_exp_string = "SINGLE_EXP" !< String to specify the opacity scheme
115character*(10), parameter :: double_exp_string = "DOUBLE_EXP" !< String to specify the opacity scheme
116
117contains
118
119!> This sets the opacity of sea water based based on one of several different schemes.
120subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
121 G, GV, US, CS, chl_2d, chl_3d)
122 type(optics_type), intent(inout) :: optics !< An optics structure that has values
123 !! set based on the opacities.
124 real, dimension(:,:), pointer :: sw_total !< Total shortwave flux into the ocean [Q R Z T-1 ~> W m-2]
125 real, dimension(:,:), pointer :: sw_vis_dir !< Visible, direct shortwave into the ocean [Q R Z T-1 ~> W m-2]
126 real, dimension(:,:), pointer :: sw_vis_dif !< Visible, diffuse shortwave into the ocean [Q R Z T-1 ~> W m-2]
127 real, dimension(:,:), pointer :: sw_nir_dir !< Near-IR, direct shortwave into the ocean [Q R Z T-1 ~> W m-2]
128 real, dimension(:,:), pointer :: sw_nir_dif !< Near-IR, diffuse shortwave into the ocean [Q R Z T-1 ~> W m-2]
129 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
130 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
131 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
132 type(opacity_cs) :: cs !< The control structure earlier set up by opacity_init.
133 real, dimension(SZI_(G),SZJ_(G)), &
134 optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentrations [mg m-3]
135 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
136 optional, intent(in) :: chl_3d !< The chlorophyll-A concentrations of each layer [mg m-3]
137
138 ! Local variables
139 integer :: i, j, k, n, is, ie, js, je, nz
140 real :: inv_sw_pen_scale ! The inverse of the e-folding scale [Z-1 ~> m-1].
141 real :: inv_nbands ! The inverse of the number of bands of penetrating
142 ! shortwave radiation [nondim]
143 real :: tmp(szi_(g),szj_(g),szk_(gv)) ! A 3-d temporary array for diagnosing opacity [Z-1 ~> m-1]
144 real :: pen_sw_tot(szi_(g),szj_(g)) ! The penetrating shortwave radiation
145 ! summed across all bands [Q R Z T-1 ~> W m-2].
146 real :: op_diag_len ! A tiny lengthscale [Z ~> m] used to remap diagnostics of opacity
147 ! from op to 1/op_diag_len * tanh(op * op_diag_len)
148 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
149
150 if (present(chl_2d) .or. present(chl_3d)) then
151 ! The optical properties are based on chlorophyll concentrations.
152 call opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
153 g, gv, us, cs, chl_2d, chl_3d)
154 else ! Use sw e-folding scale set by MOM_input
155 if (optics%nbands <= 1) then ; inv_nbands = 1.0
156 else ; inv_nbands = 1.0 / real(optics%nbands) ; endif
157
158 ! Make sure there is no division by 0.
159 inv_sw_pen_scale = 1.0 / max(cs%pen_sw_scale, 0.1*gv%Angstrom_Z, &
160 gv%dZ_subroundoff)
161 if ( cs%Opacity_scheme == double_exp ) then
162 !$OMP parallel do default(shared)
163 do k=1,nz ; do j=js,je ; do i=is,ie
164 optics%opacity_band(1,i,j,k) = inv_sw_pen_scale
165 optics%opacity_band(2,i,j,k) = 1.0 / max(cs%pen_sw_scale_2nd, &
166 0.1*gv%Angstrom_Z, gv%dZ_subroundoff)
167 enddo ; enddo ; enddo
168 if (.not.associated(sw_total) .or. (cs%pen_SW_scale <= 0.0)) then
169 !$OMP parallel do default(shared)
170 do j=js,je ; do i=is,ie ; do n=1,optics%nbands
171 optics%sw_pen_band(n,i,j) = 0.0
172 enddo ; enddo ; enddo
173 else
174 !$OMP parallel do default(shared)
175 do j=js,je ; do i=is,ie
176 optics%sw_pen_band(1,i,j) = (cs%SW_1st_EXP_RATIO) * sw_total(i,j)
177 optics%sw_pen_band(2,i,j) = (1.-cs%SW_1st_EXP_RATIO) * sw_total(i,j)
178 enddo ; enddo
179 endif
180 else
181 do k=1,nz ; do j=js,je ; do i=is,ie ; do n=1,optics%nbands
182 optics%opacity_band(n,i,j,k) = inv_sw_pen_scale
183 enddo ; enddo ; enddo ; enddo
184 if (.not.associated(sw_total) .or. (cs%pen_SW_scale <= 0.0)) then
185 !$OMP parallel do default(shared)
186 do j=js,je ; do i=is,ie ; do n=1,optics%nbands
187 optics%sw_pen_band(n,i,j) = 0.0
188 enddo ; enddo ; enddo
189 else
190 !$OMP parallel do default(shared)
191 do j=js,je ; do i=is,ie ; do n=1,optics%nbands
192 optics%sw_pen_band(n,i,j) = cs%pen_SW_frac * inv_nbands * sw_total(i,j)
193 enddo ; enddo ; enddo
194 endif
195 endif
196 endif
197
198 if (query_averaging_enabled(cs%diag)) then
199 if (cs%id_sw_pen > 0) then
200 !$OMP parallel do default(shared)
201 do j=js,je ; do i=is,ie
202 pen_sw_tot(i,j) = 0.0
203 do n=1,optics%nbands
204 pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
205 enddo
206 enddo ; enddo
207 call post_data(cs%id_sw_pen, pen_sw_tot, cs%diag)
208 endif
209 if (cs%id_sw_vis_pen > 0) then
210 if (cs%opacity_scheme == manizza_05) then
211 !$OMP parallel do default(shared)
212 do j=js,je ; do i=is,ie
213 pen_sw_tot(i,j) = 0.0
214 do n=1,min(optics%nbands,2)
215 pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
216 enddo
217 enddo ; enddo
218 else
219 !$OMP parallel do default(shared)
220 do j=js,je ; do i=is,ie
221 pen_sw_tot(i,j) = 0.0
222 do n=1,optics%nbands
223 pen_sw_tot(i,j) = pen_sw_tot(i,j) + optics%sw_pen_band(n,i,j)
224 enddo
225 enddo ; enddo
226 endif
227 call post_data(cs%id_sw_vis_pen, pen_sw_tot, cs%diag)
228 endif
229 do n=1,optics%nbands ; if (cs%id_opacity(n) > 0) then
230 op_diag_len = 1.0e-10*us%m_to_Z ! A minimal extinction depth to constrain the range of opacity [Z ~> m]
231 !$OMP parallel do default(shared)
232 do k=1,nz ; do j=js,je ; do i=is,ie
233 ! Remap opacity (op) to 1/L * tanh(op * L) where L is one Angstrom.
234 ! This gives a nearly identical value when op << 1/L but allows one to
235 ! record the values even at reduced precision when opacity is huge (i.e. opaque).
236 tmp(i,j,k) = tanh(op_diag_len * optics%opacity_band(n,i,j,k)) / op_diag_len
237 enddo ; enddo ; enddo
238 call post_data(cs%id_opacity(n), tmp, cs%diag)
239 endif ; enddo
240 endif
241
242end subroutine set_opacity
243
244
245!> This sets the "blue" band opacity based on chlorophyll A concentrations
246!! The red portion is lumped into the net heating at the surface.
247subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, &
248 G, GV, US, CS, chl_2d, chl_3d)
249 type(optics_type), intent(inout) :: optics !< An optics structure that has values
250 !! set based on the opacities.
251 real, dimension(:,:), pointer :: sw_total !< Total shortwave flux into the ocean [Q R Z T-1 ~> W m-2]
252 real, dimension(:,:), pointer :: sw_vis_dir !< Visible, direct shortwave into the ocean [Q R Z T-1 ~> W m-2]
253 real, dimension(:,:), pointer :: sw_vis_dif !< Visible, diffuse shortwave into the ocean [Q R Z T-1 ~> W m-2]
254 real, dimension(:,:), pointer :: sw_nir_dir !< Near-IR, direct shortwave into the ocean [Q R Z T-1 ~> W m-2]
255 real, dimension(:,:), pointer :: sw_nir_dif !< Near-IR, diffuse shortwave into the ocean [Q R Z T-1 ~> W m-2]
256 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
257 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
258 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
259 type(opacity_cs) :: CS !< The control structure.
260 real, dimension(SZI_(G),SZJ_(G)), &
261 optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentrations [mg m-3]
262 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
263 optional, intent(in) :: chl_3d !< A 3-d field of chlorophyll-A concentrations [mg m-3]
264
265 real :: chl_data(SZI_(G),SZJ_(G)) ! The chlorophyll A concentrations in a layer [mg m-3].
266 real :: Inv_nbands ! The inverse of the number of bands of penetrating
267 ! shortwave radiation [nondim]
268 real :: Inv_nbands_nir ! The inverse of the number of bands of penetrating
269 ! near-infrared radiation [nondim]
270 real :: SW_pen_tot ! The sum across the bands of the penetrating
271 ! shortwave radiation [Q R Z T-1 ~> W m-2].
272 real :: SW_vis_tot ! The sum across the visible bands of shortwave
273 ! radiation [Q R Z T-1 ~> W m-2].
274 real :: SW_nir_tot ! The sum across the near infrared bands of shortwave
275 ! radiation [Q R Z T-1 ~> W m-2].
276 character(len=128) :: mesg
277 integer :: i, j, k, n, is, ie, js, je, nz, nbands
278 logical :: multiband_vis_input, multiband_nir_input, total_sw_input
279
280 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
281
282! In this model, the Morel (modified) and Manizza (modified) schemes
283! use the "blue" band in the parameterizations to determine the e-folding
284! depth of the incoming shortwave attenuation. The red portion is lumped
285! into the net heating at the surface.
286! Adding Ohlmann scheme. Needs sw_total and chl as inputs. Produces 2 penetrating bands.
287! This implementation follows that in CESM-POP using a lookup table in log10(chl) space.
288! The table is initialized in subroutine init_ohlmann and the coefficients are recovered
289! with routines lookup_ohlmann_swpen and lookup_ohlmann_opacity.
290! Note that this form treats the IR solar input implicitly: the sum of partitioning
291! coefficients < 1.0. The remainder is non-penetrating and is deposited in first layer
292! irrespective of thickness. The Ohlmann (2003) paper states that the scheme is not valid
293! for vertcal grids with first layer thickness < 2.0 meters.
294!
295! Ohlmann, J.C. Ocean radiant heating in climate models. J. Climate, 16, 1337-1351, 2003.
296!
297! Morel, A., Optical modeling of the upper ocean in relation to its biogenous
298! matter content (case-i waters)., J. Geo. Res., {93}, 10,749--10,768, 1988.
299!
300! Manizza, M., C. L. Quere, A. Watson, and E. T. Buitenhuis, Bio-optical
301! feedbacks among phytoplankton, upper ocean physics and sea-ice in a
302! global model, Geophys. Res. Let., , L05,603, 2005.
303
304 nbands = optics%nbands
305
306 if (nbands <= 1) then ; inv_nbands = 1.0
307 else ; inv_nbands = 1.0 / real(nbands) ; endif
308
309 if (nbands <= 2) then ; inv_nbands_nir = 0.0
310 else ; inv_nbands_nir = 1.0 / real(nbands - 2.0) ; endif
311
312 if (.not.(associated(sw_total) .or. (associated(sw_vis_dir) .and. associated(sw_vis_dif) .and. &
313 associated(sw_nir_dir) .and. associated(sw_nir_dif)) )) then
314 if (.not.cs%warning_issued) then
315 call mom_error(warning, &
316 "opacity_from_chl called without any shortwave flux arrays allocated.\n"//&
317 "Consider setting PEN_SW_NBANDS = 0 if no shortwave fluxes are being used.")
318 endif
319 cs%warning_issued = .true.
320 endif
321
322 multiband_vis_input = (associated(sw_vis_dir) .and. associated(sw_vis_dif))
323 multiband_nir_input = (associated(sw_nir_dir) .and. associated(sw_nir_dif))
324 total_sw_input = associated(sw_total)
325
326 chl_data(:,:) = 0.0
327 if (present(chl_3d)) then
328 do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,1) ; enddo ; enddo
329 do k=1,nz ; do j=js,je ; do i=is,ie
330 if ((g%mask2dT(i,j) > 0.0) .and. (chl_3d(i,j,k) < 0.0)) then
331 write(mesg,'(" Negative chl_3d of ",(1pe12.4)," found at i,j,k = ", &
332 & 3(1x,I0), " lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
333 chl_3d(i,j,k), i, j, k, g%geoLonT(i,j), g%geoLatT(i,j)
334 call mom_error(fatal, "MOM_opacity opacity_from_chl: "//trim(mesg))
335 endif
336 enddo ; enddo ; enddo
337 elseif (present(chl_2d)) then
338 do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_2d(i,j) ; enddo ; enddo
339 do j=js,je ; do i=is,ie
340 if ((g%mask2dT(i,j) > 0.0) .and. (chl_2d(i,j) < 0.0)) then
341 write(mesg,'(" Negative chl_2d of ",(1pe12.4)," at i,j = ", &
342 & I0,", ",I0," lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
343 chl_data(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
344 call mom_error(fatal, "MOM_opacity opacity_from_chl: "//trim(mesg))
345 endif
346 enddo ; enddo
347 else
348 call mom_error(fatal, "Either chl_2d or chl_3d must be present in a call to opacity_form_chl.")
349 endif
350
351 select case (cs%opacity_scheme)
352 case (manizza_05)
353 !$OMP parallel do default(shared) private(SW_vis_tot,SW_nir_tot)
354 do j=js,je ; do i=is,ie
355 sw_vis_tot = 0.0 ; sw_nir_tot = 0.0
356 if (g%mask2dT(i,j) > 0.0) then
357 if (multiband_vis_input) then
358 sw_vis_tot = sw_vis_dir(i,j) + sw_vis_dif(i,j)
359 elseif (total_sw_input) then
360 ! Follow Manizza 05 in assuming that 42% of SW is visible.
361 sw_vis_tot = 0.42 * sw_total(i,j)
362 endif
363 if (multiband_nir_input) then
364 sw_nir_tot = sw_nir_dir(i,j) + sw_nir_dif(i,j)
365 elseif (total_sw_input) then
366 sw_nir_tot = sw_total(i,j) - sw_vis_tot
367 endif
368 endif
369
370 ! Band 1 is Manizza blue.
371 optics%sw_pen_band(1,i,j) = cs%blue_frac*sw_vis_tot
372 ! Band 2 (if used) is Manizza red.
373 if (nbands > 1) &
374 optics%sw_pen_band(2,i,j) = (1.0-cs%blue_frac)*sw_vis_tot
375 ! All remaining bands are NIR, for lack of something better to do.
376 do n=3,nbands
377 optics%sw_pen_band(n,i,j) = inv_nbands_nir * sw_nir_tot
378 enddo
379 enddo ; enddo
380 case (morel_88)
381 !$OMP parallel do default(shared) private(SW_pen_tot)
382 do j=js,je ; do i=is,ie
383 sw_pen_tot = 0.0
384 if (g%mask2dT(i,j) > 0.0) then
385 if (multiband_vis_input) then
386 sw_pen_tot = sw_pen_frac_morel(chl_data(i,j), cs) * (sw_vis_dir(i,j) + sw_vis_dif(i,j))
387 elseif (total_sw_input) then
388 sw_pen_tot = sw_pen_frac_morel(chl_data(i,j), cs) * 0.5*sw_total(i,j)
389 endif
390 endif
391
392 do n=1,nbands
393 optics%sw_pen_band(n,i,j) = inv_nbands*sw_pen_tot
394 enddo
395 enddo ; enddo
396 case (ohlmann_03)
397 ! want exactly two penetrating bands. If not, throw an error.
398 if ( nbands /= 2 ) then
399 call mom_error(fatal, "opacity_from_chl: CS%opacity_scheme requires nbands==2.")
400 endif
401 !$OMP parallel do default(shared) private(SW_vis_tot)
402 do j=js,je ; do i=is,ie
403 sw_vis_tot = 0.0 ! Ohlmann does not classify as vis/nir. Using vis to add up total
404 if (g%mask2dT(i,j) < 0.5) then
405 optics%sw_pen_band(1:2,i,j) = 0. ! Make sure there is a valid value for land points
406 else
407 if (multiband_vis_input ) then ! If multiband_vis_input is true then so is multiband_nir_input
408 sw_vis_tot = ((sw_vis_dir(i,j) + sw_vis_dif(i,j)) + sw_nir_dir(i,j)) + sw_nir_dif(i,j)
409 elseif (total_sw_input) then
410 sw_vis_tot = sw_total(i,j)
411 else
412 call mom_error(fatal, "No shortwave input was provided.")
413 endif
414
415 ! Bands 1-2 (Ohlmann factors A with coefficients for Table 1a)
416 optics%sw_pen_band(1:2,i,j) = lookup_ohlmann_swpen(chl_data(i,j),optics)*sw_vis_tot
417 endif
418 enddo ; enddo
419 case default
420 call mom_error(fatal, "opacity_from_chl: CS%opacity_scheme is not valid.")
421 end select
422
423 !$OMP parallel do default(shared) firstprivate(chl_data)
424 do k=1,nz
425 !! FOB
426 !!! I don't think this is what we want to do with Ohlmann.
427 !!! The surface CHL is used in developing the parameterization.
428 !!! Only the surface CHL is used above in setting optics%sw_pen_band for all schemes.
429 !!! Seems inconsistent to use depth dependent CHL in opacity calculation.
430 !!! Nevertheless, leaving as is for now.
431 !! FOB
432 if (present(chl_3d)) then
433 do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,k) ; enddo ; enddo
434 endif
435
436 select case (cs%opacity_scheme)
437 case (manizza_05)
438 do j=js,je ; do i=is,ie
439 if (g%mask2dT(i,j) <= 0.5) then
440 do n=1,optics%nbands
441 optics%opacity_band(n,i,j,k) = cs%opacity_land_value
442 enddo
443 else
444 do n=1,cs%chl_dep_bands
445 optics%opacity_band(n,i,j,k) = cs%opacity_coef(1,n) + &
446 cs%opacity_coef(2,n) * chl_data(i,j)**cs%chl_power(n)
447 enddo
448 do n=cs%chl_dep_bands+1,optics%nbands ! These bands do not depend on the chlorophyll.
449 ! Any nonzero values that were in opacity_coef(2,n) have been added to opacity_coef(1,n).
450 optics%opacity_band(n,i,j,k) = cs%opacity_coef(1,n)
451 enddo
452 endif
453 enddo ; enddo
454 case (morel_88)
455 do j=js,je ; do i=is,ie
456 optics%opacity_band(1,i,j,k) = cs%opacity_land_value
457 if (g%mask2dT(i,j) > 0.0) &
458 optics%opacity_band(1,i,j,k) = opacity_morel(chl_data(i,j), cs)
459
460 do n=2,optics%nbands
461 optics%opacity_band(n,i,j,k) = optics%opacity_band(1,i,j,k)
462 enddo
463 enddo ; enddo
464 case (ohlmann_03)
465 !! not testing for 2 bands since we did it above
466 do j=js,je ; do i=is,ie
467 if (g%mask2dT(i,j) <= 0.5) then
468 optics%opacity_band(1:2,i,j,k) = cs%opacity_land_value
469 else
470 ! Bands 1-2 (Ohlmann factors B with coefficients for Table 1a
471 optics%opacity_band(1:2,i,j,k) = lookup_ohlmann_opacity(chl_data(i,j),optics) * us%Z_to_m
472 endif
473 enddo ; enddo
474 case default
475 call mom_error(fatal, "opacity_from_chl: CS%opacity_scheme is not valid.")
476 end select
477 enddo
478
479end subroutine opacity_from_chl
480
481!> This sets the blue-wavelength opacity according to the scheme proposed by
482!! Morel and Antoine (1994).
483function opacity_morel(chl_data, CS)
484 real, intent(in) :: chl_data !< The chlorophyll-A concentration in [mg m-3]
485 type(opacity_cs) :: cs !< Opacity control structure
486 real :: opacity_morel !< The returned opacity [Z-1 ~> m-1]
487
488 real :: chl, chl2 ! The log10 of chl_data (in mg m-3), and Chl^2 [nondim]
489
490 chl = log10(min(max(chl_data,0.02),60.0)) ; chl2 = chl*chl
491 ! All frequency bands currently use the same opacities.
492 opacity_morel = 1.0 / ( (cs%opacity_coef(1,1) + cs%opacity_coef(2,1)*chl) + chl2 * &
493 ((cs%opacity_coef(3,1) + chl*cs%opacity_coef(4,1)) + &
494 chl2*(cs%opacity_coef(5,1) + chl*cs%opacity_coef(6,1))) )
495end function opacity_morel
496
497!> This sets the penetrating shortwave fraction according to the scheme proposed by
498!! Morel and Antoine (1994).
499function sw_pen_frac_morel(chl_data, CS)
500 real, intent(in) :: chl_data !< The chlorophyll-A concentration [mg m-3]
501 type(opacity_cs) :: cs !< Opacity control structure
502 real :: sw_pen_frac_morel !< The returned penetrating shortwave fraction [nondim]
503
504 ! The following are coefficients for the optical model taken from Morel and
505 ! Antoine (1994). These coefficients represent a non uniform distribution of
506 ! chlorophyll-a through the water column. Other approaches may be more
507 ! appropriate when using an interactive ecosystem model that predicts
508 ! three-dimensional chl-a values.
509 real :: chl, chl2 ! The log10 of chl_data in mg m-3, and Chl^2 [nondim]
510
511 chl = log10(min(max(chl_data,0.02),60.0)) ; chl2 = chl*chl
512 sw_pen_frac_morel = 1.0 - ( (cs%SW_pen_frac_coef(1) + cs%SW_pen_frac_coef(2)*chl) + chl2 * &
513 ((cs%SW_pen_frac_coef(3) + chl*cs%SW_pen_frac_coef(4)) + &
514 chl2*(cs%SW_pen_frac_coef(5) + chl*cs%SW_pen_frac_coef(6))) )
515end function sw_pen_frac_morel
516
517!> This subroutine returns a 2-d slice at constant j of fields from an optics_type, with the potential
518!! for rescaling these fields.
519subroutine extract_optics_slice(optics, j, G, GV, opacity, opacity_scale, penSW_top, penSW_scale, SpV_avg)
520 type(optics_type), intent(in) :: optics !< An optics structure that has values of opacities
521 !! and shortwave fluxes.
522 integer, intent(in) :: j !< j-index to extract
523 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
524 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
525 real, dimension(max(optics%nbands,1),SZI_(G),SZK_(GV)), &
526 optional, intent(out) :: opacity !< The opacity in each band, i-point, and layer [Z-1 ~> m-1],
527 !! but with units that can be altered by opacity_scale
528 !! and the presence of SpV_avg to change this to other
529 !! units like [H-1 ~> m-1 or m2 kg-1]
530 real, optional, intent(in) :: opacity_scale !< A factor by which to rescale the opacity [nondim] or
531 !! [Z H-1 ~> 1 or m3 kg-1]
532 real, dimension(max(optics%nbands,1),SZI_(G)), &
533 optional, intent(out) :: pensw_top !< The shortwave radiation [Q R Z T-1 ~> W m-2]
534 !! at the surface in each of the nbands bands
535 !! that penetrates beyond the surface skin layer.
536 real, optional, intent(in) :: pensw_scale !< A factor by which to rescale the shortwave flux [nondim]
537 !! or other units.
538 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
539 optional, intent(in) :: spv_avg !< The layer-averaged specific volume [R-1 ~> m3 kg-1]
540 !! that is used along with opacity_scale in non-Boussinesq
541 !! cases to change the opacity from distance based units to
542 !! mass-based units
543
544 ! Local variables
545 real :: scale_opacity ! A rescaling factor for opacity [nondim], or the same units as opacity_scale.
546 real :: scale_pensw ! A rescaling factor for the penetrating shortwave radiation [nondim] or the
547 ! same units as penSW_scale
548 integer :: i, is, ie, k, nz, n
549 is = g%isc ; ie = g%iec ; nz = gv%ke
550
551 scale_opacity = 1.0 ; if (present(opacity_scale)) scale_opacity = opacity_scale
552 scale_pensw = 1.0 ; if (present(pensw_scale)) scale_pensw = pensw_scale
553
554 if (present(opacity)) then
555 if (present(spv_avg)) then
556 do k=1,nz ; do i=is,ie ; do n=1,optics%nbands
557 opacity(n,i,k) = (scale_opacity * spv_avg(i,j,k)) * optics%opacity_band(n,i,j,k)
558 enddo ; enddo ; enddo
559 else
560 do k=1,nz ; do i=is,ie ; do n=1,optics%nbands
561 opacity(n,i,k) = scale_opacity * optics%opacity_band(n,i,j,k)
562 enddo ; enddo ; enddo
563 endif
564 endif
565
566 if (present(pensw_top)) then ; do i=is,ie ; do n=1,optics%nbands
567 pensw_top(n,i) = scale_pensw * optics%sw_pen_band(n,i,j)
568 enddo ; enddo ; endif
569
570end subroutine extract_optics_slice
571
572!> Set arguments to fields from the optics type.
573subroutine extract_optics_fields(optics, nbands)
574 type(optics_type), intent(in) :: optics !< An optics structure that has values of opacities
575 !! and shortwave fluxes.
576 integer, optional, intent(out) :: nbands !< The number of penetrating bands of SW radiation
577
578 if (present(nbands)) nbands = optics%nbands
579
580end subroutine extract_optics_fields
581
582!> Return the number of bands of penetrating shortwave radiation.
583function optics_nbands(optics)
584 type(optics_type), pointer :: optics !< An optics structure that has values of opacities
585 !! and shortwave fluxes.
586 integer :: optics_nbands !< The number of penetrating bands of SW radiation
587
588 if (associated(optics)) then
589 optics_nbands = optics%nbands
590 else
591 optics_nbands = 0
592 endif
593end function optics_nbands
594
595!> Apply shortwave heating below the boundary layer (when running with the bulk mixed layer inherited
596!! from GOLD) or throughout the water column.
597!!
598!! In addition, it causes all of the remaining SW radiation to be absorbed, provided that the total
599!! water column thickness is greater than H_limit_fluxes.
600!! For thinner water columns, the heating is scaled down proportionately, the assumption being that the
601!! remaining heating (which is left in Pen_SW) should go into an (absent for now) ocean bottom sediment layer.
602subroutine absorbremainingsw(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_limit_fluxes, &
603 adjustAbsorptionProfile, absorbAllSW, T, Pen_SW_bnd, &
604 eps, ksort, htot, Ttot, TKE, dSV_dT)
605
606 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
607 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
608 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
609 integer, intent(in) :: nsw !< Number of bands of penetrating
610 !! shortwave radiation.
611 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
612 real, dimension(max(1,nsw),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< Opacity in each band of penetrating
613 !! shortwave radiation [H-1 ~> m-1 or m2 kg-1].
614 !! The indices are band, i, k.
615 type(optics_type), intent(in) :: optics !< An optics structure that has values of
616 !! opacities and shortwave fluxes.
617 integer, intent(in) :: j !< j-index to work on.
618 real, intent(in) :: dt !< Time step [T ~> s].
619 real, intent(in) :: h_limit_fluxes !< If the total ocean depth is
620 !! less than this, they are scaled away
621 !! to avoid numerical instabilities
622 !! [H ~> m or kg m-2]. This would
623 !! not be necessary if a finite heat
624 !! capacity mud-layer were added.
625 logical, intent(in) :: adjustabsorptionprofile !< If true, apply
626 !! heating above the layers in which it
627 !! should have occurred to get the
628 !! correct mean depth (and potential
629 !! energy change) of the shortwave that
630 !! should be absorbed by each layer.
631 logical, intent(in) :: absorballsw !< If true, apply heating above the
632 !! layers in which it should have occurred
633 !! to get the correct mean depth (and
634 !! potential energy change) of the
635 !! shortwave that should be absorbed by
636 !! each layer.
637 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: t !< Layer potential/conservative
638 !! temperatures [C ~> degC]
639 real, dimension(max(1,nsw),SZI_(G)), intent(inout) :: pen_sw_bnd !< Penetrating shortwave heating in
640 !! each band that hits the bottom and will
641 !! will be redistributed through the water
642 !! column [C H ~> degC m or degC kg m-2],
643 !! size nsw x SZI_(G).
644 real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: eps !< Small thickness that must remain in
645 !! each layer, and which will not be
646 !! subject to heating [H ~> m or kg m-2]
647 integer, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: ksort !< Density-sorted k-indices.
648 real, dimension(SZI_(G)), optional, intent(in) :: htot !< Total mixed layer thickness [H ~> m or kg m-2].
649 real, dimension(SZI_(G)), optional, intent(inout) :: ttot !< Depth integrated mixed layer
650 !! temperature [C H ~> degC m or degC kg m-2]
651 real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: dsv_dt !< The partial derivative of specific volume
652 !! with temperature [R-1 C-1 ~> m3 kg-1 degC-1]
653 real, dimension(SZI_(G),SZK_(GV)), optional, intent(inout) :: tke !< The TKE sink from mixing the heating
654 !! throughout a layer [R Z3 T-2 ~> J m-2].
655
656 ! Local variables
657 real, dimension(SZI_(G),SZK_(GV)) :: &
658 t_chg_above ! A temperature change that will be applied to all the thick
659 ! layers above a given layer [C ~> degC]. This is only nonzero if
660 ! adjustAbsorptionProfile is true, in which case the net
661 ! change in the temperature of a layer is the sum of the
662 ! direct heating of that layer plus T_chg_above from all of
663 ! the layers below, plus any contribution from absorbing
664 ! radiation that hits the bottom.
665 real, dimension(SZI_(G)) :: &
666 h_heat, & ! The thickness of the water column that will be heated by
667 ! any remaining shortwave radiation [H ~> m or kg m-2].
668 t_chg, & ! The temperature change of thick layers due to the remaining
669 ! shortwave radiation and contributions from T_chg_above [C ~> degC].
670 pen_sw_rem ! The sum across all wavelength bands of the penetrating shortwave
671 ! heating that hits the bottom and will be redistributed through
672 ! the water column [C H ~> degC m or degC kg m-2]
673 real :: sw_trans ! fraction of shortwave radiation that is not
674 ! absorbed in a layer [nondim]
675 real :: unabsorbed ! fraction of the shortwave radiation that
676 ! is not absorbed because the layers are too thin [nondim]
677 real :: ih_limit ! inverse of the total depth at which the
678 ! surface fluxes start to be limited [H-1 ~> m-1 or m2 kg-1]
679 real :: h_min_heat ! minimum thickness layer that should get heated [H ~> m or kg m-2]
680 real :: opt_depth ! optical depth of a layer [nondim]
681 real :: exp_od ! exp(-opt_depth) [nondim]
682 real :: heat_bnd ! heating due to absorption in the current
683 ! layer by the current band, including any piece that
684 ! is moved upward [C H ~> degC m or degC kg m-2]
685 real :: swa ! fraction of the absorbed shortwave that is
686 ! moved to layers above with adjustAbsorptionProfile [nondim]
687 real :: coswa_frac ! The fraction of SWa that is actually moved upward [nondim]
688 real :: min_sw_heat ! A minimum remaining shortwave heating within a timestep that will be simply
689 ! absorbed in the next layer for computational efficiency, instead of
690 ! continuing to penetrate [C H ~> degC m or degC kg m-2].
691 real :: i_habs ! The inverse of the absorption length for a minimal flux [H-1 ~> m-1 or m2 kg-1]
692 real :: epsilon ! A small thickness that must remain in each
693 ! layer, and which will not be subject to heating [H ~> m or kg m-2]
694 real :: g_hconv2 ! A conversion factor for use in the TKE calculation
695 ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
696 logical :: sw_remains ! If true, some column has shortwave radiation that
697 ! was not entirely absorbed.
698 logical :: tke_calc ! If true, calculate the implications to the
699 ! TKE budget of the shortwave heating.
700 real :: c1_6, c1_60 ! Rational fractions [nondim]
701 integer :: is, ie, nz, i, k, ks, n
702
703 if (nsw < 1) return
704
705 sw_remains = .false.
706 min_sw_heat = optics%PenSW_flux_absorb * dt
707 i_habs = optics%PenSW_absorb_Invlen
708
709 h_min_heat = 2.0*gv%Angstrom_H + gv%H_subroundoff
710 is = g%isc ; ie = g%iec ; nz = gv%ke
711 c1_6 = 1.0 / 6.0 ; c1_60 = 1.0 / 60.0
712
713 tke_calc = (present(tke) .and. present(dsv_dt))
714
715 if (optics%answer_date < 20190101) then
716 g_hconv2 = (gv%g_Earth_Z_T2 * gv%H_to_RZ) * gv%H_to_RZ
717 else
718 g_hconv2 = gv%g_Earth_Z_T2 * gv%H_to_RZ**2
719 endif
720
721 h_heat(:) = 0.0
722 if (present(htot)) then ; do i=is,ie ; h_heat(i) = htot(i) ; enddo ; endif
723
724 ! Apply penetrating SW radiation to remaining parts of layers.
725 ! Excessively thin layers are not heated to avoid runaway temps.
726 do ks=1,nz ; do i=is,ie
727 k = ks
728 if (present(ksort)) then
729 if (ksort(i,ks) <= 0) cycle
730 k = ksort(i,ks)
731 endif
732 epsilon = 0.0 ; if (present(eps)) epsilon = eps(i,k)
733
734 t_chg_above(i,k) = 0.0
735
736 if (h(i,k) > 1.5*epsilon) then
737 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
738 ! SW_trans is the SW that is transmitted THROUGH the layer
739 opt_depth = h(i,k) * opacity_band(n,i,k)
740 exp_od = exp(-opt_depth)
741 sw_trans = exp_od
742
743 ! Heating at a very small rate can be absorbed by a sufficiently thick layer or several
744 ! thin layers without further penetration.
745 if (optics%answer_date < 20190101) then
746 if (nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat*min(1.0, i_habs*h(i,k)) ) sw_trans = 0.0
747 elseif ((nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat) .and. (h(i,k) > h_min_heat)) then
748 if (nsw*pen_sw_bnd(n,i) <= min_sw_heat * (i_habs*(h(i,k) - h_min_heat))) then
749 sw_trans = 0.0
750 else
751 sw_trans = min(sw_trans, &
752 1.0 - (min_sw_heat*(i_habs*(h(i,k) - h_min_heat))) / (nsw*pen_sw_bnd(n,i)))
753 endif
754 endif
755
756 heat_bnd = pen_sw_bnd(n,i) * (1.0 - sw_trans)
757 if (adjustabsorptionprofile .and. (h_heat(i) > 0.0)) then
758 ! In this case, a fraction of the heating is applied to the
759 ! overlying water so that the mean pressure at which the shortwave
760 ! heating occurs is exactly what it would have been with a careful
761 ! pressure-weighted averaging of the exponential heating profile,
762 ! hence there should be no TKE budget requirements due to this
763 ! layer. Very clever, but this is also limited so that the
764 ! water above is not heated at a faster rate than the layer
765 ! actually being heated, i.e., SWA <= h_heat / (h_heat + h(i,k))
766 ! and takes the energetics of the rest of the heating into account.
767 ! (-RWH, ~7 years later.)
768 if (opt_depth > 1e-5) then
769 swa = ((opt_depth + (opt_depth + 2.0)*exp_od) - 2.0) / &
770 ((opt_depth + opacity_band(n,i,k) * h_heat(i)) * &
771 (1.0 - exp_od))
772 else
773 ! Use Taylor series expansion of the expression above for a
774 ! more accurate form with very small layer optical depths.
775 swa = h(i,k) * (opt_depth * (1.0 - opt_depth)) / &
776 ((h_heat(i) + h(i,k)) * (6.0 - 3.0*opt_depth))
777 endif
778 coswa_frac = 0.0
779 if (swa*(h_heat(i) + h(i,k)) > h_heat(i)) then
780 coswa_frac = (swa*(h_heat(i) + h(i,k)) - h_heat(i) ) / &
781 (swa*(h_heat(i) + h(i,k)))
782 swa = h_heat(i) / (h_heat(i) + h(i,k))
783 endif
784
785 t_chg_above(i,k) = t_chg_above(i,k) + (swa * heat_bnd) / h_heat(i)
786 t(i,k) = t(i,k) + ((1.0 - swa) * heat_bnd) / h(i,k)
787 else
788 coswa_frac = 1.0
789 t(i,k) = t(i,k) + pen_sw_bnd(n,i) * (1.0 - sw_trans) / h(i,k)
790 endif
791
792 if (tke_calc) then
793 if (opt_depth > 1e-2) then
794 tke(i,k) = tke(i,k) - coswa_frac*heat_bnd*dsv_dt(i,k)* &
795 (0.5*h(i,k)*g_hconv2) * &
796 (opt_depth*(1.0+exp_od) - 2.0*(1.0-exp_od)) / (opt_depth*(1.0-exp_od))
797 else
798 ! Use Taylor series-derived approximation to the above expression
799 ! that is well behaved and more accurate when opt_depth is small.
800 tke(i,k) = tke(i,k) - coswa_frac*heat_bnd*dsv_dt(i,k)* &
801 (0.5*h(i,k)*g_hconv2) * &
802 (c1_6*opt_depth * (1.0 - c1_60*opt_depth**2))
803 endif
804 endif
805
806 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
807 endif ; enddo
808 endif
809
810 ! Add to the accumulated thickness above that could be heated.
811 ! Only layers greater than h_min_heat thick should get heated.
812 if (h(i,k) >= 2.0*h_min_heat) then
813 h_heat(i) = h_heat(i) + h(i,k)
814 elseif (h(i,k) > h_min_heat) then
815 h_heat(i) = h_heat(i) + (2.0*h(i,k) - 2.0*h_min_heat)
816 endif
817 enddo ; enddo ! i & k loops
818
819! if (.not.absorbAllSW .and. .not.adjustAbsorptionProfile) return
820
821 ! Unless modified, there is no temperature change due to fluxes from the bottom.
822 do i=is,ie ; t_chg(i) = 0.0 ; enddo
823
824 if (absorballsw) then
825 ! If there is still shortwave radiation at this point, it could go into
826 ! the bottom (with a bottom mud model), or it could be redistributed back
827 ! through the water column.
828 do i=is,ie
829 pen_sw_rem(i) = pen_sw_bnd(1,i)
830 do n=2,nsw ; pen_sw_rem(i) = pen_sw_rem(i) + pen_sw_bnd(n,i) ; enddo
831 enddo
832 do i=is,ie ; if (pen_sw_rem(i) > 0.0) sw_remains = .true. ; enddo
833
834 ih_limit = 1.0 / h_limit_fluxes
835 do i=is,ie ; if ((pen_sw_rem(i) > 0.0) .and. (h_heat(i) > 0.0)) then
836 if (h_heat(i)*ih_limit >= 1.0) then
837 t_chg(i) = pen_sw_rem(i) / h_heat(i) ; unabsorbed = 0.0
838 else
839 t_chg(i) = pen_sw_rem(i) * ih_limit
840 unabsorbed = 1.0 - h_heat(i)*ih_limit
841 endif
842 do n=1,nsw ; pen_sw_bnd(n,i) = unabsorbed * pen_sw_bnd(n,i) ; enddo
843 endif ; enddo
844 endif ! absorbAllSW
845
846 if (absorballsw .or. adjustabsorptionprofile) then
847 do ks=nz,1,-1 ; do i=is,ie
848 k = ks
849 if (present(ksort)) then
850 if (ksort(i,ks) <= 0) cycle
851 k = ksort(i,ks)
852 endif
853
854 if (t_chg(i) > 0.0) then
855 ! Only layers greater than h_min_heat thick should get heated.
856 if (h(i,k) >= 2.0*h_min_heat) then ; t(i,k) = t(i,k) + t_chg(i)
857 elseif (h(i,k) > h_min_heat) then
858 t(i,k) = t(i,k) + t_chg(i) * (2.0 - 2.0*h_min_heat/h(i,k))
859 endif
860 endif
861 ! Increase the heating for layers above.
862 t_chg(i) = t_chg(i) + t_chg_above(i,k)
863 enddo ; enddo
864 if (present(htot) .and. present(ttot)) then
865 do i=is,ie ; ttot(i) = ttot(i) + t_chg(i) * htot(i) ; enddo
866 endif
867 endif ! absorbAllSW .or. adjustAbsorptionProfile
868
869end subroutine absorbremainingsw
870
871
872!> This subroutine calculates the total shortwave heat flux integrated over
873!! bands as a function of depth. This routine is only called for computing
874!! buoyancy fluxes for use in KPP. This routine does not update the state.
875subroutine sumswoverbands(G, GV, US, h, dz, nsw, optics, j, dt, &
876 H_limit_fluxes, absorbAllSW, iPen_SW_bnd, netPen)
877 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
878 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
879 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
880 real, dimension(SZI_(G),SZK_(GV)), &
881 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
882 real, dimension(SZI_(G),SZK_(GV)), &
883 intent(in) :: dz !< Layer vertical extent [Z ~> m].
884 integer, intent(in) :: nsw !< The number of bands of penetrating shortwave
885 !! radiation, perhaps from optics_nbands(optics),
886 type(optics_type), intent(in) :: optics !< An optics structure that has values
887 !! set based on the opacities.
888 integer, intent(in) :: j !< j-index to work on.
889 real, intent(in) :: dt !< Time step [T ~> s].
890 real, intent(in) :: h_limit_fluxes !< the total depth at which the
891 !! surface fluxes start to be limited to avoid
892 !! excessive heating of a thin ocean [H ~> m or kg m-2]
893 logical, intent(in) :: absorballsw !< If true, ensure that all shortwave
894 !! radiation is absorbed in the ocean water column.
895 real, dimension(max(nsw,1),SZI_(G)), intent(in) :: ipen_sw_bnd !< The incident penetrating shortwave
896 !! in each band at the sea surface; size nsw x SZI_(G)
897 !! [C H ~> degC m or degC kg m-2].
898 real, dimension(SZI_(G),SZK_(GV)+1), &
899 intent(inout) :: netpen !< Net penetrating shortwave heat flux at each
900 !! interface, summed across all bands
901 !! [C H ~> degC m or degC kg m-2].
902 ! Local variables
903 real :: h_heat(szi_(g)) ! thickness of the water column that receives
904 ! remaining shortwave radiation [H ~> m or kg m-2].
905 real :: pen_sw_rem(szi_(g)) ! sum across all wavelength bands of the
906 ! penetrating shortwave heating that hits the bottom
907 ! and will be redistributed through the water column
908 ! [C H ~> degC m or degC kg m-2]
909
910 real, dimension(max(nsw,1),SZI_(G)) :: pen_sw_bnd ! The remaining penetrating shortwave radiation
911 ! in each band, initially iPen_SW_bnd [C H ~> degC m or degC kg m-2]
912 real :: sw_trans ! fraction of shortwave radiation not
913 ! absorbed in a layer [nondim]
914 real :: unabsorbed ! fraction of the shortwave radiation
915 ! not absorbed because the layers are too thin [nondim].
916 real :: ih_limit ! inverse of the total depth at which the
917 ! surface fluxes start to be limited [H-1 ~> m-1 or m2 kg-1]
918 real :: min_sw_heat ! A minimum remaining shortwave heating within a timestep that will be simply
919 ! absorbed in the next layer for computational efficiency, instead of
920 ! continuing to penetrate [C H ~> degC m or degC kg m-2].
921 real :: i_habs ! The inverse of the absorption length for a minimal flux [H-1 ~> m-1 or m2 kg-1]
922 real :: h_min_heat ! minimum thickness layer that should get heated [H ~> m or kg m-2]
923 real :: opt_depth ! optical depth of a layer [nondim]
924 real :: exp_od ! exp(-opt_depth) [nondim]
925 logical :: sw_remains ! If true, some column has shortwave radiation that
926 ! was not entirely absorbed.
927
928 integer :: is, ie, nz, i, k, n
929 sw_remains = .false.
930
931 i_habs = 1e3*gv%H_to_m ! optics%PenSW_absorb_Invlen
932
933 h_min_heat = 2.0*gv%Angstrom_H + gv%H_subroundoff
934 is = g%isc ; ie = g%iec ; nz = gv%ke
935
936 if (nsw < 1) then
937 netpen(:,:) = 0.0
938 return
939 endif
940
941 pen_sw_bnd(:,:) = ipen_sw_bnd(:,:)
942 do i=is,ie ; h_heat(i) = 0.0 ; enddo
943 do i=is,ie
944 netpen(i,1) = 0.
945 do n=1,max(nsw,1)
946 netpen(i,1) = netpen(i,1) + pen_sw_bnd(n,i) ! Surface interface
947 enddo
948 enddo
949
950 ! Apply penetrating SW radiation to remaining parts of layers.
951 ! Excessively thin layers are not heated to avoid runaway temps.
952 min_sw_heat = optics%PenSW_flux_absorb*dt ! Default of 2.5e-11*US%T_to_s*GV%m_to_H
953 do k=1,nz
954
955 do i=is,ie
956 netpen(i,k+1) = 0.
957
958 if (h(i,k) > 0.0) then
959 do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
960 ! SW_trans is the SW that is transmitted THROUGH the layer
961 opt_depth = dz(i,k) * optics%opacity_band(n,i,j,k)
962 exp_od = exp(-opt_depth)
963 sw_trans = exp_od
964
965 ! Heating at a very small rate can be absorbed by a sufficiently thick layer or several
966 ! thin layers without further penetration.
967 if (optics%answer_date < 20190101) then
968 if (nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat*min(1.0, i_habs*h(i,k)) ) sw_trans = 0.0
969 elseif ((nsw*pen_sw_bnd(n,i)*sw_trans < min_sw_heat) .and. (h(i,k) > h_min_heat)) then
970 if (nsw*pen_sw_bnd(n,i) <= min_sw_heat * (i_habs*(h(i,k) - h_min_heat))) then
971 sw_trans = 0.0
972 else
973 sw_trans = min(sw_trans, &
974 1.0 - (min_sw_heat*(i_habs*(h(i,k) - h_min_heat))) / (nsw*pen_sw_bnd(n,i)))
975 endif
976 endif
977
978 pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
979 netpen(i,k+1) = netpen(i,k+1) + pen_sw_bnd(n,i)
980 endif ; enddo
981 endif ! h(i,k) > 0.0
982
983 ! Add to the accumulated thickness above that could be heated.
984 ! Only layers greater than h_min_heat thick should get heated.
985 if (h(i,k) >= 2.0*h_min_heat) then
986 h_heat(i) = h_heat(i) + h(i,k)
987 elseif (h(i,k) > h_min_heat) then
988 h_heat(i) = h_heat(i) + (2.0*h(i,k) - 2.0*h_min_heat)
989 endif
990 enddo ! i loop
991 enddo ! k loop
992
993 if (absorballsw) then
994
995 ! If there is still shortwave radiation at this point, it could go into
996 ! the bottom (with a bottom mud model), or it could be redistributed back
997 ! through the water column.
998 do i=is,ie
999 pen_sw_rem(i) = pen_sw_bnd(1,i)
1000 do n=2,nsw ; pen_sw_rem(i) = pen_sw_rem(i) + pen_sw_bnd(n,i) ; enddo
1001 enddo
1002 do i=is,ie ; if (pen_sw_rem(i) > 0.0) sw_remains = .true. ; enddo
1003
1004 ih_limit = 1.0 / h_limit_fluxes
1005 do i=is,ie ; if ((pen_sw_rem(i) > 0.0) .and. (h_heat(i) > 0.0)) then
1006 if (h_heat(i)*ih_limit < 1.0) then
1007 unabsorbed = 1.0 - h_heat(i)*ih_limit
1008 else
1009 unabsorbed = 0.0
1010 endif
1011 do n=1,nsw ; pen_sw_bnd(n,i) = unabsorbed * pen_sw_bnd(n,i) ; enddo
1012 endif ; enddo
1013
1014 endif ! absorbAllSW
1015
1016end subroutine sumswoverbands
1017
1018
1019
1020!> This routine initializes the opacity module, including an optics_type.
1021subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
1022 type(time_type), target, intent(in) :: time !< The current model time.
1023 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
1024 type(verticalgrid_type), intent(in) :: gv !< model vertical grid structure
1025 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1026 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1027 !! parameters.
1028 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
1029 !! output.
1030 type(opacity_cs) :: cs !< Opacity control structure
1031 type(optics_type) :: optics !< An optics structure that has parameters
1032 !! set and arrays allocated here.
1033 ! Local variables
1034 character(len=200) :: tmpstr
1035 character(len=40) :: mdl = "MOM_opacity"
1036 character(len=40) :: bandnum, shortname
1037 character(len=200) :: longname
1038 character(len=40) :: scheme_string
1039 ! This include declares and sets the variable "version".
1040# include "version_variable.h"
1041 real :: opacity_coefs(6) ! Pairs of opacity coefficients [Z-1 ~> m-1] for blue, red and
1042 ! near-infrared radiation with parameterizations following the
1043 ! functional form from Manizza et al., GRL 2005, namely in the form
1044 ! opacity = coef_1 + coef_2 * chl**pow for each band.
1045 real :: opacity_powers(3) ! Powers of chlorophyll [nondim] for blue, red and near-infrared
1046 ! radiation bands, in expressions for opacity of the form
1047 ! opacity = coef_1 + coef_2 * chl**pow.
1048 real :: extinction_coefs(6) ! Extinction length coefficients [Z ~> m] for penetrating shortwave
1049 ! radiation in the form proposed by Morel and Antoine (1994), namely
1050 ! opacity = 1 / (sum(n=1:6, Coef(n) * log10(Chl)**(n-1)))
1051 real :: sw_pen_frac_coefs(6) ! Coefficients for the shortwave radiation fraction [nondim] in a
1052 ! fifth order polynomial fit as a funciton of log10(Chlorophyll).
1053 real :: pensw_absorb_minthick ! A thickness that is used to absorb the remaining shortwave heat
1054 ! flux when that flux drops below PEN_SW_FLUX_ABSORB [H ~> m or kg m-2]
1055 real :: pensw_minthick_dflt ! The default for PenSW_absorb_minthick [m]
1056 real :: i_nir_bands ! The inverse of the number of near-infrared bands being used [nondim]
1057 real, allocatable :: band_wavelengths(:) ! The bounding wavelengths for the penetrating shortwave
1058 ! radiation bands [nm]
1059 real, allocatable :: band_wavelen_default(:) ! The defaults for band_wavelengths [nm]
1060 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags
1061 integer :: isd, ied, jsd, jed, nz, n
1062 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1063
1064 cs%diag => diag
1065
1066 ! Read all relevant parameters and write them to the model log.
1067 call log_version(param_file, mdl, version, '')
1068
1069! parameters for CHL_A routines
1070 call get_param(param_file, mdl, "VAR_PEN_SW", cs%var_pen_sw, &
1071 "If true, use one of the CHL_A schemes specified by "//&
1072 "OPACITY_SCHEME to determine the e-folding depth of "//&
1073 "incoming short wave radiation.", default=.false.)
1074
1075 cs%opacity_scheme = no_scheme ; scheme_string = ''
1076 if (cs%var_pen_sw) then
1077 call get_param(param_file, mdl, "OPACITY_SCHEME", tmpstr, &
1078 "This character string specifies how chlorophyll "//&
1079 "concentrations are translated into opacities. Currently "//&
1080 "valid options include:\n"//&
1081 " \t\t MANIZZA_05 - Use Manizza et al., GRL, 2005. \n"//&
1082 " \t\t MOREL_88 - Use Morel, JGR, 1988. \n"//&
1083 " \t\t OHLMANN_03 - Use Ohlmann, J Clim, 2003. (only use if dz(1)>2.0m)", &
1084 default=manizza_05_string)
1085 if (len_trim(tmpstr) > 0) then
1086 tmpstr = uppercase(tmpstr)
1087 select case (tmpstr)
1088 case (manizza_05_string)
1089 cs%opacity_scheme = manizza_05 ; scheme_string = manizza_05_string
1090 case (morel_88_string)
1091 cs%opacity_scheme = morel_88 ; scheme_string = morel_88_string
1092 case (ohlmann_03_string)
1093 cs%opacity_scheme = ohlmann_03 ; scheme_string = ohlmann_03_string
1094 case default
1095 call mom_error(fatal, "opacity_init: #DEFINE OPACITY_SCHEME "//&
1096 trim(tmpstr) // "in input file is invalid.")
1097 end select
1098 call mom_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5)
1099 endif
1100 if (cs%opacity_scheme == no_scheme) then
1101 call mom_error(warning, "opacity_init: No scheme has successfully "//&
1102 "been specified for the opacity. Using the default MANIZZA_05.")
1103 cs%opacity_scheme = manizza_05 ; scheme_string = manizza_05_string
1104 endif
1105
1106 call get_param(param_file, mdl, "BLUE_FRAC_SW", cs%blue_frac, &
1107 "The fraction of the penetrating shortwave radiation "//&
1108 "that is in the blue band.", default=0.5, units="nondim")
1109 else
1110 call get_param(param_file, mdl, "EXP_OPACITY_SCHEME", tmpstr, &
1111 "This character string specifies which exponential "//&
1112 "opacity scheme to utilize. Currently "//&
1113 "valid options include:\n"//&
1114 " \t\t SINGLE_EXP - Single Exponent decay. \n"//&
1115 " \t\t DOUBLE_EXP - Double Exponent decay.", &
1116 default=single_exp_string)!New default for "else" above (non-Chl scheme)
1117 if (len_trim(tmpstr) > 0) then
1118 tmpstr = uppercase(tmpstr)
1119 select case (tmpstr)
1120 case (single_exp_string)
1121 cs%opacity_scheme = single_exp ; scheme_string = single_exp_string
1122 case (double_exp_string)
1123 cs%opacity_scheme = double_exp ; scheme_string = double_exp_string
1124 end select
1125 call mom_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5)
1126 endif
1127 call get_param(param_file, mdl, "PEN_SW_SCALE", cs%pen_sw_scale, &
1128 "The vertical absorption e-folding depth of the penetrating shortwave radiation.", &
1129 units="m", default=0.0, scale=us%m_to_Z)
1130 !BGR/ Added for opacity_scheme==double_exp read in 2nd exp-decay and fraction
1131 if (cs%Opacity_scheme == double_exp ) then
1132 call get_param(param_file, mdl, "PEN_SW_SCALE_2ND", cs%pen_sw_scale_2nd, &
1133 "The (2nd) vertical absorption e-folding depth of the "//&
1134 "penetrating shortwave radiation (use if SW_EXP_MODE==double.)", &
1135 units="m", default=0.0, scale=us%m_to_Z)
1136 call get_param(param_file, mdl, "SW_1ST_EXP_RATIO", cs%sw_1st_exp_ratio, &
1137 "The fraction of 1st vertical absorption e-folding depth "//&
1138 "penetrating shortwave radiation if SW_EXP_MODE==double.",&
1139 units="nondim", default=0.0)
1140 elseif (cs%OPACITY_SCHEME == single_exp) then
1141 !/Else disable 2nd_exp scheme
1142 cs%pen_sw_scale_2nd = 0.0
1143 cs%sw_1st_exp_ratio = 1.0
1144 endif
1145 call get_param(param_file, mdl, "PEN_SW_FRAC", cs%pen_sw_frac, &
1146 "The fraction of the shortwave radiation that penetrates "//&
1147 "below the surface.", units="nondim", default=0.0)
1148
1149 endif
1150 call get_param(param_file, mdl, "PEN_SW_NBANDS", optics%nbands, &
1151 "The number of bands of penetrating shortwave radiation.", &
1152 default=1)
1153 if (cs%Opacity_scheme == double_exp ) then
1154 if (optics%nbands /= 2) call mom_error(fatal, &
1155 "set_opacity: \Cannot use a double_exp opacity scheme with nbands!=2.")
1156 elseif (cs%Opacity_scheme == single_exp ) then
1157 if (optics%nbands /= 1) call mom_error(fatal, &
1158 "set_opacity: \Cannot use a single_exp opacity scheme with nbands!=1.")
1159 elseif (cs%Opacity_scheme == ohlmann_03 ) then
1160 if (optics%nbands /= 2) call mom_error(fatal, &
1161 "set_opacity: \OHLMANN_03 scheme requires nbands==2")
1162 endif
1163
1164 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
1165 "This sets the default value for the various _ANSWER_DATE parameters.", &
1166 default=99991231)
1167 call get_param(param_file, mdl, "OPTICS_ANSWER_DATE", optics%answer_date, &
1168 "The vintage of the order of arithmetic and expressions in the optics calculations. "//&
1169 "Values below 20190101 recover the answers from the end of 2018, while "//&
1170 "higher values use updated and more robust forms of the same expressions.", &
1171 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
1172 if (.not.gv%Boussinesq) optics%answer_date = max(optics%answer_date, 20230701)
1173
1174 call get_param(param_file, mdl, "PEN_SW_FLUX_ABSORB", optics%PenSW_flux_absorb, &
1175 "A minimum remaining shortwave heating rate that will be simply absorbed in "//&
1176 "the next sufficiently thick layers for computational efficiency, instead of "//&
1177 "continuing to penetrate. The default, 2.5e-11 degC m s-1, is about 1e-4 W m-2 "//&
1178 "or 0.08 degC m century-1, but 0 is also a valid value.", &
1179 default=2.5e-11, units="degC m s-1", scale=us%degC_to_C*gv%m_to_H*us%T_to_s)
1180
1181 if (optics%answer_date < 20190101) then ; pensw_minthick_dflt = 0.001 ; else ; pensw_minthick_dflt = 1.0 ; endif
1182 call get_param(param_file, mdl, "PEN_SW_ABSORB_MINTHICK", pensw_absorb_minthick, &
1183 "A thickness that is used to absorb the remaining penetrating shortwave heat "//&
1184 "flux when it drops below PEN_SW_FLUX_ABSORB.", &
1185 default=pensw_minthick_dflt, units="m", scale=gv%m_to_H)
1186 optics%PenSW_absorb_Invlen = 1.0 / (pensw_absorb_minthick + gv%H_subroundoff)
1187
1188 ! The defaults for the following coefficients are taken from Manizza et al., GRL, 2005.
1189 call get_param(param_file, mdl, "OPACITY_VALUES_MANIZZA", opacity_coefs, &
1190 "Pairs of opacity coefficients for blue, red and near-infrared radiation with "//&
1191 "parameterizations following the functional form from Manizza et al., GRL 2005, "//&
1192 "namely in the form opacity = coef_1 + coef_2 * chl**pow for each band. Although "//&
1193 "coefficients are set for 3 bands, more or less bands may actually be used, with "//&
1194 "extra bands following the same properties as band 3.", &
1195 units="m-1", scale=us%Z_to_m, defaults=(/0.0232, 0.074, 0.225, 0.037, 2.86, 0.0/), &
1196 do_not_log=(cs%opacity_scheme/=manizza_05))
1197 call get_param(param_file, mdl, "CHOROPHYLL_POWER_MANIZZA", opacity_powers, &
1198 "Powers of chlorophyll for blue, red and near-infrared radiation bands in "//&
1199 "expressions for opacity of the form opacity = coef_1 + coef_2 * chl**pow.", &
1200 units="nondim", defaults=(/0.674, 0.629, 0.0/), &
1201 do_not_log=(cs%opacity_scheme/=manizza_05))
1202
1203 ! The defaults for the following coefficients are taken from Morel and Antoine (1994).
1204 call get_param(param_file, mdl, "OPACITY_VALUES_MOREL", extinction_coefs, &
1205 "Shortwave extinction length coefficients for shortwave radiation in the form "//&
1206 "proposed by Morel (1988), opacity = 1 / (sum(Coef(n) * log10(Chl)**(n-1))).", &
1207 units="m", scale=us%m_to_Z, defaults=(/7.925, -6.644, 3.662, -1.815, -0.218, 0.502/), &
1208 do_not_log=(cs%opacity_scheme/=morel_88))
1209 call get_param(param_file, mdl, "SW_PEN_FRAC_COEFS_MOREL", sw_pen_frac_coefs, &
1210 "Coefficients for the shortwave radiation fraction in a fifth order polynomial "//&
1211 "fit as a function of log10(Chlorophyll).", &
1212 units="nondim", defaults=(/0.321, 0.008, 0.132, 0.038, -0.017, -0.007/), &
1213 do_not_log=(cs%opacity_scheme/=morel_88))
1214
1215 if (.not.allocated(optics%min_wavelength_band)) &
1216 allocate(optics%min_wavelength_band(optics%nbands))
1217 if (.not.allocated(optics%max_wavelength_band)) &
1218 allocate(optics%max_wavelength_band(optics%nbands))
1219
1220 ! Set the wavelengths of the opacity bands
1221 allocate(band_wavelengths(optics%nbands+1), source=0.0)
1222 allocate(band_wavelen_default(optics%nbands+1), source=0.0)
1223 if (cs%opacity_scheme == manizza_05) then
1224 if (optics%nbands >= 1) band_wavelen_default(2) = 550.0
1225 if (optics%nbands >= 2) band_wavelen_default(3) = 700.0
1226 if (optics%nbands >= 3) then
1227 i_nir_bands = 1.0 / real(optics%nbands - 2)
1228 do n=3,optics%nbands
1229 band_wavelen_default(n+1) = 2800. - (optics%nbands-n)*2100.0*i_nir_bands
1230 enddo
1231 endif
1232 endif
1233 call get_param(param_file, mdl, "OPACITY_BAND_WAVELENGTHS", band_wavelengths, &
1234 "The bounding wavelengths for the various bands of shortwave radiation, with "//&
1235 "defaults that depend on the setting for OPACITY_SCHEME.", &
1236 units="nm", defaults=band_wavelen_default, do_not_log=(optics%nbands<2))
1237 do n=1,optics%nbands
1238 optics%min_wavelength_band(n) = band_wavelengths(n)
1239 optics%max_wavelength_band(n) = band_wavelengths(n+1)
1240 enddo
1241 deallocate(band_wavelengths, band_wavelen_default)
1242
1243 ! Set opacity scheme dependent parameters.
1244
1245 if (cs%opacity_scheme == manizza_05) then
1246 allocate(cs%opacity_coef(2,optics%nbands))
1247 allocate(cs%chl_power(optics%nbands))
1248 do n=1,min(3,optics%nbands)
1249 cs%opacity_coef(1,n) = opacity_coefs(2*n-1) ; cs%opacity_coef(2,n) = opacity_coefs(2*n)
1250 cs%chl_power(n) = opacity_powers(n)
1251 enddo
1252 ! All remaining bands use the same properties as NIR, for lack of something better to do.
1253 do n=4,optics%nbands
1254 cs%opacity_coef(1,n) = cs%opacity_coef(1,n-1) ; cs%opacity_coef(2,n) = cs%opacity_coef(2,n-1)
1255 cs%chl_power(n) = cs%chl_power(n-1)
1256 enddo
1257 ! Determine the last band that is dependent on chlorophyll.
1258 cs%chl_dep_bands = optics%nbands
1259 do n=optics%nbands,1,-1
1260 if (cs%chl_power(n) /= 0.0) exit
1261 cs%chl_dep_bands = n - 1
1262 enddo
1263 do n=cs%chl_dep_bands+1,optics%nbands
1264 if (cs%opacity_coef(2,n) /= 0.0) then
1265 call mom_error(warning, "set_opacity: A non-zero value of the chlorophyll dependence in "//&
1266 "OPACITY_VALUES_MANIZZA was set for a band with zero power in its chlorophyll dependence "//&
1267 "as set by CHOROPHYLL_POWER_MANIZZA.")
1268 cs%opacity_coef(1,n) = cs%opacity_coef(1,n) + cs%opacity_coef(2,n)
1269 cs%opacity_coef(2,n) = 0.0
1270 endif
1271 enddo
1272
1273 elseif (cs%opacity_scheme == morel_88) then
1274 ! The Morel opacity scheme represents a non uniform distribution of chlorophyll-a through the
1275 ! water column. Other approaches may be more appropriate when using an interactive ecosystem
1276 ! model that predicts three-dimensional chl-a values.
1277 allocate(cs%opacity_coef(6, optics%nbands))
1278 allocate(cs%sw_pen_frac_coef(6))
1279
1280 ! As presently implemented, all frequency bands use the same opacities.
1281 do n=1,optics%nbands
1282 cs%opacity_coef(1:6,n) = extinction_coefs(1:6)
1283 enddo
1284 cs%sw_pen_frac_coef(:) = sw_pen_frac_coefs(:)
1285 endif
1286
1287 call get_param(param_file, mdl, "OPACITY_LAND_VALUE", cs%opacity_land_value, &
1288 "The value to use for opacity over land. The default is "//&
1289 "10 m-1 - a value for muddy water.", units="m-1", default=10.0, scale=us%Z_to_m)
1290
1291 cs%warning_issued = .false.
1292
1293 if (.not.allocated(optics%opacity_band)) &
1294 allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz), source=0.0)
1295 if (.not.allocated(optics%sw_pen_band)) &
1296 allocate(optics%sw_pen_band(optics%nbands,isd:ied,jsd:jed))
1297 allocate(cs%id_opacity(optics%nbands), source=-1)
1298
1299 cs%id_sw_pen = register_diag_field('ocean_model', 'SW_pen', diag%axesT1, time, &
1300 'Penetrating shortwave radiation flux into ocean', 'W m-2', conversion=us%QRZ_T_to_W_m2)
1301 cs%id_sw_vis_pen = register_diag_field('ocean_model', 'SW_vis_pen', diag%axesT1, time, &
1302 'Visible penetrating shortwave radiation flux into ocean', 'W m-2', conversion=us%QRZ_T_to_W_m2)
1303 do n=1,optics%nbands
1304 write(bandnum,'(I0)') n
1305 shortname = 'opac_'//trim(bandnum)
1306 longname = 'Opacity for shortwave radiation in band '//trim(bandnum)// &
1307 ', saved as L^-1 tanh(Opacity * L) for L = 10^-10 m'
1308 cs%id_opacity(n) = register_diag_field('ocean_model', shortname, diag%axesTL, time, &
1309 longname, 'm-1', conversion=us%m_to_Z)
1310 enddo
1311
1312 if (cs%opacity_scheme == ohlmann_03) then
1313 ! Set up the lookup table
1314 call init_ohlmann_table(optics)
1315 endif
1316
1317end subroutine opacity_init
1318
1319!> Initialize the lookup table for Ohlmann solar penetration scheme.
1320!! Step size in Chl is a constant in log-space to make lookups easy.
1321!! Step size is fine enough that nearest neighbor lookup is sufficiently
1322!! accurate.
1323subroutine init_ohlmann_table(optics)
1324
1325 implicit none
1326
1327 type(optics_type), intent(inout) :: optics
1328
1329 ! Local variables
1330
1331 !! These are the data from Ohlmann (2003) Table 1a with additional
1332 !! values provided by C. Ohlmann and implemented in CESM-POP by B. Briegleb
1333 integer, parameter :: nval_tab1a = 31
1334 real, parameter, dimension(nval_tab1a) :: &
1335 chl_tab1a = (/ &
1336 .001, .005, .01, .02, &
1337 .03, .05, .10, .15, &
1338 .20, .25, .30, .35, &
1339 .40, .45, .50, .60, &
1340 .70, .80, .90, 1.00, &
1341 1.50, 2.00, 2.50, 3.00, &
1342 4.00, 5.00, 6.00, 7.00, &
1343 8.00, 9.00, 10.00 /)
1344
1345 real, parameter, dimension(nval_tab1a) :: &
1346 a1_tab1a = (/ &
1347 0.4421, 0.4451, 0.4488, 0.4563, &
1348 0.4622, 0.4715, 0.4877, 0.4993, &
1349 0.5084, 0.5159, 0.5223, 0.5278, &
1350 0.5326, 0.5369, 0.5408, 0.5474, &
1351 0.5529, 0.5576, 0.5615, 0.5649, &
1352 0.5757, 0.5802, 0.5808, 0.5788, &
1353 0.56965, 0.55638, 0.54091, 0.52442, &
1354 0.50766, 0.49110, 0.47505 /)
1355
1356 real, parameter, dimension(nval_tab1a) :: &
1357 a2_tab1a = (/ &
1358 0.2981, 0.2963, 0.2940, 0.2894, &
1359 0.2858, 0.2800, 0.2703, 0.2628, &
1360 0.2571, 0.2523, 0.2481, 0.2444, &
1361 0.2411, 0.2382, 0.2356, 0.2309, &
1362 0.2269, 0.2235, 0.2206, 0.2181, &
1363 0.2106, 0.2089, 0.2113, 0.2167, &
1364 0.23357, 0.25504, 0.27829, 0.30274, &
1365 0.32698, 0.35056, 0.37303 /)
1366
1367 real, parameter, dimension(nval_tab1a) :: &
1368 b1_tab1a = (/ &
1369 0.0287, 0.0301, 0.0319, 0.0355, &
1370 0.0384, 0.0434, 0.0532, 0.0612, &
1371 0.0681, 0.0743, 0.0800, 0.0853, &
1372 0.0902, 0.0949, 0.0993, 0.1077, &
1373 0.1154, 0.1227, 0.1294, 0.1359, &
1374 0.1640, 0.1876, 0.2082, 0.2264, &
1375 0.25808, 0.28498, 0.30844, 0.32932, &
1376 0.34817, 0.36540, 0.38132 /)
1377
1378 real, parameter, dimension(nval_tab1a) :: &
1379 b2_tab1a = (/ &
1380 0.3192, 0.3243, 0.3306, 0.3433, &
1381 0.3537, 0.3705, 0.4031, 0.4262, &
1382 0.4456, 0.4621, 0.4763, 0.4889, &
1383 0.4999, 0.5100, 0.5191, 0.5347, &
1384 0.5477, 0.5588, 0.5682, 0.5764, &
1385 0.6042, 0.6206, 0.6324, 0.6425, &
1386 0.66172, 0.68144, 0.70086, 0.72144, &
1387 0.74178, 0.76190, 0.78155 /)
1388
1389 !! Make the table big enough so step size is smaller
1390 !! in log-space that any increment in Table 1a
1391 integer, parameter :: nval_lut=401
1392 real :: chl, log10chl_lut, w1, w2
1393 integer :: n,m,mm1,err
1394
1395 allocate(optics%a1_lut(nval_lut),optics%b1_lut(nval_lut),&
1396 & optics%a2_lut(nval_lut),optics%b2_lut(nval_lut),&
1397 & stat=err)
1398 if ( err /= 0 ) then
1399 call mom_error(fatal,"init_ohlmann: Cannot allocate lookup table")
1400 endif
1401
1402 optics%chl_min = chl_tab1a(1)
1403 optics%log10chl_min = log10(chl_tab1a(1))
1404 optics%log10chl_max = log10(chl_tab1a(nval_tab1a))
1405 optics%dlog10chl = (optics%log10chl_max - optics%log10chl_min)/(nval_lut-1)
1406
1407 ! step through the lookup table
1408 m = 2
1409 do n=1,nval_lut
1410 log10chl_lut = optics%log10chl_min + (n-1)*optics%dlog10chl
1411 chl = 10.0**log10chl_lut
1412 chl = max(chl_tab1a(1),min(chl,chl_tab1a(nval_tab1a)))
1413
1414 ! find interval in Table 1a (m-1,m]
1415 do while (chl > chl_tab1a(m))
1416 m = m + 1
1417 enddo
1418 mm1 = m-1
1419
1420 ! interpolation weights
1421 w2 = (chl - chl_tab1a(mm1))/(chl_tab1a(m) - chl_tab1a(mm1))
1422 w1 = 1. - w2
1423
1424 ! fill in the tables
1425 optics%a1_lut(n) = w1*a1_tab1a(mm1) + w2*a1_tab1a(m)
1426 optics%a2_lut(n) = w1*a2_tab1a(mm1) + w2*a2_tab1a(m)
1427 optics%b1_lut(n) = w1*b1_tab1a(mm1) + w2*b1_tab1a(m)
1428 optics%b2_lut(n) = w1*b2_tab1a(mm1) + w2*b2_tab1a(m)
1429 enddo
1430
1431 return
1432end subroutine init_ohlmann_table
1433
1434!> Get the partion of total solar into bands from Ohlmann lookup table
1435function lookup_ohlmann_swpen(chl,optics) result(A)
1436
1437 implicit none
1438
1439 real, intent(in) :: chl
1440 type(optics_type), intent(in) :: optics
1441 real, dimension(2) :: a
1442
1443 ! Local variables
1444
1445 real :: log10chl
1446 integer :: n
1447
1448 ! Make sure we are in the table
1449 if (chl > optics%chl_min) then
1450 log10chl = min(log10(chl),optics%log10chl_max)
1451 else
1452 log10chl = optics%log10chl_min
1453 endif
1454 ! Do a nearest neighbor lookup
1455 n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1
1456
1457 a(1) = optics%a1_lut(n)
1458 a(2) = optics%a2_lut(n)
1459
1460end function lookup_ohlmann_swpen
1461
1462!> Get the opacity (decay scale) from Ohlmann lookup table
1463function lookup_ohlmann_opacity(chl,optics) result(B)
1464
1465 implicit none
1466 real, intent(in) :: chl
1467 type(optics_type), intent(in) :: optics
1468 real, dimension(2) :: b
1469
1470 ! Local variables
1471 real :: log10chl
1472 integer :: n
1473
1474 ! Make sure we are in the table
1475 if (chl > optics%chl_min) then
1476 log10chl = min(log10(chl),optics%log10chl_max)
1477 else
1478 log10chl = optics%log10chl_min
1479 endif
1480 ! Do a nearest neighbor lookup
1481 n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1
1482
1483 b(1) = optics%b1_lut(n)
1484 b(2) = optics%b2_lut(n)
1485
1486 return
1487end function lookup_ohlmann_opacity
1488
1489subroutine opacity_end(CS, optics)
1490 type(opacity_cs) :: cs !< Opacity control structure
1491 type(optics_type) :: optics !< An optics type structure that should be deallocated.
1492
1493 if (allocated(cs%id_opacity)) &
1494 deallocate(cs%id_opacity)
1495 if (allocated(cs%opacity_coef)) &
1496 deallocate(cs%opacity_coef)
1497 if (allocated(cs%sw_pen_frac_coef)) &
1498 deallocate(cs%sw_pen_frac_coef)
1499 if (allocated(cs%chl_power)) &
1500 deallocate(cs%chl_power)
1501 if (allocated(optics%sw_pen_band)) &
1502 deallocate(optics%sw_pen_band)
1503 if (allocated(optics%opacity_band)) &
1504 deallocate(optics%opacity_band)
1505 if (allocated(optics%max_wavelength_band)) &
1506 deallocate(optics%max_wavelength_band)
1507 if (allocated(optics%min_wavelength_band)) &
1508 deallocate(optics%min_wavelength_band)
1509 if (allocated(optics%a1_lut)) deallocate(optics%a1_lut)
1510 if (allocated(optics%a2_lut)) deallocate(optics%a2_lut)
1511 if (allocated(optics%b1_lut)) deallocate(optics%b1_lut)
1512 if (allocated(optics%b2_lut)) deallocate(optics%b2_lut)
1513end subroutine opacity_end
1514
1515!> \namespace mom_opacity
1516!!
1517!! opacity_from_chl:
1518!! In this routine, the Morel (modified) or Manizza (modified)
1519!! schemes use the "blue" band in the parameterizations to determine
1520!! the e-folding depth of the incoming shortwave attenuation. The red
1521!! portion is lumped into the net heating at the surface.
1522!!
1523!! Morel, A., 1988: Optical modeling of the upper ocean in relation
1524!! to its biogenous matter content (case-i waters)., J. Geo. Res.,
1525!! 93, 10,749-10,768.
1526!!
1527!! Manizza, M., C. LeQuere, A. J. Watson, and E. T. Buitenhuis, 2005:
1528!! Bio-optical feedbacks among phytoplankton, upper ocean physics
1529!! and sea-ice in a global model, Geophys. Res. Let., 32, L05603,
1530!! doi:10.1029/2004GL020778.
1531
1532!! Ohlmann, J.C., 2003: Ocean radiant heating in climate models.
1533!! J. Climate, 16, 1337-1351, 2003.
1534
1535end module mom_opacity