Added more gfortran debug flag option and fixed FPEs.

src/Makefile

@@ -21,6 +21,9 @@ $(info DEBUG setting: '$(DEBUG)')
 ifeq ($(DEBUG), 1)
   ifeq ($(findstring gfortran,$(notdir $(FC))),gfortran)
     FCFLAGS := -g -Wall -Wextra -Wconversion -Og -pedantic -fcheck=bounds -fmax-errors=5 -std=f2003 -fall-intrinsics
+# More Debug Flags
+#      FCFLAGS := -g -Og -Wall -Wextra -pedantic -fcheck=all -fmax-errors=5 -std=f2003 \
+#      -fall-intrinsics -fbacktrace -ffpe-trap=invalid,zero,overflow
   else ifeq ($(FC),ifort)
     FCFLAGS := -g -warn all -check bounds -implicitnone -O0 -error-limit 5
   endif

src/canopy_calcs.F90

@@ -90,7 +90,11 @@ SUBROUTINE canopy_calcs(nn)
                 ubzref   = sqrt((uref**2.0) + (vref**2.0))
 
 ! ... get scaled canopy model profile and sub-canopy layers
-                zhc         = zk/hcmref
+                if (hcmref > 0) then
+                    zhc = zk / hcmref
+                else
+                    zhc = 0
+                end if
                 cansublays  = min(floor(hcmref/modres),modlays)
 
 ! ... check for valid model vegetation types
@@ -164,7 +168,6 @@ SUBROUTINE canopy_calcs(nn)
                             end if
 
 ! ... calculate zero-plane displacement height/hc and surface (soil+veg) roughness lengths/hc
-
                             call canopy_zpd(zhc(1:cansublays), fafraczInt(1:cansublays), &
                                 ubzref, z0ghc, lambdars, cdrag, pai, hcmref, hgtref, &
                                 z0ref, vtyperef, lu_opt, z0_opt, d_h, zo_h)

src/canopy_profile_mod.F90

@@ -376,7 +376,11 @@ SUBROUTINE CANOPY_ZPD( ZHC, FCLAI, UBZREF, Z0GHC, &
         dha =  1.0 - (cosh(qstar*nrat*FCLAI(1))/cosh(qstar*nrat))
         fafraczInt_tota = IntegrateTrapezoid( ZHC,(cosh(qstar*nrat*FCLAI)*ZHC) )
         fafraczInt_totb = IntegrateTrapezoid( ZHC, cosh(qstar*nrat*FCLAI) )
-        dhb = fafraczInt_tota/fafraczInt_totb
+        if (fafraczInt_totb > 0) then
+            dhb = fafraczInt_tota/fafraczInt_totb
+        else
+            dhb = 0
+        end if
 
         ! Final zero-plane displacement (zpd) height
         d_h = dha * dhb