From 5726d941dc709c1feadfc9b60a9dac7781b49d30 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Fri, 27 Sep 2024 15:25:57 -0600 Subject: [PATCH 1/2] Don't compute log10(chl) for small chl if chl <= chl_min, we already have log10(chl_min) stored as a parameter. Avoiding the computation of log10(chl) in those cases prevents the possibility of a floating point exception, and lets us replace a max() statement with an if check so it shouldn't affect performance. This requires adding chl_min to the optics_type (which already has log10chl_min and log10chl_max). --- src/parameterizations/vertical/MOM_opacity.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 831607d2db..9d2339a440 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -45,6 +45,7 @@ module MOM_opacity !! Lookup tables for Ohlmann solar penetration scheme !! These would naturally exist as private module variables but that is prohibited in MOM6 real :: dlog10chl !< Chl increment within lookup table + real :: chl_min !< Lower bound of Chl in lookup table real :: log10chl_min !< Lower bound of Chl in lookup table real :: log10chl_max !< Upper bound of Chl in lookup table real, allocatable, dimension(:) :: a1_lut,& !< Coefficient for band 1 @@ -1303,6 +1304,7 @@ subroutine init_ohlmann_table(optics) call MOM_error(FATAL,"init_ohlmann: Cannot allocate lookup table") endif + optics%chl_min = chl_tab1a(1) optics%log10chl_min = log10(chl_tab1a(1)) optics%log10chl_max = log10(chl_tab1a(nval_tab1a)) optics%dlog10chl = (optics%log10chl_max - optics%log10chl_min)/(nval_lut-1) @@ -1349,7 +1351,11 @@ function lookup_ohlmann_swpen(chl,optics) result(A) integer :: n ! Make sure we are in the table - log10chl = max(optics%log10chl_min,min(log10(chl),optics%log10chl_max)) + if (chl > optics%chl_min) then + log10chl = min(log10(chl),optics%log10chl_max) + else + log10chl = optics%log10chl_min + endif ! Do a nearest neighbor lookup n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1 From 8fb8aaf83b9afb8e10ea87f1084d0066f2707b47 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Tue, 1 Oct 2024 13:53:13 -0600 Subject: [PATCH 2/2] Avoid chl=0 in lookup_ohlmann_opacity() previous commit only avoiding log10(0) in lookup_ohlmann_swpen() but we have another log10(chl) in the opacity function --- src/parameterizations/vertical/MOM_opacity.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 9d2339a440..b8b8c56d21 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -1377,7 +1377,11 @@ function lookup_ohlmann_opacity(chl,optics) result(B) integer :: n ! Make sure we are in the table - log10chl = max(optics%log10chl_min,min(log10(chl),optics%log10chl_max)) + if (chl > optics%chl_min) then + log10chl = min(log10(chl),optics%log10chl_max) + else + log10chl = optics%log10chl_min + endif ! Do a nearest neighbor lookup n = nint( (log10chl - optics%log10chl_min)/optics%dlog10chl ) + 1