Changeset 14768


Ignore:
Timestamp:
04/26/13 13:57:09 (12 years ago)
Author:
adhikari
Message:

CHG: we now can evaluate solid Earth's response to steady-state ice geometry. Response to transient geomety to be implemented soon

Location:
issm/trunk-jpl/src/c/modules/GiaDeflectionCorex
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp

    r14752 r14768  
    3737
    3838struct blocky{
    39         double zhload[Ntime][N3G];
     39        double zhload[Ntime];
    4040};
    4141     
     42struct blockn{
     43        int irate;
     44};
     45
     46struct blockrad{
     47        double distrad;
     48};
     49
     50struct blocks{
     51        double aswokm;
     52        double asrpos;
     53        double swok;
     54};
    4255
    4356extern "C" {
    44         int distme_( int* pidisk,int* piedge,int* pisolt);
     57        int distme_( int* pidisk,int* piedge);
     58        int what0_( int* pidisk,int* piedge);
    4559        extern struct blockp blockp_;
    4660        extern struct blockt blockt_;
    4761        extern struct blocko blocko_;
    4862        extern struct blocky blocky_;
     63        extern struct blockn blockn_;
     64        extern struct blockrad blockrad_;
     65        extern struct blocks blocks_;
    4966}
     67
    5068/*}}}*/
    5169     
    … …  
    5876        int iedge=1; //c iedge ......... = 1 square-edged, = 2 elliptical disc x-section  (see naruse.f)
    5977        int idisk=1; // disk #
    60         int isolt=1002;
     78        int irate=0; // =0 fetch w solution (m) only; =1 dw/dt (mm/yr) only
    6179
    6280        /*intermediary: */
    … …  
    84102        /*some debug info: */
    85103        int disk_id;
    86         /*}}}*/
    87104
     105/*}}}*/
     106     
    88107        /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/
    89108        ri=arguments->ri;
    … …  
    104123        /*}}}*/
    105124
    106         /*Documentation for pset coming from naruse.f:
    107         c pset(1) ....... h (thickness of elastic lithosphere in km)
    108         c pset(2) ....... viscosity of sub-layer in Pa*s
    109         c pset(3) ....... elastic rigidity of lithosphere in dynes/cm^2
    110         c pset(4) ....... elastic rigidity of sub-lithosphere in dynes/cm^2
    111         c pset(5) ....... density of lithosphere in gm/cc
    112         c pset(6) ....... density of sub-lithosphere in gm/cc
    113         c pset(7) ....... =dset(1)
    114         c dset(1) ....... disk radius in km
     125        // Now, let's set pset from the data that we got in input to GiaDeflectionCorex:
    115126
    116         Now, let's set pset from the data that we got in input to GiaDeflectionCorex:
    117         */
    118127        blockp_.pset[0]=lithosphere_thickness;
    119128        blockp_.pset[1]=mantle_viscosity;
    … …  
    123132        blockp_.pset[5]=mantle_density;
    124133        blockp_.pset[6]=re;
    125         blocko_.rhoi=rho_ice;
     134//      blocko_.rhoi=rho_ice; // use this for non-benchmark experiments
     135//      blockp_.pset[6]=800000.0; // for testing benchmark, should be in meters
     136        blocko_.rhoi=1000.0;   // use this for benchmark not rho_ice = 917 kg m^-3
    126137     
    127138        /*loading history: */
    128         blocky_.zhload[0][0]=current_he;
    129         blocky_.zhload[1][0]=current_he;
     139        blocky_.zhload[0]=current_he;
     140        blocky_.zhload[1]=current_he;
     141//      blocky_.zhload[0]=3000; // for testing benchmark
     142//      blocky_.zhload[1]=3000; // for testing benchmark 
    130143
    131144        /*times: */
    132         blockt_.time[0]=0.1;
    133         blockt_.time[1]=2000;
    134         blockt_.time[2]=1000;
     145        blockt_.time[0]=1.e-4;    // in kyr
     146        blockt_.time[1]=2500e+0;  // in kyr
     147        blockt_.time[2]=2400e+0;  // control this for benchmark experiments
    135148
     149        /*irate: */
     150        blockn_.irate=irate;
     151
     152        /*radial distance of i-th element: */
     153        blockrad_.distrad=ri/1000.0; // in km
     154//      blockrad_.distrad=500.0; // for testing benchmark
    136155
    137156        /*Call distme driver: */
    138         distme_(&idisk,&iedge,&isolt);
     157        distme_(&idisk,&iedge);
     158       
     159        /*Call what0 driver: */
     160        what0_(&idisk,&iedge);
    139161
    140         //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]);
     162        /*this is the solution: */
     163        wi = blocks_.aswokm;
    141164
    142165        /*allocate output pointer: */
    143166        *pwi=wi;
     167
     168//      printf("wi: %g deflection: %g \n",wi,blocks_.aswokm);
     169
    144170}

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/distme.f

    r14739 r14768  
    1       subroutine distme(idisk,iedge,isolt)
     1      subroutine distme(idisk,iedge)
    22      implicit double precision (a-h,o-y)
    33      parameter (N3G = 1)
    … …  
    88      double precision time(Ntimp),dmi(Ntime),bi(Ntime),dumbt(Ntimp)
    99      double precision hload(Ntime),qpat(Ntime),qt(Ntime)
    10       real zradii(N3G),zhload(N3G,Ntime)
     10      double precision zradii(N3G),zhload(Ntime),rhoi,distrad
    1111      common /blockp/ pset
     12      common /blockrad/ distrad
    1213      common /blockt/ time,bi,dmi
     14      common /blocky/ zhload
    1315      common /blocko/ rhoi,hload
    14       common /blocky/ zhload
    1516      data g /9.832186d0/, yearco /3.15576d7/, eradm/6.371d6/
    1617      data dpi /3.1415926535897932d0/, dzero/0.0d0/
    … …  
    2324      dumbt(k) = time(k)
    2425  776 continue
     26c      write(6,*) time(1), time(2), time(3)
     27c      write(6,*) zhload(1), zhload(2), distrad, rhoi
     28c      write(6,*) pset(1), pset(2), pset(3), pset(4), pset(5), pset(6)
     29c      write(6,*) pset(7)
    2530c      call dvecpr(time,Ntime,'::::: time @ distme.f :::::',79,0,0)
    2631c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    2732      do 39 itime = 1, Ntime
    28       hload(itime) = dble( zhload(idisk,itime) )
     33      hload(itime) = dble( zhload(itime) )
    2934   39 continue
    3035c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    3843      bi(i) = hload(i-1) - ( dmi(i)*dumbt(i-1) ) 
    3944   70 continue
     45c      write(6,*) zhload(1,1), zhload(1,2)
    4046c      call dvecpr(hload,Ntime,'::::: hload @ distme.f :::::',79,0,0)
    4147c      call dvecpr(dmi,Ntime,'::::: load slope @ distme.f :::::',79,0,0)

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/freed.f

    r14739 r14768  
    66     1ac11
    77      common /blockz/ zkp
    8       common /blockn/ ngx,ngy,irate
     8      common /blockn/ irate
    99      data zero /0.0d0/, g /9.832186d0/
    1010c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/ojrule.f

    r14739 r14768  
    1       subroutine ojrule(dk,endr,bcin)
     1      subroutine ojrule(dk,bcin)
    22      implicit double precision(a-h,o-z)
    33      parameter (nhank = 1024)
    4       parameter (nrv = 91)
    54      parameter (N3G = 1)
    65      parameter (npat = N3G)
    … …  
    1110      double precision pset(7)
    1211c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    13       real swok, asrpos, aswokm(npat), dist_rad
    14       common /blockRad/ dist_rad
     12      double precision swok, asrpos, aswokm, distrad
     13      common /blockrad/ distrad
    1514c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    1615      common /blockp/ pset
    1716      common /blocki/ INTERN
    18       common /blockn/ ngx,ngy,irate
     17      common /blockn/ irate
    1918      common /blocks/ aswokm,asrpos,swok
    2019      data zero /0.0d0/, one /1.0d0/, two /2.0d0/, three /3.0d0/,
    … …  
    2221      if(INTERN.ge.npat) INTERN = 0
    2322      INTERN = INTERN + 1
    24 c     write(6,9099) INTERN
    25 c9099 format('    INTERN'/1h 1p1i10)
    2623      iprate = irate + 1
    27 c
    28 c Brute force integration for the inverse
    29 c Hankel transform.  (Actually this is Simpson's rule)
    30 c
    31 c SA: open an output file
    32 c SA:     open(unit=33,file='outputs/r',status='new')
    33 c
    3424c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    35 c here I just compute for given r
    36 c      dr = endr/dfloat(nrv)
    3725      bath = dk / three
    38 c      r = zero
    39 c      do 450 ir = 1, nrv
    40 c      r = r + dr
    41 c      rpos(ir) = r
    4226c SA :::: rpose should be normalized wrt lithosphere thickness = 100 km
    4327c SA :: give r is normalized dist_rad :: r == dist_rad / h
    44       r = dist_rad / 100.00
     28c      write(6,*) irate, distrad
     29      r = distrad / 100.00
    4530      rpos = r
    4631      ak = zero
    … …  
    9782c
    9883      subroutine wolfc(rpos,wok)
    99       parameter (nrv = 91)
    10084      parameter (N3G = 1)
    10185      parameter (npat = N3G)
    … …  
    10589      double precision pset(7)
    10690c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    107       real swok,asrpos, aswokm(npat)
     91      double precision swok,asrpos, aswokm
    10892c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    10993c     real swoki(80), sri(80), swoko(160), sro(160)
    11094      common /blockp/ pset
    11195      common /blocki/ INTERN
    112       common /blockn/ ngx,ngy,irate
     96      common /blockn/ irate
    11397      common /blocks/ aswokm,asrpos,swok
    11498c
    … …  
    121105      swok = hscale * sngl(wok)
    122106      asrpos = hsckm * sngl(rpos)
    123       aswokm(INTERN) = swok
     107      aswokm = swok
    124108c   45 continue
    125109c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    147131
    148132      subroutine rates(rpos,wok)
    149       parameter (nrv = 91)
    150133      parameter (N3G = 1)
    151134      parameter (npat = N3G)
    152135c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    153136      double precision wok, rpos
    154 c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    155137      double precision pset(7)
    156 c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    157       real swok, asrpos, aswokm(npat)
     138      double precision swok, asrpos, aswokm
    158139c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    159140      common /blockp/ pset
    160141      common /blocki/ INTERN
    161       common /blockn/ ngx,ngy,irate
     142      common /blockn/ irate
    162143      common /blocks/ aswokm,asrpos,swok
    163144      data ngyo /050201/
    … …  
    176157     1                  * ( sngl(pset(4))/ sngl(pset(2)) )
    177158      asrpos = hsckm * sngl(rpos)
    178       aswokm(INTERN) = swok
     159      aswokm = swok
    179160c   45 continue
    180161c SA :::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    184165c Like Wolf (1985), do "inside" and "outside" curves separately.
    185166c
    186       INP = INTERN - 1
    187       if(INP.gt.6) INP = 6
     167c      INP = INTERN - 1
     168c      if(INP.gt.6) INP = 6
    188169c
    189170c next call ploting .
    190171c
    191       call dvecpr(pset,7,'pset in ojrule',79,0,0)
    192       call svecpr(swok,nrv,'swok in rates - mm per yr',79,0,0)
     172c      call dvecpr(pset,7,'pset in ojrule',79,0,0)
     173c      call svecpr(swok,nrv,'swok in rates - mm per yr',79,0,0)
    193174      return
    194175      end

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/qwise.f

    r14739 r14768  
    22      implicit double precision (a-h,o-z)
    33      double precision decay(2)
    4       common /blockn/ ngx,ngy,irate
     4      common /blockn/ irate
    55      iprate = irate + 1
    66c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/stot.f

    r14739 r14768  
    33      parameter (Nafter = 1)
    44      parameter (Ntime = 2)
    5       parameter (nhank = 2048)
     5      parameter (nhank = 1024)
    66      parameter (Ntimp = Ntime + Nafter)
    77      double precision decay(2)

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f

    r14739 r14768  
    1       subroutine what0(idisk,icheck,iedge,endr)
     1      subroutine what0(idisk,iedge)
    22      implicit double precision (a-h,o-z)
    3       parameter (nhank = 2048)
     3      parameter (nhank = 1024)
    44      parameter (N = nhank/2)
    55      double precision dekay1(nhank),dekay2(nhank),amp0(nhank),
    … …  
    1010c     double precision trsl(Nrsl)
    1111      integer maxk
     12c      double precision aswokm, asrpos, swok
     13c      common /blocks/ aswokm,asrpos,swok
    1214      common /blockp/ pset
    1315      common /blockz/ zkp
    … …  
    3638      dk = endk/dfloat(nhank)
    3739c
    38       if(icheck.eq.3) go to 999
     40c      if(icheck.eq.3) go to 999
    3941      if(idisk.gt.1) go to 7001
    4042      ak = zero
    … …  
    121123c for each of N3G disks in "aswokm(nrv,N3G)" .
    122124c
    123       call ojrule(dk,endr,bcin)
     125      call ojrule(dk,bcin)
    124126      go to 999
    125127c9499 tspan = trsl(irsl)
    … …  
    136138 9000 continue
    137139c     write(6,1) bcin(nhank),zksamp(nhank)
    138       call ojrule(dk,endr,bcin)
     140      call ojrule(dk,bcin)
    139141c     call dvecpr(cinner,nhank,'cinner at time t',79,0,0)
    140142c     call dvecpr(bcin,nhank,'bcin at time t',79,0,0)
Note: See TracChangeset for help on using the changeset viewer.