ref/tools/chopmod [FIX]: avoid division by zero

src/synt/ref/tools/chopmod.f

@@ -36,10 +36,12 @@ c                   for a dominant period
 c V1.7   17/03/99   provide different density transformations
 c V1.8   22/12/11   migrated to gfortran
 c V1.9   22/12/11   support input of files with transverse isotropic model 
+c V1.10  08/04/19   set appropriate water level in computation of Q
+c                   values in order to avoid division by zero
 c
       program chopmod
       character*70 version
-      parameter(version='CHOPMOD   V1.9   spherical --> flat earth')
+      parameter(version='CHOPMOD   V1.10   spherical --> flat earth')
 c----------------------------------------------------------------------
 c 
 c STRATEGY
@@ -117,7 +119,7 @@ c dispersion
       logical disperse
       real*8 domper, q, nureforig
 c calculations
-      double precision zos, zus, zzo, zzu, zzm, rrm
+      double precision zos, zus, zzo, zzu, zzm, rrm, qmwl, qkwl
       double precision zualpha, zubeta, zurho, zsalpha, zsbeta, zsrho
 c commandline
       integer maxopt, lastarg 
@@ -375,10 +377,13 @@ c of wavelength in the full path an therefor the same damping under a given
 c quality factor.
 c
         qbeta(layer)=qm(nl)
-        qalpha(layer)=(1.d0/qm(nl))
+        qmwl=max(1.e-5,qm(nl))
+        qkwl=max(1.e-5,qk(nl))
+        qalpha(layer)=(1.d0/qmwl)
      &                  *4.d0*fbeta(layer)**2/(3.d0*falpha(layer)**2)+
-     &                (1.d0/qk(nl))
+     &                (1.d0/qkwl)
      &                  *(1.d0-(4.d0*fbeta(layer)**2/(3.d0*falpha(layer)**2)))
+        qalpha(layer)=max(1.e-10,qalpha(layer))
         qalpha(layer)=(1.d0/qalpha(layer))
       enddo
 c----------------------------------------------------------------------