Bug fix: Prevent div-by-zero in apm_phys_mod.F

This commit fixes the problem raised in geoschem/geos-chem #1100. We now (1) set YJ = 0 and (2) prevent division by YJ if YJ is very small. This only affects APM simulations in GEOS-Chem Classic. Signed-off-by: Bob Yantosca <[email protected]>
geoschem · Jan 19, 2022 · 1b55678 · 1b55678
1 parent 13e422e
commit 1b55678
Showing 1 changed file with 20 additions and 5 deletions.
diff --git a/APM/apm_phys_mod.F b/APM/apm_phys_mod.F
@@ -1344,6 +1344,10 @@ SUBROUTINE APM_PHYS( II,        JJ,      LL,       NCOAG1,
       ENDIF
 
       DTNGC = DT
+
+      ! Initialize YJ to avoid div-by-zero errors (bmy, 19 Jan 2022)
+      YJ = 0d0
+
 !      IF(CACID.GT.5.d5)THEN
       IF(CACID.GT.5.d5.or.(IFNUCL.GE.13.and.CHOM.GT.5.d5))THEN
  80      CONTINUE
@@ -1372,11 +1376,22 @@ SUBROUTINE APM_PHYS( II,        JJ,      LL,       NCOAG1,
             stop
          ENDIF
 
-         IF(DT.GT.(DNMAX/YJ)) THEN
-            DTNGC = DNMAX/YJ    ! small timestep for large J
-            DTNGC = MAX(DTNGC,1.8d2)  !set minimum DTNGC
-         ELSE
-            DTNGC = DT
+         !--------------------------------------------------------
+         ! NOTE: Rewrite code to avoid div-by-zero if YJ=0
+         !  -- Bob Yantosca (19 Jan 2022)
+         !IF(DT.GT.(DNMAX/YJ)) THEN
+         !   DTNGC = DNMAX/YJ    ! small timestep for large J
+         !   DTNGC = MAX(DTNGC,1.8d2)  !set minimum DTNGC
+         !ELSE
+         !   DTNGC = DT
+         !ENDIF
+         !--------------------------------------------------------
+         DTNGC = DT
+         IF ( ABS(YJ) > 1.0d-30 ) THEN
+            IF(DT.GT.(DNMAX/YJ)) THEN
+               DTNGC = DNMAX/YJ ! small timestep for large J
+               DTNGC = MAX(DTNGC,1.8d2) !set minimum DTNGC
+            ENDIF
          ENDIF
 
          IF(DTLEFT.GT.(1.5*DTNGC)) THEN