Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 15139) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 15140) @@ -21,8 +21,4 @@ }; -struct blockn{ - int irate; -}; - struct blockrad{ double distrad; @@ -30,5 +26,6 @@ struct blocks{ - double aswokm; + double aswokm_w; + double aswokm_dwdt; }; @@ -39,5 +36,4 @@ extern struct blockp blockp_; extern struct blocko blocko_; - extern struct blockn blockn_; extern struct blockrad blockrad_; extern struct blocks blocks_; @@ -63,5 +59,4 @@ /*gia solution parameters: */ int iedge=0; - int irate=0; /*gia inputs: */ @@ -111,5 +106,4 @@ rho_ice=arguments->rho_ice; disk_id=arguments->idisk; - irate=arguments->irate; iedge=arguments->iedge; yts=arguments->yts; @@ -118,6 +112,4 @@ /*Modify inputs to match naruse code: */ - //from our model, irate comes in with values in [1,2], which maps into [0,1] in naruse: - irate=irate-1; Ntime=numtimes; Ntimm=Ntime-1; @@ -155,7 +147,4 @@ blockt_dmi=xNew(Ntimm); - /*irate: */ - blockn_.irate=irate; - /*radial distance of i-th element: */ blockrad_.distrad=ri/1000.0; // in km @@ -169,14 +158,8 @@ /*output solution: */ - if(irate==0){ - wi = blocks_.aswokm; - *pwi=wi; - *pdwidt=0; - } - else if (irate==1){ - *pwi=0; - dwidt = blocks_.aswokm; - *pdwidt=dwidt; - } + wi = blocks_.aswokm_w; + dwidt = blocks_.aswokm_dwdt; + *pwi=wi; + *pdwidt=dwidt; /*Free ressources: */ Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f (revision 15139) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f (revision 15140) @@ -2,9 +2,8 @@ 1,amps) implicit double precision (a-h,o-z) - double precision decay(2),amps(3) + double precision decay(2),amps(5) double precision ac0,ac1,ac2,ac3,ac4,ac5,ac6,ac7,ac8,ac9,ac10, 1ac11 common /blockz/ zkp - common /blockn/ irate data zero /0.0d0/, g /9.832186d0/ c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -215,7 +214,6 @@ amps(2) = decdif*( decay(2) * ( a1 - a2*decay(2) ) - a0 ) amps(3) = a2 - if(irate.eq.0) go to 999 - amps(1) = - decay(1) * amps(1) - amps(2) = - decay(2) * amps(2) + amps(4) = - decay(1) * amps(1) + amps(5) = - decay(2) * amps(2) go to 999 9990 write(6,998) idgen Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 15139) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 15140) @@ -1,21 +1,22 @@ - subroutine ojrule(dk,bcin) + subroutine ojrule(dk,bcin_w,bcin_dwdt) implicit double precision(a-h,o-z) parameter (nhank = 1024) - double precision yvalue(nhank),bcin(nhank) - double precision wok, rpos + double precision yvalue_w(nhank),yvalue_dwdt(nhank) + double precision bcin_w(nhank),bcin_dwdt(nhank) + double precision wok_w,wok_dwdt,rpos + double precision swok_w,swok_dwdt double precision pset(7) - double precision aswokm, distrad + double precision aswokm_w,aswokm_dwdt,distrad c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: common /blockrad/ distrad common /blockp/ pset - common /blockn/ irate - common /blocks/ aswokm + common /blocks/ aswokm_w,aswokm_dwdt data zero /0.0d0/, one /1.0d0/, two /2.0d0/, three /3.0d0/, 1rescal/ 1.0d0/ + data yearco /3.15576d7/ c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - iprate = irate + 1 bath = dk / three c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c rpose should be normalized wrt lithosphere thickness +c rpos should be normalized wrt lithosphere thickness c give r is normalized dist_rad :: r == dist_rad / h c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -30,5 +31,6 @@ rak = ak * r rarg = dbesj0( rak ) - yvalue(ik) = bcin(ik) * rarg + yvalue_w(ik) = bcin_w(ik) * rarg + yvalue_dwdt(ik) = bcin_dwdt(ik) * rarg 425 continue c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -37,65 +39,25 @@ c find the area under the curve using the Simpson's rule formulas c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - sumde = zero + sumde_w = zero + sumde_dwdt = zero do 300 int = 1, nhank intp1 = int + 1 ide = 2 + ( (-1)**intp1 + 1 ) fide = dfloat(ide) - sumde = ( fide * yvalue(int) ) + sumde + sumde_w = ( fide * yvalue_w(int) ) + sumde_w + sumde_dwdt = ( fide * yvalue_dwdt(int) ) + sumde_dwdt 300 continue - wok = bath * sumde -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - go to (1000,2000), iprate - 1000 call wolfc(rpos,wok) - go to 999 - 2000 call rates(rpos,wok) - 999 return - - end -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - subroutine wolfc(rpos,wok) - double precision wok, rpos - double precision pset(7) - double precision swok,asrpos, aswokm -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - common /blockp/ pset - common /blocks/ aswokm -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c make single prec. and return to dimensional units. + wok_w = bath * sumde_w + wok_dwdt = bath * sumde_dwdt + c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: hscale = sngl(pset(1)) hsckm = hscale / 1.0e3 - swok = hscale * sngl(wok) - asrpos = hsckm * sngl(rpos) - aswokm = swok + swok_w = hscale * sngl(wok_w) + aswokm_w = swok_w + swok_dwdt = (hscale * yearco * 1.0e3 * sngl(wok_dwdt)) + 1 * ( sngl(pset(4))/ sngl(pset(2)) ) + aswokm_dwdt = swok_dwdt + return - return end -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - subroutine rates(rpos,wok) - double precision wok, rpos - double precision pset(7) - double precision swok, asrpos, aswokm -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - common /blockp/ pset - common /blocks/ aswokm -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - data ngyo /050201/ - data ngyii /050201/ - data yearco /3.15576d7/ -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c make single prec. and return to dimensional units. -c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - hscale = sngl(pset(1)) - hsckm = hscale / 1.0e3 - swok = (hscale * yearco * 1.0e3 * sngl(wok)) - 1 * ( sngl(pset(4))/ sngl(pset(2)) ) - asrpos = hsckm * sngl(rpos) - aswokm = swok - - return - end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f (revision 15139) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f (revision 15140) @@ -1,3 +1,4 @@ - subroutine pwise(t,ta,tb,xi1,xi2,slope,ycept,decay,bhaq) + subroutine pwise(t,ta,tb,xi1,xi2,xi3,xi4,slope,ycept,decay, + 1bhaq_w,bhaq_dwdt) implicit double precision (a-h,o-z) double precision decay(2) @@ -19,14 +20,14 @@ eb2 = dexp(gbt2) c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c define xi1 term: +c define xit1 term: c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - xi1t =(ycept/decay(1)) * (eb1 - ea1) - + xit1 =(ycept/decay(1)) * (eb1 - ea1) - 1(slope/(decay(1)*decay(1))) * 2 ( (1.0d0 - tb*decay(1))*eb1 3 - (1.0d0 - ta*decay(1))*ea1 ) c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c define xi2 term: +c define xit2 term: c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - xi2t =(ycept/decay(2)) * (eb2 - ea2) - + xit2 =(ycept/decay(2)) * (eb2 - ea2) - 1(slope/(decay(2)*decay(2))) * 2 ( (1.0d0 - tb*decay(2))*eb2 @@ -36,5 +37,6 @@ c ABOVE IS THE NON-DEGENERATE CASE c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - bhaq = (xi1 * xi1t) + (xi2 * xi2t) + bhaq_w = (xi1 * xit1) + (xi2 * xit2) + bhaq_dwdt = (xi3 * xit1) + (xi4 * xit2) return end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f (revision 15139) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f (revision 15140) @@ -1,7 +1,6 @@ - subroutine qwise(t,ta,qjadon,xi0,xi1,xi2,slope,ycept,decay,bhaq) + subroutine qwise(t,ta,qjadon,xi0,xi1,xi2,xi3,xi4,slope,ycept, + 1decay,bhaq_w,bhaq_dwdt) implicit double precision (a-h,o-z) double precision decay(2) - common /blockn/ irate - iprate = irate + 1 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c This subroutine retrieves the convolution for the J-th linear piece-wise @@ -18,15 +17,16 @@ xg1 = xi1/(decay(1)*decay(1)) xg2 = xi2/(decay(2)*decay(2)) + xg3 = xi3/(decay(1)*decay(1)) + xg4 = xi4/(decay(2)*decay(2)) gb1 = decay(1)*ycept gb2 = decay(2)*ycept - go to ( 66, 67 ), iprate c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c define xi0t term: +c define xit0 term: c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - 66 xi0t = (xi0 + qjadon) * ( ( slope * t ) + ycept ) + xit0_w = (xi0 + qjadon) * ( ( slope * t ) + ycept ) c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c define xi1t term: +c define xit1 term: c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - xi1t = xg1 * ( + xit1_w = xg1 * ( 1 gb1 + slope * ( ( t * decay(1) ) - 1.0d0 ) 2 - ( gb1 + slope * ( ( ta * decay(1) ) - 1.0d0 )) @@ -34,12 +34,11 @@ 4 ) c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: -c define xi2t term: +c define xit2 term: c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - xi2t = xg2 * ( + xit2_w = xg2 * ( 1 gb2 + slope * ( ( t * decay(2) ) - 1.0d0 ) 2 - ( gb2 + slope * ( ( ta * decay(2) ) - 1.0d0 ) ) 3 * dexp( decay(2) * (ta - t) ) 4 ) - go to 90 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c And the rate equivalents: @@ -47,11 +46,11 @@ c having corrected in x1t, x2t pass). c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - 67 xi0t = (xi0 + qjadon) * slope - xi1t =-xg1 * ( + xit0_dwdt = (xi0 + qjadon) * slope + xit1_dwdt =-xg3 * ( 1 slope 2 + ( gb1 + slope * ( ( ta * decay(1) ) - 1.0d0 )) 3 * dexp( decay(1) * (ta - t) ) 4 ) - xi2t =-xg2 * ( + xit2_dwdt =-xg4 * ( 1 slope 2 + ( gb2 + slope * ( ( ta * decay(2) ) - 1.0d0 )) @@ -61,5 +60,6 @@ c add terms for the J-th (and final) interval contribution. c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - 90 bhaq = xi0t + xi1t + xi2t + bhaq_w = xit0_w + xit1_w + xit2_w + bhaq_dwdt = xit0_dwdt + xit1_dwdt + xit2_dwdt return end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f (revision 15139) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f (revision 15140) @@ -1,3 +1,4 @@ - subroutine stot(ikval,qjadon,fltng,Ntimp,Ntimm,time,bi,dmi) + subroutine stot(ikval,qjadon,fltng_w,fltng_dwdt,Ntimp,Ntimm, + 1time,bi,dmi) implicit double precision (a-h,o-z) integer Ntimp,Ntimm @@ -8,7 +9,7 @@ double precision time(Ntimp),bi(Ntimm),dmi(Ntimm) double precision dekay1(nhank),dekay2(nhank),amp0(nhank), - 1amp1(nhank),amp2(nhank) + 1amp1(nhank),amp2(nhank),amp3(nhank),amp4(nhank) c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - common /blockm/ dekay1,dekay2,amp0,amp1,amp2 + common /blockm/ dekay1,dekay2,amp0,amp1,amp2,amp3,amp4 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c This subroutine returns the inverse Laplace transform to the @@ -36,4 +37,6 @@ xi1 = amp1(ikval) xi2 = amp2(ikval) + xi3 = amp3(ikval) + xi4 = amp4(ikval) c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -41,5 +44,6 @@ c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: t = time(Ntimp) - sumb = 0.0d0 + sumb_w = 0.0d0 + sumb_dwdt = 0.0d0 ta = time(1) tb = time(2) @@ -56,5 +60,6 @@ c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: if( t . gt . ta . and . t . le . tb) go to 38 - call pwise(t,ta,tb,xi1,xi2,slope,ycept,decay,bhaq) + call pwise(t,ta,tb,xi1,xi2,xi3,xi4,slope,ycept,decay, + 1bhaq_w,bhaq_dwdt) go to 39 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -62,10 +67,13 @@ c evaluation still has to consider the load itself c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - 38 call qwise(t,ta,qjadon,xi0,xi1,xi2,slope,ycept,decay,bhaq) - 39 sumb = bhaq + sumb + 38 call qwise(t,ta,qjadon,xi0,xi1,xi2,xi3,xi4,slope,ycept,decay, + 1bhaq_w,bhaq_dwdt) + 39 sumb_w = bhaq_w + sumb_w + sumb_dwdt = bhaq_dwdt + sumb_dwdt ta = time(i + 1) tb = time(i + 2) 97 continue - fltng = sumb + fltng_w = sumb_w + fltng_dwdt = sumb_dwdt return end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f (revision 15139) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f (revision 15140) @@ -5,8 +5,10 @@ parameter (N = nhank/2) double precision dekay1(nhank),dekay2(nhank),amp0(nhank), - 1amp1(nhank),amp2(nhank),zksam(nhank),zksamp(nhank) - double precision decay(2),pset(7),amps(3), + 1amp1(nhank),amp2(nhank),amp3(nhank),amp4(nhank), + 1zksam(nhank),zksamp(nhank) + double precision decay(2),pset(7),amps(5), 1decta(2),dyri1(nhank),dyri2(nhank),sna(nhank) - double precision cinner(nhank),bcin(nhank) + double precision cinner_w(nhank),cinner_dwdt(nhank) + double precision bcin_w(nhank),bcin_dwdt(nhank) double precision time(Ntimp),bi(Ntimm),dmi(Ntimm) integer maxk @@ -14,5 +16,5 @@ common /blockp/ pset common /blockz/ zkp - common /blockm/ dekay1,dekay2,amp0,amp1,amp2 + common /blockm/ dekay1,dekay2,amp0,amp1,amp2,amp3,amp4 data yearco /3.15576d7/, pi /3.1415926535897932384d0/ data g /9.832186d0/, four /4.d0/, two /2.0d0/, @@ -69,4 +71,6 @@ amp1(ik) = amps(1) amp2(ik) = amps(2) + amp3(ik) = amps(4) + amp4(ik) = amps(5) zksam(ik) = zkd zksamp(ik) = zkp @@ -89,7 +93,9 @@ pref = diku / ( diku + rghm ) qjadon = one / ( four * zksamp(ik) * urat ) - call stot(ik,qjadon,fltng,Ntimp,Ntimm,time,bi,dmi) - cinner(ik) = - fltng * pref * twoap - bcin(ik) = cinner(ik) * dbesj1(xakap) + call stot(ik,qjadon,fltng_w,fltng_dwdt,Ntimp,Ntimm,time,bi,dmi) + cinner_w(ik) = - fltng_w * pref * twoap + cinner_dwdt(ik) = - fltng_dwdt * pref * twoap + bcin_w(ik) = cinner_w(ik) * dbesj1(xakap) + bcin_dwdt(ik) = cinner_dwdt(ik) * dbesj1(xakap) 8000 continue c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -99,5 +105,5 @@ c for each of N3G disks in "aswokm(nrv,N3G)" . c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - call ojrule(dk,bcin) + call ojrule(dk,bcin_w,bcin_dwdt) go to 999 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -108,10 +114,13 @@ pref = diku / ( diku + rghm ) qjadon = one / ( four * zksamp(ik) * urat ) - call stot(ik,qjadon,fltng,Ntimp,Ntimm,time,bi,dmi) - cinner(ik) = - fltng * pref * twoap - bcin(ik) = cinner(ik) * oxakap * ( dsin(xakap) * oxakap + call stot(ik,qjadon,fltng_w,fltng_dwdt,Ntimp,Ntimm,time,bi,dmi) + cinner_w(ik) = - fltng_w * pref * twoap + cinner_dwdt(ik) = - fltng_dwdt * pref * twoap + bcin_w(ik) = cinner_w(ik) * oxakap * ( dsin(xakap) * oxakap + 1 - dcos(xakap) ) + bcin_dwdt(ik) = cinner_dwdt(ik) * oxakap * ( dsin(xakap) * oxakap 1 - dcos(xakap) ) 9000 continue - call ojrule(dk,bcin) + call ojrule(dk,bcin_w,bcin_dwdt) 999 return