getparms.f

linux环境下利用fortran语言开发

                if(iso .EQ. 'y') then
                    niso = 0
                    write(temp1,*) niso,'     layer#',i,' is isotropic'
                    if(choice .EQ. 'L' .OR. choice .EQ. 'l') then
                    write(stderr,"(/,' enter lambda, mu (Pa) and 
     &rho (g/cm3): ',$)")
                    read(stdin,*) lambda,mu,rho

*   isotropic layer is a special case of transversely isotropic layer

                    write(temp1,*) lambda,
     &                             '     Lame constant lambda (Pa)'
                    write(temp1,*) mu,'     shear rigidity mu (Pa)'
                    write(temp1,*) rho,'     density (g/cm3)'
            rho = rho*1.e3
                    param(ii+1) = lambda+2.*mu
                    param(ii+2) = param(ii+1)
                    param(ii+3) = mu
                    param(ii+4) = mu
                    param(ii+5) = lambda
                    param(ii+6) = rho

                    else
                write(stderr,
     & "(/,' enter Vp (m/s), Vs (m/s), and rho (g/cm3): ',$)")
                    read(stdin,*) Vp, Vs, rho

                    write(temp1,*) Vp,
     &                        '     compressional velocity (m/s)'
                    write(temp1,*) Vs,'     shear velocity (m/s)'
                    write(temp1,*) rho,'     density (g/cm3)'
            rho = rho*1.e3
                    param(ii+1) = rho*Vp*Vp
                    param(ii+2) = param(ii+1)
                    param(ii+3) = rho*Vs*Vs
                    param(ii+4) = param(ii+3)
                    param(ii+5) = param(ii+1)-2.*param(ii+3)
                    param(ii+6) = rho

                    endif
*   axis of anisotropy is taken horizontal so that rays slowness 
*   and axis are in practice never aligned (Cf the way displacement
*   polarizations are computed in polar.f)

                    param(ii+7) = 1.
                    param(ii+8) = 0.
                    param(ii+9) = 0.
                else
                    niso = 1
                    write(temp1,*) niso,'      layer#',i,
     &                             ' is anisotropic'
                    if(choice .EQ. 'L' .OR. choice .EQ. 'l') then
                    write(stderr,"(/,' enter 5 elastic constants 
     &A,C,N,L,F (Pa), and rho (g/cm3): ',$)")
                    read(stdin,*) param(ii+1),param(ii+2),param(ii+3),
     &              param(ii+4),param(ii+5),rho
            param(ii+6) = rho*1.e3
              write(temp1,*) param(ii+1),'     A (C33) (Pa)'
              write(temp1,*) param(ii+2),'     C (C11) (Pa)'
              write(temp1,*) param(ii+3),'     N (C66) (Pa)'
              write(temp1,*) param(ii+4),'     L (C44) (Pa)'
              write(temp1,*) param(ii+5),'     F (C13) (Pa)'
        write(temp1,*) rho,'     density (g/cm3)'
                    else
              write(stderr,
     & "(/,' enter Thomsen parameters Vp (m/s), Vs (m/s),'
     &,/,/,5x,'gamma, delta, epsilon, rho (g/cm3): ',$)")
                    read(stdin,*) alpha,beta,gamma,delta,epsilon,rho
        write(temp1,*) alpha,'     Vp (slow P velocity) (m/s)'
        write(temp1,*) beta,'     Vs (slow S velocity) (m/s)'
        write(temp1,*) gamma,'     gamma (S velocity anisotropy)'
        write(temp1,*) delta,'     delta (anellipticity parameter)'
        write(temp1,*) epsilon,'     epsilon (P velocity anisotropy)'
        write(temp1,*) rho,'     density (g/cm3)'
            rho = rho*1.e3
                    param(ii+2) = rho*alpha*alpha
                    param(ii+4) = rho*beta*beta
                    param(ii+1) = param(ii+2)*(1.e0+2.e0*epsilon)
                    param(ii+3) = param(ii+4)*(1.e0+2.e0*gamma)
                    param(ii+5) = param(ii+2)-param(ii+4)
                    param(ii+5) = param(ii+5)*(param(ii+5)+
     &                      2.e0*delta*param(ii+2))
                    param(ii+5) = sqrt(param(ii+5)) - param(ii+4)
                    param(ii+6) = rho

                    endif
                    write(stderr,"(/,' enter azimuth and dip in degrees 
     &of anisotropy axis: ',$)")
                    read(stdin,*) azi,dip
            write(temp1,*) azi,dip,
     &'    azimuth and dip of anisotropy axis (degrees)'
                    azi = azi*deg2rad
                    dip = dip*deg2rad
                    param(ii+7) = cos(azi)*cos(dip)
                    param(ii+8) = sin(azi)*cos(dip)
                    param(ii+9) = sin(dip)
                endif

                if(i .NE. nlayer) then
50          write(stderr,"(/,' enter thickness of layer (m): ',$)")
                    read(stdin,*) thick
            geom(i) = geom(i-1)+thick
        if(thick .LE. 0.e0) then
        write(stderr,"(/,' *** invalid geometry, please try again')")
        goto 50
        endif
                    write(temp1,*) thick, '     thickness of layer#',i
                endif

        triso(i) = niso

                ii = ii+9

20          continue
          
          close(temp1)
          return
      endif

1000   write(stderr,*) 'ERROR in opening ',parfile,' iostat=',ioerror,
     & '. routine: getmodel'
      return
      end



      subroutine getoutput()

*
*   gets output parameters from standard input
*
*   OUTPUT VARIABLES:
*
*       block "output"  graphics/verbose parameters:
*
*                       verbose		selects verbose option
*                       trace		selects graphics option
*                       scale		scale for graphics
*                       xmin, ymin      lower left corner of plot
*
*       block "io"      input/output logical unit numbers:
*
*                       stdin		standard input logical unit #
*                       stdout		standard output logical unit #
*                       stderr		standard error logical unit #
*                       temp1		temporary logical unit number
*                       temp2		temporary logical unit number
*                       temp3		temporary logical unit number
*

      integer verbose, trace
      integer stdin, stdout, stderr, temp1, temp2, temp3
      real scale, xmin, ymin

      common /output/ verbose, trace, scale, xmin, ymin
      common /io/ stdin, stdout, stderr, temp1, temp2, temp3

      character itrace, iverb

      write(stderr,"(/,' *** TRISO: Ray Tracing over ',$)") 
      write(stderr,"(' Tabular Transversely Isotropic Media. ***',/)")
      write(stderr,"(' do you want a ray tracing display? (y/n) ',$)")
      read(stdin,'(A)') itrace

      if(itrace .EQ. 'y') then
        trace = 1
        write(stderr,"(/,' enter plot size in inches: ',$)")
        read(stdin,*) scale
        write(stderr,"(/,' do you want some verbose during run time? 
     &(y/n) ',$)")
        read(stdin,'(A)') iverb
        if(iverb .EQ. 'y') then
                verbose = 1
        else
                verbose = 0
        endif
      else
        write(stderr,"(/,' do you want some verbose during run time? 
     &(y/n) ',$)")
        read(stdin,'(A)') iverb
        if(iverb .EQ. 'y') then
                verbose = 1
        else
                verbose = 0
        endif
      endif

      return
      end



      subroutine getpath()

*
*   reads up to 200 raypaths
*
*   OUTPUT VARIABLES:
*
*       block "path"    raypath parameters:
*
*                       mode		sequence of modes along 
*                                       raypaths.
*                       interface	sequence of interfaces 
*                                       along raypaths.
*                       nsegment	number of segments in 
*                                       raypaths.
*                       npath		number of raypaths.
*
*   GLOBAL VARIABLES:
*
*       block "io"      input/output logical unit numbers
*

      integer nsegment(200), i, ioerror, j, nlayer, triso(100)
      real param(900), geom(0:99)
      integer stdin, stdout, stderr, temp1, temp2, temp3
      integer mode(200,100),interface(200,0:100), npath, ilayer
      common /path/ mode, interface, nsegment, npath
      common /model/ param, geom, nlayer, triso
      common /io/ stdin, stdout, stderr, temp1, temp2, temp3

      character*80 parfile
      character ifile
      character*2 char2(100), char1

      write(stderr,
     & "(/,' *** ATTENTION: NOW READING RAYPATH PARAMETERS')")

      write(stderr,"(/,' do you want to read the path from an ',$)")
      write(stderr,"(' existing file? (y/n) ',$)")
      read(stdin,'(A)') ifile

      if(ifile .EQ. 'y') then
        write(stderr,"(/,' enter name of input parameter file: ',$)")
        read(stdin,100) parfile
100     format(A80)
        open(unit=temp2,file=parfile,status='unknown',iostat=ioerror,
     &          err=1000)
          rewind(temp2)

        read(temp2,*) npath
        if(npath .GT. 200) then
        write(stderr,"(/,' *** too many paths: number set to 200.')")
        npath = 200
        endif

*   all rays must start from the surface

        do 150 j=1,npath

        interface(j,0) = 0

*   read the path from parfile

        read(temp2,*) nsegment(j)
        if(nsegment(j) .GT. 100) then
        write(stderr,
     & "(/,' *** #of segments in path#',i3,' too large,')") j
        write(stderr,"(/,'    recommend job interruption.')")
        nsegment(j) = 100
        endif
        
        do10 i=1,nsegment(j)
                read(temp2,'(1X,I3,3X,A)') interface(j,i), char1
                if( (abs(interface(j,i)-interface(j,i-1)).NE.
     &   1) .OR. (interface(j,i).EQ.0 .AND. i.NE.nsegment(j)) .OR. 
     &   (interface(j,i).GE.nlayer) ) then
                        write(stderr,"(/,' *** invalid ray path#',i3,
     &   'segment#',i3,': recommend job interruption...')") j,i
                endif

                ilayer = max0(interface(j,i),interface(j,i-1))
                if(triso(ilayer) .EQ. 0) then
                if(char1 .EQ. 'S' .OR. char1 .EQ. 's') then
                        mode(j,i) = 1
                else if(char1 .EQ. 'P' .OR. char1 .EQ. 'p') then
                        mode(j,i) = 2
                else
                        mode(j,i) = 2
                        write(stderr,"(/,' *** invalid mode: path#',
     &   i3,'segment#',i3,'; recommend job interruption...')") j,i
                endif
                else
                if(char1 .EQ. 'SP' .OR. char1 .EQ. 'sp') then
                        mode(j,i) = 1
                else if(char1 .EQ. 'QP' .OR. char1 .EQ. 'qp') then
                        mode(j,i) = 2
                else if(char1 .EQ. 'QS' .OR. char1 .EQ. 'qs') then
                        mode(j,i) = 3
                else
                        mode(j,i) = 2
                        write(stderr,"(/,' *** invalid mode: path#',
     &   i3,'segment#',i3,'; recommend job interruption...')") j,i
                endif

                endif
10         continue

        if(interface(j,nsegment(j)) .NE. 0) then
                write(stderr,"(/,' *** invalid raypath#',i3,
     &   ': recommend job interruption')") j
        endif
150        continue

        close(temp2)
        return
      else
        write(stderr,"(/,' enter name of output parameter file: ',$)")
        read(stdin,200) parfile
200     format(A80)
        open(unit=temp2,file=parfile,status='unknown',iostat=ioerror,
     &          err=1000)
        rewind(temp2)
222     write(stderr,
     & "(/,' enter number of events to be traced (max 200): ',
     & $)")
        read(stdin,*) npath
        if(npath .GT. 200) then
        write(stderr,"(/,' *** too many paths, please reenter:')")
        goto 222
        endif
        write(temp2,*) npath, '     # of paths'

        do250 j=1,npath

        write(stderr,"(/,' NOW READING PATH #',i3,'\n')") j
333     write(stderr,"(/,' enter number of ray segments in the ',$)")
        write(stderr,"(' path (max 100): ',$)")
        read(stdin,*) nsegment(j)
        if(nsegment(j) .GT. 100) then
        write(stderr,"(/,' *** too many segments, please reenter:')")
        goto 333
        endif 
        write(temp2,*) nsegment(j), '     # of segments in path#',j
        
*   all rays must start from the surface

        interface(j,0) = 0

*   read the path from standard input and save it in parfile

50      do20 i=1,nsegment(j)    
60           write(stderr,"(/,' enter interface index for ',$)")
             write(stderr,"(' end of segment #',i3,': ',$)") i
                read(stdin,*) interface(j,i)
                if( (abs(interface(j,i)-interface(j,i-1)).NE.
     &   1) .OR. (interface(j,i).EQ.0 .AND. i.NE.nsegment(j)) .OR. 
     &   (interface(j,i).GE.nlayer) ) then
                        write(stderr,"(/,' *** invalid ray path, please
     & reenter:')")
                        goto 60
                endif
                ilayer = max0(interface(j,i),interface(j,i-1))
70              if(triso(ilayer) .EQ. 0) then
           write(stderr,"(/,' enter propagation mode of segment #',i3 
     &,' (P or S): ',$)") i
                read(stdin,'(A)') char2(i)
                if(char2(i) .EQ. 'S' .OR. char2(i) .EQ. 's') then
                        mode(j,i) = 1
                else if(char2(i) .EQ. 'P' .OR. char2(i) .EQ. 'p') then
                        mode(j,i) = 2
                else
                        write(stderr,"(/,' *** invalid mode, please
     & reenter:')")
                        goto 70
                endif
                else
           write(stderr,"(/,' enter propagation mode of segment #',i3 
     &,' (SP, QP or QS): ',$)") i
                read(stdin,'(A)') char2(i)
                if(char2(i) .EQ. 'SP' .OR. char2(i) .EQ. 'sp') then
                        mode(j,i) = 1
                else if(char2(i) .EQ. 'QP' .OR. char2(i) .EQ. 'qp') then
                        mode(j,i) = 2
                else if(char2(i) .EQ. 'QS' .OR. char2(i) .EQ. 'qs') then
                        mode(j,i) = 3
                else
                        write(stderr,"(/,' *** invalid mode, please
     & reenter:')")
                        goto 70
                endif

                endif
20         continue
        
        if(interface(j,nsegment(j)) .NE. 0) then
                write(stderr,"(/,' *** invalid raypath (must end at 
     & surface), please reenter:')")
                goto 50
        endif

        do30 i=1,nsegment(j)

                write(temp2,"(1X,I3,3X,2A,I4)") interface(j,i),char2(i),
     &   '      destination and mode of segment#',i 

30         continue

250        continue

        close(temp2)

        return
      endif

1000   write(stderr,*) 'ERROR in opening ',parfile,' iostat=',ioerror,
     &   ' routine getpath'
      return
      end