portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Routines used to calculate the opacity of the ocean. | |
| 6 | module MOM_opacity | |
| 7 | ||
| 8 | use MOM_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data | |
| 9 | use MOM_diag_mediator, only : query_averaging_enabled, register_diag_field | |
| 10 | use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING | |
| 11 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 12 | use MOM_string_functions, only : uppercase | |
| 13 | use MOM_grid, only : ocean_grid_type | |
| 14 | use MOM_unit_scaling, only : unit_scale_type | |
| 15 | use MOM_variables, only : thermo_var_ptrs | |
| 16 | use MOM_verticalGrid, only : verticalGrid_type | |
| 17 | ||
| 18 | implicit none ; private | |
| 19 | ||
| 20 | #include <MOM_memory.h> | |
| 21 | ||
| 22 | public set_opacity, opacity_init, opacity_end | |
| 23 | public extract_optics_slice, extract_optics_fields, optics_nbands | |
| 24 | public absorbRemainingSW, sumSWoverBands | |
| 25 | ||
| 26 | !> This type is used to store information about ocean optical properties | |
| 27 | type, 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 | ||
| 63 | end type optics_type | |
| 64 | ||
| 65 | !> The control structure with parameters for the MOM_opacity module | |
| 66 | type, 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 | !>@} | |
| 104 | end type opacity_CS | |
| 105 | ||
| 106 | !>@{ Coded integers to specify the opacity scheme | |
| 107 | integer, parameter :: NO_SCHEME = 0, MANIZZA_05 = 1, MOREL_88 = 2, SINGLE_EXP = 3, DOUBLE_EXP = 4,& | |
| 108 | & OHLMANN_03 = 5 | |
| 109 | !>@} | |
| 110 | ||
| 111 | character*(10), parameter :: MANIZZA_05_STRING = "MANIZZA_05" !< String to specify the opacity scheme | |
| 112 | character*(10), parameter :: MOREL_88_STRING = "MOREL_88" !< String to specify the opacity scheme | |
| 113 | character*(10), parameter :: OHLMANN_03_STRING = "OHLMANN_03" !< String to specify the opacity scheme | |
| 114 | character*(10), parameter :: SINGLE_EXP_STRING = "SINGLE_EXP" !< String to specify the opacity scheme | |
| 115 | character*(10), parameter :: DOUBLE_EXP_STRING = "DOUBLE_EXP" !< String to specify the opacity scheme | |
| 116 | ||
| 117 | contains | |
| 118 | ||
| 119 | !> This sets the opacity of sea water based based on one of several different schemes. | |
| 120 | 12 | subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, & |
| 121 | 0 | 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 | 24 | real :: tmp(SZI_(G),SZJ_(G),SZK_(GV)) ! A 3-d temporary array for diagnosing opacity [Z-1 ~> m-1] |
| 144 | 24 | 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 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 149 | ||
| 150 | 12 | 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 | 0 | G, GV, US, CS, chl_2d, chl_3d) |
| 154 | else ! Use sw e-folding scale set by MOM_input | |
| 155 | 12 | if (optics%nbands <= 1) then ; Inv_nbands = 1.0 |
| 156 | 0 | 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 | 12 | GV%dZ_subroundoff) |
| 161 | 12 | if ( CS%Opacity_scheme == DOUBLE_EXP ) then |
| 162 | !$OMP parallel do default(shared) | |
| 163 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 164 | 0 | 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 | 0.1*GV%Angstrom_Z, GV%dZ_subroundoff) |
| 167 | enddo ; enddo ; enddo | |
| 168 | 0 | if (.not.associated(sw_total) .or. (CS%pen_SW_scale <= 0.0)) then |
| 169 | !$OMP parallel do default(shared) | |
| 170 | 0 | do j=js,je ; do i=is,ie ; do n=1,optics%nbands |
| 171 | 0 | optics%sw_pen_band(n,i,j) = 0.0 |
| 172 | enddo ; enddo ; enddo | |
| 173 | else | |
| 174 | !$OMP parallel do default(shared) | |
| 175 | 0 | do j=js,je ; do i=is,ie |
| 176 | 0 | optics%sw_pen_band(1,i,j) = (CS%SW_1st_EXP_RATIO) * sw_total(i,j) |
| 177 | 0 | 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 | 13014912 | do k=1,nz ; do j=js,je ; do i=is,ie ; do n=1,optics%nbands |
| 182 | 12960000 | optics%opacity_band(n,i,j,k) = inv_sw_pen_scale |
| 183 | enddo ; enddo ; enddo ; enddo | |
| 184 | 12 | if (.not.associated(sw_total) .or. (CS%pen_SW_scale <= 0.0)) then |
| 185 | !$OMP parallel do default(shared) | |
| 186 | 173532 | do j=js,je ; do i=is,ie ; do n=1,optics%nbands |
| 187 | 172800 | optics%sw_pen_band(n,i,j) = 0.0 |
| 188 | enddo ; enddo ; enddo | |
| 189 | else | |
| 190 | !$OMP parallel do default(shared) | |
| 191 | 0 | do j=js,je ; do i=is,ie ; do n=1,optics%nbands |
| 192 | 0 | 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 | 12 | if (query_averaging_enabled(CS%diag)) then |
| 199 | 12 | if (CS%id_sw_pen > 0) then |
| 200 | !$OMP parallel do default(shared) | |
| 201 | 0 | do j=js,je ; do i=is,ie |
| 202 | 0 | Pen_SW_tot(i,j) = 0.0 |
| 203 | 0 | do n=1,optics%nbands |
| 204 | 0 | Pen_SW_tot(i,j) = Pen_SW_tot(i,j) + optics%sw_pen_band(n,i,j) |
| 205 | enddo | |
| 206 | enddo ; enddo | |
| 207 | 0 | call post_data(CS%id_sw_pen, Pen_SW_tot, CS%diag) |
| 208 | endif | |
| 209 | 12 | if (CS%id_sw_vis_pen > 0) then |
| 210 | 0 | if (CS%opacity_scheme == MANIZZA_05) then |
| 211 | !$OMP parallel do default(shared) | |
| 212 | 0 | do j=js,je ; do i=is,ie |
| 213 | 0 | Pen_SW_tot(i,j) = 0.0 |
| 214 | 0 | do n=1,min(optics%nbands,2) |
| 215 | 0 | 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 | 0 | do j=js,je ; do i=is,ie |
| 221 | 0 | Pen_SW_tot(i,j) = 0.0 |
| 222 | 0 | do n=1,optics%nbands |
| 223 | 0 | 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 | 0 | call post_data(CS%id_sw_vis_pen, Pen_SW_tot, CS%diag) |
| 228 | endif | |
| 229 | 24 | do n=1,optics%nbands ; if (CS%id_opacity(n) > 0) then |
| 230 | 0 | 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 | 0 | 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 | 0 | tmp(i,j,k) = tanh(op_diag_len * optics%opacity_band(n,i,j,k)) / op_diag_len |
| 237 | enddo ; enddo ; enddo | |
| 238 | 0 | call post_data(CS%id_opacity(n), tmp, CS%diag) |
| 239 | endif ; enddo | |
| 240 | endif | |
| 241 | ||
| 242 | 12 | end 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. | |
| 247 | 0 | subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, & |
| 248 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | nbands = optics%nbands |
| 305 | ||
| 306 | 0 | if (nbands <= 1) then ; Inv_nbands = 1.0 |
| 307 | 0 | else ; Inv_nbands = 1.0 / real(nbands) ; endif |
| 308 | ||
| 309 | 0 | if (nbands <= 2) then ; Inv_nbands_nir = 0.0 |
| 310 | 0 | else ; Inv_nbands_nir = 1.0 / real(nbands - 2.0) ; endif |
| 311 | ||
| 312 | 0 | 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 | 0 | if (.not.CS%warning_issued) then |
| 315 | call MOM_error(WARNING, & | |
| 316 | "opacity_from_chl called without any shortwave flux arrays allocated.\n"//& | |
| 317 | 0 | "Consider setting PEN_SW_NBANDS = 0 if no shortwave fluxes are being used.") |
| 318 | endif | |
| 319 | 0 | CS%warning_issued = .true. |
| 320 | endif | |
| 321 | ||
| 322 | 0 | multiband_vis_input = (associated(sw_vis_dir) .and. associated(sw_vis_dif)) |
| 323 | 0 | multiband_nir_input = (associated(sw_nir_dir) .and. associated(sw_nir_dif)) |
| 324 | 0 | total_sw_input = associated(sw_total) |
| 325 | ||
| 326 | 0 | chl_data(:,:) = 0.0 |
| 327 | 0 | if (present(chl_3d)) then |
| 328 | 0 | do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,1) ; enddo ; enddo |
| 329 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 330 | 0 | 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 | 0 | chl_3d(i,j,k), i, j, k, G%geoLonT(i,j), G%geoLatT(i,j) |
| 334 | 0 | call MOM_error(FATAL, "MOM_opacity opacity_from_chl: "//trim(mesg)) |
| 335 | endif | |
| 336 | enddo ; enddo ; enddo | |
| 337 | 0 | elseif (present(chl_2d)) then |
| 338 | 0 | do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_2d(i,j) ; enddo ; enddo |
| 339 | 0 | do j=js,je ; do i=is,ie |
| 340 | 0 | 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 | 0 | chl_data(i,j), i, j, G%geoLonT(i,j), G%geoLatT(i,j) |
| 344 | 0 | call MOM_error(FATAL, "MOM_opacity opacity_from_chl: "//trim(mesg)) |
| 345 | endif | |
| 346 | enddo ; enddo | |
| 347 | else | |
| 348 | 0 | call MOM_error(FATAL, "Either chl_2d or chl_3d must be present in a call to opacity_form_chl.") |
| 349 | endif | |
| 350 | ||
| 351 | 0 | select case (CS%opacity_scheme) |
| 352 | case (MANIZZA_05) | |
| 353 | !$OMP parallel do default(shared) private(SW_vis_tot,SW_nir_tot) | |
| 354 | 0 | do j=js,je ; do i=is,ie |
| 355 | 0 | SW_vis_tot = 0.0 ; SW_nir_tot = 0.0 |
| 356 | 0 | if (G%mask2dT(i,j) > 0.0) then |
| 357 | 0 | if (multiband_vis_input) then |
| 358 | 0 | SW_vis_tot = sw_vis_dir(i,j) + sw_vis_dif(i,j) |
| 359 | 0 | elseif (total_sw_input) then |
| 360 | ! Follow Manizza 05 in assuming that 42% of SW is visible. | |
| 361 | 0 | SW_vis_tot = 0.42 * sw_total(i,j) |
| 362 | endif | |
| 363 | 0 | if (multiband_nir_input) then |
| 364 | 0 | SW_nir_tot = sw_nir_dir(i,j) + sw_nir_dif(i,j) |
| 365 | 0 | elseif (total_sw_input) then |
| 366 | 0 | SW_nir_tot = sw_total(i,j) - SW_vis_tot |
| 367 | endif | |
| 368 | endif | |
| 369 | ||
| 370 | ! Band 1 is Manizza blue. | |
| 371 | 0 | optics%sw_pen_band(1,i,j) = CS%blue_frac*sw_vis_tot |
| 372 | ! Band 2 (if used) is Manizza red. | |
| 373 | 0 | if (nbands > 1) & |
| 374 | 0 | 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 | 0 | do n=3,nbands |
| 377 | 0 | 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 | 0 | do j=js,je ; do i=is,ie |
| 383 | 0 | SW_pen_tot = 0.0 |
| 384 | 0 | if (G%mask2dT(i,j) > 0.0) then |
| 385 | 0 | if (multiband_vis_input) then |
| 386 | 0 | SW_pen_tot = SW_pen_frac_morel(chl_data(i,j), CS) * (sw_vis_dir(i,j) + sw_vis_dif(i,j)) |
| 387 | 0 | elseif (total_sw_input) then |
| 388 | 0 | SW_pen_tot = SW_pen_frac_morel(chl_data(i,j), CS) * 0.5*sw_total(i,j) |
| 389 | endif | |
| 390 | endif | |
| 391 | ||
| 392 | 0 | do n=1,nbands |
| 393 | 0 | 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 | 0 | if ( nbands /= 2 ) then |
| 399 | 0 | 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 | 0 | do j=js,je ; do i=is,ie |
| 403 | 0 | SW_vis_tot = 0.0 ! Ohlmann does not classify as vis/nir. Using vis to add up total |
| 404 | 0 | if (G%mask2dT(i,j) < 0.5) then |
| 405 | 0 | optics%sw_pen_band(1:2,i,j) = 0. ! Make sure there is a valid value for land points |
| 406 | else | |
| 407 | 0 | if (multiband_vis_input ) then ! If multiband_vis_input is true then so is multiband_nir_input |
| 408 | 0 | SW_vis_tot = ((sw_vis_dir(i,j) + sw_vis_dif(i,j)) + sw_nir_dir(i,j)) + sw_nir_dif(i,j) |
| 409 | 0 | elseif (total_sw_input) then |
| 410 | 0 | SW_vis_tot = sw_total(i,j) |
| 411 | else | |
| 412 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | if (present(chl_3d)) then |
| 433 | 0 | do j=js,je ; do i=is,ie ; chl_data(i,j) = chl_3d(i,j,k) ; enddo ; enddo |
| 434 | endif | |
| 435 | ||
| 436 | 0 | select case (CS%opacity_scheme) |
| 437 | case (MANIZZA_05) | |
| 438 | 0 | do j=js,je ; do i=is,ie |
| 439 | 0 | if (G%mask2dT(i,j) <= 0.5) then |
| 440 | 0 | do n=1,optics%nbands |
| 441 | 0 | optics%opacity_band(n,i,j,k) = CS%opacity_land_value |
| 442 | enddo | |
| 443 | else | |
| 444 | 0 | do n=1,CS%chl_dep_bands |
| 445 | optics%opacity_band(n,i,j,k) = CS%opacity_coef(1,n) + & | |
| 446 | 0 | CS%opacity_coef(2,n) * chl_data(i,j)**CS%chl_power(n) |
| 447 | enddo | |
| 448 | 0 | 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 | 0 | optics%opacity_band(n,i,j,k) = CS%opacity_coef(1,n) |
| 451 | enddo | |
| 452 | endif | |
| 453 | enddo ; enddo | |
| 454 | case (MOREL_88) | |
| 455 | 0 | do j=js,je ; do i=is,ie |
| 456 | 0 | optics%opacity_band(1,i,j,k) = CS%opacity_land_value |
| 457 | 0 | if (G%mask2dT(i,j) > 0.0) & |
| 458 | 0 | optics%opacity_band(1,i,j,k) = opacity_morel(chl_data(i,j), CS) |
| 459 | ||
| 460 | 0 | do n=2,optics%nbands |
| 461 | 0 | 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 | 0 | do j=js,je ; do i=is,ie |
| 467 | 0 | if (G%mask2dT(i,j) <= 0.5) then |
| 468 | 0 | 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 | 0 | 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 | 0 | call MOM_error(FATAL, "opacity_from_chl: CS%opacity_scheme is not valid.") |
| 476 | end select | |
| 477 | enddo | |
| 478 | ||
| 479 | 0 | end subroutine opacity_from_chl |
| 480 | ||
| 481 | !> This sets the blue-wavelength opacity according to the scheme proposed by | |
| 482 | !! Morel and Antoine (1994). | |
| 483 | 0 | function 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 | 0 | 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 | 0 | Chl2*(CS%opacity_coef(5,1) + Chl*CS%opacity_coef(6,1))) ) |
| 495 | 0 | end function opacity_morel |
| 496 | ||
| 497 | !> This sets the penetrating shortwave fraction according to the scheme proposed by | |
| 498 | !! Morel and Antoine (1994). | |
| 499 | 0 | function 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 | 0 | 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 | 0 | Chl2*(CS%SW_pen_frac_coef(5) + Chl*CS%SW_pen_frac_coef(6))) ) |
| 515 | 0 | end 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. | |
| 519 | 2160 | subroutine 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 | 2160 | is = G%isc ; ie = G%iec ; nz = GV%ke |
| 550 | ||
| 551 | 2160 | scale_opacity = 1.0 ; if (present(opacity_scale)) scale_opacity = opacity_scale |
| 552 | 2160 | scale_penSW = 1.0 ; if (present(penSW_scale)) scale_penSW = penSW_scale |
| 553 | ||
| 554 | 2160 | if (present(opacity)) then |
| 555 | 720 | if (present(SpV_avg)) then |
| 556 | 0 | do k=1,nz ; do i=is,ie ; do n=1,optics%nbands |
| 557 | 0 | 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 | 13014720 | do k=1,nz ; do i=is,ie ; do n=1,optics%nbands |
| 561 | 12960000 | opacity(n,i,k) = scale_opacity * optics%opacity_band(n,i,j,k) |
| 562 | enddo ; enddo ; enddo | |
| 563 | endif | |
| 564 | endif | |
| 565 | ||
| 566 | 347760 | if (present(penSW_top)) then ; do i=is,ie ; do n=1,optics%nbands |
| 567 | 345600 | penSW_top(n,i) = scale_penSW * optics%sw_pen_band(n,i,j) |
| 568 | enddo ; enddo ; endif | |
| 569 | ||
| 570 | 2880 | end subroutine extract_optics_slice |
| 571 | ||
| 572 | !> Set arguments to fields from the optics type. | |
| 573 | 0 | subroutine 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 | 0 | if (present(nbands)) nbands = optics%nbands |
| 579 | ||
| 580 | 0 | end subroutine extract_optics_fields |
| 581 | ||
| 582 | !> Return the number of bands of penetrating shortwave radiation. | |
| 583 | 732 | function 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 | 732 | if (associated(optics)) then |
| 589 | 732 | optics_nbands = optics%nbands |
| 590 | else | |
| 591 | 0 | optics_nbands = 0 |
| 592 | endif | |
| 593 | 732 | end 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. | |
| 602 | 1440 | subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_limit_fluxes, & |
| 603 | 720 | adjustAbsorptionProfile, absorbAllSW, T, Pen_SW_bnd, & |
| 604 | 1440 | 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 | 1440 | 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 | 720 | 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 | 1440 | T_chg, & ! The temperature change of thick layers due to the remaining |
| 669 | ! shortwave radiation and contributions from T_chg_above [C ~> degC]. | |
| 670 | 720 | 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 | 720 | if (nsw < 1) return |
| 704 | ||
| 705 | 720 | SW_Remains = .false. |
| 706 | 720 | min_SW_heat = optics%PenSW_flux_absorb * dt |
| 707 | 720 | I_Habs = optics%PenSW_absorb_Invlen |
| 708 | ||
| 709 | 720 | h_min_heat = 2.0*GV%Angstrom_H + GV%H_subroundoff |
| 710 | 720 | is = G%isc ; ie = G%iec ; nz = GV%ke |
| 711 | 720 | C1_6 = 1.0 / 6.0 ; C1_60 = 1.0 / 60.0 |
| 712 | ||
| 713 | 720 | TKE_calc = (present(TKE) .and. present(dSV_dT)) |
| 714 | ||
| 715 | 720 | if (optics%answer_date < 20190101) then |
| 716 | 0 | g_Hconv2 = (GV%g_Earth_Z_T2 * GV%H_to_RZ) * GV%H_to_RZ |
| 717 | else | |
| 718 | 720 | g_Hconv2 = GV%g_Earth_Z_T2 * GV%H_to_RZ**2 |
| 719 | endif | |
| 720 | ||
| 721 | 92880 | h_heat(:) = 0.0 |
| 722 | 720 | 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 | 6534720 | do ks=1,nz ; do i=is,ie |
| 727 | 6480000 | k = ks |
| 728 | 6480000 | if (present(ksort)) then |
| 729 | 0 | if (ksort(i,ks) <= 0) cycle |
| 730 | 0 | k = ksort(i,ks) |
| 731 | endif | |
| 732 | 6480000 | epsilon = 0.0 ; if (present(eps)) epsilon = eps(i,k) |
| 733 | ||
| 734 | 6480000 | T_chg_above(i,k) = 0.0 |
| 735 | ||
| 736 | 6480000 | if (h(i,k) > 1.5*epsilon) then |
| 737 | 12960000 | 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 | 0 | opt_depth = h(i,k) * opacity_band(n,i,k) |
| 740 | 0 | exp_OD = exp(-opt_depth) |
| 741 | 0 | 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 | 0 | if (optics%answer_date < 20190101) then |
| 746 | 0 | 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 | 0 | elseif ((nsw*Pen_SW_bnd(n,i)*SW_trans < min_SW_heat) .and. (h(i,k) > h_min_heat)) then |
| 748 | 0 | if (nsw*Pen_SW_bnd(n,i) <= min_SW_heat * (I_Habs*(h(i,k) - h_min_heat))) then |
| 749 | 0 | SW_trans = 0.0 |
| 750 | else | |
| 751 | SW_trans = min(SW_trans, & | |
| 752 | 0 | 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 | 0 | Heat_bnd = Pen_SW_bnd(n,i) * (1.0 - SW_trans) |
| 757 | 0 | 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 | 0 | 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 | 0 | (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 | 0 | ((h_heat(i) + h(i,k)) * (6.0 - 3.0*opt_depth)) |
| 777 | endif | |
| 778 | 0 | coSWa_frac = 0.0 |
| 779 | 0 | 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 | 0 | (SWa*(h_heat(i) + h(i,k))) |
| 782 | 0 | SWa = h_heat(i) / (h_heat(i) + h(i,k)) |
| 783 | endif | |
| 784 | ||
| 785 | 0 | T_chg_above(i,k) = T_chg_above(i,k) + (SWa * Heat_bnd) / h_heat(i) |
| 786 | 0 | T(i,k) = T(i,k) + ((1.0 - SWa) * Heat_bnd) / h(i,k) |
| 787 | else | |
| 788 | 0 | coSWa_frac = 1.0 |
| 789 | 0 | T(i,k) = T(i,k) + Pen_SW_bnd(n,i) * (1.0 - SW_trans) / h(i,k) |
| 790 | endif | |
| 791 | ||
| 792 | 0 | if (TKE_calc) then |
| 793 | 0 | 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 | 0 | (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 | 0 | (C1_6*opt_depth * (1.0 - C1_60*opt_depth**2)) |
| 803 | endif | |
| 804 | endif | |
| 805 | ||
| 806 | 0 | 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 | 6534000 | if (h(i,k) >= 2.0*h_min_heat) then |
| 813 | 6480000 | h_heat(i) = h_heat(i) + h(i,k) |
| 814 | 0 | elseif (h(i,k) > h_min_heat) then |
| 815 | 0 | 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 | 87120 | do i=is,ie ; T_chg(i) = 0.0 ; enddo |
| 823 | ||
| 824 | 720 | 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 | 87120 | do i=is,ie |
| 829 | 86400 | Pen_SW_rem(i) = Pen_SW_bnd(1,i) |
| 830 | 87120 | do n=2,nsw ; Pen_SW_rem(i) = Pen_SW_rem(i) + Pen_SW_bnd(n,i) ; enddo |
| 831 | enddo | |
| 832 | 87120 | do i=is,ie ; if (Pen_SW_rem(i) > 0.0) SW_Remains = .true. ; enddo |
| 833 | ||
| 834 | 720 | Ih_limit = 1.0 / H_limit_fluxes |
| 835 | 87120 | do i=is,ie ; if ((Pen_SW_rem(i) > 0.0) .and. (h_heat(i) > 0.0)) then |
| 836 | 0 | if (h_heat(i)*Ih_limit >= 1.0) then |
| 837 | 0 | T_chg(i) = Pen_SW_rem(i) / h_heat(i) ; unabsorbed = 0.0 |
| 838 | else | |
| 839 | 0 | T_chg(i) = Pen_SW_rem(i) * Ih_limit |
| 840 | 0 | unabsorbed = 1.0 - h_heat(i)*Ih_limit |
| 841 | endif | |
| 842 | 0 | do n=1,nsw ; Pen_SW_bnd(n,i) = unabsorbed * Pen_SW_bnd(n,i) ; enddo |
| 843 | endif ; enddo | |
| 844 | endif ! absorbAllSW | |
| 845 | ||
| 846 | 720 | if (absorbAllSW .or. adjustAbsorptionProfile) then |
| 847 | 6534720 | do ks=nz,1,-1 ; do i=is,ie |
| 848 | 6480000 | k = ks |
| 849 | 6480000 | if (present(ksort)) then |
| 850 | 0 | if (ksort(i,ks) <= 0) cycle |
| 851 | 0 | k = ksort(i,ks) |
| 852 | endif | |
| 853 | ||
| 854 | 6480000 | if (T_chg(i) > 0.0) then |
| 855 | ! Only layers greater than h_min_heat thick should get heated. | |
| 856 | 0 | if (h(i,k) >= 2.0*h_min_heat) then ; T(i,k) = T(i,k) + T_chg(i) |
| 857 | 0 | elseif (h(i,k) > h_min_heat) then |
| 858 | 0 | 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 | 6534000 | T_chg(i) = T_chg(i) + T_chg_above(i,k) |
| 863 | enddo ; enddo | |
| 864 | 720 | if (present(htot) .and. present(Ttot)) then |
| 865 | 0 | do i=is,ie ; Ttot(i) = Ttot(i) + T_chg(i) * htot(i) ; enddo |
| 866 | endif | |
| 867 | endif ! absorbAllSW .or. adjustAbsorptionProfile | |
| 868 | ||
| 869 | 2160 | end 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. | |
| 875 | 0 | subroutine sumSWoverBands(G, GV, US, h, dz, nsw, optics, j, dt, & |
| 876 | 0 | 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 | 0 | real :: h_heat(SZI_(G)) ! thickness of the water column that receives |
| 904 | ! remaining shortwave radiation [H ~> m or kg m-2]. | |
| 905 | 0 | 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 | 0 | 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 | 0 | SW_Remains = .false. |
| 930 | ||
| 931 | 0 | I_Habs = 1e3*GV%H_to_m ! optics%PenSW_absorb_Invlen |
| 932 | ||
| 933 | 0 | h_min_heat = 2.0*GV%Angstrom_H + GV%H_subroundoff |
| 934 | 0 | is = G%isc ; ie = G%iec ; nz = GV%ke |
| 935 | ||
| 936 | 0 | if (nsw < 1) then |
| 937 | 0 | netPen(:,:) = 0.0 |
| 938 | 0 | return |
| 939 | endif | |
| 940 | ||
| 941 | 0 | pen_SW_bnd(:,:) = iPen_SW_bnd(:,:) |
| 942 | 0 | do i=is,ie ; h_heat(i) = 0.0 ; enddo |
| 943 | 0 | do i=is,ie |
| 944 | 0 | netPen(i,1) = 0. |
| 945 | 0 | do n=1,max(nsw,1) |
| 946 | 0 | 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 | 0 | min_SW_heat = optics%PenSW_flux_absorb*dt ! Default of 2.5e-11*US%T_to_s*GV%m_to_H |
| 953 | 0 | do k=1,nz |
| 954 | ||
| 955 | 0 | do i=is,ie |
| 956 | 0 | netPen(i,k+1) = 0. |
| 957 | ||
| 958 | 0 | if (h(i,k) > 0.0) then |
| 959 | 0 | 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 | 0 | opt_depth = dz(i,k) * optics%opacity_band(n,i,j,k) |
| 962 | 0 | exp_OD = exp(-opt_depth) |
| 963 | 0 | 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 | 0 | if (optics%answer_date < 20190101) then |
| 968 | 0 | 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 | 0 | elseif ((nsw*Pen_SW_bnd(n,i)*SW_trans < min_SW_heat) .and. (h(i,k) > h_min_heat)) then |
| 970 | 0 | if (nsw*Pen_SW_bnd(n,i) <= min_SW_heat * (I_Habs*(h(i,k) - h_min_heat))) then |
| 971 | 0 | SW_trans = 0.0 |
| 972 | else | |
| 973 | SW_trans = min(SW_trans, & | |
| 974 | 0 | 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 | 0 | Pen_SW_bnd(n,i) = Pen_SW_bnd(n,i) * SW_trans |
| 979 | 0 | 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 | 0 | if (h(i,k) >= 2.0*h_min_heat) then |
| 986 | 0 | h_heat(i) = h_heat(i) + h(i,k) |
| 987 | 0 | elseif (h(i,k) > h_min_heat) then |
| 988 | 0 | 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 | 0 | 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 | 0 | do i=is,ie |
| 999 | 0 | Pen_SW_rem(i) = Pen_SW_bnd(1,i) |
| 1000 | 0 | do n=2,nsw ; Pen_SW_rem(i) = Pen_SW_rem(i) + Pen_SW_bnd(n,i) ; enddo |
| 1001 | enddo | |
| 1002 | 0 | do i=is,ie ; if (Pen_SW_rem(i) > 0.0) SW_Remains = .true. ; enddo |
| 1003 | ||
| 1004 | 0 | Ih_limit = 1.0 / H_limit_fluxes |
| 1005 | 0 | do i=is,ie ; if ((Pen_SW_rem(i) > 0.0) .and. (h_heat(i) > 0.0)) then |
| 1006 | 0 | if (h_heat(i)*Ih_limit < 1.0) then |
| 1007 | 0 | unabsorbed = 1.0 - h_heat(i)*Ih_limit |
| 1008 | else | |
| 1009 | 0 | unabsorbed = 0.0 |
| 1010 | endif | |
| 1011 | 0 | do n=1,nsw ; Pen_SW_bnd(n,i) = unabsorbed * Pen_SW_bnd(n,i) ; enddo |
| 1012 | endif ; enddo | |
| 1013 | ||
| 1014 | endif ! absorbAllSW | |
| 1015 | ||
| 1016 | end subroutine sumSWoverBands | |
| 1017 | ||
| 1018 | ||
| 1019 | ||
| 1020 | !> This routine initializes the opacity module, including an optics_type. | |
| 1021 | 1 | subroutine 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 | 1 | real, allocatable :: band_wavelengths(:) ! The bounding wavelengths for the penetrating shortwave |
| 1058 | ! radiation bands [nm] | |
| 1059 | 1 | 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 | 1 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke |
| 1063 | ||
| 1064 | 1 | CS%diag => diag |
| 1065 | ||
| 1066 | ! Read all relevant parameters and write them to the model log. | |
| 1067 | 1 | 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 | 1 | "incoming short wave radiation.", default=.false.) |
| 1074 | ||
| 1075 | 1 | CS%opacity_scheme = NO_SCHEME ; scheme_string = '' |
| 1076 | 1 | 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 | 0 | default=MANIZZA_05_STRING) |
| 1085 | 0 | if (len_trim(tmpstr) > 0) then |
| 1086 | 0 | tmpstr = uppercase(tmpstr) |
| 1087 | 0 | select case (tmpstr) |
| 1088 | case (MANIZZA_05_STRING) | |
| 1089 | 0 | CS%opacity_scheme = MANIZZA_05 ; scheme_string = MANIZZA_05_STRING |
| 1090 | case (MOREL_88_STRING) | |
| 1091 | 0 | CS%opacity_scheme = MOREL_88 ; scheme_string = MOREL_88_STRING |
| 1092 | case (OHLMANN_03_STRING) | |
| 1093 | 0 | CS%opacity_scheme = OHLMANN_03 ; scheme_string = OHLMANN_03_STRING |
| 1094 | case default | |
| 1095 | call MOM_error(FATAL, "opacity_init: #DEFINE OPACITY_SCHEME "//& | |
| 1096 | 0 | trim(tmpstr) // "in input file is invalid.") |
| 1097 | end select | |
| 1098 | 0 | call MOM_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5) |
| 1099 | endif | |
| 1100 | 0 | if (CS%opacity_scheme == NO_SCHEME) then |
| 1101 | call MOM_error(WARNING, "opacity_init: No scheme has successfully "//& | |
| 1102 | 0 | "been specified for the opacity. Using the default MANIZZA_05.") |
| 1103 | 0 | 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 | 0 | "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 | 1 | default=Single_Exp_String)!New default for "else" above (non-Chl scheme) |
| 1117 | 1 | if (len_trim(tmpstr) > 0) then |
| 1118 | 1 | tmpstr = uppercase(tmpstr) |
| 1119 | 1 | select case (tmpstr) |
| 1120 | case (SINGLE_EXP_STRING) | |
| 1121 | 1 | CS%opacity_scheme = SINGLE_EXP ; scheme_string = SINGLE_EXP_STRING |
| 1122 | case (DOUBLE_EXP_STRING) | |
| 1123 | 1 | CS%opacity_scheme = DOUBLE_EXP ; scheme_string = DOUBLE_EXP_STRING |
| 1124 | end select | |
| 1125 | 1 | 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 | 1 | 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 | 1 | 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 | 0 | 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 | 0 | units="nondim", default=0.0) |
| 1140 | 1 | elseif (CS%OPACITY_SCHEME == Single_Exp) then |
| 1141 | !/Else disable 2nd_exp scheme | |
| 1142 | 1 | CS%pen_sw_scale_2nd = 0.0 |
| 1143 | 1 | 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 | 1 | "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 | 1 | default=1) |
| 1153 | 1 | if (CS%Opacity_scheme == DOUBLE_EXP ) then |
| 1154 | 0 | if (optics%nbands /= 2) call MOM_error(FATAL, & |
| 1155 | 0 | "set_opacity: \Cannot use a double_exp opacity scheme with nbands!=2.") |
| 1156 | 1 | elseif (CS%Opacity_scheme == SINGLE_EXP ) then |
| 1157 | 1 | if (optics%nbands /= 1) call MOM_error(FATAL, & |
| 1158 | 0 | "set_opacity: \Cannot use a single_exp opacity scheme with nbands!=1.") |
| 1159 | 0 | elseif (CS%Opacity_scheme == OHLMANN_03 ) then |
| 1160 | 0 | if (optics%nbands /= 2) call MOM_error(FATAL, & |
| 1161 | 0 | "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 | 1 | 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 | 1 | default=default_answer_date, do_not_log=.not.GV%Boussinesq) |
| 1172 | 1 | 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 | 1 | default=2.5e-11, units="degC m s-1", scale=US%degC_to_C*GV%m_to_H*US%T_to_s) |
| 1180 | ||
| 1181 | 1 | 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 | 1 | default=PenSW_minthick_dflt, units="m", scale=GV%m_to_H) |
| 1186 | 1 | 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 | 1 | 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 | 1 | 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 | 1 | 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 | 1 | do_not_log=(CS%opacity_scheme/=MOREL_88)) |
| 1214 | ||
| 1215 | 1 | if (.not.allocated(optics%min_wavelength_band)) & |
| 1216 | 1 | allocate(optics%min_wavelength_band(optics%nbands)) |
| 1217 | 1 | if (.not.allocated(optics%max_wavelength_band)) & |
| 1218 | 1 | allocate(optics%max_wavelength_band(optics%nbands)) |
| 1219 | ||
| 1220 | ! Set the wavelengths of the opacity bands | |
| 1221 | 3 | allocate(band_wavelengths(optics%nbands+1), source=0.0) |
| 1222 | 3 | allocate(band_wavelen_default(optics%nbands+1), source=0.0) |
| 1223 | 1 | if (CS%opacity_scheme == MANIZZA_05) then |
| 1224 | 0 | if (optics%nbands >= 1) band_wavelen_default(2) = 550.0 |
| 1225 | 0 | if (optics%nbands >= 2) band_wavelen_default(3) = 700.0 |
| 1226 | 0 | if (optics%nbands >= 3) then |
| 1227 | 0 | I_NIR_bands = 1.0 / real(optics%nbands - 2) |
| 1228 | 0 | do n=3,optics%nbands |
| 1229 | 0 | 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 | 1 | units="nm", defaults=band_wavelen_default, do_not_log=(optics%nbands<2)) |
| 1237 | 2 | do n=1,optics%nbands |
| 1238 | 1 | optics%min_wavelength_band(n) = band_wavelengths(n) |
| 1239 | 2 | optics%max_wavelength_band(n) = band_wavelengths(n+1) |
| 1240 | enddo | |
| 1241 | 1 | deallocate(band_wavelengths, band_wavelen_default) |
| 1242 | ||
| 1243 | ! Set opacity scheme dependent parameters. | |
| 1244 | ||
| 1245 | 1 | if (CS%opacity_scheme == MANIZZA_05) then |
| 1246 | 0 | allocate(CS%opacity_coef(2,optics%nbands)) |
| 1247 | 0 | allocate(CS%chl_power(optics%nbands)) |
| 1248 | 0 | do n=1,min(3,optics%nbands) |
| 1249 | 0 | CS%opacity_coef(1,n) = opacity_coefs(2*n-1) ; CS%opacity_coef(2,n) = opacity_coefs(2*n) |
| 1250 | 0 | 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 | 0 | do n=4,optics%nbands |
| 1254 | 0 | CS%opacity_coef(1,n) = CS%opacity_coef(1,n-1) ; CS%opacity_coef(2,n) = CS%opacity_coef(2,n-1) |
| 1255 | 0 | CS%chl_power(n) = CS%chl_power(n-1) |
| 1256 | enddo | |
| 1257 | ! Determine the last band that is dependent on chlorophyll. | |
| 1258 | 0 | CS%chl_dep_bands = optics%nbands |
| 1259 | 0 | do n=optics%nbands,1,-1 |
| 1260 | 0 | if (CS%chl_power(n) /= 0.0) exit |
| 1261 | 0 | CS%chl_dep_bands = n - 1 |
| 1262 | enddo | |
| 1263 | 0 | do n=CS%chl_dep_bands+1,optics%nbands |
| 1264 | 0 | 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 | 0 | "as set by CHOROPHYLL_POWER_MANIZZA.") |
| 1268 | 0 | CS%opacity_coef(1,n) = CS%opacity_coef(1,n) + CS%opacity_coef(2,n) |
| 1269 | 0 | CS%opacity_coef(2,n) = 0.0 |
| 1270 | endif | |
| 1271 | enddo | |
| 1272 | ||
| 1273 | 1 | 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 | 0 | allocate(CS%opacity_coef(6, optics%nbands)) |
| 1278 | 0 | allocate(CS%sw_pen_frac_coef(6)) |
| 1279 | ||
| 1280 | ! As presently implemented, all frequency bands use the same opacities. | |
| 1281 | 0 | do n=1,optics%nbands |
| 1282 | 0 | CS%opacity_coef(1:6,n) = extinction_coefs(1:6) |
| 1283 | enddo | |
| 1284 | 0 | 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 | 1 | "10 m-1 - a value for muddy water.", units="m-1", default=10.0, scale=US%Z_to_m) |
| 1290 | ||
| 1291 | 1 | CS%warning_issued = .false. |
| 1292 | ||
| 1293 | 1 | if (.not.allocated(optics%opacity_band)) & |
| 1294 | 1310776 | allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz), source=0.0) |
| 1295 | 1 | if (.not.allocated(optics%sw_pen_band)) & |
| 1296 | 1 | allocate(optics%sw_pen_band(optics%nbands,isd:ied,jsd:jed)) |
| 1297 | 2 | 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 | 1 | '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 | 1 | 'Visible penetrating shortwave radiation flux into ocean', 'W m-2', conversion=US%QRZ_T_to_W_m2) |
| 1303 | 2 | do n=1,optics%nbands |
| 1304 | 1 | write(bandnum,'(I0)') n |
| 1305 | 1 | shortname = 'opac_'//trim(bandnum) |
| 1306 | longname = 'Opacity for shortwave radiation in band '//trim(bandnum)// & | |
| 1307 | 1 | ', 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 | 2 | longname, 'm-1', conversion=US%m_to_Z) |
| 1310 | enddo | |
| 1311 | ||
| 1312 | 1 | if (CS%opacity_scheme == OHLMANN_03) then |
| 1313 | ! Set up the lookup table | |
| 1314 | 0 | call init_ohlmann_table(optics) |
| 1315 | endif | |
| 1316 | ||
| 1317 | 1 | end 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. | |
| 1323 | 0 | subroutine 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 | 0 | & stat=err) |
| 1398 | 0 | if ( err /= 0 ) then |
| 1399 | 0 | call MOM_error(FATAL,"init_ohlmann: Cannot allocate lookup table") |
| 1400 | endif | |
| 1401 | ||
| 1402 | 0 | optics%chl_min = chl_tab1a(1) |
| 1403 | 0 | optics%log10chl_min = log10(chl_tab1a(1)) |
| 1404 | 0 | optics%log10chl_max = log10(chl_tab1a(nval_tab1a)) |
| 1405 | 0 | optics%dlog10chl = (optics%log10chl_max - optics%log10chl_min)/(nval_lut-1) |
| 1406 | ||
| 1407 | ! step through the lookup table | |
| 1408 | 0 | m = 2 |
| 1409 | 0 | do n=1,nval_lut |
| 1410 | 0 | log10chl_lut = optics%log10chl_min + (n-1)*optics%dlog10chl |
| 1411 | 0 | chl = 10.0**log10chl_lut |
| 1412 | 0 | chl = max(chl_tab1a(1),min(chl,chl_tab1a(nval_tab1a))) |
| 1413 | ||
| 1414 | ! find interval in Table 1a (m-1,m] | |
| 1415 | 0 | do while (chl > chl_tab1a(m)) |
| 1416 | 0 | m = m + 1 |
| 1417 | enddo | |
| 1418 | 0 | mm1 = m-1 |
| 1419 | ||
| 1420 | ! interpolation weights | |
| 1421 | 0 | w2 = (chl - chl_tab1a(mm1))/(chl_tab1a(m) - chl_tab1a(mm1)) |
| 1422 | 0 | w1 = 1. - w2 |
| 1423 | ||
| 1424 | ! fill in the tables | |
| 1425 | 0 | optics%a1_lut(n) = w1*a1_tab1a(mm1) + w2*a1_tab1a(m) |
| 1426 | 0 | optics%a2_lut(n) = w1*a2_tab1a(mm1) + w2*a2_tab1a(m) |
| 1427 | 0 | optics%b1_lut(n) = w1*b1_tab1a(mm1) + w2*b1_tab1a(m) |
| 1428 | 0 | optics%b2_lut(n) = w1*b2_tab1a(mm1) + w2*b2_tab1a(m) |
| 1429 | enddo | |
| 1430 | ||
| 1431 | 0 | return |
| 1432 | end subroutine init_ohlmann_table | |
| 1433 | ||
| 1434 | !> Get the partion of total solar into bands from Ohlmann lookup table | |
| 1435 | 0 | function 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 | 0 | if (chl > optics%chl_min) then |
| 1450 | 0 | log10chl = min(log10(chl),optics%log10chl_max) |
| 1451 | else | |
| 1452 | 0 | log10chl = optics%log10chl_min |
| 1453 | endif | |
| 1454 | ! Do a nearest neighbor lookup | |
| 1455 | 0 | n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1 |
| 1456 | ||
| 1457 | 0 | A(1) = optics%a1_lut(n) |
| 1458 | 0 | A(2) = optics%a2_lut(n) |
| 1459 | ||
| 1460 | 0 | end function lookup_ohlmann_swpen |
| 1461 | ||
| 1462 | !> Get the opacity (decay scale) from Ohlmann lookup table | |
| 1463 | 0 | function 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 | 0 | if (chl > optics%chl_min) then |
| 1476 | 0 | log10chl = min(log10(chl),optics%log10chl_max) |
| 1477 | else | |
| 1478 | 0 | log10chl = optics%log10chl_min |
| 1479 | endif | |
| 1480 | ! Do a nearest neighbor lookup | |
| 1481 | 0 | n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1 |
| 1482 | ||
| 1483 | 0 | B(1) = optics%b1_lut(n) |
| 1484 | 0 | B(2) = optics%b2_lut(n) |
| 1485 | ||
| 1486 | 0 | return |
| 1487 | end function lookup_ohlmann_opacity | |
| 1488 | ||
| 1489 | 1 | subroutine 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 | 1 | if (allocated(CS%id_opacity)) & |
| 1494 | 1 | deallocate(CS%id_opacity) |
| 1495 | 1 | if (allocated(CS%opacity_coef)) & |
| 1496 | 0 | deallocate(CS%opacity_coef) |
| 1497 | 1 | if (allocated(CS%sw_pen_frac_coef)) & |
| 1498 | 0 | deallocate(CS%sw_pen_frac_coef) |
| 1499 | 1 | if (allocated(CS%chl_power)) & |
| 1500 | 0 | deallocate(CS%chl_power) |
| 1501 | 1 | if (allocated(optics%sw_pen_band)) & |
| 1502 | 1 | deallocate(optics%sw_pen_band) |
| 1503 | 1 | if (allocated(optics%opacity_band)) & |
| 1504 | 1 | deallocate(optics%opacity_band) |
| 1505 | 1 | if (allocated(optics%max_wavelength_band)) & |
| 1506 | 1 | deallocate(optics%max_wavelength_band) |
| 1507 | 1 | if (allocated(optics%min_wavelength_band)) & |
| 1508 | 1 | deallocate(optics%min_wavelength_band) |
| 1509 | 1 | if (allocated(optics%a1_lut)) deallocate(optics%a1_lut) |
| 1510 | 1 | if (allocated(optics%a2_lut)) deallocate(optics%a2_lut) |
| 1511 | 1 | if (allocated(optics%b1_lut)) deallocate(optics%b1_lut) |
| 1512 | 1 | if (allocated(optics%b2_lut)) deallocate(optics%b2_lut) |
| 1513 | 1 | end 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 | ||
| 1535 | 0 | end module MOM_opacity |