Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f
===================================================================
--- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f	(revision 14840)
+++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f	(revision 14841)
@@ -8,15 +8,13 @@
      1decta(2),dyri1(nhank),dyri2(nhank),sna(nhank)
       double precision cinner(nhank),bcin(nhank)
-c     double precision trsl(Nrsl)
       integer maxk
-c      double precision aswokm, asrpos, swok
-c      common /blocks/ aswokm,asrpos,swok
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       common /blockp/ pset
       common /blockz/ zkp
       common /blockm/ dekay1,dekay2,amp0,amp1,amp2
-c     common /blockk/ trsl
       data yearco /3.15576d7/, pi /3.1415926535897932384d0/
       data g /9.832186d0/, four /4.d0/, two /2.0d0/,
      1 one /1.0d0/, zero/0.0d0/ , maxk/64/
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       twopi = two * pi
       r2 = pset(6)
@@ -26,7 +24,8 @@
       h  = pset(1)
       urat = u1/u2
-c SA: ==================================================
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+c  alphap is dimensionless disk radius
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       alphap = pset(7)/pset(1)
-c     alphap is dimensionless disk radius
       twoap = two * alphap
       rghm = ( r1 * g * h * alphap ) / (two * u2)
@@ -38,5 +37,4 @@
       dk = endk/dfloat(nhank)
 c
-c      if(icheck.eq.3) go to 999
       if(idisk.gt.1) go to 7001
       ak = zero
@@ -52,6 +50,7 @@
       e2 = dexp(zkp2)
       e4 = dexp(zkp4)
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       call freed(r2,u2,r1,u1,h,zkd,e1,e2,e4,b0,b1,a2,a1,a0,decay,amps)
-c
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       decta(1) = decay(1)/tmxyr
       decta(2) = decay(2)/tmxyr
@@ -59,9 +58,9 @@
       dyri2(ik) = decta(2)
       sna(ik) = pikn
-c                            Form vectors for full construction in pwise.f
-c                            & stot.f
-c                            Note that freed will produce decay spectra
-c                            defined as positive, ie. negative decay must
-c                            reinsert a minus sign.
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+c Form vectors for full construction in pwise.f and stot.f
+c Note that freed will produce decay spectra defined as positive, 
+c ie. negative decay must reinsert a minus sign.
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       dekay1(ik) = decay(1)
       dekay2(ik) = decay(2)
@@ -73,19 +72,5 @@
  7000 continue
  7001 continue
-c     call dvecpr(zksam,nhank,'k sampled like Wolf',79,0,0)
-c     call dvecpr(zksamp,nhank,'k - dimensionless',79,0,0)
-c     call dvecpr(amp0,nhank,'amp0 of k from freed',79,0,0)
-c     call dvecpr(amp1,nhank,'amp1 of k from freed',79,0,0)
-c     call dvecpr(amp2,nhank,'amp2 of k from freed',79,0,0)
-c     call dvecpr(dekay1,nhank,'dekay1 of k from freed',79,0,0)
-c     call dvecpr(dekay2,nhank,'dekay2 of k from freed',79,0,0)
-c
-c A call below perfomed plots of the two decay times in years
-c sup {-1} vs. wave number k
-c
-c     if(icheck.gt.2) go to 49
-c     call wplot(dyri1,dyri2,sna)
-c     go to 999
-c
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 c The following looped call sets up the free solution convolved with the
 c load function q hat.  Note that the returned vector set "cinner" is the
@@ -97,9 +82,6 @@
 c elliptical cross section.  Note loops 8500,8000 and 9500,9000 for the two
 c cases, respectively.
-c
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    49 go to (8499,9499), iedge
-c8499 tspan = trsl(irsl)
-c                                     trsl was nondimesionalized
-c                                     in the routine distme.f
  8499 do 8000 ik = 1, nhank
       xakap = zksamp(ik)*alphap
@@ -111,19 +93,13 @@
       bcin(ik) = cinner(ik) * dbesj1(xakap)
  8000 continue
-c      write(6,133) diku,xakap,pset(1),pset(7),alphap,urat
-c  133 format('      last diku  xakap in what0'/1h ,1p6d16.8)
-c      write(6,1) bcin(nhank),zksamp(nhank),fltng
-c    1 format('      last bcin  zksamp and fltng in what0'/1h ,1p3d16.8)
-c      write(6,2) cinner(nhank),cinner(1),qjadon
-c    2 format(' last and first cinner in what0 and qjadon'/1h ,1p3d16.8)
-c
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 c "ojrule.f" computes the inverse Hankel trasform with a simple
 c Simpson's rule.  The routine "ojrule" is buliding a set of solutions stored
 c in common "blocks" in r or "asrpos(nrv) ", and computed rate or displacement
 c for each of N3G disks in "aswokm(nrv,N3G)" . 
-c
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       call ojrule(dk,bcin)
       go to 999
-c9499 tspan = trsl(irsl)
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  9499 do 9000 ik = 1, nhank
       xakap = zksamp(ik)*alphap
@@ -137,8 +113,6 @@
      1 - dcos(xakap) )
  9000 continue
-c     write(6,1) bcin(nhank),zksamp(nhank)
       call ojrule(dk,bcin)
-c     call dvecpr(cinner,nhank,'cinner at time t',79,0,0)
-c     call dvecpr(bcin,nhank,'bcin at time t',79,0,0)
   999 return
+
       end
