C=================================================================
C-----     MEDGEN - MEDIUM GENERATOR
C==================================================================
      integer*2 n
      REAL cr(10)                  ! surfer water
      REAL Z1(801 ,5  ),Y(5  ,3)
      REAL Z2(801 ,5  )
!      REAL*8 RF(8),RG(8)
      real*4 ziu,zil,zfu,zfl
      REAL Z3(801 ,5  )
      INTEGER NG(5  )
!      COMPLEX*16 RT(16)
      COMPLEX W1(801 ,5  )
      COMPLEX W2(801 ,5  )
      COMPLEX W3(801 ,5  )
      COMMON/HDa/ NG,Y
      COMMON/HDb/ Z1,Z2,Z3
      COMMON/HDc/ W1,W2,W3
      real*8 coefx,coefz,coefc,coefr
      INTEGER IEOF(5), ifrwrd
      character*2 chdnm,sfrflnm
      character*4 chlxt
      character*80 flnm,mixflnm,wtrflnm,rlfflnm,polflnm
      DATA MPTR/130/
!      DATA  RF     /0.536594727634303,
!     &              0.673553876466301,
!     &              0.806390346541482,
!     &              0.946626691482681,
!     &              1.10456582640887,
!     &              1.29085359436889,
!     &              1.52319634917456,
!     &              2.11559898850822/
!      DATA   RG    /0.125863025116237,
!     &              0.218733072447135,
!     &              0.297038764222607,
!     &              0.346159504241035,
!     &              0.360780665842672,
!     &              0.330883622941456,
!     &              0.235373306621273,
!     &              0.949957005556266/
      DATA NRT/8/,CY/8./,RIPL,RIPR/0.64,1.2/,LMSTEP/4/
      LAR=801
      MGL=50
      NTP=5
      iwss=9         ! water sound speed
      n=nargs()
      if(n.le.1) then
        write(*,*)'Usage: FRMEDGEN headdatafilename PROTOCOLFLNM'
        stop
      end if
      call getarg(1,flnm,n)
      open(unit=10,file=flnm,mode='read')
      OPEN(UNIT=12,FILE='A12.wpr',FORM='UNFORMATTED')
      OPEN(UNIT=13,FILE='A13.wpr',FORM='UNFORMATTED')
      OPEN(UNIT=14,FILE='A14.wpr',FORM='UNFORMATTED')
      OPEN(UNIT=15,FILE='A15.wpr',FORM='UNFORMATTED')
      OPEN(UNIT=16,FILE='A16.wpr',FORM='UNFORMATTED')
!      OPEN(UNIT=20,FILE='A20.wpr',FORM='UNFORMATTED')
!      OPEN(UNIT=21,FILE='A21.wpr',FORM='UNFORMATTED')
      OPEN(UNIT=20,FILE='A20.wpr',FORM='binary')
      OPEN(UNIT=21,FILE='A21.wpr',FORM='binary')
      IDBPR=0
      n=nargs()
      if (n.ge.3) then
        call getarg(2,flnm,n)
        iprtcl=43
          open(unit=iprtcl,file=flnm,mode='write')
        idbpr=1
        write(iprtcl,*)' ----- Reading phase -----'
      end if
c-----------------------
      I=0
C-----                       FREQUENCY
      I=I+1
      cmfr=1.0
      cmfr0=1.0
      READ (10,*,ERR=1100,END=1100) FR
C                            DEPTH GRID STEP
      I=I+1
      READ (10,*,ERR=1100,END=1100) DZ
C                            REFERENCE SOUND SPEED
      I=I+1
      READ (10,*,ERR=1100,END=1100) CRF
C                            FINAL HORIZONTAL
      I=I+1
      READ (10,*,ERR=1100,END=1100) YFN
C                            INITIAL NODE FOR PATTERN
      I=I+1
      READ (10,*,ERR=1100,END=1100) MPI
C                 AVERAGE INTERVAL IN NODES FOR PATTERN
      I=I+1
      READ (10,*,ERR=1100,END=1100) MPA
C                 COMPUTATION STEPS PER ONE PATTERN STRING
      I=I+1
      READ (10,*,ERR=1100,END=1100) ISTEP
