Changeset 15140


Ignore:
Timestamp:
05/29/13 13:53:49 (12 years ago)
Author:
adhikari
Message:

CHG: now one can compute both W and dWdt solutions at the same time

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

Legend:

Unmodified
Added
Removed

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

    r15037 r15140  
    2121};
    2222
    23 struct blockn{
    24         int irate;
    25 };
    26 
    2723struct blockrad{
    2824        double distrad;
    … …  
    3026
    3127struct blocks{
    32         double aswokm;
     28        double aswokm_w;
     29        double aswokm_dwdt;
    3330};
    3431
    … …  
    3936        extern struct blockp blockp_;
    4037        extern struct blocko blocko_;
    41         extern struct blockn blockn_;
    4238        extern struct blockrad blockrad_;
    4339        extern struct blocks blocks_;
    … …  
    6359        /*gia solution parameters: */
    6460        int iedge=0;
    65         int irate=0;
    6661
    6762        /*gia inputs: */
    … …  
    111106        rho_ice=arguments->rho_ice;
    112107        disk_id=arguments->idisk;
    113         irate=arguments->irate;
    114108        iedge=arguments->iedge;
    115109        yts=arguments->yts;
    … …  
    118112
    119113        /*Modify inputs to match naruse code: */
    120         //from our model, irate comes in with values in [1,2], which maps into [0,1] in naruse:
    121         irate=irate-1;
    122114        Ntime=numtimes;
    123115        Ntimm=Ntime-1;
    … …  
    155147        blockt_dmi=xNew<IssmDouble>(Ntimm);
    156148
    157         /*irate: */
    158         blockn_.irate=irate;
    159 
    160149        /*radial distance of i-th element: */
    161150        blockrad_.distrad=ri/1000.0; // in km
    … …  
    169158
    170159        /*output solution: */
    171         if(irate==0){
    172                 wi = blocks_.aswokm;
    173                 *pwi=wi;
    174                 *pdwidt=0;
    175         }
    176         else if (irate==1){
    177                 *pwi=0;
    178                 dwidt = blocks_.aswokm;
    179                 *pdwidt=dwidt;
    180         }
     160        wi = blocks_.aswokm_w;
     161        dwidt = blocks_.aswokm_dwdt;
     162        *pwi=wi;
     163        *pdwidt=dwidt;
    181164
    182165        /*Free ressources: */

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

    r14768 r15140  
    22     1,amps)       
    33      implicit double precision (a-h,o-z)
    4       double precision decay(2),amps(3)
     4      double precision decay(2),amps(5)
    55      double precision ac0,ac1,ac2,ac3,ac4,ac5,ac6,ac7,ac8,ac9,ac10,
    66     1ac11
    77      common /blockz/ zkp
    8       common /blockn/ irate
    98      data zero /0.0d0/, g /9.832186d0/
    109c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    215214      amps(2) =  decdif*( decay(2) * ( a1 - a2*decay(2) ) - a0 )
    216215      amps(3) = a2
    217       if(irate.eq.0) go to 999
    218       amps(1) = - decay(1) * amps(1)
    219       amps(2) = - decay(2) * amps(2)
     216      amps(4) = - decay(1) * amps(1)
     217      amps(5) = - decay(2) * amps(2)
    220218      go to 999
    221219 9990 write(6,998) idgen

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

    r15024 r15140  
    1       subroutine ojrule(dk,bcin)
     1      subroutine ojrule(dk,bcin_w,bcin_dwdt)
    22      implicit double precision(a-h,o-z)
    33      parameter (nhank = 1024)
    4       double precision yvalue(nhank),bcin(nhank)
    5       double precision wok, rpos
     4      double precision yvalue_w(nhank),yvalue_dwdt(nhank)
     5      double precision bcin_w(nhank),bcin_dwdt(nhank)
     6      double precision wok_w,wok_dwdt,rpos
     7      double precision swok_w,swok_dwdt
    68      double precision pset(7)
    7       double precision aswokm, distrad
     9      double precision aswokm_w,aswokm_dwdt,distrad
    810c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    911      common /blockrad/ distrad
    1012      common /blockp/ pset
    11       common /blockn/ irate
    12       common /blocks/ aswokm
     13      common /blocks/ aswokm_w,aswokm_dwdt
    1314      data zero /0.0d0/, one /1.0d0/, two /2.0d0/, three /3.0d0/,
    1415     1rescal/ 1.0d0/
     16      data yearco /3.15576d7/
    1517c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    16       iprate = irate + 1
    1718      bath = dk / three
    1819c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    19 c rpose should be normalized wrt lithosphere thickness
     20c rpos should be normalized wrt lithosphere thickness
    2021c give r is normalized dist_rad :: r == dist_rad / h
    2122c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    3031      rak = ak * r
    3132      rarg = dbesj0( rak )
    32       yvalue(ik) = bcin(ik) * rarg
     33      yvalue_w(ik) = bcin_w(ik) * rarg
     34      yvalue_dwdt(ik) = bcin_dwdt(ik) * rarg
    3335  425 continue
    3436c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    3739c find the area under the curve using the Simpson's rule formulas
    3840c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    39       sumde = zero
     41      sumde_w = zero
     42      sumde_dwdt = zero
    4043      do 300 int = 1, nhank
    4144      intp1 = int + 1
    4245      ide = 2 + ( (-1)**intp1 + 1 )
    4346      fide = dfloat(ide)
    44       sumde = ( fide * yvalue(int) ) + sumde
     47      sumde_w = ( fide * yvalue_w(int) ) + sumde_w
     48      sumde_dwdt = ( fide * yvalue_dwdt(int) ) + sumde_dwdt
    4549  300 continue
    46       wok = bath * sumde
    47 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    48       go to (1000,2000), iprate
    49  1000 call wolfc(rpos,wok)
    50       go to 999
    51  2000 call rates(rpos,wok)
    52   999 return
    53 
    54       end
    55 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    56 
    57 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    58       subroutine wolfc(rpos,wok)
    59       double precision wok, rpos
    60       double precision pset(7)
    61       double precision swok,asrpos, aswokm
    62 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    63       common /blockp/ pset
    64       common /blocks/ aswokm
    65 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    66 c make single prec. and return to dimensional units.
     50      wok_w = bath * sumde_w
     51      wok_dwdt = bath * sumde_dwdt
     52     
    6753c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    6854      hscale = sngl(pset(1))
    6955      hsckm = hscale / 1.0e3
    70       swok = hscale * sngl(wok)
    71       asrpos = hsckm * sngl(rpos)
    72       aswokm = swok
     56      swok_w = hscale * sngl(wok_w)
     57      aswokm_w = swok_w
     58      swok_dwdt = (hscale * yearco * 1.0e3 * sngl(wok_dwdt))
     59     1                  * ( sngl(pset(4))/ sngl(pset(2)) )
     60      aswokm_dwdt = swok_dwdt
     61      return
    7362
    74       return
    7563      end
    76 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    77 
    78 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    79       subroutine rates(rpos,wok)
    80       double precision wok, rpos
    81       double precision pset(7)
    82       double precision swok, asrpos, aswokm
    83 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    84       common /blockp/ pset
    85       common /blocks/ aswokm
    86 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    87       data ngyo /050201/
    88       data ngyii /050201/
    89       data yearco /3.15576d7/
    90 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    91 c make single prec. and return to dimensional units.
    92 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    93       hscale = sngl(pset(1))
    94       hsckm = hscale / 1.0e3
    95       swok = (hscale * yearco * 1.0e3 * sngl(wok))
    96      1                  * ( sngl(pset(4))/ sngl(pset(2)) )
    97       asrpos = hsckm * sngl(rpos)
    98       aswokm = swok
    99 
    100       return
    101       end

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

    r14739 r15140  
    1       subroutine pwise(t,ta,tb,xi1,xi2,slope,ycept,decay,bhaq)
     1      subroutine pwise(t,ta,tb,xi1,xi2,xi3,xi4,slope,ycept,decay,
     2     1bhaq_w,bhaq_dwdt)
    23      implicit double precision (a-h,o-z)
    34      double precision decay(2)
    … …  
    1920      eb2 = dexp(gbt2)
    2021c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    21 c define xi1 term:
     22c define xit1 term:
    2223c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    23       xi1t =(ycept/decay(1)) * (eb1 - ea1) -
     24      xit1 =(ycept/decay(1)) * (eb1 - ea1) -
    2425     1(slope/(decay(1)*decay(1))) *
    2526     2                            ( (1.0d0 - tb*decay(1))*eb1
    2627     3                            - (1.0d0 - ta*decay(1))*ea1 )
    2728c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    28 c define xi2 term:
     29c define xit2 term:
    2930c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    30       xi2t =(ycept/decay(2)) * (eb2 - ea2) -
     31      xit2 =(ycept/decay(2)) * (eb2 - ea2) -
    3132     1(slope/(decay(2)*decay(2))) *
    3233     2                            ( (1.0d0 - tb*decay(2))*eb2
    … …  
    3637c ABOVE IS THE NON-DEGENERATE CASE
    3738c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    38       bhaq = (xi1 * xi1t) + (xi2 * xi2t)   
     39      bhaq_w = (xi1 * xit1) + (xi2 * xit2)   
     40      bhaq_dwdt = (xi3 * xit1) + (xi4 * xit2)   
    3941      return
    4042      end

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

    r14768 r15140  
    1       subroutine qwise(t,ta,qjadon,xi0,xi1,xi2,slope,ycept,decay,bhaq)
     1      subroutine qwise(t,ta,qjadon,xi0,xi1,xi2,xi3,xi4,slope,ycept,
     2     1decay,bhaq_w,bhaq_dwdt)
    23      implicit double precision (a-h,o-z)
    34      double precision decay(2)
    4       common /blockn/ irate
    5       iprate = irate + 1
    65c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    76c This subroutine retrieves the convolution for the J-th linear piece-wise
    … …  
    1817      xg1 = xi1/(decay(1)*decay(1))
    1918      xg2 = xi2/(decay(2)*decay(2))
     19      xg3 = xi3/(decay(1)*decay(1))
     20      xg4 = xi4/(decay(2)*decay(2))
    2021      gb1 = decay(1)*ycept
    2122      gb2 = decay(2)*ycept
    22       go to ( 66, 67 ), iprate
    2323c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    24 c define xi0t term:
     24c define xit0 term:
    2525c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    26    66 xi0t = (xi0 + qjadon) * ( ( slope * t ) + ycept )
     26      xit0_w = (xi0 + qjadon) * ( ( slope * t ) + ycept )
    2727c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    28 c define xi1t term:
     28c define xit1 term:
    2929c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    30       xi1t = xg1 * (
     30      xit1_w = xg1 * (
    3131     1              gb1 + slope * ( ( t * decay(1) ) - 1.0d0 )
    3232     2          - ( gb1 + slope * ( ( ta * decay(1) ) - 1.0d0 ))
    … …  
    3434     4                       )
    3535c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    36 c define xi2t term:
     36c define xit2 term:
    3737c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    38       xi2t = xg2 * (
     38      xit2_w = xg2 * (
    3939     1              gb2 + slope * ( ( t * decay(2) ) - 1.0d0 )
    4040     2          - ( gb2 + slope * ( ( ta * decay(2) ) - 1.0d0 ) )
    4141     3                                   * dexp( decay(2) * (ta - t) )
    4242     4                       )
    43       go to 90
    4443c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    4544c And the rate equivalents:
    … …  
    4746c having corrected in x1t, x2t pass).
    4847c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    49    67 xi0t = (xi0 + qjadon) * slope
    50       xi1t =-xg1 * (
     48      xit0_dwdt = (xi0 + qjadon) * slope
     49      xit1_dwdt =-xg3 * (
    5150     1              slope 
    5251     2     + ( gb1 + slope * ( ( ta * decay(1) ) - 1.0d0 ))
    5352     3                                   * dexp( decay(1) * (ta - t) )
    5453     4                       )
    55       xi2t =-xg2 * (
     54      xit2_dwdt =-xg4 * (
    5655     1              slope
    5756     2     + ( gb2 + slope * ( ( ta * decay(2) ) - 1.0d0 ))
    … …  
    6160c add terms for the J-th (and final) interval contribution.
    6261c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    63    90 bhaq = xi0t + xi1t + xi2t
     62      bhaq_w = xit0_w + xit1_w + xit2_w
     63      bhaq_dwdt = xit0_dwdt + xit1_dwdt + xit2_dwdt
    6464      return
    6565      end

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

    r15024 r15140  
    1       subroutine stot(ikval,qjadon,fltng,Ntimp,Ntimm,time,bi,dmi)
     1      subroutine stot(ikval,qjadon,fltng_w,fltng_dwdt,Ntimp,Ntimm,
     2     1time,bi,dmi)
    23      implicit double precision (a-h,o-z)
    34      integer Ntimp,Ntimm
    … …  
    89      double precision time(Ntimp),bi(Ntimm),dmi(Ntimm)
    910      double precision dekay1(nhank),dekay2(nhank),amp0(nhank),
    10      1amp1(nhank),amp2(nhank)
     11     1amp1(nhank),amp2(nhank),amp3(nhank),amp4(nhank)
    1112c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    12       common /blockm/ dekay1,dekay2,amp0,amp1,amp2
     13      common /blockm/ dekay1,dekay2,amp0,amp1,amp2,amp3,amp4
    1314c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    1415c  This subroutine returns the inverse Laplace transform to the
    … …  
    3637      xi1 = amp1(ikval)
    3738      xi2 = amp2(ikval)
     39      xi3 = amp3(ikval)
     40      xi4 = amp4(ikval)
    3841
    3942c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    4144c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    4245      t = time(Ntimp)
    43       sumb = 0.0d0
     46      sumb_w = 0.0d0
     47      sumb_dwdt = 0.0d0
    4448      ta = time(1)
    4549      tb = time(2)
    … …  
    5660c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    5761      if( t . gt . ta . and . t . le . tb) go to 38
    58       call pwise(t,ta,tb,xi1,xi2,slope,ycept,decay,bhaq)
     62      call pwise(t,ta,tb,xi1,xi2,xi3,xi4,slope,ycept,decay,
     63     1bhaq_w,bhaq_dwdt)
    5964      go to 39
    6065c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    6267c evaluation still has to consider the load itself
    6368c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    64    38 call qwise(t,ta,qjadon,xi0,xi1,xi2,slope,ycept,decay,bhaq)
    65    39 sumb = bhaq + sumb
     69   38 call qwise(t,ta,qjadon,xi0,xi1,xi2,xi3,xi4,slope,ycept,decay,
     70     1bhaq_w,bhaq_dwdt)
     71   39 sumb_w = bhaq_w + sumb_w
     72      sumb_dwdt = bhaq_dwdt + sumb_dwdt
    6673      ta = time(i + 1)
    6774      tb = time(i + 2)
    6875   97 continue
    69       fltng = sumb
     76      fltng_w = sumb_w
     77      fltng_dwdt = sumb_dwdt
    7078      return
    7179      end

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

    r15025 r15140  
    55      parameter (N = nhank/2)
    66      double precision dekay1(nhank),dekay2(nhank),amp0(nhank),
    7      1amp1(nhank),amp2(nhank),zksam(nhank),zksamp(nhank)
    8       double precision decay(2),pset(7),amps(3),
     7     1amp1(nhank),amp2(nhank),amp3(nhank),amp4(nhank),
     8     1zksam(nhank),zksamp(nhank)
     9      double precision decay(2),pset(7),amps(5),
    910     1decta(2),dyri1(nhank),dyri2(nhank),sna(nhank)
    10       double precision cinner(nhank),bcin(nhank)
     11      double precision cinner_w(nhank),cinner_dwdt(nhank)
     12      double precision bcin_w(nhank),bcin_dwdt(nhank)
    1113      double precision time(Ntimp),bi(Ntimm),dmi(Ntimm)
    1214      integer maxk
    … …  
    1416      common /blockp/ pset
    1517      common /blockz/ zkp
    16       common /blockm/ dekay1,dekay2,amp0,amp1,amp2
     18      common /blockm/ dekay1,dekay2,amp0,amp1,amp2,amp3,amp4
    1719      data yearco /3.15576d7/, pi /3.1415926535897932384d0/
    1820      data g /9.832186d0/, four /4.d0/, two /2.0d0/,
    … …  
    6971      amp1(ik) = amps(1) 
    7072      amp2(ik) = amps(2) 
     73      amp3(ik) = amps(4) 
     74      amp4(ik) = amps(5) 
    7175      zksam(ik) = zkd
    7276      zksamp(ik) = zkp
    … …  
    8993      pref = diku / ( diku + rghm )
    9094      qjadon = one / ( four * zksamp(ik) * urat )
    91       call stot(ik,qjadon,fltng,Ntimp,Ntimm,time,bi,dmi)
    92       cinner(ik) = - fltng * pref * twoap
    93       bcin(ik) = cinner(ik) * dbesj1(xakap)
     95      call stot(ik,qjadon,fltng_w,fltng_dwdt,Ntimp,Ntimm,time,bi,dmi)
     96      cinner_w(ik) = - fltng_w * pref * twoap
     97      cinner_dwdt(ik) = - fltng_dwdt * pref * twoap
     98      bcin_w(ik) = cinner_w(ik) * dbesj1(xakap)
     99      bcin_dwdt(ik) = cinner_dwdt(ik) * dbesj1(xakap)
    94100 8000 continue
    95101c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    99105c for each of N3G disks in "aswokm(nrv,N3G)" .
    100106c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    101       call ojrule(dk,bcin)
     107      call ojrule(dk,bcin_w,bcin_dwdt)
    102108      go to 999
    103109c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    … …  
    108114      pref = diku / ( diku + rghm )
    109115      qjadon = one / ( four * zksamp(ik) * urat )
    110       call stot(ik,qjadon,fltng,Ntimp,Ntimm,time,bi,dmi)
    111       cinner(ik) = - fltng * pref * twoap
    112       bcin(ik) = cinner(ik) * oxakap * ( dsin(xakap) * oxakap
     116      call stot(ik,qjadon,fltng_w,fltng_dwdt,Ntimp,Ntimm,time,bi,dmi)
     117      cinner_w(ik) = - fltng_w * pref * twoap
     118      cinner_dwdt(ik) = - fltng_dwdt * pref * twoap
     119      bcin_w(ik) = cinner_w(ik) * oxakap * ( dsin(xakap) * oxakap
     120     1 - dcos(xakap) )
     121      bcin_dwdt(ik) = cinner_dwdt(ik) * oxakap * ( dsin(xakap) * oxakap
    113122     1 - dcos(xakap) )
    114123 9000 continue
    115       call ojrule(dk,bcin)
     124      call ojrule(dk,bcin_w,bcin_dwdt)
    116125  999 return
    117126
Note: See TracChangeset for help on using the changeset viewer.