Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (revision 21577) @@ -0,0 +1,171 @@ +/*!\file GiaDeflectionCorex + * \brief: GIA solution from Erik Ivins. + * Compute deflection wi from a single disk of radius re, load history hes for + * numtimes time steps. + */ + +#include "./GiaDeflectionCorex.h" + +#include "../../classes/classes.h" +#include "../../shared/shared.h" +#include "../../toolkits/toolkits.h" +#include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h" + +/*External blocks: {{{*/ +struct blockp{ + double pset[7]; +}; + +struct blocko{ + double rhoi; +}; + +struct blockrad{ + double distrad; +}; + +struct blocks{ + double aswokm_w; + double aswokm_dwdt; +}; + +extern "C" { + int distme_(int* pNtime,int* pNtimp,int* pNtimm,double* time,double* bi,double* dmi,double* zhload); + + int what0_(int* piedge,int* pNtimp,int* pNtimm,double* time,double* bi,double* dmi); + extern struct blockp blockp_; + extern struct blocko blocko_; + extern struct blockrad blockrad_; + extern struct blocks blocks_; +} + +/*}}}*/ + +void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){ + + /*intermediary: */ + int i; + + /*output: */ + IssmDouble wi=0; + IssmDouble dwidt=0; + + /*inputs: {{{*/ + /*constant: */ + int idisk=1; // disk # + IssmDouble yts; + + /*coming from the model structure, runtime configurable:*/ + /*gia solution parameters: */ + int iedge=0; + + /*gia inputs: */ + IssmDouble ri; //radial distance from center of disk to vertex i + IssmDouble re; //radius of disk + IssmDouble* hes; //loading history (in ice thickness) + IssmDouble* times; //loading history times + int numtimes; //loading history length + IssmDouble currenttime; + int Ntime; // number of times with load history + int Ntimm; // Ntime-1 : for slope/y-cept of load segments + int Ntimp; // Ntime+1 : for evaluation time + IssmDouble* blockt_time=NULL; + IssmDouble* blockt_bi=NULL; + IssmDouble* blockt_dmi=NULL; + IssmDouble* blocky_zhload=NULL; + + /*gia material parameters: */ + IssmDouble lithosphere_shear_modulus; + IssmDouble lithosphere_density; + IssmDouble mantle_shear_modulus; + IssmDouble mantle_viscosity; + IssmDouble mantle_density; + IssmDouble lithosphere_thickness; + + /*ice properties: */ + IssmDouble rho_ice; + + /*some debug info: */ + int disk_id; + +/*}}}*/ + + /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/ + ri = arguments->ri; + re = arguments->re; + hes = arguments->hes; + times = arguments->times; + numtimes = arguments->numtimes; + currenttime = arguments->currenttime; + lithosphere_shear_modulus = arguments->lithosphere_shear_modulus; + lithosphere_density = arguments->lithosphere_density; + mantle_shear_modulus = arguments->mantle_shear_modulus; + mantle_viscosity = arguments->mantle_viscosity; + mantle_density = arguments->mantle_density; + lithosphere_thickness = arguments->lithosphere_thickness; + rho_ice = arguments->rho_ice; + disk_id = arguments->idisk; + iedge = arguments->iedge; + yts = arguments->yts; + + /*}}}*/ + + /*Modify inputs to match naruse code: */ + Ntime=numtimes; + Ntimm=Ntime-1; + Ntimp=Ntime+1; + + /*Prepare block inputs for fortran distme and what0 routines of the naruse code: {{{*/ + /*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; + blockp_.pset[2]=lithosphere_shear_modulus; + blockp_.pset[3]=mantle_shear_modulus; + blockp_.pset[4]=lithosphere_density; + blockp_.pset[5]=mantle_density; + blockp_.pset[6]=re; + blocko_.rhoi=rho_ice; + + /*loading history: */ + blocky_zhload=xNew(Ntime); + for(i=0;i(Ntimp); + for (i=0;i(Ntimm); + + /*dmi: */ + blockt_dmi=xNew(Ntimm); + + /*radial distance of i-th element: */ + blockrad_.distrad=ri/1000.0; // in km + /*}}}*/ + + /*Call distme driver: */ + distme_(&Ntime,&Ntimp,&Ntimm,blockt_time,blockt_bi,blockt_dmi,blocky_zhload); + + /*Call what0 driver: */ + what0_(&iedge,&Ntimp,&Ntimm,blockt_time,blockt_bi,blockt_dmi); + + /*output solution: */ + wi = blocks_.aswokm_w; + dwidt = blocks_.aswokm_dwdt; + *pwi=wi; + *pdwidt=dwidt; + + /*Free ressources: */ + xDelete(blockt_time); + xDelete(blockt_bi); + xDelete(blockt_dmi); + xDelete(blocky_zhload); + +} Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (revision 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (revision 21577) @@ -0,0 +1,13 @@ +/*!\file: GiaDeflectionCorex.h + * \brief header file for ... + */ + +#ifndef _GIADEFLECTIONCOREX_H +#define _GIADEFLECTIONCOREX_H + +#include "../../classes/classes.h" + +/* local prototypes: */ +void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments); + +#endif /* _GIADEFLECTIONCOREX_H */ Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 21577) @@ -0,0 +1,73 @@ + subroutine distme(Ntime,Ntimp,Ntimm,time,bi,dmi,zhload) + implicit double precision (a-h,o-y) + integer Ntime,Ntimp,Ntimm + parameter (Nafter=1) + double precision pset(7) + double precision time(Ntimp),dmi(Ntimm),bi(Ntimm),dumbt(Ntimp) + double precision hload(Ntime),qpat(Ntime),qt(Ntime) + double precision zhload(Ntime),rhoi,distrad +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + common /blockp/ pset + common /blockrad/ distrad + common /blocko/ rhoi + data g /9.832186d0/, yearco /3.15576d7/, eradm/6.371d6/ + data dpi /3.1415926535897932d0/, dzero/0.0d0/ +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c The units of time(Ntimp) are ka and the height of the load in meters. +c The slope, then for example, is in units of meters per ka. +c Note that "dumbt( )" is designed to perserve the initial "time( )" variable. +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 776 k = 1, Ntimp + dumbt(k) = time(k) + 776 continue +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 39 itime = 1, Ntime + hload(itime) = dble( zhload(itime) ) + 39 continue +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c now set up a piece-wise history: bi() = y-intercept +c dmi() = slope +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 70 i = 2, Ntime + 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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +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. +c (tfact is in seconds) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + tfact = pset(2)/pset(4) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c get all times as dimensionless +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 20 jt = 1, Nafter + time(Ntime + jt) = ( dumbt(Ntime + jt) * yearco * 1.0d3 ) / tfact + 20 continue + do 75 ind = 1, Ntimm + dmi(ind) = dmi(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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 80 iq = 1, Ntime + qpat(iq) = hload(iq)*rhoi*g + qt(iq) = qpat(iq) / pset(4) + 80 continue +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c As the final step in this routine, create a dimensionless stress load from +c qp. Here we'll use bi( ) and dmi( ) vectors with dimensionless time. Then +c qp (and it's piece-wise decomposition) is ready for the direct dimensionless +c integrals for the inverse Laplace transform and inverse Hankel transform +c without further mutiplicative factors. +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 85 i = 2, Ntime + 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 Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f (revision 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f (revision 21577) @@ -0,0 +1,222 @@ + subroutine freed(r2,u2,r1,u1,h,zk,e1,e2,e4,b0,b1,a2,a1,a0,decay + 1,amps) + implicit double precision (a-h,o-z) + double precision decay(2),amps(5) + double precision ac0,ac1,ac2,ac3,ac4,ac5,ac6,ac7,ac8,ac9,ac10, + 1ac11 + common /blockz/ zkp + data zero /0.0d0/, g /9.832186d0/ +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c Given the inputs to this subroutine(r2 through zk in the call +c statement above), the outputs are coefficients of "s" that +c are crucial to the Laplace transform inversion. From b0 and b1 +c we can compute the decay poles (or eigenvalues). +c +c This is NOT true in our case, though. - SA +c NOTE IN THE CODE THAT A CALL TO THIS SUBROUTINE NEED NOT +c BE MADE AT EACH TIME STEP --- BUT WILL HAVE TO BE CALLED +c IN THE NUMERICAL INTEGRATION FOR COMPUTING THE INVERSE HANKEL +c TRANSFORM **** +c +c Each term should be returned as dimensionless +c h => length u2 => stress taumx2 => time +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + zkp2 = zkp*zkp + ur = u1/u2 + ghu2 = (g*h) / u2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac0 dimensional units are stress times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac0 = 4.0d0*ur*zkp2*( 1.0d0 + e4 + + 1 2.0d0*e2*(1.0d0 + 2.0d0*zkp2) ) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac1 dimensional units are stress times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac1 = 2.0d0*r1*ghu2*zkp*(1.0d0 - e4 + 4.0d0*zkp*e2) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac2 dimensional units are stress^2 times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac2 = 8.0d0*ur*ur*(-1.0d0 + e1)* + 1 (1.0d0 + e1)*(1.0d0 + e2)*zkp2 +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac3 dimensional units are stress^2 times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac3 = + 1 2.0d0*zkp*ghu2*ur*((r1 + r2)*(1.0d0 + e4) + + 2 2.0d0*(r2 - r1)*e2*( 1.0d0 + 2.0d0*zkp2 )) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac4 dimensional units are stress^2 times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac4 = ghu2*ghu2*r1*(r2 - r1)* + 1 (1.0d0 - e4 + 4.0d0*zkp*e2) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac5 dimensional units are stress^3 times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac5 = + 1 4.0d0*zkp2*ur*ur*ur*(1.0d0 - e2 - 2.0d0*e1*zkp)* + 2 (1.0d0 - e2 + 2.0d0*e1*zkp) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac6 dimensional units are stress^3 times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac6 = + 1 2.0d0*zkp*ur*ur*(1.0d0 - e4 + 2 - 4.0d0*e2*zkp)*ghu2*r2 +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac7 dimensional units are stress^3 times l^-2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac7 = + 1 ur * ( ( (1.0d0 - e1)*(1.0d0 + e1) )**2)*r1*(r2 - r1) + 2 * ( ghu2*ghu2 ) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac8 dimensional units are stress^0 times l^-1 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac8 = -2.0d0*zkp*(1.0d0 + e2*(1.0d0 + 2.0d0*zkp*(1.0d0 + zkp))) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac9 dimensional units are stress^1 times l^-1 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac9 = + 1 ( 4.0d0*zkp*u1 - + 2 g*h*(r2 - r1)*(1.0d0 + e2*(1.0d0 + 2.0d0*zkp*(1.0d0 + zkp))) + 3 ) / u2 +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac10 dimensional units are stress^2 times l^-1 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac10 = + 1 -2.0d0*zkp*ur*ur*( 1.0d0 - e2 + 2 - 2.0d0*zkp*e2*(1.0d0 + zkp) ) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c ac11 dimensional units are stress^2 times l^-1 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ac11 = + 1 ghu2*ur*(r2 - r1)*(1.0d0 - e2*(1.0d0 + 2.0d0*zkp)) +c DIMESIONLESS +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c Other functions may be found in file "apcw0.record" +c (Nov. 9 1996) +c The following is a Mathematica version of the isolation of the +c coefficeints of the L transform variable s in the denominator. +c Here is where the set-up is performed to obtain the "free decay" +c times (with the Hankel transform variable "zk" embedded. Note that +c some greater efficency could be achieved by further simplifying the +c combinations of "acn" functions which are now a series of function +c subroutines in the fortran code. The corresponding Mathematica +c session is "twolayer.Linversion" dated Nov. 23, 1996. +c +c In[59]:= +c Together[%] +c Out[59]= +c ac2 + ac3 - ac4 + 2 ac5 - 2 ac6 + 2 ac7 +c --------------------------------------------- +c ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7 +c In[61]:= +c Simplify[Coefficient[els,s^2]] +c Out[61]= +c 1 +c In[65]:= +c eslnos = +c ac5/(ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7) - +c ac6/(ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7) + +c +c ac7/(ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7) +c +c Out[65]= +c ac5 +c --------------------------------------------- - +c ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7 +c +c ac6 +c --------------------------------------------- + +c ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7 +c +c ac7 +c --------------------------------------------- +c ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7 +c In[66]:= +c Simplify[%] +c Out[66]= +c ac5 - ac6 + ac7 +c --------------------------------------------- +c ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7 +c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c a common denominator factor is: bc +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + bc = ( ac0 - ac1 + + 1 ac2 + ac3 + 2 - ac4 + ac5 - + 3 ac6 + ac7 ) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c b1: Denominator coefficent of s: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + b1 = + 1 ( ac2 + ac3 + 2 - ac4 + ( 2.0d0 * ac5 ) + 3 - ( 2.0d0 * ac6 ) + 4 + ( 2.0d0 * ac7 ) ) / bc +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c b0: Denominator coefficent of s^0: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + b0 = + 1 ( ac5 - ac6 + + 2 ac7 ) / bc +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c the eigenvaules are just the solution of the quadratic in s: +c so return as "decay" +c *** Note that the decay times are defined as positive *** +c if a negative inverse decay time is returned there is an error! +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + sb1 = b1*b1 + fb0 = 4.0d0*b0 + diff = sb1 - fb0 + if(diff.le.zero) go to 25 + rs = dsqrt( diff ) + decay(1) = -( - b1 - rs ) / 2.0d0 + decay(2) = -( - b1 + rs ) / 2.0d0 + go to 26 + 25 idgen = 100 + go to 9990 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c and for the numerator part of the quadratic s dependence +c the Mathematica session is: +c +c Out[14]= +c 2 +c ac10 + ac11 + (2 ac10 + 2 ac11 + ac9) s + (ac10 + ac11 + ac8 + ac9) s +c---------------------------------------------------------------------- +c ac0 - ac1 + ac2 + ac3 - ac4 + ac5 - ac6 + ac7 +c +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 26 a0 = ( ac10 + ac11 ) / bc +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + a1 = ( 2.0d0*( ac10 + ac11 ) + 1 + ac9 ) / bc +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + a2 = (ac10 + ac11 + 1 + ac8 + ac9) / bc +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c and the following terms are the amplitudes of the inverse Laplace +c transform solution for the non-q part. (See the boxed equation on +c page 4 of the Nov. 23 1996 notes.) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + decdif = 1.0d0/(decay(2) - decay(1)) + amps(1) = -decdif*( decay(1) * ( a1 - a2*decay(1) ) - a0 ) + amps(2) = decdif*( decay(2) * ( a1 - a2*decay(2) ) - a0 ) + amps(3) = a2 + amps(4) = - decay(1) * amps(1) + amps(5) = - decay(2) * amps(2) + go to 999 + 9990 write(6,998) idgen + 998 format(' idgen val ** fatal error ** degenerate e.v.'/1h ,1p,1i12) + 999 return + end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 21577) @@ -0,0 +1,63 @@ + subroutine ojrule(dk,bcin_w,bcin_dwdt) + implicit double precision(a-h,o-z) + parameter (nhank = 1024) + 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_w,aswokm_dwdt,distrad +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + common /blockrad/ distrad + common /blockp/ pset + 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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + bath = dk / three +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c rpos 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 form the yvalue's for the Simpson's rule formulas +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 425 ik = 1, nhank + ak = ak + dk + rak = ak * r + rarg = dbesj0( rak ) + yvalue_w(ik) = bcin_w(ik) * rarg + yvalue_dwdt(ik) = bcin_dwdt(ik) * rarg + 425 continue +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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + sumde_w = zero + sumde_dwdt = zero + do 300 int = 1, nhank + intp1 = int + 1 + ide = 2 + ( (-1)**intp1 + 1 ) + fide = dfloat(ide) + sumde_w = ( fide * yvalue_w(int) ) + sumde_w + sumde_dwdt = ( fide * yvalue_dwdt(int) ) + sumde_dwdt + 300 continue + wok_w = bath * sumde_w + wok_dwdt = bath * sumde_dwdt + +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + hscale = sngl(pset(1)) + hsckm = hscale / 1.0e3 + 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 + + end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f (revision 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f (revision 21577) @@ -0,0 +1,42 @@ + 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) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c This subroutine retrieves the convolution for the ith linear piece-wise +c q hat function (the load shape or Bessel function part having +c been removed) with the free-decay solution. (see notes of +c 12-31-96 "Convolution in time"). The convolution is returned as "bhaq". +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + tbt = tb - t + tat = ta - t + gat1 = tat * decay(1) + gat2 = tat * decay(2) + gbt1 = tbt * decay(1) + gbt2 = tbt * decay(2) + ea1 = dexp(gat1) + ea2 = dexp(gat2) + eb1 = dexp(gbt1) + eb2 = dexp(gbt2) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c define xit1 term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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 xit2 term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + xit2 =(ycept/decay(2)) * (eb2 - ea2) - + 1(slope/(decay(2)*decay(2))) * + 2 ( (1.0d0 - tb*decay(2))*eb2 + 3 - (1.0d0 - ta*decay(2))*ea2 ) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c add terms for the i-th interval contribution. +c ABOVE IS THE NON-DEGENERATE CASE +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f (revision 21577) @@ -0,0 +1,65 @@ + 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) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c This subroutine retrieves the convolution for the J-th linear piece-wise +c q hat function (the load shape or Bessel function part having been removed) +c with the free-decay solution. (see notes of 3-27-97 "convo.ice" Mathematica +c session). The convolution is returned as "bhaq". +c +c THIS ROUTINE REPLACES pwise.f ONLY FOR t < time(Ntime) * +c (such that the load is still in place at time t). * +c +c Note irate = 1 case has to be applied to the linear term only (freed.f applies +c this correction to exponential terms) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c define xit0 term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + xit0_w = (xi0 + qjadon) * ( ( slope * t ) + ycept ) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c define xit1 term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + xit1_w = xg1 * ( + 1 gb1 + slope * ( ( t * decay(1) ) - 1.0d0 ) + 2 - ( gb1 + slope * ( ( ta * decay(1) ) - 1.0d0 )) + 3 * dexp( decay(1) * (ta - t) ) + 4 ) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c define xit2 term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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 ) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c And the rate equivalents: +c (sign switch due to freed.f already +c having corrected in x1t, x2t pass). +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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 ) + xit2_dwdt =-xg4 * ( + 1 slope + 2 + ( gb2 + slope * ( ( ta * decay(2) ) - 1.0d0 )) + 3 * dexp( decay(2) * (ta - t) ) + 4 ) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c add terms for the J-th (and final) interval contribution. +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f (revision 21577) @@ -0,0 +1,79 @@ + subroutine stot(ikval,qjadon,fltng_w,fltng_dwdt,Ntimp,Ntimm, + 1time,bi,dmi) + implicit double precision (a-h,o-z) + integer Ntimp,Ntimm + parameter (Nafter = 1) + parameter (nhank = 1024) + double precision decay(2) + double precision pset(7) + double precision time(Ntimp),bi(Ntimm),dmi(Ntimm) + double precision dekay1(nhank),dekay2(nhank),amp0(nhank), + 1amp1(nhank),amp2(nhank),amp3(nhank),amp4(nhank) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + common /blockm/ dekay1,dekay2,amp0,amp1,amp2,amp3,amp4 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c This subroutine returns the inverse Laplace transform to the +c time-domain for the vertical displacement at time t for Hankel wavenumber +c ikval. (In general this routine needs to be called nhank times). +c The main derivation uses the Faltung theorem of Laplace transforms. +c (1-1-97) NEW CASE OF 3-27-97 IS FOR t(Ntime + i) < t(Ntime) OR IN OTHER +c WORDS, THE LOAD STILL IN PLACE AT t. OPTION CALL to qwise.f +c PERFORMS THIS. +c +c A theory for the degenerate case was worked out but has been removed +c as an option from this code. +c +c Definition of tspan: nondimensional time span backwards form present +c when this routine is first called the dimensional +c equivalent might be say tspan = 12 ka, then 11 and +c then finally tspan = 0. +c Additional note for r.s.l calculations: the routines qwise and pwise +c are identical to the previous case for computations of present-day only +c vertical deformation field. +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + decay(1) = dekay1(ikval) + decay(2) = dekay2(ikval) + xi0 = amp0(ikval) + xi1 = amp1(ikval) + xi2 = amp2(ikval) + xi3 = amp3(ikval) + xi4 = amp4(ikval) + +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c note that tspan must be updated in the calling routine "what0.f" +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + t = time(Ntimp) + sumb_w = 0.0d0 + sumb_dwdt = 0.0d0 + ta = time(1) + tb = time(2) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 97 i = 1,Ntimm +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c note that this "if" prevents adding load +c segments of "future" times when computing +c an r.s.l. history (10-06-98). +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + if(t.lt.ta) go to 97 + slope=dmi(i) + ycept=bi(i) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + if( t . gt . ta . and . t . le . tb) go to 38 + call pwise(t,ta,tb,xi1,xi2,xi3,xi4,slope,ycept,decay, + 1bhaq_w,bhaq_dwdt) + go to 39 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c Note that qwise is employed only for the J-th Q hat term when t for +c evaluation still has to consider the load itself +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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_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 21577) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f (revision 21577) @@ -0,0 +1,127 @@ + subroutine what0(iedge,Ntimp,Ntimm,time,bi,dmi) + implicit double precision (a-h,o-z) + integer Ntimp,Ntimm + parameter (nhank = 1024) + parameter (N = nhank/2) + double precision dekay1(nhank),dekay2(nhank),amp0(nhank), + 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_w(nhank),cinner_dwdt(nhank) + double precision bcin_w(nhank),bcin_dwdt(nhank) + double precision time(Ntimp),bi(Ntimm),dmi(Ntimm) + integer maxk +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + common /blockp/ pset + common /blockz/ zkp + 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/, + 1 one /1.0d0/, zero/0.0d0/ , maxk/64/ +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + twopi = two * pi + r2 = pset(6) + u2 = pset(4) + r1 = pset(5) + u1 = pset(3) + h = pset(1) + urat = u1/u2 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c alphap is dimensionless disk radius +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + alphap = pset(7)/pset(1) + twoap = two * alphap + rghm = ( r1 * g * h * alphap ) / (two * u2) + taumx = pset(2)/pset(4) + tmxyr = taumx / yearco +c + dfac = dfloat(nhank)/dfloat(maxk) + endk = dfloat(nhank)/dfac + dk = endk/dfloat(nhank) +c + ak = zero + do 7000 ik = 1,nhank + ak = ak + dk + zkp = ak + pikn = (6.371d6 * zkp) / h + zkd = pikn / 6.371d6 +c + zkp2 = 2.0d0 * zkp + zkp4 = 4.0d0 * zkp + e1 = dexp(zkp) + e2 = dexp(zkp2) + e4 = dexp(zkp4) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + call freed(r2,u2,r1,u1,h,zkd,e1,e2,e4,b0,b1,a2,a1,a0,decay,amps) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + decta(1) = decay(1)/tmxyr + decta(2) = decay(2)/tmxyr + dyri1(ik) = decta(1) + dyri2(ik) = decta(2) + sna(ik) = pikn +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c Form vectors for full construction in pwise.f and stot.f +c Note that freed will produce decay spectra defined as positive, +c ie. negative decay must reinsert a minus sign. +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + dekay1(ik) = decay(1) + dekay2(ik) = decay(2) + amp0(ik) = amps(3) + amp1(ik) = amps(1) + amp2(ik) = amps(2) + amp3(ik) = amps(4) + amp4(ik) = amps(5) + zksam(ik) = zkd + zksamp(ik) = zkp + 7000 continue +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c The following looped call sets up the free solution convolved with the +c load function q hat. Note that the returned vector set "cinner" is the +c inner-most part of the arguement of the inverse Hankel trans. integral. +c It is time-dependent and the loop is for the k-dependancy. The time for +c calculation is given in the vector "time(Ntimp)" in the routine stot.f that is +c called below. Note that the sign on cinner(ik) below is for a load directed +c downward. ** For iedge = 1 assume sq. edge load and for iedge = 2 assume an +c elliptical cross section. Note loops 8500,8000 and 9500,9000 for the two +c cases, respectively. +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 49 go to (8499,9499), iedge + 8499 do 8000 ik = 1, nhank + xakap = zksamp(ik)*alphap + diku = xakap * urat + pref = diku / ( diku + rghm ) + qjadon = one / ( four * zksamp(ik) * urat ) + 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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c "ojrule.f" computes the inverse Hankel trasform with a simple +c Simpson's rule. The routine "ojrule" is buliding a set of solutions stored +c in common "blocks" in r or "asrpos(nrv) ", and computed rate or displacement +c for each of N3G disks in "aswokm(nrv,N3G)" . +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + call ojrule(dk,bcin_w,bcin_dwdt) + go to 999 +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 9499 do 9000 ik = 1, nhank + xakap = zksamp(ik)*alphap + oxakap = one/xakap + diku = xakap * urat + pref = diku / ( diku + rghm ) + qjadon = one / ( four * zksamp(ik) * urat ) + 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_w,bcin_dwdt) + 999 return + + end Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21576) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21577) @@ -77,5 +77,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.transient.isthermal",TransientIsthermalEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgroundingline",TransientIsgroundinglineEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgiaivins",TransientIsgiaivinsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgia",TransientIsgiaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.isesa",TransientIsesaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.transient.isdamageevolution",TransientIsdamageevolutionEnum)); @@ -85,5 +85,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.transient.iscoupler",TransientIscouplerEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.giaivins.cross_section_shape",GiaIvinsCrossSectionShapeEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum)); /*For stress balance only*/ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 21576) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 21577) @@ -29,5 +29,5 @@ #include "./GetVectorFromInputsx/GetVectorFromInputsx.h" #include "./GetVectorFromControlInputsx/GetVectorFromControlInputsx.h" -#include "./GiaIvinsDeflectionCorex/GiaIvinsDeflectionCorex.h" +#include "./GiaDeflectionCorex/GiaDeflectionCorex.h" #include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h" #include "./SetActiveNodesLSMx/SetActiveNodesLSMx.h"