C                      DB VALUE PER ONE POSITION IN PATTERN
      I=I+1
      READ (10,*,ERR=1100,END=1100) SCL
C                  COMP. STEPS PER APPROXIMATION of medium
      I=I+1
      READ (10,*,ERR=1100,END=1100) LMSTEP
c---------------------pole range intervals
        read(10,*)irt
        write(20) irt
        if(idbpr.ne.0)
     +   write(iprtcl,*)'Number of pole regions:',irt
        read(10,*)chdnm
        write(20)chdnm
        if(idbpr.ne.0)
     +   write(iprtcl,*)'Poles denominator char:',chdnm
        do i=1,irt
          read(10,*)xi,xf,chlxt
          write(20) xi,xf,chlxt
        if(idbpr.ne.0)
     +   write(iprtcl,*)'Poles region limits/xtnsn:',xi,xf,chlxt
        end do
      VRF=2.*3.141592653589793*FR/CRF
      CW=(DZ*VRF)**2/12.
C--------------EXTERNAL MEMORY FIELD WRITING REGIONS
      READ(10,*,ERR=1200,END=1200)NXRG
      WRITE(20) NXRG
        if(idbpr.ne.0)
     +   write(iprtcl,*)'Number of results regions:',nxrg
      IF(NXRG.gt.0) then
      DO NN=1,NXRG
      READ(10,*,ERR=1200,END=1200) XI,XF,ZiU,ZiL,zfu,zfl
      WRITE(20,ERR=1200) XI,XF,ZiU,ZiL,zfu,zfl
        if(idbpr.ne.0) then
          write(iprtcl,*)'       Left/right abscissa limits:',xi,xf
          write(iprtcl,*) 'Left upper/lower ordinata limits:',ziu,zil
          write(iprtcl,*)'Right upper/lower ordinata limits:',zfu,zfl
        end if
      end do
      end if
C-----HYDROGRAPHY
      read(10,*)mixflnm
      read(10,*)wtrflnm,coefz,coefx,coefc
      read(10,*)rlfflnm,coefr
      read(10,*)polflnm
      read(10,*)sfrflnm
      open(unit=11,file=mixflnm,mode='read')
      cmn=100.0
      IU=11
      READ (11,*,END=1090)NG
c----------------------------surfer water file
      if (len_trim(wtrflnm).gt.0) then
      open(unit=iwss,file=wtrflnm,iostat=nios
     +    ,mode='read',status='old') ! water sound speed file
      if (nios.eq.0) then
        read(iwss,'(a6)')flnm
        read(iwss,*)nzw,nyw
        read(iwss,*)zmnw,zmxw
        read(iwss,*)ymnw,ymxw
        read(iwss,*)cmn,cmx
        if(idbpr.ne.0) then
          write(iprtcl,*)'Z and X dimensions:',nzw,nyw
          write(iprtcl,*)'Lo/Hi on Z:',zmnw,zmxw
          write(iprtcl,*)'Lo/Hi on Y:',ymnw,ymxw
          write(iprtcl,*)'Lo/Hi on C:',cmn,cmx
        end if
        ng(1)=nzw+2
        IT=1
        zmnw=coefz*zmnw
        zmxw=coefz*zmxw
        ymnw=coefx*ymnw
        ymxw=coefx*ymxw
        cmn=coefc*cmn
        cmx=coefc*cmx
        dyw=(ymxw-ymnw)/(nyw-1)
        do iy=1,nyw
          yw=ymnw+(iy-1)*dyw
          if(yw.eq.ymnw) then
            dzw=(zmxw-zmnw)/(nzw-1)
            z1(1,1)=zmnw   !  2 addititional hrzts
            z1(2,1)=zmnw
            do iz=1,nzw
              z1(iz+2,1)=zmnw+(iz-1)*dzw
            end do
          end if
          izo=3
          do iz=1,nzw,10
            j=nzw-iz+1
            if(j.gt.10)j=10
            read(iwss,*)(cr(i),i=1,j)
            if(iz.eq.1) then
              W1(1,1)=(CRF/CMPLX(coefc*cr(1),0.0))**2
              W1(2,1)=(CRF/CMPLX(coefc*cr(1),0.0))**2
            end if
            do i=1,j
              W1(izo,1)=(CRF/CMPLX(coefc*cr(I),0.0))**2
              if(cmn.gt.cr(i))cmn=cr(i)
              izo=izo+1
            end do
          end do
          write(*,*)' surf yw,j,izo:',yw,j,izo-1
          read(iwss,'(a6)')skip  ! blank line
          L=NG(IT)
          IU=IT+11
          WRITE (IU)yw,(Z1(I,IT),W1(I,IT),I=1,L)
          ML=L*2+1
        end do
        close(unit=iwss)
      end if
      end if
