App-SeismicUnixGui
view release on metacpan or search on metacpan
lib/App/SeismicUnixGui/fortran/src/mmodpg4L_SU_Aug27_20_emilio/mmodpg.for view on Meta::CPAN
endif
c
call read_par_r4('MINIMUN DISTANCE (KM)?',xmin)
call read_par_r4('MAXIMUN DISTANCE (KM)?',xmax)
call read_par_r4('MINIMUN TIME (SEC)?',tmin)
call read_par_r4('MAXIMUN TIME (SEC)?',tmax)
c **** RE-CHECK PARAMETERS ****
ID=0
call read_par_i4('1- RECHECK PARAMETERS,0- NO',ID)
if(ID.eq.1) go to 15
c
c Transformation Matrix between data array and world coordinates
c
rvinvd = 0.0
if(idred.eq.0) rvinvd = rvinv
c *******************
tr(1) = datax1 - datadx
tr(2) = datadx
tr(3) = 0.0
tr(4) = 0.0
tr(4) = datat1 - datadt - tr(1)*rvinvd
! tr(5) = 0.0
tr(5) = -datadx*rvinvd
! tr(6) = datadt
c
c ***** READ VELOCITY DEPTH MODEL ******
c
call READMMOD(VT,VB,DZ,VST,VSB,RHOT,RHOB,nl)
c
call read_par_i4('Working layer number ?',lch)
if(lch.lt.1) lch = 1
if(lch.gt.nl) lch = nl
c
c ***************** Open X-T Window ****************
c
call pgbegin(0,'?',1,1)
c call pgpaper(13.5,0.5)
c call pgpaper(12.0,8.5/11.0)
call pgpaper(16.0,0.55)
c call pgpaper(10.0,0.5)
call pgask(flag)
************************************
c
10 continue
c
c Check for vt-vb too small.
c
A3 = 0.001
c
do 20 i = 1, nl
20 if(ABS(VT(I)-VB(I)).le.A3) vb(i) = vt(i)
c
c Reflections at the bottom of layers decided automatically
c based on velocity discontinuities
c
do 34 I=1,nl-1
IR(I)=0
34 if(ABS(VT(I+1)-VB(I)).GT.A3) IR(I)=1
IR(nl) = 0 !** No reflection at the bottom of model
c
c *** COMPUTATIONS ***
c
DZ1TEM=DZ(1)
DZ(1)=DZ(1)-(SDEPTH+RDEPTH)/2. ! Correct only if 1rst layer
c is a constant vel. layer.
call txpr(nl,VT,VB,DZ,PMIN,PMAX,DP,IR,multin,ntp,ILA,P,X,T)
DZ(1)=DZ1TEM
c write(*,*) 'total # of computed points ',ntp
c
c *** PLOTTING ***
c
55 continue
c call pgvport(0.1,0.9,0.1,0.9)
call pgvport(0.075,1.0,0.08,0.925)
call pgwindow(xmin,xmax,tmax,tmin)
call pgbox('BCTN',0.0,0,'BCTN',0.0,0)
call pglabel('X(km)','Tred (sec)',
+ 'Primary P-wave T-X Curves, 1-D Model')
c
if(idrdtr.eq.1) then
call pggray(f,ntrmax,nsmax,1,ntr,1,ns,cmax,cmin,tr)
endif
c
do 140 j=1,lch
nplj = 0
do 120 i=1,ntp
c ** SELECT LAYER J **
if(ABS(ILA(I)).NE.J) go to 120
nplj = nplj+1
xa1(nplj) = X(I)
xa5(nplj) = -X(I)
xa2(nplj) = T(I) - X(I) * rvinv !** RED. TIME
xa3(nplj) = P(I)
xa4(nplj) = T(I) - X(I) * P(I) !** TAU
120 continue
if(nplj.eq.0) go to 140
c
c ** PLOT LAYER J **
c
icolor = icolor + 1
if(icolor.gt.15) icolor = 1
call pgsci(icolor)
c call pgsci(3)
c
c ** X - T Plot **
c
call pgline(nplj,xa1, xa2)
call pgline(nplj,xa5, xa2)
c
c ** TAU - P Plot **
c
c call gks$polyline(nplj, xa3, xa4)
c
140 continue
c
c Draw digitized X-T data
c
if(idrxy.eq.1) then
c Positive offsets branch
if(ndxyp.gt.0) then
do ixy = 1,ndxyp
lib/App/SeismicUnixGui/fortran/src/mmodpg4L_SU_Aug27_20_emilio/mmodpg.for view on Meta::CPAN
call read_par_r4('Multiplicative constant ??',vmf)
do i = 1,nl
VT(i) = vmf * VT(i)
VB(i) = vmf * VB(i)
enddo
endif
c
c **************************
c
icolor = 1
call pgpage
call pgsci(1)
go to 10
c
255 continue
c
c write modified model to terminal
c
write(*,*) ' '
write(*,*) 'MODIFIED MODEL:'
call WRIMOD2(nl,VT,VB,DZ,VST,VSB,RHOT,RHOB)
c
c write modified model to file mmodpg.out
c
OPEN(UNIT=IOUT,FILE='mmodpg.out',STATUS='UNKNOWN',
+ FORM='UNFORMATTED')
do K=1,NL+1
write(IOUT) VT(K),VB(K),DZ(K),
+ VST(K),VSB(K),RHOT(K),RHOB(K)
enddo
CLOSE(UNIT=IOUT)
c
write(*,*) 'This model has been written to file:'
write(*,*) '*** mmodpg.out ***'
c
call pgend
500 format(' 1- VTOP = ',f6.3,', 2- VBOT = ',f6.3,' (km
+/s)')
502 format(' 3- DZ = ',f6.3,' (km)',', 4- VTOP and VBOT')
503 format(' 7- Increment = ',f7.4,' (km or km/s)')
504 format(' 9- Clip for data, 10- DP = ',f9.6,
+' (s/km)')
end
c ****************************************************
c
SUBROUTINE rdata(f, m, n, ntr, ns, datadt, fmin, fmax)
INTEGER m,n,ntr,ns,idtusec
REAL f(m,m), fmin, fmax
CHARACTER*40 FIN
LOGICAL EX
c
c Read data parameters
c
iin=26
OPEN(UNIT=iin,FILE='parmmod',STATUS='OLD')
read(iin,*) ntr,ns,idtusec
close(iin)
datadt = float(idtusec) * 1e-6
c
c115 write(*,*) 'INPUT DATA FILE NAME ?? '
c READ(5,'(A)') FIN
c INQUIRE(FILE=FIN,EXIST=EX)
c if(.NOT.EX) then
c write(*,*)'FILE DOES NOT EXIST, TRY AGAIN WITH A NEW NAME'
c go to 115
c endif
c
c OPEN(UNIT=iin,FILE=FIN,STATUS='OLD',FORM='UNFORMATTED')
c
c Read data File
c
OPEN(UNIT=iin,FILE='datammod',STATUS='OLD',FORM='UNFORMATTED')
k=1
120 READ(iin) (f(k,i), i=1,ns)
if(k.GE.ntr) go to 125
k=k+1
go to 120
125 CLOSE(UNIT=iin)
c
fmin = 1e30
fmax = -1e30
do 20 i=1,ntr
do 10 j=1,ns
fmin = min(f(i,j),fmin)
fmax = max(f(i,j),fmax)
10 continue
20 continue
c
write(*,*) 'Data min, max = ',fmin,fmax
write(*,*)
c
END
c
c ************************************************
c
subroutine read_datxy(x,y,n,fin,iin,iwrit)
c
c Reads file containing (x,y) pairs and an arbitrary number of comment lines
c in between. The (x,y) pairs are returned in arrays x and y. The
c comment lines are written on the terminal.
c
c n = number of (x,y) pairs in the file (output).
c
c fin = Default input file name (input)
c iin = reading input unit (input).
c iwrit = if iwrit.eq.1, (x,y) pairs are written on screen (input).
c
dimension x(*),y(*)
character*40 fin
character*40 comment
LOGICAL EX
c
c ********* read input file *********
c
c Check for default file "fin". If file does not exist, then ask for
c a file name
c
c write(*,*) fin
go to 117
115 write(*,*) 'Input File Name ?? '
( run in 1.046 second using v1.01-cache-2.11-cpan-8f98c5d2c55 )