Index: ../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f =================================================================== --- ../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 14834) +++ ../trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 14835) @@ -1,13 +1,15 @@ 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,16 +38,11 @@ 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,10 +55,12 @@ 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,11 +75,9 @@ 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