source: issm/oecreview/Archive/14312-15392/ISSM-14834-14835.diff

../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f

@@ -1 +1 @@
       subroutine distme(idisk,iedge)
       implicit double precision (a-h,o-y)
       parameter (N3G = 1)
-      parameter (Ntime = 2)
+      parameter (Ntime = 5)
+      parameter (Ntimm = Ntime-1)
       parameter (Nafter = 1)
       parameter (Ntimp = Ntime + Nafter)
       double precision pset(7)
-      double precision time(Ntimp),dmi(Ntime),bi(Ntime),dumbt(Ntimp)
+      double precision time(Ntimp),dmi(Ntimm),bi(Ntimm),dumbt(Ntimp)
       double precision hload(Ntime),qpat(Ntime),qt(Ntime)
       double precision zradii(N3G),zhload(Ntime),rhoi,distrad
+c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
       common /blockp/ pset
       common /blockrad/ distrad 
       common /blockt/ time,bi,dmi
@@ -36 +38 @@
 c now set up a piece-wise history: bi() = y-intercept 
 c                                 dmi() = slope 
 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
-      bi(1) = 0.0d0
-      dmi(1) = hload(1) / dumbt(1) 
       do 70 i = 2, Ntime
-      dmi(i) = ( hload(i) - hload(i-1) )/( dumbt(i)  - dumbt(i-1) )
-      bi(i) = hload(i-1) - ( dmi(i)*dumbt(i-1) )  
+      dmi(i-1) = ( hload(i) - hload(i-1) )/( dumbt(i)  - dumbt(i-1) )
+      bi(i-1) = hload(i-1) - ( dmi(i)*dumbt(i-1) )  
    70 continue
 c      write(6,*) zhload(1,1), zhload(1,2) 
-c      call dvecpr(hload,Ntime,'::::: hload @ distme.f :::::',79,0,0)
-c      call dvecpr(dmi,Ntime,'::::: load slope @ distme.f :::::',79,0,0)
-c      call dvecpr(bi,Ntime,'::::: load y-cept @ distme.f :::::',79,0,0)
 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 c With pset(6) in mks units, lets convert the piecewise linear formulas
 c for the time-dependent ice load heights to dimensionless values w.r.t. time.
@@ -58 +55 @@
       do 20 jt = 1, Nafter
       time(Ntime + jt) = ( dumbt(Ntime + jt) * yearco * 1.0d3 ) / tfact
    20 continue
-      do 75 ind = 1, Ntime 
+      do 75 ind = 1, Ntimm 
       dmi(ind) =  dmi(ind) / (( yearco * 1.0d3 ) / tfact )
-      time(ind) = ( dumbt(ind) * yearco * 1.0d3 ) / tfact
    75 continue
+      do 77 j = 1, Ntime 
+      time(j) = ( dumbt(j) * yearco * 1.0d3 ) / tfact
+   77 continue
 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 c create an incremental load in Pa and non-dimensionalized:
 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
@@ -76 +75 @@
 c integrals for the inverse Laplace transform and inverse Hankel transform
 c without further mutiplicative factors.
 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
-      bi(1) = 0.0d0
-      dmi(1) = qt(1) / time(1) 
       do 85 i = 2, Ntime
-      dmi(i) = ( qt(i) - qt(i-1) )/( time(i)  - time(i-1) )
-      bi(i) = qt(i-1) - ( dmi(i)*time(i-1) )  
+      dmi(i-1) = ( qt(i) - qt(i-1) )/( time(i)  - time(i-1) )
+      bi(i-1) = qt(i-1) - ( dmi(i-1)*time(i-1) )  
    85 continue
   999 return
       end