view release on metacpan or  search on metacpan

lib/App/SeismicUnixGui/fortran/src/immodpg.for  view on Meta::CPAN

! Open X-T Window ****************
!       print *,' 1. start '
!       call pgbegin(0,'?',1,1)  !ask what device to use

	   call pgbegin(0,' ',1,1)  ! use default device

!  Overall terminal size 
!  in inches width_in,aspect 11 " wide and 8/11 high
	   call pgpaper(10.75,0.75)
!	   call pgpaper(16.0,0.55) 
!  XLEFT_bot, XRIGHT_TOP, YBOT, YTOP)
!	   call pgvport(-5.,1.,0.,1.)
!      call pgsvp(0.0,0.5,0.1,1.0)
       call pgask(flag)
       
       a2_prior = -1.0
       a1_prior =  0.
       do 8  i = 1,current_layer_number
              k = 2*i - 1
              va_prior(k) = vt(i)
              if(vt(i).gt.a2_prior) a2=vt(i)
              za_prior(k) = a1_prior
              va_prior(k+1) = vb(i)
              if(vb(i).gt.a2) a2=vb(i)
              a1_prior = a1_prior + dz(i)
              za_prior(k+1) = a1_prior
8      continue
! End of setup
! user sees an EMPTY BLACK SCREEN
!       print *,' 1. start of principal input loop'
10	   continue

************************************
! START OF ALL INTERACTIONS WITH THE USER
!       print *,' L511 start of interaction with the user'

! Check for Vtop-Vbottom too small.-km/s
	    A3 = 0.00001
!       print*, '1D. immodpg.for,made it, current_layer_number=',
!     +  current_layer_number

       do i = 1, current_layer_number

	   if(ABS(VT(i)-VB(i)).le.A3) vb(i) = vt(i)

       end do

! Reflections at the bottom of layers decided automatically
! based on velocity discontinuities
	   do i=1,nl-1
         IR(i)=0 
    	 if(ABS(VT(i+1)-VB(i)).GT.A3) IR(i)=1     
        end do 
        
! No reflection at the bottom of model
	   IR(nl) = 0  
	
!	print*,'reflection at layer i,=',i,IR(i)

!      COMPUTATIONS ***
	   DZ1TEM=DZ(1)

!      Correct only if 1st layer is a constant vel. layer.
	   DZ(1)=DZ(1)-(SDEPTH+RDEPTH)/2.

	   call txpr(nl,VT,VB,DZ,PMIN,PMAX,DP,IR,multin,ntp,ILA,P,X,T)
	   DZ(1)=DZ1TEM

55	   continue

! PLOTTING and REPLOTTING when correct option turns
       call pgpage ! clear screen
! for pgvport
! left_bottom_X, right_top_X,left_bottom_Y,right_top_Y
! 0.7= Xfraction of black screen occupied by seismic data
       call pgvport(0.15,0.75,0.15,0.9)
! increase vertical exaggeration by decreasing the 
! second number -- original at 0.75-- to a smaller value
       call pgvport(0.15,0.6,0.15,0.9)
!      call pgvport(0.1,0.9,0.1,0.9)
!      character size, in proportion to 1.0
       call pgsch(1.)
       call pgwindow(xmin,xmax,tmax,tmin)
       call pgbox('BCTN',0.0,0,'BCTN',0.0,0)
       call pglabel('X(km)','Tred (sec)',
     + 'Body wave, T-X Curves, 1-D Model')
!
! plot seismic image in a window
       if(idrdtr.eq.1) then
!        print *, 'L561, clip_min,clip_max=',clip_min,clip_max
!        print *, 'L562, ns,ntr=',ns,ntr
!        print *, 'L563, ntrmax,nsmax=',ntrmax,nsmax  
!        clip_min = -1.00000
!        clip_max = 1.00000
	   call pggray(Amp,ntrmax,nsmax,1,ntr,1,ns,clip_max,clip_min,tr)
!         print *, 'L543, plot seismic image in window;ns=',ns 

       endif
!
! plot modelled arrivals by layer
       do 140 j=1,current_layer_number
         nplj = 0
         array_ntp(j) = ntp

! plot first breaks in a layer
	   do 120 i=1,ntp

! SELECT LAYER J **
! nplj - number of points in layer j
	    if(ABS(ILA(i)).NE.J) go to 120
	       nplj = nplj+1
	       xa1(nplj)  =   X(i)
	       xa2(nplj)  =   T(i) - X(i) * rvinv_kmps !** RED. TIME
	       xa3(nplj)  =   P(i)
              xa4(nplj)  =   T(i) - X(i) *  P(i) !** TAU
!             get ready to output all the data
              xout(nplj,j) = X(i)
              tout(nplj,j) = T(i)

!        write(*,*)'*** test output *****************************'