Copy of 2e05202, Added a div0 protection to phosphorus cycling

rgknox · Jan 31, 2023 · 85e6083 · 85e6083
1 parent d09dfe8
commit 85e6083
Showing 1 changed file with 27 additions and 11 deletions.
diff --git a/components/elm/src/biogeochem/PhosphorusStateUpdate3Mod.F90 b/components/elm/src/biogeochem/PhosphorusStateUpdate3Mod.F90
@@ -183,18 +183,34 @@ subroutine PhosphorusStateUpdate3(bounds,num_soilc, filter_soilc, num_soilp, fil
 
                 if (temp_solutionp(c,j) < 0.0_r8) then
 
-                    col_pf%labilep_to_secondp_vr(c,j) = max(col_pf%labilep_to_secondp_vr(c,j),1.e-15)/ &
-                            (max(col_pf%labilep_to_secondp_vr(c,j),1.e-15)+col_pf%sminp_leached_vr(c,j))* &
-                            (temp_solutionp(c,j) + col_pf%labilep_to_secondp_vr(c,j)*dt + &
+                   if( abs(col_pf%labilep_to_secondp_vr(c,j)+col_pf%sminp_leached_vr(c,j)) >1.e-20_r8 )then
+
+                      col_pf%labilep_to_secondp_vr(c,j) = col_pf%labilep_to_secondp_vr(c,j)/ &
+                           (col_pf%labilep_to_secondp_vr(c,j)+col_pf%sminp_leached_vr(c,j))* &
+                           (temp_solutionp(c,j) + col_pf%labilep_to_secondp_vr(c,j)*dt + &
                             col_pf%sminp_leached_vr(c,j)*dt) /dt
-                    col_pf%sminp_leached_vr(c,j) = col_pf%sminp_leached_vr(c,j)/ &
-                            (col_pf%labilep_to_secondp_vr(c,j)+col_pf%sminp_leached_vr(c,j))* &
-                            (temp_solutionp(c,j) + col_pf%labilep_to_secondp_vr(c,j)*dt + &
-                            col_pf%sminp_leached_vr(c,j)*dt) /dt
-                       temp_solutionp(c,j) = 0.0_r8
-                       col_ps%solutionp_vr(c,j) = 0.0_r8
-                       col_ps%labilep_vr(c,j) = 0.0_r8
-                 else
+
+                      col_pf%sminp_leached_vr(c,j) = col_pf%sminp_leached_vr(c,j)/ &
+                           (col_pf%labilep_to_secondp_vr(c,j)+col_pf%sminp_leached_vr(c,j))* &
+                           (temp_solutionp(c,j) + col_pf%labilep_to_secondp_vr(c,j)*dt + &
+                           col_pf%sminp_leached_vr(c,j)*dt) /dt
+                   else
+                      ! If there is nothing there to drive proportions, just split it
+                      col_pf%labilep_to_secondp_vr(c,j) = 0.5_r8 * &
+                           (temp_solutionp(c,j) + col_pf%labilep_to_secondp_vr(c,j)*dt + &
+                           col_pf%sminp_leached_vr(c,j)*dt) /dt
+
+                      col_pf%sminp_leached_vr(c,j) = 0.5_r8 * &
+                           (temp_solutionp(c,j) + col_pf%labilep_to_secondp_vr(c,j)*dt + &
+                           col_pf%sminp_leached_vr(c,j)*dt) /dt
+
+                   end if
+
+                   temp_solutionp(c,j) = 0.0_r8
+                   col_ps%solutionp_vr(c,j) = 0.0_r8
+                   col_ps%labilep_vr(c,j) = 0.0_r8
+
+                else
                        ! sorbp = smax*solutionp/(ks+solutionp)
                        ! sorbp + solutionp = smax*solutionp/(ks+solutionp) + solutionp = total p pool after competition
                        ! solve quadratic function to get equilibrium solutionp and adsorbp pools