c----------------------------surfer relief file
      if (len_trim(rlfflnm).gt.0) then
      open(unit=iwss,file=rlfflnm,iostat=nios
     +    ,mode='read',status='old') ! relief file
      if (nios.eq.0) then
        IT=2
        ng(2)=2
        read(iwss,*)j,i,dzw
        dzw=coefr*dzw
        write(*,*)' medgen - rlf trans dpth:',dzw
        do k=1,j
          read(iwss,*)z,yw
          z=coefr*z
          write(*,*)' relief',yw,z
          z1(1,2)=z
          z1(2,2)=z+dzw
          if(k.eq.1) then
            W1(1,2)=(CRF/CMPLX(1.5,0.0))**2
            W1(2,2)=(CRF/CMPLX(1.5,0.0))**2
          end if
          L=NG(IT)
          IU=IT+11
          WRITE (IU)yw,(Z1(I,IT),W1(I,IT),I=1,L)
          WRITE (*,*)'medgen rlf:',yw,(Z1(I,IT),W1(I,IT),I=1,L)
          ML=L*2+1
        end do
        close(unit=iwss)
      end if   !--------relief surfer file ended
      end if
   31 CONTINUE
C-----
      READ(11,*,END=51)IT
      IF(IDBPR.GT.0)PRINT *,' it=',IT
      IF(IT.LE.0)GOTO 51
      READ(11,*,END=51)YS
      IF(IDBPR.GT.0)PRINT *,' y=',YS
      READ(11,*,ERR=1070,END=1070)NI,NF
      IF(IDBPR.GT.0)PRINT *,' ni,nf:',NI,NF
      IF(NI.GT.NF)GOTO 33
      DO 32 I=NI,NF
         K=I
         READ (11,*,END=1080)F
         IF(IDBPR.GT.0)PRINT *,' depth=',F
         Z1(K,IT)=F
   32 CONTINUE
   33 CONTINUE
      READ(11,*,ERR=1070,END=1070)NI,NF
      IF(IDBPR.GT.0)PRINT *,' ni,nf:',NI,NF
      IF(NI.GT.NF)GOTO 43
      DO I=NI,NF
         K=I
         READ(11,*,END=1060)G
         write(*,*)' real=',g
         IF(IT.le.3) then
           READ(11,*,END=1060)H
           write(*,*)' imag=',h
           W1(K,IT)=(CRF/CMPLX(G,H))**2
           if (cmn.gt.g)cmn=g
         else
           W1(K,IT)=G
         end if
      end do ! i
   43 CONTINUE
      L=NG(IT)
      IU=IT+11
      WRITE (IU)YS,(Z1(I,IT),W1(I,IT),I=1,L)
      ML=L*2+1
      GOTO 31
51    CONTINUE
C-----CLOSE AND REWIND
      write(*,*)' input files finished'
       END FILE 12
       REWIND 12
       END FILE 13
       REWIND 13
       END FILE 14
       REWIND 14
       END FILE 15
       REWIND 15
       END FILE 16
       REWIND 16
      AMP=1./MPA
      KPTR=MPTR*MPA+MPI-1
      NDPW=NG(1)+NG(2)+NG(3)
      NSRW=NG(1)+NG(2)
      NPRW=NG(1)+1
      NRDP=NG(4)+NG(5)
      NRPR=1
      NRSR=NG(4)
C-----------------------------SOURCE DEFINITIONS
      NFL1=0
      NFL2=0
      NFL3=0
      ifrwrd=0
231   CONTINUE
C  -
      READ(10,*,END=295)IT
