Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 14839) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 14840) @@ -5,12 +5,9 @@ parameter (npat = N3G) double precision yvalue(nhank),bcin(nhank) -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision wok, rpos -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision pset(7) -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision swok, asrpos, aswokm, distrad common /blockrad/ distrad -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: common /blockp/ pset common /blocki/ INTERN @@ -19,18 +16,17 @@ data zero /0.0d0/, one /1.0d0/, two /2.0d0/, three /3.0d0/, 1rescal/ 1.0d0/ - if(INTERN.ge.npat) INTERN = 0 - INTERN = INTERN + 1 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: iprate = irate + 1 -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: bath = dk / three -c SA :::: rpose should be normalized wrt lithosphere thickness = 100 km -c SA :: give r is normalized dist_rad :: r == dist_rad / h -c write(6,*) irate, distrad - r = distrad / 100.00 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c rpose should be normalized wrt lithosphere thickness +c give r is normalized dist_rad :: r == dist_rad / h +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + r = distrad / (pset(1) / 1.0d3) rpos = r ak = zero -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c form the yvalue's for the Simpson's rule formulas -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: do 425 ik = 1, nhank ak = ak + dk @@ -39,10 +35,9 @@ yvalue(ik) = bcin(ik) * rarg 425 continue -c correct to end point val. in Simp. -c Rule - yvalue(nhank) = bcin(nhank) * rarg / two -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c correct to end point val. in Simp. Rule +c yvalue(nhank) = bcin(nhank) * rarg / two c find the area under the curve using the Simpson's rule formulas -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: sumde = zero do 300 int = 1, nhank @@ -52,25 +47,6 @@ sumde = ( fide * yvalue(int) ) + sumde 300 continue -c wok(ir) = bath * sumde wok = bath * sumde -c 450 continue -c -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: -c SA close(33) -c SA -c call dvecpr(wok,nrv,'wok',79,0,0) -c SA write(6,71) ak,bcin(nhank) -c SA 71 format(' OJRULE ak last bcin'/1h ,1p2d16.8) -c SA write(6,72) r, bath -c SA 72 format(' OJRULE r bath'/1h ,1p2d16.8) -c SA write(6,70) fide,sumde -c SA 70 format(' OJRULE fide sumde'/1h ,1p2d16.8) -c -c double loop concluded and now for some single prec. IO: -c -c wolfc.f does w vs. r for 1 values of t -c rates.f does w dot vs. r for 1 value of t -c -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: go to (1000,2000), iprate 1000 call wolfc(rpos,wok) @@ -78,98 +54,59 @@ 2000 call rates(rpos,wok) 999 return -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: + end -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine wolfc(rpos,wok) parameter (N3G = 1) parameter (npat = N3G) -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision wok, rpos -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision pset(7) -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision swok,asrpos, aswokm -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: -c real swoki(80), sri(80), swoko(160), sro(160) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: common /blockp/ pset common /blocki/ INTERN common /blockn/ irate common /blocks/ aswokm,asrpos,swok -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c make single prec. and return to dimensional units. -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: hscale = sngl(pset(1)) hsckm = hscale / 1.0e3 -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: -c do 45 i = 1, nrv swok = hscale * sngl(wok) asrpos = hsckm * sngl(rpos) aswokm = swok -c 45 continue -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: -c -c Like Wolf (1985), do "inside" and "outside" curves seperately. -c -c do 70 i = 1,80 -c swoki(i) = swok(i) -c sri(i) = asrpos(i) -c 70 continue -c do 80 i = 81, nrv -c ic2 = i - 80 -c swoko(ic2) = swok(i) -c sro(ic2) = asrpos(i) -c 80 continue -c -c and the second "outside" plot -c -c call svecpr(swoki,80,'swoki in meters',79,0,0) -c call svecpr(swoko,160,'swoko in meters',79,0,0) -c call svecpr(sri,80,'sri in km',79,0,0) -c call svecpr(sro,160,'sro in km',79,0,0) + return end +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine rates(rpos,wok) parameter (N3G = 1) parameter (npat = N3G) -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision wok, rpos double precision pset(7) double precision swok, asrpos, aswokm -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: common /blockp/ pset common /blocki/ INTERN common /blockn/ irate common /blocks/ aswokm,asrpos,swok +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: data ngyo /050201/ data ngyii /050201/ data yearco /3.15576d7/ -c write(6,41) irate -c 41 format(' irate in rates'/1h 1p1i10) -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c make single prec. and return to dimensional units. -c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: hscale = sngl(pset(1)) hsckm = hscale / 1.0e3 -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: -c do 45 i = 1, nrv swok = (hscale * yearco * 1.0e3 * sngl(wok)) 1 * ( sngl(pset(4))/ sngl(pset(2)) ) asrpos = hsckm * sngl(rpos) aswokm = swok -c 45 continue -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: -c SA: write displacement (w) in file #39 -c SA: call svecpr(asrpos,nrv,'asrpos in rates - km',79,0,0) -c -c Like Wolf (1985), do "inside" and "outside" curves separately. -c -c INP = INTERN - 1 -c if(INP.gt.6) INP = 6 -c -c next call ploting . -c -c call dvecpr(pset,7,'pset in ojrule',79,0,0) -c call svecpr(swok,nrv,'swok in rates - mm per yr',79,0,0) + return end