Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 14767) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 14768) @@ -37,15 +37,33 @@ struct blocky{ - double zhload[Ntime][N3G]; + double zhload[Ntime]; }; +struct blockn{ + int irate; +}; + +struct blockrad{ + double distrad; +}; + +struct blocks{ + double aswokm; + double asrpos; + double swok; +}; extern "C" { - int distme_( int* pidisk,int* piedge,int* pisolt); + int distme_( int* pidisk,int* piedge); + int what0_( int* pidisk,int* piedge); extern struct blockp blockp_; extern struct blockt blockt_; extern struct blocko blocko_; extern struct blocky blocky_; + extern struct blockn blockn_; + extern struct blockrad blockrad_; + extern struct blocks blocks_; } + /*}}}*/ @@ -58,5 +76,5 @@ int iedge=1; //c iedge ......... = 1 square-edged, = 2 elliptical disc x-section (see naruse.f) int idisk=1; // disk # - int isolt=1002; + int irate=0; // =0 fetch w solution (m) only; =1 dw/dt (mm/yr) only /*intermediary: */ @@ -84,6 +102,7 @@ /*some debug info: */ int disk_id; - /*}}}*/ +/*}}}*/ + /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/ ri=arguments->ri; @@ -104,16 +123,6 @@ /*}}}*/ - /*Documentation for pset coming from naruse.f: - c pset(1) ....... h (thickness of elastic lithosphere in km) - c pset(2) ....... viscosity of sub-layer in Pa*s - c pset(3) ....... elastic rigidity of lithosphere in dynes/cm^2 - c pset(4) ....... elastic rigidity of sub-lithosphere in dynes/cm^2 - c pset(5) ....... density of lithosphere in gm/cc - c pset(6) ....... density of sub-lithosphere in gm/cc - c pset(7) ....... =dset(1) - c dset(1) ....... disk radius in km + // Now, let's set pset from the data that we got in input to GiaDeflectionCorex: - Now, let's set pset from the data that we got in input to GiaDeflectionCorex: - */ blockp_.pset[0]=lithosphere_thickness; blockp_.pset[1]=mantle_viscosity; @@ -123,22 +132,39 @@ blockp_.pset[5]=mantle_density; blockp_.pset[6]=re; - blocko_.rhoi=rho_ice; +// blocko_.rhoi=rho_ice; // use this for non-benchmark experiments +// blockp_.pset[6]=800000.0; // for testing benchmark, should be in meters + blocko_.rhoi=1000.0; // use this for benchmark not rho_ice = 917 kg m^-3 /*loading history: */ - blocky_.zhload[0][0]=current_he; - blocky_.zhload[1][0]=current_he; + blocky_.zhload[0]=current_he; + blocky_.zhload[1]=current_he; +// blocky_.zhload[0]=3000; // for testing benchmark +// blocky_.zhload[1]=3000; // for testing benchmark /*times: */ - blockt_.time[0]=0.1; - blockt_.time[1]=2000; - blockt_.time[2]=1000; + blockt_.time[0]=1.e-4; // in kyr + blockt_.time[1]=2500e+0; // in kyr + blockt_.time[2]=2400e+0; // control this for benchmark experiments + /*irate: */ + blockn_.irate=irate; + + /*radial distance of i-th element: */ + blockrad_.distrad=ri/1000.0; // in km +// blockrad_.distrad=500.0; // for testing benchmark /*Call distme driver: */ - distme_(&idisk,&iedge,&isolt); + distme_(&idisk,&iedge); + + /*Call what0 driver: */ + what0_(&idisk,&iedge); - //printf("disk id: %i bi: [%g,%g] dmi: [%g,%g] \n",disk_id,blockt_.bi[0], blockt_.bi[1],blockt_.dmi[0],blockt_.dmi[1]); + /*this is the solution: */ + wi = blocks_.aswokm; /*allocate output pointer: */ *pwi=wi; + +// printf("wi: %g deflection: %g \n",wi,blocks_.aswokm); + } Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 14767) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 14768) @@ -1,3 +1,3 @@ - subroutine distme(idisk,iedge,isolt) + subroutine distme(idisk,iedge) implicit double precision (a-h,o-y) parameter (N3G = 1) @@ -8,9 +8,10 @@ double precision time(Ntimp),dmi(Ntime),bi(Ntime),dumbt(Ntimp) double precision hload(Ntime),qpat(Ntime),qt(Ntime) - real zradii(N3G),zhload(N3G,Ntime) + double precision zradii(N3G),zhload(Ntime),rhoi,distrad common /blockp/ pset + common /blockrad/ distrad common /blockt/ time,bi,dmi + common /blocky/ zhload common /blocko/ rhoi,hload - common /blocky/ zhload data g /9.832186d0/, yearco /3.15576d7/, eradm/6.371d6/ data dpi /3.1415926535897932d0/, dzero/0.0d0/ @@ -23,8 +24,12 @@ dumbt(k) = time(k) 776 continue +c write(6,*) time(1), time(2), time(3) +c write(6,*) zhload(1), zhload(2), distrad, rhoi +c write(6,*) pset(1), pset(2), pset(3), pset(4), pset(5), pset(6) +c write(6,*) pset(7) c call dvecpr(time,Ntime,'::::: time @ distme.f :::::',79,0,0) c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: do 39 itime = 1, Ntime - hload(itime) = dble( zhload(idisk,itime) ) + hload(itime) = dble( zhload(itime) ) 39 continue c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -38,4 +43,5 @@ bi(i) = 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) Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f (revision 14767) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f (revision 14768) @@ -6,5 +6,5 @@ 1ac11 common /blockz/ zkp - common /blockn/ ngx,ngy,irate + common /blockn/ irate data zero /0.0d0/, g /9.832186d0/ c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 14767) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 14768) @@ -1,6 +1,5 @@ - subroutine ojrule(dk,endr,bcin) + subroutine ojrule(dk,bcin) implicit double precision(a-h,o-z) parameter (nhank = 1024) - parameter (nrv = 91) parameter (N3G = 1) parameter (npat = N3G) @@ -11,10 +10,10 @@ double precision pset(7) c SA ::::::::::::::::::::::::::::::::::::::::::::::::: - real swok, asrpos, aswokm(npat), dist_rad - common /blockRad/ dist_rad + double precision swok, asrpos, aswokm, distrad + common /blockrad/ distrad c SA ::::::::::::::::::::::::::::::::::::::::::::::::: common /blockp/ pset common /blocki/ INTERN - common /blockn/ ngx,ngy,irate + common /blockn/ irate common /blocks/ aswokm,asrpos,swok data zero /0.0d0/, one /1.0d0/, two /2.0d0/, three /3.0d0/, @@ -22,25 +21,11 @@ if(INTERN.ge.npat) INTERN = 0 INTERN = INTERN + 1 -c write(6,9099) INTERN -c9099 format(' INTERN'/1h 1p1i10) iprate = irate + 1 -c -c Brute force integration for the inverse -c Hankel transform. (Actually this is Simpson's rule) -c -c SA: open an output file -c SA: open(unit=33,file='outputs/r',status='new') -c c SA ::::::::::::::::::::::::::::::::::::::::::::::::: -c here I just compute for given r -c dr = endr/dfloat(nrv) bath = dk / three -c r = zero -c do 450 ir = 1, nrv -c r = r + dr -c rpos(ir) = r c SA :::: rpose should be normalized wrt lithosphere thickness = 100 km c SA :: give r is normalized dist_rad :: r == dist_rad / h - r = dist_rad / 100.00 +c write(6,*) irate, distrad + r = distrad / 100.00 rpos = r ak = zero @@ -97,5 +82,4 @@ c subroutine wolfc(rpos,wok) - parameter (nrv = 91) parameter (N3G = 1) parameter (npat = N3G) @@ -105,10 +89,10 @@ double precision pset(7) c SA ::::::::::::::::::::::::::::::::::::::::::::::::: - real swok,asrpos, aswokm(npat) + double precision swok,asrpos, aswokm c SA ::::::::::::::::::::::::::::::::::::::::::::::::: c real swoki(80), sri(80), swoko(160), sro(160) common /blockp/ pset common /blocki/ INTERN - common /blockn/ ngx,ngy,irate + common /blockn/ irate common /blocks/ aswokm,asrpos,swok c @@ -121,5 +105,5 @@ swok = hscale * sngl(wok) asrpos = hsckm * sngl(rpos) - aswokm(INTERN) = swok + aswokm = swok c 45 continue c SA ::::::::::::::::::::::::::::::::::::::::::::::::: @@ -147,17 +131,14 @@ subroutine rates(rpos,wok) - parameter (nrv = 91) parameter (N3G = 1) parameter (npat = N3G) c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision wok, rpos -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: double precision pset(7) -c SA ::::::::::::::::::::::::::::::::::::::::::::::::: - real swok, asrpos, aswokm(npat) + double precision swok, asrpos, aswokm c SA ::::::::::::::::::::::::::::::::::::::::::::::::: common /blockp/ pset common /blocki/ INTERN - common /blockn/ ngx,ngy,irate + common /blockn/ irate common /blocks/ aswokm,asrpos,swok data ngyo /050201/ @@ -176,5 +157,5 @@ 1 * ( sngl(pset(4))/ sngl(pset(2)) ) asrpos = hsckm * sngl(rpos) - aswokm(INTERN) = swok + aswokm = swok c 45 continue c SA ::::::::::::::::::::::::::::::::::::::::::::::::: @@ -184,11 +165,11 @@ c Like Wolf (1985), do "inside" and "outside" curves separately. c - INP = INTERN - 1 - if(INP.gt.6) INP = 6 +c INP = INTERN - 1 +c if(INP.gt.6) INP = 6 c c next call ploting . c - call dvecpr(pset,7,'pset in ojrule',79,0,0) - call svecpr(swok,nrv,'swok in rates - mm per yr',79,0,0) +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 Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f (revision 14767) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f (revision 14768) @@ -2,5 +2,5 @@ implicit double precision (a-h,o-z) double precision decay(2) - common /blockn/ ngx,ngy,irate + common /blockn/ irate iprate = irate + 1 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f (revision 14767) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f (revision 14768) @@ -3,5 +3,5 @@ parameter (Nafter = 1) parameter (Ntime = 2) - parameter (nhank = 2048) + parameter (nhank = 1024) parameter (Ntimp = Ntime + Nafter) double precision decay(2) Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f (revision 14767) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f (revision 14768) @@ -1,5 +1,5 @@ - subroutine what0(idisk,icheck,iedge,endr) + subroutine what0(idisk,iedge) implicit double precision (a-h,o-z) - parameter (nhank = 2048) + parameter (nhank = 1024) parameter (N = nhank/2) double precision dekay1(nhank),dekay2(nhank),amp0(nhank), @@ -10,4 +10,6 @@ c double precision trsl(Nrsl) integer maxk +c double precision aswokm, asrpos, swok +c common /blocks/ aswokm,asrpos,swok common /blockp/ pset common /blockz/ zkp @@ -36,5 +38,5 @@ dk = endk/dfloat(nhank) c - if(icheck.eq.3) go to 999 +c if(icheck.eq.3) go to 999 if(idisk.gt.1) go to 7001 ak = zero @@ -121,5 +123,5 @@ c for each of N3G disks in "aswokm(nrv,N3G)" . c - call ojrule(dk,endr,bcin) + call ojrule(dk,bcin) go to 999 c9499 tspan = trsl(irsl) @@ -136,5 +138,5 @@ 9000 continue c write(6,1) bcin(nhank),zksamp(nhank) - call ojrule(dk,endr,bcin) + 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)