Changeset 14841
- Timestamp:
- 05/01/13 16:52:43 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/what0.f
r14768 r14841 8 8 1decta(2),dyri1(nhank),dyri2(nhank),sna(nhank) 9 9 double precision cinner(nhank),bcin(nhank) 10 c double precision trsl(Nrsl)11 10 integer maxk 12 c double precision aswokm, asrpos, swok 13 c common /blocks/ aswokm,asrpos,swok 11 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 14 12 common /blockp/ pset 15 13 common /blockz/ zkp 16 14 common /blockm/ dekay1,dekay2,amp0,amp1,amp2 17 c common /blockk/ trsl18 15 data yearco /3.15576d7/, pi /3.1415926535897932384d0/ 19 16 data g /9.832186d0/, four /4.d0/, two /2.0d0/, 20 17 1 one /1.0d0/, zero/0.0d0/ , maxk/64/ 18 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 21 19 twopi = two * pi 22 20 r2 = pset(6) … … 26 24 h = pset(1) 27 25 urat = u1/u2 28 c SA: ================================================== 26 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 27 c alphap is dimensionless disk radius 28 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 29 29 alphap = pset(7)/pset(1) 30 c alphap is dimensionless disk radius31 30 twoap = two * alphap 32 31 rghm = ( r1 * g * h * alphap ) / (two * u2) … … 38 37 dk = endk/dfloat(nhank) 39 38 c 40 c if(icheck.eq.3) go to 99941 39 if(idisk.gt.1) go to 7001 42 40 ak = zero … … 52 50 e2 = dexp(zkp2) 53 51 e4 = dexp(zkp4) 52 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 54 53 call freed(r2,u2,r1,u1,h,zkd,e1,e2,e4,b0,b1,a2,a1,a0,decay,amps) 55 c 54 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 56 55 decta(1) = decay(1)/tmxyr 57 56 decta(2) = decay(2)/tmxyr … … 59 58 dyri2(ik) = decta(2) 60 59 sna(ik) = pikn 61 c Form vectors for full construction in pwise.f62 c &stot.f63 c Note that freed will produce decay spectra64 c defined as positive, ie. negative decay must65 c reinsert a minus sign.60 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 61 c Form vectors for full construction in pwise.f and stot.f 62 c Note that freed will produce decay spectra defined as positive, 63 c ie. negative decay must reinsert a minus sign. 64 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 66 65 dekay1(ik) = decay(1) 67 66 dekay2(ik) = decay(2) … … 73 72 7000 continue 74 73 7001 continue 75 c call dvecpr(zksam,nhank,'k sampled like Wolf',79,0,0) 76 c call dvecpr(zksamp,nhank,'k - dimensionless',79,0,0) 77 c call dvecpr(amp0,nhank,'amp0 of k from freed',79,0,0) 78 c call dvecpr(amp1,nhank,'amp1 of k from freed',79,0,0) 79 c call dvecpr(amp2,nhank,'amp2 of k from freed',79,0,0) 80 c call dvecpr(dekay1,nhank,'dekay1 of k from freed',79,0,0) 81 c call dvecpr(dekay2,nhank,'dekay2 of k from freed',79,0,0) 82 c 83 c A call below perfomed plots of the two decay times in years 84 c sup {-1} vs. wave number k 85 c 86 c if(icheck.gt.2) go to 49 87 c call wplot(dyri1,dyri2,sna) 88 c go to 999 89 c 74 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 90 75 c The following looped call sets up the free solution convolved with the 91 76 c load function q hat. Note that the returned vector set "cinner" is the … … 97 82 c elliptical cross section. Note loops 8500,8000 and 9500,9000 for the two 98 83 c cases, respectively. 99 c 84 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 100 85 49 go to (8499,9499), iedge 101 c8499 tspan = trsl(irsl)102 c trsl was nondimesionalized103 c in the routine distme.f104 86 8499 do 8000 ik = 1, nhank 105 87 xakap = zksamp(ik)*alphap … … 111 93 bcin(ik) = cinner(ik) * dbesj1(xakap) 112 94 8000 continue 113 c write(6,133) diku,xakap,pset(1),pset(7),alphap,urat 114 c 133 format(' last diku xakap in what0'/1h ,1p6d16.8) 115 c write(6,1) bcin(nhank),zksamp(nhank),fltng 116 c 1 format(' last bcin zksamp and fltng in what0'/1h ,1p3d16.8) 117 c write(6,2) cinner(nhank),cinner(1),qjadon 118 c 2 format(' last and first cinner in what0 and qjadon'/1h ,1p3d16.8) 119 c 95 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 120 96 c "ojrule.f" computes the inverse Hankel trasform with a simple 121 97 c Simpson's rule. The routine "ojrule" is buliding a set of solutions stored 122 98 c in common "blocks" in r or "asrpos(nrv) ", and computed rate or displacement 123 99 c for each of N3G disks in "aswokm(nrv,N3G)" . 124 c 100 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 125 101 call ojrule(dk,bcin) 126 102 go to 999 127 c 9499 tspan = trsl(irsl)103 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 128 104 9499 do 9000 ik = 1, nhank 129 105 xakap = zksamp(ik)*alphap … … 137 113 1 - dcos(xakap) ) 138 114 9000 continue 139 c write(6,1) bcin(nhank),zksamp(nhank)140 115 call ojrule(dk,bcin) 141 c call dvecpr(cinner,nhank,'cinner at time t',79,0,0)142 c call dvecpr(bcin,nhank,'bcin at time t',79,0,0)143 116 999 return 117 144 118 end
Note:
See TracChangeset
for help on using the changeset viewer.