C     PRINT *,IT
      IF(IT.LE.0)GOTO 295
      IF(IT.EQ.2)NFL2=NFL2+1
      IF(IT.EQ.3.AND.NFL3.EQ.0)NFL3=1
C  -
      READ(10,*,END=1307)YS
      kkk=ifrwrd
      iii=ifrwrd
      IF(ifrwrd.EQ.0)YSR=YS
      write(*,*)'ys,yfn,ifrwrd,iii',YS,yfn,ifrwrd,iii
      IF(ifrwrd.EQ.0) then
        if(yfn.ge.ys) then
          ifrwrd=1
         else
          ifrwrd=-1
        end if
      end if
      iii=ifrwrd
      write(*,*)'ysr,yfn,ifrwrd,iii',YSr,yfn,ifrwrd,iii
      !pause
      iii=ifrwrd    ! misterious bug of assigning value in comparison expression
      IF(iii.gt.0.AND.YS.GT.YSR) then
        YSR=YS
      end if
      write(*,*)'ys,yfn,ifrwrd,kkk',YS,yfn,ifrwrd,kkk
      !pause
      iii=ifrwrd
      IF(iii.lt.0.AND.YS.LT.YSR)YSR=YS
      write(*,*)'ys,yfn',YS,yfn
      write(*,*)' ifrwrd,iii',ifrwrd,iii
      !pause
C  -
      READ (10,*,END=1307) cmx
C  -
      READ (10,*,END=1307) NSRAR
      IF(NSRAR.EQ.0)NSRAR=1
      IF(NSRAR.LT.0)GOTO 271
C  -
      READ(10,*,END=1307)NI,NF
      IF(NI.GT.NF)GOTO 247
      DO 239 I=NI,NF
       READ (10,*,END=1307)F
       Z1(I,1)=F
C      Z1(I,1)=VRF*F
  239 CONTINUE
  247 CONTINUE
      READ (10,*,END=1307)NI,NF
      IF(NI.GT.NF)GOTO 263
      DO 255 I=NI,NF
       READ(10,*,END=1307)G,H
       W1(I,1)=CMPLX(G,H)
  255 CONTINUE
  263 CONTINUE
      GOTO 287
C  -  EQUIDISTANT SOURCES ARRAY
  271 CONTINUE
      NSRAR=-NSRAR
      READ (10,*,END=1307) G
C      Z1(1,1)=VRF*G
      Z1(1,1)=G
      READ (10,*,END=1307) G
C      G=VRF*G
      READ (10,*,END=1307) E,F
      W1(1,1)=CMPLX(E,F)
      READ (10,*,END=1307) E,F
      W1(1,2)=CMPLX(E,F)
      IF(NSRAR.LT.2)GOTO 287
      DO 279 I=2,NSRAR
C      Z1(I,1)=Z1(I-1,1)+VRF*G
      Z1(I,1)=Z1(I-1,1)+G
      W1(I,1)=W1(I-1,1)*W1(1,2)
  279 CONTINUE
  287 CONTINUE
      WRITE(21)IT
      L=NSRAR
      dir=0.0 !---------------obsolete
      WRITE (21)YS,cmn,cmx,crf,DIR,NSRAR,(Z1(I,1),W1(I,1),I=1,L)
      WRITE (*,*)' medgen-source definition results, ysrc',YS
      WRITE (*,*)' medgen-cmn,cmx,crf',cmn,cmx,crf
      WRITE (*,*)' medgen-dir,nsrar',DIR,NSRAR
      WRITE (*,*)' medgen-depth,amp',(Z1(I,1),W1(I,1),I=1,L)
      GOTO 231
295   CONTINUE
      L=-1
      WRITE(21)L
      END FILE 21
      REWIND 21
C  -
      DY=SIGN(CY*CRF/FR,YFN-YSR)
C      DZ=VRF*DZ
      K=2*NRT+33
      F=0.
