-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathEXSIM12.for
4282 lines (3603 loc) · 184 KB
/
EXSIM12.for
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
program EXSIM12 !! Released in May 3, 2012, Modified version of EXSIM_dmb by G. Atkinson and K. Assatourians
c EXSIM is a stochastic finite-fault program, following the methodology of:
c Motazedian,Dariush and Gail M. Atkinson, (2005)
c "Stochastic Finite-Fault Modeling Based on a Dynamic Corner Frequency",
c Bulletin of the Seismological Society of America,95,995-1010.
c EXSIM has been modified several times since its original development.
c This version was modified from "EXSIM_dmb", where
c EXSIM_dmb is a version of EXSIM that was modified by D. Boore as described in Boore (2009 BSSA)
c Boore's modifications addressed some consistency of terminology issues, and improved low-frequencies.
c
c This EXSIM is designed to work for a list of observation sites for a
c user-defined fault geometry and location. Reads the input parameters and
c outputs acceleration time series, and the average Fourier
c Spectra and PSA over trials. It also generates averages.
c
c An example input file of program parameters is given is Example_EXSIM11.par
c Other input files for crustal amplifications, site amplifications and empirical filter may also be needed.
* Dates: Original EXSIM 2003-2005: based on modifications to FINSIM
* FINSIM written by Beresnev and Atkinson (1997 SRL; also BSSA)
* FINSIM modified by EXSIM by Motazedian and Atkinson (2005 BSSA)
* EXSIM modified by DM Boore as per log file below
* 08/17/08 - Modifications by D. M. Boore
* 08/17/08 - Corrected typo "ISlipWeigth" and "falg", and allow all
* unit slip weights if flag -1.
* 08/18/08 - Read in y (=vrup/beta), write R in PSA_FA output
* Revise column headers of PSA_FA output.
* Lots of changes to output, including writing
* separate PSA_FA files for each site (to make it
* easier to import into a graphics program).
* Allow the site coordinates to be entered as
* isitecoordflag coords
* 1 lat,long
* 2 R, Az
* 3 N, E
* 08/19/08 - Reverse order of fault input, and allow use of Wells
* and Coppersmith if FL=0 or FW=0, put magnitude input before
* FL, FW (because the later needs amag). Also get fault
* type (needed for Wells and Coppermsith).
* 08/20/08 - Change averaging algorithm (log avg for PSA, E(FAS^2)^1/2 for FAS)
* Return E(FAS^2)^1/2, not log of the value, out of the
* binning routine SAMPLE, called by compute FACCN.
*
* Add input variable iflagscalefactor to determine whether to use
* scale factor based on the following:
* iflagscalefactor integrand of integral
* 1 velocity spectrum squared
* 2 acceleration spectrum squared
* 3 D. Boore's factor based on acc. sq. high-frequency asymptotes
* Add input variable to determine if use arithmetic or
* geometric average for PSA
* 08/21/08 - Added scalefactor (h) to computeStochasticWave output and print out
* the scalefactor and corner frequency for each subfault for
* the first simulation. Also include an option for tapering the
* scaling factor.
* 08/22/08 - Write out the time series for the first realization only.
* 08/25/08 - Read in qmin
* In checkFFTpoints, increased upper index so that 2**n=32768.
* Also, doubled dimensions of wave, totalWave, and subwave.
* In many subroutines replaced dimensions of arrays passed
* in through the parameter list with "(*)".
* 08/27/08 - Added dur_sub to common /sub/, in order to print it out at some point.
* 11/06/08 - Add option to compute RMS of PSA, and also use same options for computing the average of FAS
* (for study purposes).
* 11/07/08 - Changed way that averages are computed. Now the accumulated sum is computed
* in computeAverageResponse and computeAverageFA and this is used in the main program
* to compute an average. But a modification was needed to allow
* for negative values of log10(y) for geometric means. I set the average to -9.99
* if accum_fas = 0.0 or accum_psa = 0.0.
* 11/08/08 - I may be having problems with the rms average of FAS, so compute this differently now.
* 11/15/08 - Increase maximum number of nl, nw to 200.
* 11/25/08 - I am concerned about the calls to ran1. Note that the distributed version of
* exsim uses urand, not ran1. Apparently I susbtituted ran1 for urand, but I think
* I made a mistake. If used to
* generate a series of "random" numbers, ran1 should be called the
* first time with a negative integer (idum) as a seed.
* Ran1 performs some initializations, including assigning
* values to Saved variables iv and iy, as well as changing
* idum. In exsim ran1 is called for a number of things,
* possibly including locateHypocentreRandomly,
* createRandomWeight, computeStochasticWave, and
* computing a random delay for the time series from
* each subsource. I would think that the idum, iv,
* and iy values would need to be in sync with one another,
* but this will not be the case if ran1 is called with
* a different value for idum that is not the same as in
* previous calls. In the previous version of exsim, iseed
* is used as idum in the calls in all calls except for
* the call for computing the random delay. But iseed was
* not constrained to be less than 0.0 initially, so the ran1
* initialization was probably not done. And also there was
* one call to ran1 that used islipweight for idum, which
* could have values less than, equal, to, or grater than 0.0
* (in my runs it was lsess than 0.0). When less than 0.0, as
* it was in all of my runs, the iv and iy would have been
* initialzed and would have been used in the next call to
* ran1, but the idum would not have been the reset value
* of islipweight. I do not know the consequences of all of
* this. I have changed exsim to contrain iseed to be less
* than 0 initially and have used iseed in all calls to ran1.
* 11/27/08 - Added duration calculation
* 11/28/08 - Use 1/f0 for risetime
* 11/28/08 - Change input file to have a comment line before the name of the scale
* factor file.
* 11/30/08 - Move computation of random slipweights
* if Islipweight .eq. 1.0 from main program to getInputParameters
* (a more logical place because the slipweights array is filled
* in the same place for all values of ISlipWeight).
* 11/30/08 - Move "findSubFaultMoment" out of loop over isite
* 12/01/08 - Use "if then" in loop over subsources in
* finding NumberofActiveSubs
* 12/01/08 - Remove "NumberofActivesSubs" from the arguments passed
* to subroutines computeStochasticWave and sourceSpectra, as
* the variable is not used.
* 12/01/08 - Added writing of istart, istop, nptsTotalWave, maxPointsOfWave
* in WritePar (I'm trying to figure out if findindex is needed and if
* code setting istop is correct).
* 11/28/08 - Use 1/firstElementF0 for risetime
* 12/01/08 - Use original exsim risetime
* 12/01/08 - Add option to choose one of three ways of obtaining risetime
* 12/06/08 - Allow iteration of random hypocenters for each station
* 12/08/08 - Read initial seed.
* 12/09/08 - Add scaling_coefficient, following Dariush's suggestion.
* 12/09/08 - Increase length of output file names.
* 12/17/08 - Move computations that do not depend on the loop over
* nsims out of the loop. This required defining f0, etc as
* arrays (one value per subsource). This is computationally efficient.
* Another advantage is that
* I can calculate the maximum duration for the total wave before doing
* the subsource wave calculation, and therefore I can use dynamic
* allocation for the time series array dimension. This required
* changing common sub not to include arrays of variables that
* depend on the subsource, passing these things through argument lists.
* I renamed common /sub/ to common /params/. I also removed dur_sub
* from this common block.
* I also compute dur_sub outside of computeStochasticWave and pass the
* value into that subroutine through the argument list.
* 12/17/08 - Major revision involving the time series (use tpads and front and back, do not truncate
* the motions from each subsource; use dynamic allocations). NOTE:
* Using Mavroeidis and Papageorgiou has not been tested, and it will probably not work
* because I have not taken into account the added npadl. I comment out a section of code
* before the call to M&P that may have to do with determining the index for which the pulse
* starts. It would be easy to fix this later using t_arrive(i,j).
* 12/22/08 - If nl = nw = 1, over-ride the computation of a random hypocenter and
* set n_hypocenters, i0, j0 = 1, 1, 1
* 12/31/08 - Write sqrt(nl*nw)*fas to help in deciding if that factor can correct for
* the incoherent summation.
* 12/31/08 - replace "np" for the taper correction at low frequencies by "ftaper_f0",
* the ratio of ftaper and f0 read from the input file. Also, there is no
* need for "taper_scalefactor", because ftaper_f0 .eq. 0.0 implies
* no tapering toward low frequencies.
* 01/04/09 - Apply a new taper that attempts to correct for the incoherent summation,
* the decay of the subsource spectra for frequencies less than the subsource
* corner frequency (following Frankel, 1995), and the improper long-period spectral
* level.
* 01/06/09 - Remove "scaling_coefficient" from input and program.
* - Put computation of c4taper, fc4taper outside of the sourcespectra function.
* - Revise format of input, using a header line before each line containing input.
* - Include a low-cut filter.
* - Compute average pga and pgv.
* 01/10/09 - Increase max number of frequencies to 500, without changing anything else.
* 01/13/09 - "apply_taper" was read in from the parameter file, but it has not been used since
* I added the Frankel-like scale factor (on 01/04/09). Thus I have formally removed it from this
* version of the program.
* - Add SD, PSV to FAS, PSA output file and change order of FAS and PSA.
* 02/06/09 - Compute and write out husid time series and 95%-5% duration
* - Change input to be closer to smsim input
* - Read only output file stem name
* 02/14/09 - Correct error in setting i0, j0 if n_hypocenters = 1 and nsites > 1
* 02/16/09 - Change input, use stress scaling to adjust Wells and Coppersmith FL, FW
* - Compute RJB, RCD, print all distances for each site.
* - Include an option to move the sites to positions relative to
* the midpoint of the fault or to the tip of the fault
* if move_site .eq. .true. and isitecoordflag .eq. 2 (r and epi)
* (reset isitecoordflag to 3 in this case)
* or if move_site .eq. .true. and isitecoordflag .eq. 3 and one of the
* site locations = 0.0
* - Moved writePar before loops, and just after call to getInputParameters.
* 02/17/09 - Modify findDistanceAndAzimuth subroutine
* - Modify findSubfaultDistance subroutine
* - Change FaultStrikeDeg to FaultStrike
* - Delete convertDegreeToRadian subroutine
* - Delete findDistance subroutine
* - Delete shortestSubfaultDistance subroutine
* - Read in dl, dw, and compute nl, nw
* 02/18/09 - Replace calls to get_time with a call to the system routine date_time
* 02/23/09 - Hardwire in choice of 75% rather than 95% in determining the duration.
* 75% is recommended by Ou and Herrmann (SRL) and is consistent with a few
* Husid plots I made.
* 05/28/09 - Specify a hypocenter location by distances along the fault and down the dip.
* - Add period column to psa_fs output.
* 12/05/09 - Allow for up to a six-segment path duration (Atkinson and Boore
* use a four-segment function).
* The form of the input matches that in the smsim programs.
* 12/06/09 - Pass path duration parameters through calling arguments rather than common.
* - Obtain number of frequencies in the crustal and site amp files from the
* first line of the files (not from the params file).
* - Move the slip weight params before the site coord params and obtain slip weights
* from a text file (if slip weights are to be specified).
* - Modify order of site coord params so that it is easy to edit if do not want
* to do sims at all sites in the list (just change nsites; the list of sites need
* not be edited if the first nsites are those for which simulations are desired).
* - changed "y" to "vrup_beta" and "v" to "vrup" (the former names makes it very
* difficult to search for those variables).
*
* May/2012 - D.M. Boore version modified by Assatourians to clean up, with some changes in input conventions
* These were mostly changes to simplify inputs and outputs, provide additional flexibility for outputs,
* allow use of an empirical filter with variability to modify motions; changes noted throughout
* in a commented version of this code that we have (EXSIM11) but not all such comments are retained below.
* use dflib !Modified/added by KA, Aug.2010
real wave(:), totalWave(:), vel_total(:), husid(:)
allocatable :: wave, totalWave, vel_total, husid
dimension slipWeights(200,200),
: slat(100),slon(100),
: dist(100),
: accum_fas(500), accum_psa(500),
: averagelogPSA(500),averagePSA(500),
: averagesqFA(500),averageFA(500),nnFA(500)
dimension weightedMoment(200,200)
dimension freq(500),psa(500),fa(500),fh(500),hd1(500),hd2(500)
real subfaultMoment(200,200), f0(200,200), risetime(200,200)
integer NumberOfActiveSubs(200,200)
real subfaultDistance(200,200), actualSlip(200,200),
: dur_sub(200,200), dur_path(200,200),
: t_arrive(200,200), delay(200,200), freq_out(10),
: freq4intrp(500), psa4intrp(500)
real dur_75_05, arias
real alogFASsum(500), alogPSAsum(500)
integer nFASsum(500), nPSAsum(500)
real avgavgFAS(500), avgavgPSA(500), avgavgPSA_out(10)
character f4head(10)*5
character *3 chr3, chr3_
real amag, r_ref, rlow(10), a_s(10), b_s(10), m_s(10)
real stutter
dimension afreq1(500), amp1(500), afreq2(500), amp2(500)
dimension inormal(10), xrandom1(10), xrandom2(10),
* xrmin(10), xrmax(10)
dimension siteLocation(300,2)
character f_crustal_amps*30,f_site_amps*30,f_empirical_amps*30 !Modified/added by KA, Aug.2010
character InputFileName*256
character f_stem*120
character fpar*120, facc*120, fpsa*120, f_h_f0*120, fhusid*120,
: f_times*120, f_dur_75_05*120, f_dist_psa*120
character ftmp*120 !Modified/added by KA, APR.2012
logical f_exist, rmv_trend, specify_length, specify_width,
: move_site, write_site_files
character fault_type*10
character datx*8, time_start*10, time_stop*10
character*30 fparrandom
real rpathdur(6), pathdur(6)
integer(2) n_ofargs !Modified/added by KA, Aug.2010
integer(4) i_system !Modified/added by KA, Feb.2012
common /params/
: rho,beta,iKapa,fmax,
: Qmin,Q0,Qeta, !Modified/added by KA, Mar.2012. Return to Q0*F**eta definition
! : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, !Modified/added by KA, Mar.2012. Return to Q0*F**eta definition
: iwind,eps,eta, !Modified/added by KA, Mar.2012, eps/eta added
: n_crustal_amps,n_site_amps,n_empirical_amps,totalMoment, !Modified/added by KA, Aug.2010
: c4taper, fc4taper,
: nsubs, flocut, nslope, pi, twopi
common /par/ iFFTsub, dt, npadl, npadt
common /seed/ iseed
common /fnames/ f_crustal_amps,f_site_amps,f_empirical_amps !Modified/added by KA, Aug.2010
integer soFarNumberOfSubs
integer nptsTotalWave
c nptsTotalWave is the total number of points in the final timeseries.
nu_ctl = 99
ioPar = nu_ctl - 1
ioPSA = nu_ctl - 2
ioHUS = nu_ctl - 3
ioACC = nu_ctl - 4
ioresp1 = nu_ctl - 5
ioresp2 = nu_ctl - 6
ioa = nu_ctl - 7
nu_h = nu_ctl - 8
nu_i0j0 = nu_ctl - 9
nu_times = nu_ctl - 10
nu_dur_75_05 = nu_ctl - 11
nu_dist_psa = nu_ctl - 20
iotmp=nu_ctl-25 !Added by Karen for getting tmp output in one file (20104030)
ioaccall=nu_ctl-30 !Added by Karen for getting acc output in one file (20101015)
open(unit=ioaccall,file='ACC.out',status='unknown') !Added by Karen for getting acc output in one file (20101015)
* Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME
call DATE_AND_TIME( datx, time_start )
* Date is returned as 'CCYYMMDD'
* Time is returned as 'hhmmss.sss'
pi=4.0*atan(1.0)
twopi = 2.0 * pi
n_of_args=nargs() !Modified/added by KA, Aug.2010
f_exist = .false.
do while (.not. f_exist)
InputFileName = ' '
if(n_of_args.gt.1)then !Modified/added by KA, Aug.2010
call getarg(1,InputFileName) !Modified/added by KA, Aug.2010
else !Modified/added by KA, Aug.2010
n_of_args=0 !Modified/added by KA, Aug.2010
write(*, '(a\)')
: ' Enter name of input parameter file '//
: '(cr = EXSIM12.params): ' !Modified/added by KA, Apr.2012
read(*, '(a)') InputFileName
endif !Modified/added by KA, Aug.2010
if (InputFileName(1:4) .eq. ' ') !Modified/added by KA, Apr.2012
: InputFileName = 'EXSIM12.params' !Modified/added by KA, Apr.2012
call trim_c(InputFileName, nc_f_in)
inquire(file=InputFileName(1:nc_f_in), exist=f_exist)
if (.not. f_exist) then
write(*,'(a)') ' ******* FILE DOES NOT EXIST ******* '
end if
end do
call check_directories() !Modified/added by KA, Mar.2012: Check/creat directories for outputs
open(unit=nu_ctl, file=InputFileName(1:nc_f_in), status='unknown')
print *, ' Control file: '//InputFileName(1:nc_f_in)
c Read (and echo) the input parameters that are common to all grid points.
call getInputParameters(nu_ctl,
: FaultStrike,FaultDip,h,
: FaultLat,FaultLon, move_site, write_site_files,
: siteLocation,numberOfSites,isitecoordflag,
: fault_type,
: FaultLength, FaultWidth,
: dl, dw, nl,nw, nsubs,
: specify_length, specify_width, stress_ref,
: vrup_beta,
: hyp_loc_fl, hyp_loc_fw, i0_in,j0_in,n_hypocenters,
: tpadl, tpadt, dt, beta, rho, amag,
: stress,iKapa,fmax,
! : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, !Modified/added by KA, Mar.2012. Return to Q0*F**eta definition
: Qmin,Q0,Qeta, !Modified/added by KA, Mar.2012. Return to Q0*F**eta definition
: r_ref, nsprd_segs, rlow, a_s, b_s, m_s,
: rpathdur, pathdur, durslope, ndur_hinges,
: iwind,eps,eta, !Modified/added by KA, Mar.2012, eps/eta added
: flocut, nslope,
: nfreq,freq1,freq2,
: nfout, freq_out,
: f_stem,
: f_crustal_amps,
: f_site_amps,
: f_empirical_amps, !Modified/added by KA, Aug.2010
: iseed, nsims,damp,
: islipweight, slipWeights,
: iRow,jColumn,pulsingPercent,
: iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak,
: iDynamicFlag,
: iflagscalefactor,
: iflagfas_avg, iflagpsa_avg,
: i_rise_time)
close(nu_ctl)
vrup=vrup_beta * beta
subFaultRadius=sqrt((dl*dw)/pi)
!%%%%%%%%%%
riseTime_original=subFaultRadius/vrup
!%%%%%%%%%%
call trim_c(f_stem, nc_f_stem)
fpar = ' '
fpar = f_stem(1:nc_f_stem)//'_parameters.out'
call trim_c(fpar, nc_fpar)
* Initialize q_f (freq-dependent Q model)
dummy = q_f_setup(Qmin,Q0,Qeta) !Modified/added by KA, Mar.2012. Return to Q0*F**eta definition
! dummy = q_f_setup(fr1, qr1, s1, ft1, ft2, !Modified/added by KA, Mar.2012. Return to Q0*F**eta definition
! : fr2, qr2, s2, c_q) !Modified/added by KA, Mar.2012. Return to Q0*F**eta definition
* Initialize gsprd (geometric sprading)
numsource = 1
dummy = gsprd_f_setup(r_ref,nsprd_segs,rlow,
: a_s,b_s,m_s,
: numsource,amag)
npadl = tpadl/dt
npadt = tpadt/dt
c Generate some logarithmically-spaced frequencies
freq(1)=freq1
if (nfreq.ge.2) then
frinc=alog(freq2/freq1)/(nfreq-1)
do k=2,nfreq
freq(k)=freq1*exp((k-1)*frinc)
enddo
endif
c End of Generate some logarithmically-spaced frequencies
iseed = -iabs(iseed)
if(i0_in.eq.0.or.j0_in.eq.0)
: call locateHypocentreRandomly(i0,j0,nl,nw)
NoOfEffectiveSubfaults=float(nl)*(pulsingPercent/100.0)
NoOfEffectiveSubfaults=NoOfEffectiveSubfaults/2. !!double side propogation
if(NoOfEffectiveSubfaults.le.1) NoOfEffectiveSubfaults=1
totalMoment=10.**(1.5*amag+16.05)
avgSubFaultMoment=totalMoment/float(nw*nl)
firstElementF0=
: 4.9e+6*beta*(stress/avgSubFaultMoment)**(1.0/3.0)
F0main=4.9e+6*beta*(stress/totalMoment)**(1.0/3.0)
call findSubFaultMoment(nl,nw,slipWeights
: ,totalMoment,weightedMoment)
f4head = '00000'
do i = 1, nfout
if (freq_out(i) .lt. 0.0 .or. freq_out(i) .ge. 10.0) then
write(f4head(i)(1:5),'(f5.2)') freq_out(i)
else
write(f4head(i)(2:5),'(f4.1)') freq_out(i)
end if
end do
f_dur_75_05 = ' '
f_dur_75_05 = f_stem(1:nc_f_stem)//'_dur_75_05_dist_psa.out'
call trim_c(f_dur_75_05, nc_f_dur_75_05)
open(nu_dur_75_05, file='./other/'//f_dur_75_05(1:nc_f_dur_75_05), !Modified/added by KA, Feb.2012
: status='unknown')
write(nu_dur_75_05,'(2x,a,1x,a,2x,a,2x,a,1x,a,1x,a,
: 1x,a, 1x,a,1x,a,1x,a,1x,a,5x,a, 3x,a,10(7x,a4, 3x,a8)
: )')
: 'isite','ihyp','i0','j0','dur_75_05(sec)',
: 'arias(cm/s)','sitecoord(1)', 'sitecoord(2)',
: 'isitecoordflag','site_lat_degrees',
: 'site_lon_degrees', 'd_jb', 'd_cd2f',
: ('freq', 'psa'//f4head(i), i = 1, nfout)
DO isite=1,numberOfSites ! Loop on Sites around the Fault
SiteLat=siteLocation(isite,1)
SiteLon= siteLocation(isite,2)
print *, "Working on Site #",isite,SiteLat, SiteLon
* Compute some distances:
call site_coords_in_degrees(
: FaultLat, FaultLon,
: SiteLat, SiteLon, isitecoordflag,
: site_lat_degrees, site_lon_degrees )
h_min_c = 3.0 ! Campbell depth to seismogenic region
w1 = 0.0
w2 = faultwidth
s1 = 0.0
s2 = faultlength
call dist_3df(
: site_lat_degrees, site_lon_degrees,
: FaultLat, FaultLon, h, FaultStrike, FaultDip,
: w1, w2, s1, s2,
: h_min_c, d_jb, az_jb, d_cd2f, az_cd2f, d_c, az_c,
: d_sta_n, d_sta_e, irgn_cd2f, irgn_c, irgn_jb)
nnFA=0
totalDuration=0
fi2=0
call findDistanceAndAzimuth(FaultLat,FaultLon,SiteLat,
: SiteLon,R,fi2,isitecoordflag)
c Note that R, azm fi2 are the distance, azm. w.r.t. origin (not epi),
c Note also I must input faultstrike in degrees and its converted to rad.
alogFASsum = 0.0
alogPSAsum = 0.0
nFASsum = 0
nPSAsum = 0
alogPGAsum = 0.0
alogPGVsum = 0.0
DO ihypo = 1, n_hypocenters
if (n_hypocenters .gt. 1) then
call locateHypocentreRandomly(i0,j0,nl,nw)
else if (i0_in. eq. 0 .or. j0_in .eq. 0) then
call locateHypocentreRandomly(i0,j0,nl,nw)
else
i0 = i0_in
j0 = j0_in
end if
nsubfaults = 0
t_arrive_min=10000.
t_end_max=0.
risetime_max = 0.0
DO i=1,nl
DO j=1,nw
nsubfaults = nsubfaults + 1
subFaultMoment(i,j)=weightedMoment(i,j)
if(iDynamicFlag.eq.0) then
NumberOfActiveSubs(i,j)=1
else
NumberOfActiveSubs(i,j)=NumberOfPulsingSubs(i0,j0,i,j
: ,nl,nw,NoOfEffectiveSubfaults)
if(NumberOfActiveSubs(i,j).eq.0) then
NumberOfActiveSubs(i,j)=1
end if
end if
f0(i,j)=
: firstElementF0*(NumberOfActiveSubs(i,j)**(-1.0/3.0))
if (i_rise_time .eq. 1) then
risetime(i,j) = riseTime_original
else if (i_rise_time .eq. 2) then
riseTime(i,j) = 1.0/f0(i,j)
else
print *,' ERROR: invalid value of i_rise_time, = ',
: i_rise_time
stop
end if
if(risetime(i,j) .gt. risetime_max) then
risetime_max = risetime(i,j)
end if
actualSlip(i,j)=1.e-22*subFaultMoment(i,j)/
: (rho*beta**2.*dl*dw)
subfaultDistance(i,j)=
: findSubfaultDistance(R,h,FaultStrike,
: fi2,FaultDip,dl,dw,i,j)
c calculate duration of subsource time history at distance
c subfaultDistance
call dur_path_cmp(subfaultDistance(i,j),
: rpathdur, pathdur, durslope, ndur_hinges,
: dur_path(i,j)) ! computes path part of duration
dur_sub(i,j) = dur_path(i,j) + risetime(i,j)
if ((i-i0) .ne. 0 .or. (j-j0) .ne. 0) then
delay(i,j)=sqrt((dl*(i-i0))**2.+(dw*(j-j0))**2.)/vrup
else
delay(i,j)=0. ! delay from hypocenter to itself
end if
t_arrive(i,j) = delay(i,j) + subfaultDistance(i,j)/beta
if (t_arrive(i,j) .lt. t_arrive_min)
: t_arrive_min = t_arrive(i,j)
t_end= t_arrive(i,j) + dur_sub(i,j)
if (t_end .gt. t_end_max) then
imax = i
jmax = j
t_end_max = t_end
end if
END DO ! loop over nw to calculate quantities not dependent on nsims
END DO ! loop over nl to calculate quantities not dependent on nsims
nptsTotalWave=(t_end_max - t_arrive_min + tpadl + tpadt +
: risetime_max)/dt ! need risetime_max because the stutter
! delay can shift the subsource time series
! corresponding to t_end_max by an amount up to
! risetime of that subsource; to be safe I
! use risetimemax rather than risetime(imax,jmax),
! where imax, jmax are the indices of the subsource
! corresponding to t_end_max.
signnpw2 = +1.0
call get_npw2(nptsTotalWave,signnpw2,npw2TotalWave)
if (write_site_files) then
if (ihypo .eq. 1) then
f_times = ' '
f_times = f_stem(1:nc_f_stem)//'_times_s000.out'
write(f_times(nc_f_stem+9:nc_f_stem+11),'(i3.3)') isite
call trim_c(f_times, nc_f_times)
end if
end if
HypoDistance=findSubfaultDistance(R,h,FaultStrike,
: fi2,FaultDip,dl,dw,i0,j0)
c now do actual generation and summation of subfault time series ***
! Initialize arrays (index over frequency) containing cumulative sums for averages
averagelogPSA=0.0
averagePSA=0.0
averageFA=0.0
accum_fas = 0.0
accum_psa = 0.0
nnFA = 0
accum_logpga = 0.0
accum_logpgv = 0.0
accum_dur_75_05 = 0.0
accum_arias = 0.0
DO isim=1,nsims
print *, ' Site isite, Hypocenter iteration ihypo, '//
: 'Simulation number isim = ', isite, ihypo, isim
allocate (totalWave(npw2TotalWave))
allocate (husid(npw2TotalWave))
totalWave = 0.0
DO i=1,nl
DO j=1,nw
nsubsource = npadl + dur_sub(i,j)/dt + npadt
signnpw2 = +1.0
call get_npw2(nsubsource,signnpw2,iFFTsub)
allocate (wave(iFFTsub))
wave=0
call computeStochasticWave(
: wave,
: subFaultMoment(i,j), subfaultDistance(i,j),
: F0Main, f0(i,j), risetime(i,j), dur_sub(i,j),
: nl,nw,
: iflagscalefactor, hij) ! dmb
c subsource accelerogram is randomly delayed to
c simulate complexity in slip process
c subtract minimum delay to eliminate static time shift
stutter = Ran1(iseed) * riseTime(i,j)
ishift = ((t_arrive(i,j) - t_arrive_min) + stutter)/dt
if (npadl+ishift+dur_sub(i,j)+npadt .gt.
: npw2TotalWave) then
print *,'npadl+ishift+dur_sub(i,j)+'//
: 'npadt.gt.npw2TotalWave'
print *, npadl+ishift+dur_sub(i,j)+npadt,
: npw2TotalWave
end if
! No need to call shift---just figure out the shifted index and use it
! to fill the totalWave vector.
c add current accelerogram to total wave field
nhigh = min(iFFTsub, npw2TotalWave - ishift)
DO k = 1, nhigh
totalWave(k+ishift)=totalWave(k+ishift)+wave(k)
END DO
deallocate (wave)
END DO ! loop over nw
END DO ! loop over nl
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c include Mavroeidis and Papageorgiou, 2003 approach
if( int(iPapaFlag).eq.1)then
call MavroPapa(totalWave,npw2TotalWave,
* iPapaFlag,PapaGama,PapaNu,PapaT0,amag,PapaPeak)
endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c ************************** end of summation **********************
c accelerogram and model parameters are saved only once
pga=findPeakValue(npw2TotalWave,totalWave)
accum_logpga = accum_logpga + alog10(pga)
allocate (vel_total(npw2TotalWave))
vel_total = 0.0
call acc2v(totalWave, npw2TotalWave, dt, .false., vel_total)
pgv=findPeakValue(npw2TotalWave,vel_total)
accum_logpgv = accum_logpgv + alog10(pgv)
deallocate(vel_total)
rmv_trend = .false.
husid = 0.0
call accsqint(totalWave, npw2TotalWave, dt, rmv_trend,
: husid)
husid_end = husid(npw2TotalWave)
do i = 1, npw2TotalWave
husid(i) = husid(i)/husid_end
end do
call locate(husid,npw2TotalWave,0.05,j05)
call locate(husid,npw2TotalWave,0.75,j75)
accum_dur_75_05 = accum_dur_75_05 +
: alog10(float(j75-j05)*dt) ! for geometric mean
accum_arias = accum_arias +
: alog10((pi/(2.0*981.0))*husid_end) ! for geometric mean
if (write_site_files) then !Added by Karen for getting acc output of individual simulations (20120415)
facc = ' '
write(chr3,'(i3.3)')isite !Added by Karen for getting acc output of individual simulations (20120415)
write(chr3_,'(i3.3)')isim !Added by Karen for getting acc output of individual simulations (20120415)
facc = trim(f_stem)//'S'//chr3//'iter'//chr3_//'.acc' !Added by Karen for getting acc output of individual simulations (20120415)
call trim_c(facc,nc_facc)
open (ioAcc,file='./ACC/'//facc(1:nc_facc), !Added by Karen for getting acc output of individual simulations (20120415)
: status='unknown')
call writeAcc(ioAcc, totalWave, npw2TotalWave, pga, !Added by Karen for getting acc output of individual simulations (20120415)
: isim,isite,HypoDistance,d_cd2f, fpar, !Added by Karen for getting acc output of individual simulations (20120415)
: amag,siteLat,siteLon,h,r,fi2*180./pi) !Added by Karen for getting acc output of individual simulations (20120415)
close(ioAcc) !Added by Karen for getting acc output of individual simulations (20120415)
endif
! if (ihypo.eq.1 .and. isim .eq. 1) then !Modified/added by KA, Mar.2012 to include all husid time series
if (ihypo.eq.1) then !Modified/added by KA, Mar.2012 to include all husid time series
if(isim.eq.1)then
call writeAcc(ioaccall, totalWave, npw2TotalWave, pga, !Added by Karen for getting acc output in one file (20101015)
: isim,isite,HypoDistance,d_cd2f, fpar, !Added by Karen for getting acc output in one file (20101015)
: amag,siteLat,siteLon,h,r,fi2*180./pi) !Added by Karen for getting acc output in one file (20101015)
endif
if (write_site_files) then
fhusid = ' '
fhusid = f_stem(1:nc_f_stem)//'_husid_s000iter000.out'
write(fhusid(nc_f_stem+9:nc_f_stem+11),'(i3.3)') isite
write(fhusid(nc_f_stem+16:nc_f_stem+18),'(i3.3)') isim
call trim_c(fhusid, nc_fhusid)
open (ioHus,file='./other/'//fhusid(1:nc_fhusid), !Modified/added by KA to direct output in folder "other", Feb.2012
1 status='unknown')
dur_75_05_sim1 = float(j75-j05)*dt
arias_sim1 = (pi/(2.0*981.0))*husid_end
call writeHUS(ioHus, husid, npw2TotalWave,
: arias_sim1, dur_75_05_sim1,
: isim,isite,HypoDistance,d_cd2f,fpar,
: amag,siteLat,siteLon,h,r,fi2*180./pi) !Added by Karen for reporting earthquake parameters in husid files, May.2012
close(ioHUS)
end if
end if
dur=float(npw2TotalWave-1)*dt-5.0
psa=0.0
call computeResponse(totalWave,dt,dur,damp,nfreq,freq,psa)
call computeAverageResponse(accum_psa,psa,nfreq,nsims,
: iflagpsa_avg) ! dmb
FA=0.0
call computeFACCN(totalWave,npw2TotalWave,dt,nfreq,freq,FA)
call computeAverageFA(accum_fas,FA,nfreq,nsims,nnFA,
: iflagfas_avg)
write (*,"(' *** finished with trial',i4)") isim
deallocate(totalWave)
deallocate(husid)
if (write_site_files) then
ftmp = ' '
ftmp = trim(f_stem)//'S'//chr3//'iter'//chr3_//'.tmp' !Added by Karen for getting spectral output of every simulation (20120430)
call trim_c(ftmp,nc_ftmp) !Added by Karen for getting spectral output of every simulation (20120430)
open (iotmp,file='./PSA/'//ftmp(1:nc_ftmp),status='unknown') !Added by Karen for getting spectral output of every simulation (20120430)
!Added by Karen for getting spectral output of every simulation (20120430)
call writePSA_FA(iotmp, freq,PGA,PGV, !Added by Karen for getting spectral output of every simulation (20120430)
: psa,nfreq,isim,0,damp,amag, !Added by Karen for getting spectral output of every simulation (20120430)
: siteLat,siteLon,h,r,fi2*180./pi,d_cd2f,FA, !Added by Karen for getting spectral output of every simulation (20120430)
: isite, HypoDistance,fPar, nl, nw) !Added by Karen for getting spectral output of every simulation (20120430)
close(iotmp) !Added by Karen for getting spectral output of every simulation (20120430)
endif
END DO ! End loop over nsims
averagePGA = 10.0**(accum_logpga/float(nsims))
alogPGAsum = alogPGAsum + alog10(averagePGA)
averagePGV = 10.0**(accum_logpgv/float(nsims))
alogPGVsum = alogPGVsum + alog10(averagePGV)
DO i = 1, nfreq
if (iflagpsa_avg .eq. 2 .and. accum_psa(i) .eq. 0.0) then
averagePSA(i) = -9.99
else
if (iflagpsa_avg .eq. 1) then ! arithmetic
averagePSA(i) = accum_psa(i)/float(nsims)
else if (iflagpsa_avg .eq. 2) then ! geometric
averagePSA(i) = 10.0**(accum_psa(i)/float(nsims))
else if (iflagpsa_avg .eq. 3) then ! rms
averagePSA(i) = sqrt(accum_psa(i)/float(nsims))
else
print *,
: ' ERROR, iflagpsa_avg not 1, 2, or 3; = ', iflagpsa_avg
stop
end if
end if
if (iflagfas_avg .eq. 2 .and. accum_fas(i) .eq. 0.0) then
averageFA(i) = -9.99
else
if (iflagfas_avg .eq. 1) then ! arithmetic
averageFA(i) = accum_fas(i)/float(nnFA(i))
else if (iflagfas_avg .eq. 2) then ! geometric
averageFA(i) = 10.0**(accum_fas(i)/float(nnFA(i)))
else if (iflagfas_avg .eq. 3) then ! rms
averageFA(i) = sqrt(accum_fas(i))
else
print *,
: ' ERROR, iflagfas_avg not 1, 2, or 3; = ', iflagfas_avg
stop
end if
end if
if (averagePSA(i) .gt. 0.0) then
alogPSAsum(i) = alogPSAsum(i) + alog10(averagePSA(i))
nPSAsum(i) = nPSAsum(i) + 1
end if
if (averageFA(i) .gt. 0.0) then
alogFASsum(i) = alogFASsum(i) + alog10(averageFA(i))
nFASsum(i) = nFASsum(i) + 1
end if
END DO ! loop over nfreq
dur_75_05 = 10.0**(accum_dur_75_05/float(nsims))
arias = 10.0**(accum_arias/float(nsims))
if(ihypo==n_hypocenters)then
write(nu_dur_75_05,'(4x,i3,
: 2x,i3,
: 1x,i3, 1x,i3,
: 6x,f9.2, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: 1p, 1x,e10.3)',advance='no') !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: isite,
: ihypo,
: i0, j0,
: dur_75_05,
: arias
else
write(nu_dur_75_05,'(4x,i3,
: 2x,i3,
: 1x,i3, 1x,i3,
: 6x,f9.2, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: 1p, 1x,e10.3)') !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: isite,
: ihypo,
: i0, j0,
: dur_75_05,
: arias
endif
END DO !loop over random hypocenters
* Compute average of average FAS and PSA
avgavgPGA = 10.0**(alogPGAsum/float(n_hypocenters))
avgavgPGV = 10.0**(alogPGVsum/float(n_hypocenters))
DO i = 1, nfreq
avgavgPSA(i) = 10.0**(alogPSAsum(i)/float(nPSAsum(i)))
avgavgFAS(i) = 10.0**(alogFASsum(i)/float(nFASsum(i)))
END DO ! end loop over nfreq
nfreq4intrp = nfreq
do i = 1, nfreq4intrp
freq4intrp(i) = freq(i)
psa4intrp(i) = avgavgPSA(i)
end do
do i = 1, nfout
if (freq_out(i) .eq. -1.0) then ! pgv
avgavgPSA_out(i) = avgavgPGV
else if (freq_out(i) .ge. 99.0) then ! pga
avgavgPSA_out(i) = avgavgPGA
else
avgavgPSA_out(i) =
: yintrf( freq_out(i), freq4intrp, psa4intrp, nfreq4intrp)
end if
end do
fpsa = ' '
fpsa = f_stem(1:nc_f_stem)//'_psa_fs_s000.out'
write(fpsa(nc_f_stem+10:nc_f_stem+12),'(i3.3)') isite
call trim_c(fpsa, nc_fpsa)
c open (ioPsa,file=fpsa(1:nc_fpsa),status='unknown') !Modified/added by KA, Mar.2012 to include average spectra in PSA folder
open (ioPsa,file='./PSA/'//fpsa(1:nc_fpsa),status='unknown') !Modified/added by KA, Mar.2012 to include average spectra in PSA folder
call writePSA_FA(ioPSA, freq,avgavgPGA,avgavgPGV,
: avgavgPSA,nfreq,nsims,1,damp,amag, !Modified/added by KA, Mar.2012 to include average spectra in PSA folder
: siteLat,siteLon,h,r,fi2*180./pi,d_cd2f,avgavgFAS,
: isite, HypoDistance,fPar, nl, nw)
close(ioPSA)
write(nu_dur_75_05, '( !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: 5x,f9.3, 4x,f9.3, 14x,i1, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: 8x,f9.3, 8x,f9.3, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: 1x,f8.2, 1x,f8.2, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: 1p,10(1x,e10.3, 1x,e10.3))') !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: SiteLat, SiteLon, isitecoordflag, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: site_lat_degrees, site_lon_degrees, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: d_jb, d_cd2f, !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
: (freq_out(i), avgavgPSA_out(i), i= 1, nfout) !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
END DO ! End loop over sites !Modified/added by KA, Apr.2012 to join summary PSA with duration/arias intensity this feature will work correctly if non-random hypocenter is selected
c
c
c
open (ioPar,file='./other/'//fpar(1:nc_fpar),status='unknown') !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
call writePar(ioPar, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: FaultStrike,FaultDip,h, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: FaultLat, FaultLon, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: siteLocation,numberOfSites, move_site, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: fault_type, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: FaultLength, FaultWidth,nl,nw,dl,dw, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: specify_length, specify_width, stress_ref, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: vrup_beta, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: hyp_loc_fl, hyp_loc_fw, i0_in,j0_in,n_hypocenters, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: iseed, nsims, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: amag, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: stress, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: pulsingPercent,iDynamicFlag, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: i_rise_time, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: iflagscalefactor, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: iflagfas_avg, iflagpsa_avg, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: tpadl, tpadt, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: r_ref, nsprd_segs, rlow, a_s, b_s, m_s, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: rpathdur, pathdur, durslope, ndur_hinges, !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: actualSlip !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
: ) !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
!Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
close(ioPar) !Karen&Katsu2010/11/10-Moved to this section to allow all parameters to be reported properly
c
c
c
print *
* Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME
call DATE_AND_TIME( datx, time_stop )
* Date is returned as 'CCYYMMDD'