Index: /issm/trunk-jpl/configs/config-macosx64-larour-gia.sh =================================================================== --- /issm/trunk-jpl/configs/config-macosx64-larour-gia.sh (revision 14739) +++ /issm/trunk-jpl/configs/config-macosx64-larour-gia.sh (revision 14739) @@ -0,0 +1,31 @@ +#!/bin/sh + +#petsc 3.2 +#mpich 1.4 + +pythonversion=2.7 + +./configure \ + --prefix=$ISSM_DIR \ + --with-matlab-dir=$MATLAB_DIR \ + --with-python-dir=$ISSM_DIR/externalpackages/python/install\ + --with-python-numpy-dir=$ISSM_DIR/externalpackages/python/install/Library/Frameworks/Python.framework/Versions/$pythonversion/lib/python$pythonversion/site-packages/numpy/core/include/numpy\ + --with-triangle-dir=$ISSM_DIR/externalpackages/triangle/install \ + --with-metis-dir=$ISSM_DIR/externalpackages/metis/install \ + --with-mpi-include=$ISSM_DIR/externalpackages/mpich2/install/include \ + --with-mpi-lib=" $ISSM_DIR/externalpackages/mpich2/install/lib/libpmpich.a $ISSM_DIR/externalpackages/mpich2/install/lib/libmpich.a $ISSM_DIR/externalpackages/mpich2/install/lib/libmpl.a " \ + --with-chaco-dir=$ISSM_DIR/externalpackages/chaco/install \ + --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a" \ + --with-math-lib="/usr/lib/libm.dylib" \ + --with-graphics-lib="/usr/X11/lib/libX11.dylib" \ + --with-cxxoptflags=" -fno-common -no-cpp-precomp -fexceptions -arch x86_64 -mmacosx-version-min=10.5 -O3 -DNDEBUG -w "\ + --with-numthreads=8 \ + --with-petsc-dir=$ISSM_DIR/externalpackages/petsc/install \ + --with-petsc-arch=$ISSM_ARCH \ + --with-blas-lapack-dir=$ISSM_DIR/externalpackages/petsc/install \ + --with-blacs-dir=$ISSM_DIR/externalpackages/petsc/install/ \ + --with-scalapack-dir=$ISSM_DIR/externalpackages/petsc/install/ \ + --with-mumps-dir=$ISSM_DIR/externalpackages/petsc/install/ \ + --with-dakota-dir=$ISSM_DIR/externalpackages/dakota/install\ + --with-math77-dir=$ISSM_DIR/externalpackages/math77/install\ + --with-gia Index: /issm/trunk-jpl/m4/issm_options.m4 =================================================================== --- /issm/trunk-jpl/m4/issm_options.m4 (revision 14738) +++ /issm/trunk-jpl/m4/issm_options.m4 (revision 14739) @@ -1381,16 +1381,39 @@ AC_MSG_RESULT($HAVE_RIFTS) dnl }}} + dnl math77{{{ + AC_ARG_WITH([math77-dir], + AS_HELP_STRING([--with-math77-dir=DIR], [math77 root directory.]), + [MATH77_ROOT=$withval],[MATH77_ROOT=""]) + AC_MSG_CHECKING(for math77) + + if test -d "$MATH77_ROOT"; then + + HAVE_MATH77=yes + MATH77LIB="-L$MATH77_ROOT/ -lmath77" + + AC_DEFINE([_HAVE_MATH77_],[1],[with math77 in ISSM src]) + AC_SUBST([MATH77LIB]) + else + HAVE_MATH77=no + fi + AC_MSG_RESULT($HAVE_MATH77) + dnl }}} dnl with-gia{{{ AC_ARG_WITH([gia], AS_HELP_STRING([--with-gia = YES], [compile with gia capabilities (default is yes)]), - [GIA=$withval],[GIA=yes]) + [GIA=$withval],[GIA=no]) AC_MSG_CHECKING(for gia capability compilation) if test "x$GIA" = "xyes"; then - + + if test "x$HAVE_MATH77" = "xno"; then + HAVE_GIA=no + AC_MSG_ERROR([gia requires compilation of math77 library! Reconfigure with --with-math77 option on]); + else dnl defaults HAVE_GIA=yes - AC_DEFINE([_HAVE_GIA_],[1],[with gia capability]) + fi + else HAVE_GIA=no Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 14738) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 14739) @@ -510,4 +510,11 @@ ./modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp\ ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h\ + ./modules/GiaDeflectionCorex/distme.f\ + ./modules/GiaDeflectionCorex/freed.f\ + ./modules/GiaDeflectionCorex/ojrule.f\ + ./modules/GiaDeflectionCorex/pwise.f\ + ./modules/GiaDeflectionCorex/qwise.f\ + ./modules/GiaDeflectionCorex/stot.f\ + ./modules/GiaDeflectionCorex/what0.f\ ./classes/GiaDeflectionCoreArgs.h @@ -918,4 +925,5 @@ libISSMCore_a_SOURCES = $(issm_sources) libISSMCore_a_CXXFLAGS = $(ALLCXXFLAGS) +libISSMCore_a_FFLAGS = $(AM_FFLAGS) if SHAREDLIBS @@ -925,4 +933,5 @@ libISSMCore_la_SOURCES = $(issm_sources) libISSMCore_la_LIBADD = $(PETSCLIB) $(TAOLIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB) $(ADOLCLIB) +libISSMCore_la_FFLAGS = $(AM_FFLAGS) libISSMOverload_la_SOURCES = ./shared/String/stricmp.c Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 14739) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f (revision 14739) @@ -0,0 +1,80 @@ + subroutine distme(idisk,iedge,isolt) + implicit double precision (a-h,o-y) + parameter (N3G = 1) + parameter (Ntime = 2) + parameter (Nafter = 1) + parameter (Ntimp = Ntime + Nafter) + double precision pset(7) + double precision time(Ntimp),dmi(Ntime),bi(Ntime),dumbt(Ntimp) + double precision hload(Ntime),qpat(Ntime),qt(Ntime) + real zradii(N3G),zhload(N3G,Ntime) + common /blockp/ pset + common /blockt/ time,bi,dmi + common /blocko/ rhoi,hload + common /blocky/ zhload + 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 call dvecpr(time,Ntime,'::::: time @ distme.f :::::',79,0,0) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 39 itime = 1, Ntime + hload(itime) = dble( zhload(idisk,itime) ) + 39 continue +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +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) ) + 70 continue +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. +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, Ntime + dmi(ind) = dmi(ind) / (( yearco * 1.0d3 ) / tfact ) + time(ind) = ( dumbt(ind) * yearco * 1.0d3 ) / tfact + 75 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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 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) ) + 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 14739) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f (revision 14739) @@ -0,0 +1,224 @@ + 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(3) + double precision ac0,ac1,ac2,ac3,ac4,ac5,ac6,ac7,ac8,ac9,ac10, + 1ac11 + common /blockz/ zkp + common /blockn/ ngx,ngy,irate + 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 + if(irate.eq.0) go to 999 + amps(1) = - decay(1) * amps(1) + amps(2) = - 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 14739) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f (revision 14739) @@ -0,0 +1,194 @@ + subroutine ojrule(dk,endr,bcin) + implicit double precision(a-h,o-z) + parameter (nhank = 1024) + parameter (nrv = 91) + parameter (N3G = 1) + parameter (npat = N3G) + double precision yvalue(nhank),bcin(nhank) +c SA ::::::::::::::::::::::::::::::::::::::::::::::::: + double precision wok, rpos +c SA ::::::::::::::::::::::::::::::::::::::::::::::::: + double precision pset(7) +c SA ::::::::::::::::::::::::::::::::::::::::::::::::: + real swok, asrpos, aswokm(npat), dist_rad + common /blockRad/ dist_rad +c SA ::::::::::::::::::::::::::::::::::::::::::::::::: + common /blockp/ pset + common /blocki/ INTERN + common /blockn/ ngx,ngy,irate + common /blocks/ aswokm,asrpos,swok + 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 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 + 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(ik) = bcin(ik) * rarg + 425 continue +c correct to end point val. in Simp. +c Rule + yvalue(nhank) = bcin(nhank) * rarg / two +c +c find the area under the curve using the Simpson's rule formulas +c + sumde = zero + do 300 int = 1, nhank + intp1 = int + 1 + ide = 2 + ( (-1)**intp1 + 1 ) + fide = dfloat(ide) + 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 ::::::::::::::::::::::::::::::::::::::::::::::::: + go to (1000,2000), iprate + 1000 call wolfc(rpos,wok) + go to 999 + 2000 call rates(rpos,wok) + 999 return +c SA ::::::::::::::::::::::::::::::::::::::::::::::::: + end +c + subroutine wolfc(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) +c SA ::::::::::::::::::::::::::::::::::::::::::::::::: +c real swoki(80), sri(80), swoko(160), sro(160) + common /blockp/ pset + common /blocki/ INTERN + common /blockn/ ngx,ngy,irate + common /blocks/ aswokm,asrpos,swok +c +c make single prec. and return to dimensional units. +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(INTERN) = 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 + + 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) +c SA ::::::::::::::::::::::::::::::::::::::::::::::::: + common /blockp/ pset + common /blocki/ INTERN + common /blockn/ ngx,ngy,irate + common /blocks/ aswokm,asrpos,swok + 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 make single prec. and return to dimensional units. +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(INTERN) = 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 + INP = INTERN - 1 + 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) + return + end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f (revision 14739) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/pwise.f (revision 14739) @@ -0,0 +1,40 @@ + subroutine pwise(t,ta,tb,xi1,xi2,slope,ycept,decay,bhaq) + 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 xi1 term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + xi1t =(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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + xi2t =(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 = (xi1 * xi1t) + (xi2 * xi2t) + return + end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f (revision 14739) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f (revision 14739) @@ -0,0 +1,65 @@ + subroutine qwise(t,ta,qjadon,xi0,xi1,xi2,slope,ycept,decay,bhaq) + implicit double precision (a-h,o-z) + double precision decay(2) + common /blockn/ ngx,ngy,irate + iprate = irate + 1 +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)) + gb1 = decay(1)*ycept + gb2 = decay(2)*ycept + go to ( 66, 67 ), iprate +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c define xi0t term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 66 xi0t = (xi0 + qjadon) * ( ( slope * t ) + ycept ) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c define xi1t term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + xi1t = 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 xi2t term: +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + xi2t = 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: +c (sign switch due to freed.f already +c having corrected in x1t, x2t pass). +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 67 xi0t = (xi0 + qjadon) * slope + xi1t =-xg1 * ( + 1 slope + 2 + ( gb1 + slope * ( ( ta * decay(1) ) - 1.0d0 )) + 3 * dexp( decay(1) * (ta - t) ) + 4 ) + xi2t =-xg2 * ( + 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 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + 90 bhaq = xi0t + xi1t + xi2t + return + end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f (revision 14739) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f (revision 14739) @@ -0,0 +1,71 @@ + subroutine stot(ikval,qjadon,fltng) + implicit double precision (a-h,o-z) + parameter (Nafter = 1) + parameter (Ntime = 2) + parameter (nhank = 2048) + parameter (Ntimp = Ntime + Nafter) + double precision decay(2) + double precision pset(7) + double precision time(Ntimp),bi(Ntime),dmi(Ntime) + double precision dekay1(nhank),dekay2(nhank),amp0(nhank), + 1amp1(nhank),amp2(nhank) + common /blockm/ dekay1,dekay2,amp0,amp1,amp2 + common /blockt/ time,bi,dmi +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) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: +c note that tspan must be updated in the calling routine "what0.f" +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + t = time(Ntimp) + sumb = 0.0d0 + ta = 0.0d0 + tb = time(1) +c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + do 97 i = 1,Ntime +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,slope,ycept,decay,bhaq) + 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,slope,ycept,decay,bhaq) + 39 sumb = bhaq + sumb + ta = time(i) + tb = time(i + 1) + 97 continue + fltng = sumb + return + end Index: /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f =================================================================== --- /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f (revision 14739) +++ /issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f (revision 14739) @@ -0,0 +1,142 @@ + subroutine what0(idisk,icheck,iedge,endr) + implicit double precision (a-h,o-z) + parameter (nhank = 2048) + 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), + 1decta(2),dyri1(nhank),dyri2(nhank),sna(nhank) + double precision cinner(nhank),bcin(nhank) +c double precision trsl(Nrsl) + integer maxk + common /blockp/ pset + common /blockz/ zkp + common /blockm/ dekay1,dekay2,amp0,amp1,amp2 +c common /blockk/ trsl + 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/ + twopi = two * pi + r2 = pset(6) + u2 = pset(4) + r1 = pset(5) + u1 = pset(3) + h = pset(1) + urat = u1/u2 +c SA: ================================================== + alphap = pset(7)/pset(1) +c alphap is dimensionless disk radius + 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 + if(icheck.eq.3) go to 999 + if(idisk.gt.1) go to 7001 + 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) + 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 Form vectors for full construction in pwise.f +c & stot.f +c Note that freed will produce decay spectra +c defined as positive, ie. negative decay must +c reinsert a minus sign. + dekay1(ik) = decay(1) + dekay2(ik) = decay(2) + amp0(ik) = amps(3) + amp1(ik) = amps(1) + amp2(ik) = amps(2) + zksam(ik) = zkd + zksamp(ik) = zkp + 7000 continue + 7001 continue +c call dvecpr(zksam,nhank,'k sampled like Wolf',79,0,0) +c call dvecpr(zksamp,nhank,'k - dimensionless',79,0,0) +c call dvecpr(amp0,nhank,'amp0 of k from freed',79,0,0) +c call dvecpr(amp1,nhank,'amp1 of k from freed',79,0,0) +c call dvecpr(amp2,nhank,'amp2 of k from freed',79,0,0) +c call dvecpr(dekay1,nhank,'dekay1 of k from freed',79,0,0) +c call dvecpr(dekay2,nhank,'dekay2 of k from freed',79,0,0) +c +c A call below perfomed plots of the two decay times in years +c sup {-1} vs. wave number k +c +c if(icheck.gt.2) go to 49 +c call wplot(dyri1,dyri2,sna) +c go to 999 +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 +c8499 tspan = trsl(irsl) +c trsl was nondimesionalized +c in the routine distme.f + 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) + cinner(ik) = - fltng * pref * twoap + bcin(ik) = cinner(ik) * dbesj1(xakap) + 8000 continue +c write(6,133) diku,xakap,pset(1),pset(7),alphap,urat +c 133 format(' last diku xakap in what0'/1h ,1p6d16.8) +c write(6,1) bcin(nhank),zksamp(nhank),fltng +c 1 format(' last bcin zksamp and fltng in what0'/1h ,1p3d16.8) +c write(6,2) cinner(nhank),cinner(1),qjadon +c 2 format(' last and first cinner in what0 and qjadon'/1h ,1p3d16.8) +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,endr,bcin) + go to 999 +c9499 tspan = trsl(irsl) + 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) + cinner(ik) = - fltng * pref * twoap + bcin(ik) = cinner(ik) * oxakap * ( dsin(xakap) * oxakap + 1 - dcos(xakap) ) + 9000 continue +c write(6,1) bcin(nhank),zksamp(nhank) + call ojrule(dk,endr,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) + 999 return + end