From d7337145ec3ae4aef8b4a9d6030672a8c1ed336c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 10 Dec 2021 10:09:26 -0500 Subject: [PATCH 1/4] (*)Fix extract_diabatic_member Return the diabatic_aux_CSp from extract_diabatic_member it is present as an optional argument. Somehow this was omitted when this routine was created, but without this correction the offline tracer mode returns a segmentation fault. Also, added the proper conversion factor in the register_diag_field call for e_predia, and internally calculate the interface heights in units of [Z ~> m] for dimensional consistency testing. All answers are bitwise identical in cases that ran before. --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 77ec87b230..1b68cf8211 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -294,7 +294,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & - eta ! Interface heights before diapycnal mixing [m]. + eta ! Interface heights before diapycnal mixing [Z ~> m] real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: & cn_IGW ! baroclinic internal gravity wave speeds [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: temp_diag ! Previous temperature for diagnostics [degC] @@ -326,7 +326,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_T_predia > 0) call post_data(CS%id_T_predia, tv%T, CS%diag) if (CS%id_S_predia > 0) call post_data(CS%id_S_predia, tv%S, CS%diag) if (CS%id_e_predia > 0) then - call find_eta(h, tv, G, GV, US, eta, eta_to_m=1.0, dZref=G%Z_ref) + call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref) call post_data(CS%id_e_predia, eta, CS%diag) endif @@ -2536,6 +2536,7 @@ subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, if (present(optics_CSp)) optics_CSp => CS%optics if (present(KPP_CSp)) KPP_CSp => CS%KPP_CSp if (present(energetic_PBL_CSp)) energetic_PBL_CSp => CS%energetic_PBL + if (present(diabatic_aux_CSp)) diabatic_aux_CSp => CS%diabatic_aux_CSp ! Constants within diabatic_CS if (present(evap_CFL_limit)) evap_CFL_limit = CS%evap_CFL_limit @@ -3175,7 +3176,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Layer Thickness before diabatic forcing', & trim(thickness_units), conversion=GV%H_to_MKS, v_extensive=.true.) CS%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, Time, & - 'Interface Heights before diabatic forcing', 'm') + 'Interface Heights before diabatic forcing', 'm', conversion=US%Z_to_m) if (use_temperature) then CS%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, Time, & 'Potential Temperature', 'degC') From 5172c495c3505a7565448d91ff48c487148f5500 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 10 Dec 2021 10:10:16 -0500 Subject: [PATCH 2/4] Warn if opacity_from_chl is called without fluxes Issue a warning with a helpful message if opacity_from_chl is called with no shortwave fluxes, and added logical tests to avoid a segmentation fault later in this routine. This should not happen, as it makes no sense, but it was occurring with the offline tracer code, and can be avoided by setting PEN_SW_NBANDS=0 if there are no shortwave fluxes to penetrate. Also turned the real dimensional parameter op_diag_len into a variable and set it immediately before where it is used. Many spelling errors were also corrected in MOM_opacity.F90. All answers are identical in cases that ran before. --- .../vertical/MOM_opacity.F90 | 89 +++++++++++-------- 1 file changed, 52 insertions(+), 37 deletions(-) diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 02d49d024d..a99524060b 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -47,7 +47,7 @@ module MOM_opacity end type optics_type -!> The control structure with paramters for the MOM_opacity module +!> The control structure with parameters for the MOM_opacity module type, public :: opacity_CS ; private logical :: var_pen_sw !< If true, use one of the CHL_A schemes (specified by OPACITY_SCHEME) to !! determine the e-folding depth of incoming shortwave radiation. @@ -67,6 +67,7 @@ module MOM_opacity !! The default is 10 m-1 - a value for muddy water. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. + logical :: warning_issued !< A flag that is used to avoid repetative warnings. !>@{ Diagnostic IDs integer :: id_sw_pen = -1, id_sw_vis_pen = -1 @@ -83,9 +84,6 @@ module MOM_opacity character*(10), parameter :: SINGLE_EXP_STRING = "SINGLE_EXP" !< String to specify the opacity scheme character*(10), parameter :: DOUBLE_EXP_STRING = "DOUBLE_EXP" !< String to specify the opacity scheme -real, parameter :: op_diag_len = 1e-10 !< Lengthscale L used to remap opacity - !! from op to 1/L * tanh(op * L) - contains !> This sets the opacity of sea water based based on one of several different schemes. @@ -103,24 +101,26 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_ type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(opacity_CS) :: CS !< The control structure earlier set up by opacity_init. real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentractions [mg m-3] + optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentrations [mg m-3] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - optional, intent(in) :: chl_3d !< The chlorophyll-A concentractions of each layer [mg m-3] + optional, intent(in) :: chl_3d !< The chlorophyll-A concentrations of each layer [mg m-3] ! Local variables integer :: i, j, k, n, is, ie, js, je, nz real :: inv_sw_pen_scale ! The inverse of the e-folding scale [m-1]. real :: Inv_nbands ! The inverse of the number of bands of penetrating - ! shortwave radiation. + ! shortwave radiation [nondim] logical :: call_for_surface ! if horizontal slice is the surface layer real :: tmp(SZI_(G),SZJ_(G),SZK_(GV)) ! A 3-d temporary array. real :: chl(SZI_(G),SZJ_(G),SZK_(GV)) ! The concentration of chlorophyll-A [mg m-3]. real :: Pen_SW_tot(SZI_(G),SZJ_(G)) ! The penetrating shortwave radiation ! summed across all bands [Q R Z T-1 ~> W m-2]. + real :: op_diag_len ! A tiny lengthscale [m] used to remap opacity + ! from op to 1/op_diag_len * tanh(op * op_diag_len) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (present(chl_2d) .or. present(chl_3d)) then - ! The optical properties are based on cholophyll concentrations. + ! The optical properties are based on chlorophyll concentrations. call opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, & G, GV, US, CS, chl_2d, chl_3d) else ! Use sw e-folding scale set by MOM_input @@ -199,11 +199,12 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_ call post_data(CS%id_sw_vis_pen, Pen_SW_tot, CS%diag) endif do n=1,optics%nbands ; if (CS%id_opacity(n) > 0) then + op_diag_len = 1e-10 ! A dimensional depth to constrain the range of opacity [m] !$OMP parallel do default(shared) do k=1,nz ; do j=js,je ; do i=is,ie ! Remap opacity (op) to 1/L * tanh(op * L) where L is one Angstrom. ! This gives a nearly identical value when op << 1/L but allows one to - ! store the values when opacity is divergent (i.e. opaque). + ! record the values even at reduced precision when opacity is huge (i.e. opaque). tmp(i,j,k) = tanh(op_diag_len * optics%opacity_band(n,i,j,k)) / op_diag_len enddo ; enddo ; enddo call post_data(CS%id_opacity(n), tmp, CS%diag) @@ -213,12 +214,12 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_ end subroutine set_opacity -!> This sets the "blue" band opacity based on chloophyll A concencentrations +!> This sets the "blue" band opacity based on chlorophyll A concentrations !! The red portion is lumped into the net heating at the surface. subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_nir_dif, & G, GV, US, CS, chl_2d, chl_3d) type(optics_type), intent(inout) :: optics !< An optics structure that has values - !! set based on the opacities. + !! set based on the opacities. real, dimension(:,:), pointer :: sw_total !< Total shortwave flux into the ocean [Q R Z T-1 ~> W m-2] real, dimension(:,:), pointer :: sw_vis_dir !< Visible, direct shortwave into the ocean [Q R Z T-1 ~> W m-2] real, dimension(:,:), pointer :: sw_vis_dif !< Visible, diffuse shortwave into the ocean [Q R Z T-1 ~> W m-2] @@ -229,15 +230,15 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(opacity_CS) :: CS !< The control structure. real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentractions [mg m-3] + optional, intent(in) :: chl_2d !< Vertically uniform chlorophyll-A concentrations [mg m-3] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - optional, intent(in) :: chl_3d !< A 3-d field of chlorophyll-A concentractions [mg m-3] + optional, intent(in) :: chl_3d !< A 3-d field of chlorophyll-A concentrations [mg m-3] real :: chl_data(SZI_(G),SZJ_(G)) ! The chlorophyll A concentrations in a layer [mg m-3]. real :: Inv_nbands ! The inverse of the number of bands of penetrating - ! shortwave radiation. + ! shortwave radiation [nondim] real :: Inv_nbands_nir ! The inverse of the number of bands of penetrating - ! near-infrafed radiation. + ! near-infrared radiation [nondim] real :: SW_pen_tot ! The sum across the bands of the penetrating ! shortwave radiation [Q R Z T-1 ~> W m-2]. real :: SW_vis_tot ! The sum across the visible bands of shortwave @@ -247,7 +248,7 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir type(time_type) :: day character(len=128) :: mesg integer :: i, j, k, n, is, ie, js, je, nz, nbands - logical :: multiband_vis_input, multiband_nir_input + logical :: multiband_vis_input, multiband_nir_input, total_sw_input is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -257,9 +258,9 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir ! into the net heating at the surface. ! ! Morel, A., Optical modeling of the upper ocean in relation to its biogenous -! matter content (case-i waters).,J. Geo. Res., {93}, 10,749--10,768, 1988. +! matter content (case-i waters)., J. Geo. Res., {93}, 10,749--10,768, 1988. ! -! Manizza, M., C.~L. Quere, A.~Watson, and E.~T. Buitenhuis, Bio-optical +! Manizza, M., C. L. Quere, A. Watson, and E. T. Buitenhuis, Bio-optical ! feedbacks among phytoplankton, upper ocean physics and sea-ice in a ! global model, Geophys. Res. Let., , L05,603, 2005. @@ -271,10 +272,19 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir if (nbands <= 2) then ; Inv_nbands_nir = 0.0 else ; Inv_nbands_nir = 1.0 / real(nbands - 2.0) ; endif - multiband_vis_input = (associated(sw_vis_dir) .and. & - associated(sw_vis_dif)) - multiband_nir_input = (associated(sw_nir_dir) .and. & - associated(sw_nir_dif)) + if (.not.(associated(sw_total) .or. (associated(sw_vis_dir) .and. associated(sw_vis_dif) .and. & + associated(sw_nir_dir) .and. associated(sw_nir_dif)) )) then + if (.not.CS%warning_issued) then + call MOM_error(WARNING, & + "opacity_from_chl called without any shortwave flux arrays allocated.\n"//& + "Consider setting PEN_SW_NBANDS = 0 if no shortwave fluxes are being used.") + endif + CS%warning_issued = .true. + endif + + multiband_vis_input = (associated(sw_vis_dir) .and. associated(sw_vis_dif)) + multiband_nir_input = (associated(sw_nir_dir) .and. associated(sw_nir_dif)) + total_sw_input = associated(sw_total) chl_data(:,:) = 0.0 if (present(chl_3d)) then @@ -298,7 +308,7 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir endif enddo ; enddo else - call MOM_error(FATAL, "Either chl_2d or chl_3d must be preesnt in a call to opacity_form_chl.") + call MOM_error(FATAL, "Either chl_2d or chl_3d must be present in a call to opacity_form_chl.") endif select case (CS%opacity_scheme) @@ -309,12 +319,13 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir if (G%mask2dT(i,j) > 0.5) then if (multiband_vis_input) then SW_vis_tot = sw_vis_dir(i,j) + sw_vis_dif(i,j) - else ! Follow Manizza 05 in assuming that 42% of SW is visible. + elseif (total_sw_input) then + ! Follow Manizza 05 in assuming that 42% of SW is visible. SW_vis_tot = 0.42 * sw_total(i,j) endif if (multiband_nir_input) then SW_nir_tot = sw_nir_dir(i,j) + sw_nir_dif(i,j) - else + elseif (total_sw_input) then SW_nir_tot = sw_total(i,j) - SW_vis_tot endif endif @@ -333,11 +344,13 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir !$OMP parallel do default(shared) private(SW_pen_tot) do j=js,je ; do i=is,ie SW_pen_tot = 0.0 - if (G%mask2dT(i,j) > 0.5) then ; if (multiband_vis_input) then + if (G%mask2dT(i,j) > 0.5) then + if (multiband_vis_input) then SW_pen_tot = SW_pen_frac_morel(chl_data(i,j)) * (sw_vis_dir(i,j) + sw_vis_dif(i,j)) - else + elseif (total_sw_input) then SW_pen_tot = SW_pen_frac_morel(chl_data(i,j)) * 0.5*sw_total(i,j) - endif ; endif + endif + endif do n=1,nbands optics%sw_pen_band(n,i,j) = Inv_nbands*sw_pen_tot @@ -395,7 +408,7 @@ function opacity_morel(chl_data) real :: opacity_morel !< The returned opacity [m-1] ! The following are coefficients for the optical model taken from Morel and - ! Antoine (1994). These coeficients represent a non uniform distribution of + ! Antoine (1994). These coefficients represent a non uniform distribution of ! chlorophyll-a through the water column. Other approaches may be more ! appropriate when using an interactive ecosystem model that predicts ! three-dimensional chl-a values. @@ -415,7 +428,7 @@ function SW_pen_frac_morel(chl_data) real :: SW_pen_frac_morel !< The returned penetrating shortwave fraction [nondim] ! The following are coefficients for the optical model taken from Morel and - ! Antoine (1994). These coeficients represent a non uniform distribution of + ! Antoine (1994). These coefficients represent a non uniform distribution of ! chlorophyll-a through the water column. Other approaches may be more ! appropriate when using an interactive ecosystem model that predicts ! three-dimensional chl-a values. @@ -496,7 +509,7 @@ function optics_nbands(optics) optics_nbands = optics%nbands end function optics_nbands -!> Apply shortwave heating below the boundary layer (when running with the bulk mixed layer inhereted +!> Apply shortwave heating below the boundary layer (when running with the bulk mixed layer inherited !! from GOLD) or throughout the water column. !! !! In addition, it causes all of the remaining SW radiation to be absorbed, provided that the total @@ -515,7 +528,7 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l real, dimension(SZI_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. real, dimension(max(1,nsw),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< Opacity in each band of penetrating !! shortwave radiation [H-1 ~> m-1 or m2 kg-1]. - !! The indicies are band, i, k. + !! The indices are band, i, k. type(optics_type), intent(in) :: optics !< An optics structure that has values of !! opacities and shortwave fluxes. integer, intent(in) :: j !< j-index to work on. @@ -548,7 +561,7 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l real, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: eps !< Small thickness that must remain in !! each layer, and which will not be !! subject to heating [H ~> m or kg m-2] - integer, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: ksort !< Density-sorted k-indicies. + integer, dimension(SZI_(G),SZK_(GV)), optional, intent(in) :: ksort !< Density-sorted k-indices. real, dimension(SZI_(G)), optional, intent(in) :: htot !< Total mixed layer thickness [H ~> m or kg m-2]. real, dimension(SZI_(G)), optional, intent(inout) :: Ttot !< Depth integrated mixed layer !! temperature [degC H ~> degC m or degC kg m-2] @@ -912,7 +925,7 @@ end subroutine sumSWoverBands -!> This routine initalizes the opacity module, including an optics_type. +!> This routine initializes the opacity module, including an optics_type. subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) type(time_type), target, intent(in) :: Time !< The current model time. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -922,7 +935,7 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) !! parameters. type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic !! output. - type(opacity_CS) :: CS !< Opacity control struct + type(opacity_CS) :: CS !< Opacity control structure type(optics_type) :: optics !< An optics structure that has parameters !! set and arrays allocated here. ! Local variables @@ -1083,6 +1096,8 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) "The value to use for opacity over land. The default is "//& "10 m-1 - a value for muddy water.", units="m-1", default=10.0) + CS%warning_issued = .false. + if (.not.allocated(optics%opacity_band)) & allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz)) if (.not.allocated(optics%sw_pen_band)) & @@ -1106,7 +1121,7 @@ end subroutine opacity_init subroutine opacity_end(CS, optics) - type(opacity_CS) :: CS !< Opacity control struct + type(opacity_CS) :: CS !< Opacity control structure type(optics_type) :: optics !< An optics type structure that should be deallocated. if (allocated(CS%id_opacity)) & @@ -1125,7 +1140,7 @@ end subroutine opacity_end !! !! opacity_from_chl: !! In this routine, the Morel (modified) or Manizza (modified) -!! schemes use the "blue" band in the paramterizations to determine +!! schemes use the "blue" band in the parameterizations to determine !! the e-folding depth of the incoming shortwave attenuation. The red !! portion is lumped into the net heating at the surface. !! From 049241ce1622355fe7f1639294fc06deafd061b1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 16 Dec 2021 05:22:10 -0500 Subject: [PATCH 3/4] +Rescaled optics%opacity_band Rescaled the units of optics%opacity_band from [m-1] to [Z-1 ~> m-1], with the opacity values returned from extract_optics_slice also rescaled by the same factor, which can be offset by compensating changes to the opacity_scale argument. Also rescaled 4 other internal variables and documented the units on 3 more. One uncommon parameter (SW_1ST_EXP_RATIO) listed the wrong units in its get_param call, and this was corrected, but turned out not to have been logged for any of the MOM6-examples test cases. Some compensating changes were also made in the MOM_generic_tracer module, which directly accesses the contents of the optics_type (thereby preventing it from being opaque). All answers and output in the MOM6-examples test suite are bitwise identical. --- .../vertical/MOM_bulk_mixed_layer.F90 | 2 +- .../vertical/MOM_diabatic_aux.F90 | 4 +- .../vertical/MOM_opacity.F90 | 56 +++++++++---------- src/tracer/MOM_generic_tracer.F90 | 23 ++++---- 4 files changed, 42 insertions(+), 43 deletions(-) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 046329523d..926eaaa013 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -424,7 +424,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C eps(i,k) = 0.0 ; if (k > nkmb) eps(i,k) = GV%Angstrom_H T(i,k) = tv%T(i,j,k) ; S(i,k) = tv%S(i,j,k) enddo ; enddo - if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacity_band, opacity_scale=GV%H_to_m) + if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacity_band, opacity_scale=GV%H_to_Z) do k=1,nz ; do i=is,ie d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 312d114dde..b0ca5a931e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -1105,7 +1105,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! To accommodate vanishing upper layers, we need to allow for an instantaneous ! distribution of forcing over some finite vertical extent. The bulk mixed layer ! code handles this issue properly. - H_limit_fluxes = max(GV%Angstrom_H, 1.E-30*GV%m_to_H) + H_limit_fluxes = max(GV%Angstrom_H, 1.0e-30*GV%m_to_H) ! diagnostic to see if need to create mass to avoid grounding if (CS%id_createdH>0) CS%createdH(:,:) = 0. @@ -1160,7 +1160,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! Nothing more is done on this j-slice if there is no buoyancy forcing. if (.not.associated(fluxes%sw)) cycle - if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H)) + if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%Z_to_H)) ! The surface forcing is contained in the fluxes type. ! We aggregate the thermodynamic forcing for a time step into the following: diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index a99524060b..658170beda 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -25,7 +25,7 @@ module MOM_opacity type, public :: optics_type integer :: nbands !< The number of penetrating bands of SW radiation - real, allocatable :: opacity_band(:,:,:,:) !< SW optical depth per unit thickness [m-1] + real, allocatable :: opacity_band(:,:,:,:) !< SW optical depth per unit thickness [Z-1 ~> m-1] !! The number of radiation bands is most rapidly varying (first) index. real, allocatable :: sw_pen_band(:,:,:) !< shortwave radiation [Q R Z T-1 ~> W m-2] @@ -55,15 +55,15 @@ module MOM_opacity !! water properties into the opacity (i.e., the e-folding depth) and !! (perhaps) the number of bands of penetrating shortwave radiation to use. real :: pen_sw_scale !< The vertical absorption e-folding depth of the - !! penetrating shortwave radiation [m]. + !! penetrating shortwave radiation [Z ~> m]. real :: pen_sw_scale_2nd !< The vertical absorption e-folding depth of the - !! (2nd) penetrating shortwave radiation [m]. - real :: SW_1ST_EXP_RATIO !< Ratio for 1st exp decay in Two Exp decay opacity + !! (2nd) penetrating shortwave radiation [Z ~> m]. + real :: SW_1ST_EXP_RATIO !< Ratio for 1st exp decay in Two Exp decay opacity [nondim] real :: pen_sw_frac !< The fraction of shortwave radiation that is - !! penetrating with a constant e-folding approach. + !! penetrating with a constant e-folding approach [nondim] real :: blue_frac !< The fraction of the penetrating shortwave !! radiation that is in the blue band [nondim]. - real :: opacity_land_value !< The value to use for opacity over land [m-1]. + real :: opacity_land_value !< The value to use for opacity over land [Z-1 ~> m-1]. !! The default is 10 m-1 - a value for muddy water. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to !! regulate the timing of diagnostic output. @@ -107,15 +107,15 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_ ! Local variables integer :: i, j, k, n, is, ie, js, je, nz - real :: inv_sw_pen_scale ! The inverse of the e-folding scale [m-1]. + real :: inv_sw_pen_scale ! The inverse of the e-folding scale [Z-1 ~> m-1]. real :: Inv_nbands ! The inverse of the number of bands of penetrating ! shortwave radiation [nondim] logical :: call_for_surface ! if horizontal slice is the surface layer - real :: tmp(SZI_(G),SZJ_(G),SZK_(GV)) ! A 3-d temporary array. + real :: tmp(SZI_(G),SZJ_(G),SZK_(GV)) ! A 3-d temporary array for diagnosing opacity [Z-1 ~> m-1] real :: chl(SZI_(G),SZJ_(G),SZK_(GV)) ! The concentration of chlorophyll-A [mg m-3]. real :: Pen_SW_tot(SZI_(G),SZJ_(G)) ! The penetrating shortwave radiation ! summed across all bands [Q R Z T-1 ~> W m-2]. - real :: op_diag_len ! A tiny lengthscale [m] used to remap opacity + real :: op_diag_len ! A tiny lengthscale [Z ~> m] used to remap diagnostics of opacity ! from op to 1/op_diag_len * tanh(op * op_diag_len) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -128,14 +128,14 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_ else ; Inv_nbands = 1.0 / real(optics%nbands) ; endif ! Make sure there is no division by 0. - inv_sw_pen_scale = 1.0 / max(CS%pen_sw_scale, 0.1*GV%Angstrom_m, & - GV%H_to_m*GV%H_subroundoff) + inv_sw_pen_scale = 1.0 / max(CS%pen_sw_scale, 0.1*GV%Angstrom_Z, & + GV%H_to_Z*GV%H_subroundoff) if ( CS%Opacity_scheme == DOUBLE_EXP ) then !$OMP parallel do default(shared) do k=1,nz ; do j=js,je ; do i=is,ie optics%opacity_band(1,i,j,k) = inv_sw_pen_scale optics%opacity_band(2,i,j,k) = 1.0 / max(CS%pen_sw_scale_2nd, & - 0.1*GV%Angstrom_m,GV%H_to_m*GV%H_subroundoff) + 0.1*GV%Angstrom_Z, GV%H_to_Z*GV%H_subroundoff) enddo ; enddo ; enddo if (.not.associated(sw_total) .or. (CS%pen_SW_scale <= 0.0)) then !$OMP parallel do default(shared) @@ -199,7 +199,7 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_ call post_data(CS%id_sw_vis_pen, Pen_SW_tot, CS%diag) endif do n=1,optics%nbands ; if (CS%id_opacity(n) > 0) then - op_diag_len = 1e-10 ! A dimensional depth to constrain the range of opacity [m] + op_diag_len = 1.0e-10*US%m_to_Z ! A minimal extinction depth to constrain the range of opacity [Z ~> m] !$OMP parallel do default(shared) do k=1,nz ; do j=js,je ; do i=is,ie ! Remap opacity (op) to 1/L * tanh(op * L) where L is one Angstrom. @@ -375,18 +375,18 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir enddo else ! Band 1 is Manizza blue. - optics%opacity_band(1,i,j,k) = 0.0232 + 0.074*chl_data(i,j)**0.674 + optics%opacity_band(1,i,j,k) = (0.0232 + 0.074*chl_data(i,j)**0.674) * US%Z_to_m if (nbands >= 2) & ! Band 2 is Manizza red. - optics%opacity_band(2,i,j,k) = 0.225 + 0.037*chl_data(i,j)**0.629 + optics%opacity_band(2,i,j,k) = (0.225 + 0.037*chl_data(i,j)**0.629) * US%Z_to_m ! All remaining bands are NIR, for lack of something better to do. - do n=3,nbands ; optics%opacity_band(n,i,j,k) = 2.86 ; enddo + do n=3,nbands ; optics%opacity_band(n,i,j,k) = 2.86*US%Z_to_m ; enddo endif enddo ; enddo case (MOREL_88) do j=js,je ; do i=is,ie optics%opacity_band(1,i,j,k) = CS%opacity_land_value if (G%mask2dT(i,j) > 0.5) & - optics%opacity_band(1,i,j,k) = opacity_morel(chl_data(i,j)) + optics%opacity_band(1,i,j,k) = US%Z_to_m * opacity_morel(chl_data(i,j)) do n=2,optics%nbands optics%opacity_band(n,i,j,k) = optics%opacity_band(1,i,j,k) @@ -460,7 +460,8 @@ subroutine extract_optics_slice(optics, j, G, GV, opacity, opacity_scale, penSW_ type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(max(optics%nbands,1),SZI_(G),SZK_(GV)), & - optional, intent(out) :: opacity !< The opacity in each band, i-point, and layer + optional, intent(out) :: opacity !< The opacity in each band, i-point, and layer [Z-1 ~> m-1], + !! but with units that can be altered by opacity_scale. real, optional, intent(in) :: opacity_scale !< A factor by which to rescale the opacity. real, dimension(max(optics%nbands,1),SZI_(G)), & optional, intent(out) :: penSW_top !< The shortwave radiation [Q R Z T-1 ~> W m-2] @@ -866,7 +867,7 @@ subroutine sumSWoverBands(G, GV, US, h, nsw, optics, j, dt, & if (h(i,k) > 0.0) then do n=1,nsw ; if (Pen_SW_bnd(n,i) > 0.0) then ! SW_trans is the SW that is transmitted THROUGH the layer - opt_depth = h(i,k)*GV%H_to_m * optics%opacity_band(n,i,j,k) + opt_depth = h(i,k)*GV%H_to_Z * optics%opacity_band(n,i,j,k) exp_OD = exp(-opt_depth) SW_trans = exp_OD @@ -1015,19 +1016,18 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) call MOM_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5) endif call get_param(param_file, mdl, "PEN_SW_SCALE", CS%pen_sw_scale, & - "The vertical absorption e-folding depth of the "//& - "penetrating shortwave radiation.", units="m", default=0.0) + "The vertical absorption e-folding depth of the penetrating shortwave radiation.", & + units="m", default=0.0, scale=US%m_to_Z) !BGR/ Added for opacity_scheme==double_exp read in 2nd exp-decay and fraction if (CS%Opacity_scheme == DOUBLE_EXP ) then call get_param(param_file, mdl, "PEN_SW_SCALE_2ND", CS%pen_sw_scale_2nd, & "The (2nd) vertical absorption e-folding depth of the "//& - "penetrating shortwave radiation "//& - "(use if SW_EXP_MODE==double.)",& - units="m", default=0.0) + "penetrating shortwave radiation (use if SW_EXP_MODE==double.)", & + units="m", default=0.0, scale=US%m_to_Z) call get_param(param_file, mdl, "SW_1ST_EXP_RATIO", CS%sw_1st_exp_ratio, & "The fraction of 1st vertical absorption e-folding depth "//& "penetrating shortwave radiation if SW_EXP_MODE==double.",& - units="m", default=0.0) + units="nondim", default=0.0) elseif (CS%OPACITY_SCHEME == Single_Exp) then !/Else disable 2nd_exp scheme CS%pen_sw_scale_2nd = 0.0 @@ -1094,12 +1094,12 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) call get_param(param_file, mdl, "OPACITY_LAND_VALUE", CS%opacity_land_value, & "The value to use for opacity over land. The default is "//& - "10 m-1 - a value for muddy water.", units="m-1", default=10.0) + "10 m-1 - a value for muddy water.", units="m-1", default=10.0, scale=US%Z_to_m) CS%warning_issued = .false. if (.not.allocated(optics%opacity_band)) & - allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz)) + allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz), source=0.0) if (.not.allocated(optics%sw_pen_band)) & allocate(optics%sw_pen_band(optics%nbands,isd:ied,jsd:jed)) allocate(CS%id_opacity(optics%nbands), source=-1) @@ -1114,7 +1114,7 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics) longname = 'Opacity for shortwave radiation in band '//trim(adjustl(bandnum)) & // ', saved as L^-1 tanh(Opacity * L) for L = 10^-10 m' CS%id_opacity(n) = register_diag_field('ocean_model', shortname, diag%axesTL, Time, & - longname, 'm-1') + longname, 'm-1', conversion=US%m_to_Z) enddo end subroutine opacity_init diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index 4627d0ec80..fbde0dc04e 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -54,7 +54,7 @@ module MOM_generic_tracer implicit none ; private - !> An state hidden in module data that is very much not allowed in MOM6 + !> A state hidden in module data that is very much not allowed in MOM6 ! ### This needs to be fixed logical :: g_registered = .false. @@ -83,13 +83,8 @@ module MOM_generic_tracer !> Pointer to the first element of the linked list of generic tracers. type(g_tracer_type), pointer :: g_tracer_list => NULL() - integer :: H_to_m !< Auxiliary to access GV%H_to_m in routines that do not have access to GV - end type MOM_generic_tracer_CS -! This include declares and sets the variable "version". -#include "version_variable.h" - contains !> Initializes the generic tracer packages and adds their tracers to the list @@ -104,9 +99,12 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) !! advection and diffusion module. type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct -! Local variables + ! Local variables logical :: register_MOM_generic_tracer + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer' character(len=200) :: inputdir ! The directory where NetCDF input files are. ! These can be overridden later in via the field manager? @@ -381,8 +379,6 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, call g_tracer_set_csdiag(CS%diag) #endif - CS%H_to_m = GV%H_to_m - end subroutine initialize_MOM_generic_tracer !> Column physics for generic tracers. @@ -503,7 +499,8 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! !Calculate tendencies (i.e., field changes at dt) from the sources / sinks ! - if ((G%US%L_to_m == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. (G%US%s_to_T == 1.0)) then + if ((G%US%L_to_m == 1.0) .and. (G%US%s_to_T == 1.0) .and. (G%US%Z_to_m == 1.0) .and. & + (G%US%Q_to_J_kg == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0)) then ! Avoid unnecessary copies when no unit conversion is needed. call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & G%areaT, get_diag_time_end(CS%diag), & @@ -512,7 +509,9 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, else call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & - optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, & + optics%nbands, optics%max_wavelength_band, & + sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & + opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & internal_heat=G%US%RZ_to_kg_m2*tv%internal_heat(:,:), & frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga) endif @@ -864,7 +863,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS) !nnz: fake rho0 rho0=1.0 - dzt(:,:,:) = CS%H_to_m * h(:,:,:) + dzt(:,:,:) = GV%H_to_m * h(:,:,:) sosga = global_area_mean(sfc_state%SSS, G) From 0544f9f2e31eec2e14f98720700a993907cf5ab8 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 16 Dec 2021 06:50:29 -0500 Subject: [PATCH 4/4] +(*)Avoid segmentation faults if PEN_SW_NBANDS = 0 Modified three routines in MOM_opacity to avoid segmentation faults if PEN_SW_NBANDS = 0, and hence if the optics type is not being allocated. In the case of optics_nbands(), this involved formally changing an optics_type argument into a pointer to an optics_type. (Pointers to an optics_type were already been used as the argument in all calls to optics_nbands(), but it was not always associated.) In two other routines, the change is simply to add a return call if there are 0 bands of shortwave radiation. With these changes, the single column test cases with no penetrating shortwave radiation now successfully run if PEN_SW_NBANDS = 0 and give answers that are identical to those obtained with PEN_SW_NBANDS = 1. All answers and output in cases that ran previously are bitwise identical, but there is a subtle change in a public interface. --- .../vertical/MOM_opacity.F90 | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 658170beda..9aa8fafd14 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -503,11 +503,15 @@ end subroutine extract_optics_fields !> Return the number of bands of penetrating shortwave radiation. function optics_nbands(optics) - type(optics_type), intent(in) :: optics !< An optics structure that has values of opacities + type(optics_type), pointer :: optics !< An optics structure that has values of opacities !! and shortwave fluxes. integer :: optics_nbands !< The number of penetrating bands of SW radiation - optics_nbands = optics%nbands + if (associated(optics)) then + optics_nbands = optics%nbands + else + optics_nbands = 0 + endif end function optics_nbands !> Apply shortwave heating below the boundary layer (when running with the bulk mixed layer inherited @@ -617,8 +621,10 @@ subroutine absorbRemainingSW(G, GV, US, h, opacity_band, nsw, optics, j, dt, H_l ! TKE budget of the shortwave heating. real :: C1_6, C1_60 integer :: is, ie, nz, i, k, ks, n - SW_Remains = .false. + if (nsw < 1) return + + SW_Remains = .false. min_SW_heat = optics%PenSW_flux_absorb * dt I_Habs = optics%PenSW_absorb_Invlen @@ -842,12 +848,16 @@ subroutine sumSWoverBands(G, GV, US, h, nsw, optics, j, dt, & integer :: is, ie, nz, i, k, ks, n SW_Remains = .false. - min_SW_heat = optics%PenSW_flux_absorb*dt ! Default of 2.5e-11*US%T_to_s*GV%m_to_H I_Habs = 1e3*GV%H_to_m ! optics%PenSW_absorb_Invlen h_min_heat = 2.0*GV%Angstrom_H + GV%H_subroundoff is = G%isc ; ie = G%iec ; nz = GV%ke + if (nsw < 1) then + netPen(:,:) = 0.0 + return + endif + pen_SW_bnd(:,:) = iPen_SW_bnd(:,:) do i=is,ie ; h_heat(i) = 0.0 ; enddo do i=is,ie @@ -859,6 +869,7 @@ subroutine sumSWoverBands(G, GV, US, h, nsw, optics, j, dt, & ! Apply penetrating SW radiation to remaining parts of layers. ! Excessively thin layers are not heated to avoid runaway temps. + min_SW_heat = optics%PenSW_flux_absorb*dt ! Default of 2.5e-11*US%T_to_s*GV%m_to_H do k=1,nz do i=is,ie