From 714b4b9262332b3263a82edc0e87f30f0dcaca00 Mon Sep 17 00:00:00 2001 From: Michael Levy Date: Fri, 27 Sep 2024 08:32:34 -0600 Subject: [PATCH] Refactor computation of K Paulot et al define K = (1/kg + H/kw)^-1 rewriting this expression as K = ((kw + kg*H) / (kg*kw))^-1 = (kg*kw) / (kw + kg*H) allows for situations where either kg or kw (but not both) is zero. It's also computationally more efficient, replacing three divisions and an addition with two multiplications, an addition, and a division (effectively replacing two divisions with multiplications) --- src/marbl_nhx_surface_emis_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/marbl_nhx_surface_emis_mod.F90 b/src/marbl_nhx_surface_emis_mod.F90 index 2f986ba1..1e19d6fb 100644 --- a/src/marbl_nhx_surface_emis_mod.F90 +++ b/src/marbl_nhx_surface_emis_mod.F90 @@ -72,7 +72,7 @@ subroutine marbl_nhx_surface_emis_compute( & call marbl_comp_Hstar_nhx(num_elements, ph, sst, sss, Hstar_nhx) - K(:) = c1 / (c1 / kg_nh3(:) + Hstar_nhx / kw_nh3(:)) + K(:) = (kg_nh3(:) * kw_nh3(:)) / (kw_nh3(:) + Hstar_nhx * kg_nh3(:)) nhx_surface_emis(:) = (c1 - ifrac(:)) * K(:) * Hstar_nhx(:) * max(nh4(:),c0) end subroutine marbl_nhx_surface_emis_compute