c----------------------------------------------
      write(*,*)' building resulting file,ifrwrd',ifrwrd
      !pause
      IF(ifrwrd.EQ.0) then
        if(yfn.ge.ysr) then
          ifrwrd=1
         else
          ifrwrd=-1
        end if
      end if
      WRITE(20)DZ,DY,LMSTEP,YSR,YFN,ifrwrd,CRF,VRF
     +       ,NG,NDPW,NPRW,NSRW,NRDP,NRPR,NRSR
     +       ,KPTR,MPI,MPA,AMP,ISTEP,SCL
     +       ,NRT,RIPL,RIPR,fr,polflnm,sfrflnm
      if (idbpr.ne.0) then
        write(iprtcl,*)
     +' ----- Read and written as in head data file: -----'
        write(iprtcl,*)'                      Frequency:',fr
        write(iprtcl,*)'                     Depth step:',dz
        write(iprtcl,*)'                Reference speed:',crf
        write(iprtcl,*)'                 Final abscissa:',yfn
        write(iprtcl,*)'               Pattern ini node:',mpi
        write(iprtcl,*)'     Pattern averaging interval:',mpa
        write(iprtcl,*)'      Comp.steps/pattern string:',istep
        write(iprtcl,*)'              Pattern DB rating:',scl
        write(iprtcl,*)'Comp.steps/medium approximation:',lmstep

      end if
      PRINT *,DZ,DY,LMSTEP,YSR,YFN,ifrwrd,CRF,VRF
     +       ,NG,NDPW,NPRW,NSRW,NRDP,NRPR,NRSR
     +       ,KPTR,MPI,MPA,AMP,ISTEP,SCL
     +       ,NRT,RIPL,RIPR,sfrflnm
      PRINT *,'DZ,DY,LMSTEP,YSR,YFN,ifrwrd,CRF,VRF'
     +       ,'NG,NDPW,NPRW,NSRW,NRDP,NRPR,NRSR'
     +       ,'KPTR,MPI,MPA,AMP,ISTEP,SCL'
     +       ,'NRT,RIPL,RIPR'
      !PRINT *,(RT(I),I=1,NRT)
C-----REARRANGE AND INTERPOLATE
      DO 73 J=1,NTP
         IEOF(J)=0
         DO 73 I=1,2
   73    CALL GTSL(J,IEOF(J))
      IF(Y(1,1).NE.Y(2,1).OR.
     &   Y(2,1).NE.Y(3,1).OR.
     &   Y(3,1).NE.Y(1,1))STOP 07
      CALL PTSL(Y(1,1))
   77 YS=AMIN1(Y(1,2),Y(2,2),Y(3,2),Y(4,2),Y(5,2))
      CALL PTSL(YS)
      DO 81 I=1,NTP
         IF(YS.EQ.Y(I,2))CALL GTSL(I,IEOF(I))
   81 CONTINUE
      IF(IEOF(1)*IEOF(2)*IEOF(3)*IEOF(4)*IEOF(5).NE.1)GOTO 77
      IF(Y(1,2).NE.Y(2,2).OR.
     &Y(2,2).NE.Y(3,2).OR.
     &   Y(3,2).NE.Y(1,2))STOP 06
      END FILE 20
      REWIND 20
      PRINT 2011
      PRINT 1011
      PRINT 2011
1011  FORMAT(24X,' *  MEDGEN SUCCESSFULLY FINISHED *')
2011  FORMAT(24X,' *********************************')
      STOP
C              ERRORS DURING READINGS
1060  PRINT 1061
1061  FORMAT(64X,' ERROR READING  G,H')
      STOP 16
1070  PRINT 1071
1071  FORMAT(64X,' ERROR READING  NI,NF')
      STOP 16
1080  PRINT *,F
      STOP 16
1090  PRINT 1091
1091  FORMAT(64X,' ERROR READING   NG(1:5)')
      STOP 16
1100  PRINT 1101,I
1101  FORMAT(64X,'ERROR READING FR,DZ,E.T. AT ',I2)
      PRINT *,FR,DZ,CRF,YFN,MPI,MPA,ISTEP,SCL,LMSTEP,IRT
      STOP 16
1200  PRINT 1201
1201  FORMAT(64X,'ERROR READING    NXRG,XI,XF,ZU,ZL')
      PRINT *,NXRG,XI,XF,ZiU,ZiL
      STOP 16
1307  PRINT 1310
1310  FORMAT(64X,'EOF OR ERROR AT SORCE DEFINITIONS READING')
      STOP 16
      END
C $INCLUDE:'PTSL.FOR'
C $INCLUDE:'GTSL.FOR'