      program inv_grav2d
C 
C     THE PROGRAM INVERTS GRAVITYC DATA TO PRODUCE A STEPY MODEL 
C     AS THE CAUSE OF THE FIELD OBSERVATION
C
      integer idebug
      common /debug/ idebug
c      
      external forwdr
      integer*4  status
      character*1 coord,res
      character*10 p,TITLE(8)
      character*15 plot_text
      character*80 ctext
      CHARACTER*80 FILENAME
      real Y(500),lbls(30)
      DIMENSION WSP(5000),VAR(20),VAR0(20)
      COMMON/PROF/x(500),xSTART,DELX,AZIM
      COMMON/MODEL/DELRO,Xb(20),zb(20)
      COMMON/DATA/DATA(500) 
      COMMON/VAR/NVAR
      real DATA2(500) 
      LOGICAL UNOBSCURED ,auto_refresh
      COMMON/PMODEL/MODELX,MODELY,MPTS,NBEDS
      INTEGER NBEDS,MPTS(10)
      REAL MODELX(10,10),MODELY(10,10)
      COMMON/FIT/KFIT
      COMMON/DELT/DUM(20),IPRNTT 
      COMMON/REF/XREF,TREF,AMPD,IGEOM
c
      integer nchars, jd
      character filen*80, file_log*80, crep*1
      logical lexist
c
      DATA RAD/0.017453/

C 
c
c
      write(*,1234) 
 1234 format(////,'PROGRAM INV_GRAV2D - INVERTS GRAVITY 2d DATA',
     +    /,'Version 21st February 1994 ')
c
      write(*,*) ' '
      write(*,*) ' '
      WRITE(*,*) 'Inversion of data for a gravity 2D model'
      write(*,*) '----------------------------------------'
      write(*,*) ' '
c
c
      write(*,1235)
 1235 format(//,'Do you want a LOG file (debug mode) [n] ?? >',$)
      crep = 'n'
      call creadd (crep, 2, 'yn')
c
      idebug = 0
      if (crep .eq. 'y')  idebug = 1  
c
c
c
c
c     open file for log of program performance
c     note if more than one user using this program it will create
c     seperate files
c     Construct the log name
c                       
      filen='inv_grav2d'
      Call clengt (filen, nchars)
      if (idebug .eq. 0)  go to 1238
c 
      Do 1237 jd=1,99 
        Write(file_log,1236) filen(1:nchars), jd
 1236   Format(a, '.log', i2.2 )
        call fexist (file_log, lexist)
	if (lexist) go to 1237
c
        OPEN (13, file=file_log, status='NEW', iostat=status)
        if (status.eq.0) go to 1238
        if(jd.eq.99) then
           print *, 'Cannot create (>99) .log file'
           go to 9999
        endif
 1237 continue
c
c
c
 1238 write(*,*) ' '
      write(*,*) ' '
      if (idebug .eq. 1) write(*,1239) file_log
 1239 format('Your .log file name is: ',a)
c
c

      WRITE(*,101)
 101  format(//,'Enter the name of your CONTROL file > ',$)
      READ(*,'(A)') FILENAME
      OPEN(10,FILE=FILENAME,status='old')
      write(*,105)
 105  format('Enter the name of a file for output >',$)
      READ(*,'(A)') FILENAME
      OPEN(11,FILE=FILENAME)
C-=--------------------------------------------------
      KFIDEL=0 
      IPRNTT=0
      READ(10, 1) (TITLE(I),I=1,8) 
 199  WRITE(11,15) (TITLE(I),I=1,8) 
      IF(KFIDEL.EQ.1) GO TO 200 
 15   FORMAT(1H1,10X,8A10)
 1    FORMAT(8A10)
      READ(10,*) NPTS,KFIT,ITERS,NVAR
 2    FORMAT(8I5) 
      READ(10,*) XREF,TREF,AMPD,azim 
      write(*,1001) xref,tref,ampd,azim
 1001 format(4f10.1)
 3    FORMAT(8F10.1)
c     TRY TO READ NEXT LINE AS A FILEN NAME (X,Y, Z FILE), ELSE READ IN
C     DATA VALUES FROM INPUT FILE.
      READ(10,'(A)',err=5) FILENAME
      open(23,file=filename,status='OLD',err=2001)
      write(*,'(/,a29,a)') 'Reading observed data from input file ',
     1  filename
      n=1
 4    read(23,*,end=45) x(n),y(n),data(n)
      n=n+1
      goto 4
 45   close (23)
      npts=n-1
      write(*,'(/,a10,i3,a36 )') 'There are ',npts,
     1 ' stations in the observed data line.'
      write(*,47)
 47   format('Which coord represents station (x or y) >' ,$)
      read(*,'(a)') coord  
      if ((coord.eq.'y').or.(coord.eq.'Y') ) then
        do 6 i=1,npts
 6        x(i)=y(i)
        endif
      GO TO 7
 5    BACKSPACE 10
      READ(10,*) (DATA(I),I=1,NPTS)
      read(10,*) (x(i),i=1,npts)
C
 7    continue 
 200  READ(10,*) (VAR(I),I=1,NVAR)
      do 210 i=1,nvar
 210     var0(i)=var(i)
      READ(10,*) (DUM(I),I=1,NVAR)
      P='   ABSOLUTE'
      IF(KFIT.EQ.1) P='  RELATIVE'
      WRITE(11,801) P 
 801  FORMAT(/,5X,A10,/' ITNO        RMS ERR         DELRO          XOR' 
     1'           DEPTH      DEXTENT          DIP') 
 10   FORMAT(8F10.4)
C 
c     all data now read in
      write(*,11) 
 11   format(/,'All Data Now  Read In - Calculating Step Model',/)
C     SET THE PARAMETRS UP FOR NLSQ2 
      if (idebug .eq. 1)  write(13,*) 'calling PARDIR aaa..',nvar
      CALL PARDIR(VAR,NVAR) 
      if (idebug .eq. 1)  write(13,*) 'after calling PARDIR aaa..',nvar
C 
C     SET UP THE ROUTINE FOR DIFERS
C     CALL FORSET(Z,DUM,Z,MAX,Z,NVAR,Z,0)
C 
C     SET UP THE DIF ROUTINE

      WSP(1)=XREF 
      WSP(2)=TREF 
      WSP(3)=AMPD 
      IGEOM = 1	                                   ! akak
c
C 
      if (idebug .eq. 1)  write(13,*) 'calling DIFSET..'
      CALL DIFSET(WSP,KFIT,1) 
      if (idebug .eq. 1)  write(13,*) 'after calling DIFSET..'
C 
      NIT=ITERS 
      IPRNT=0 
      DU=DUM(1)
      DUM(1)=0.0
C 
C     CALL THE MAIN INVERSION ROUTINE AND BEGIN THE GEOMETRIC INVERSION
      CALL NLSQ2(WSP,VAR,NPTS,NVAR,IPRNT,FORWDR,NIT)
C 
C     SET UP FOR THE DELRO INVERSION 
      CALL FORSET(Z,DUM,Z,MAX,Z,NVAR,Z,0)
C 
      WSP(1)=XREF 
      WSP(2)=TREF 
      WSP(3)=AMPD
      IGEOM = 0                                    ! akak
C 
      CALL DIFSET(WSP,KFIT,0) 
C 
C     CALL THE MAIN ROUTINE TO START THE INVERSION PROCESS FOR DELTA RO 
      NIT=ITERS  
      DUM(1)=DU
      CALL NLSQ2(WSP,VAR,NPTS,NVAR,IPRNT,FORWDR,NIT)
C 
C     CALL FORWD TO OBTAIN THE DIFFERENCE VECTOR 
      MAX=MAX0(NPTS,NVAR) 
      CALL FORWDR(WSP,VAR,WSP,MAX,NPTS,NVAR,RMS,2)
C 
C     OUTPUT TABLE TO START 
      RMS=SQRT(RMS/FLOAT(NPTS)) 
C 
C     RECOVER FINAL PARAMETERS 
      CALL PARINV(VAR,NVAR) 
C 
      WRITE(11,802)  
 802  FORMAT(//,1X,' FINAL MODEL:',5X,' DRO',10x,' XB(I)',10X,'ZB(
     1I)',//)
      WRITE(11,888) VAR(1),VAR(2),VAR(3)
 888  FORMAT(19X,F4.1,3X,F16.1,3X,F15.1)
      ns=(NVAR-1)/2
      DO 889 I=2,NS
 889  WRITE(11,890) VAR(2*I),VAR(2*I+1) 
 890  FORMAT(26X,F16.1,3X,F15.1)
 
      WRITE(11, 803) 
 803  FORMAT(/,9X'X',5X,'OBSERVED',5X,'CALCULATED',5X'ERROR'/)
C 
C     RECOVER THE CALCULATED
      DO 804 I=1,NPTS 
 804    WSP(NPTS+I)=DATA(I)-WSP(I)
C     
c      x=xstart-delx  

      DO 805 I=1,NPTS 
c      x=x+delx
 805  WRITE(11,806) X(i),DATA(I),WSP(NPTS+I),WSP(I)
C 
 806  FORMAT(1X,F9.1,F13.4,F15.4,F10.4) 
C 
C 
C     GENERATE THE PROFILE ARRAY 
      WSP(2*NPTS+1)=XSTART
      DO 809 I=2,NPTS 
 809    WSP(2*NPTS+I)=WSP(2*NPTS+I-1)+DELX
C 
C     TEST FOR DATA FINISH - OR ANOTHER INVERSION
      READ(10,1) (TITLE(I),I=1,8) 
      KFIDEL=0
      IF(TITLE(1).EQ.'FINISH    '   ) GO TO 3000 
      IF(TITLE(1).EQ.'REPEAT    ') KFIDEL=1
      GO TO 199 

C     PROGRAM ERROR MESSAGES
 2001 print *,'error opening file as old file',filename
 3000 continue
c
 9999 write(*,*) ' '
      stop 'End of program INV_GRAV2D..'
C
      END 


      BLOCK DATA
C     BLOCK DATA SETS THE LOCAL CONSTANTS 
C     ETS SHOULD BE CHANGED FOR EACH MACHINA
C     ETA IS MACHINE PRECISION ZERO 
C     FOR CYBER ETA=1.E-14
C     TOL IS THE SMALLEST POSITIVE NUMBER DIVIDED BY ETA
C     FOR CYBER TOL=1.0E-278
C     EPS IS THE LARGEST VALUE EXP ROUTINE WILL TAKE
C     FOR CYBER EPS=146.
C 
      COMMON/DECEXP/DECR,EXPND,TOL1,TOL2,EP1,BND,MODE,IEPS,ICORL
      COMMON/MACHEP/ETA,TOL 
      COMMON/LOCAL/EPS
      DATA DECR,EXPND,TOL1,TOL2,EP1,BND,MODE,IEPS,ICORL/
     1 .5,2.,.01,.02,.01,.05,1,0,0/ 
      DATA ETA/1.0E-14/,TOL/1.0E-38/ 
      DATA EPS/38./
      END 


      SUBROUTINE PARDIR(PAR,NPAR) 
C 
C     DIRECT TRANSFORMATION OF INVERSION PARAMETERS
C 
C
      integer idebug
      common /debug/ idebug
      DIMENSION PAR(NPAR)
c  
      if (idebug .eq. 1)  write(13,*) 'in PARDIR..',npar
      DO 10 N=2,NPAR
 10   PAR(N)=ALOG(PAR(N)) 
      if (idebug .eq. 1)  write(13,*) 'leaving PARDIR..',npar
      RETURN
      end
C 
      SUBROUTINE PARINV(PAR,NPAR) 
      integer idebug
      common /debug/ idebug
      DIMENSION PAR(NPAR)
c  
      if (idebug .eq. 1)  write(13,*) 'in PARINV..',npar
      write(*,*) npar
      write(*,*) par(2)
      DO 20 N=2,NPAR 
 20   PAR(N)=EXP(PAR(N))
      if (idebug .eq. 1)  write(13,*) 'leaving PARINV..'
      RETURN
      END 
c
      SUBROUTINE DIFERS(VAR,DIF,ND) 
C 
C     CALCULATES THE ERROR VECTOR IN DIF(I),I=1,ND FOR
C     CURRENT VARIABLES IN VAR AND ND OBSERVED VALUES 
C     DATA VALUES IN DAT(1)-(ND) IN /DATA/
C 
C     ENTRY POINT DIFSET IS USED TO INITIALISE THE LOCATION 
C 
C     OF THE PIVOT POINT (DATUM),THE VALUE OF THE OBSERVED
C     DATA AT THE PIVOT POINT,THE ESTIMATED ANOMALY AMPLITUDE 
C     THE TYPE OF RESIDUAL REQUIRED AND THE "GEOM " FACTOR
C     WHICH EFFECTS THE TYPE OF FIT (AMPLITUDE OR GEOMETRIC)
C     THET IS REQUIRED
C
      integer idebug
      common /debug/ idebug
c  
C
      DIMENSION VAR(20),DIF(ND) 
      DIMENSION CAC(500)
      COMMON/DATA/DATA(500) 
      COMMON/PROF/x(500),XSTART,DELX,AZIM
      COMMON/MODEL/DELRO,XB(20),ZB(20) 
      COMMON/FIT/IREL 
c
      COMMON/REF/XR,TR,AMP,IGEOM
      COMMON/VAR/NVAR
      DATA RAD/0.017453/
C 
      DATA DMIN/1.0/
C 
C     SET THE MODEL PARAMETERS
C 
      if (idebug .eq. 1)  write(13,*) 'in DIFERS..'
      DELRO=VAR(1) 
      NS=(NVAR-1)/2  
      DO 25 I=1,NS
      XB(I)=VAR(2*I)
 25   ZB(I)=VAR(2*I+1)
      
C 
C     CALCULATE THE ANOMALY AT PIVOT POINT
C 
      XPIV=XR
      CALL GRAV2D(XPIV,NS,AN) 
C 
      CALCR=AN 
C 
C     CHECK FOR GEOMETRIC OPTION
      AMPD=1.0
      AMPC=1. 
      IF(IGEOM.NE.1) GO TO 10 
c
C     CALCULATE THE THEORETICAL AMPLITUDE 
      DO 45 I=1,ND
      CALL GRAV2D(X(I),NS,AN) 
C     WRITE(*,127) X(I),NS,AN
 45   CAC(I)=AN 
127   FORMAT(F10.1,I5,F10.2)

      AN1=CAC(1) 

      DO 46 I=1,ND
      IF(AN1.GE.CAC(I)) GO TO 46
      AN1=CAC(I)
 46   CONTINUE   
ccc      write(*,125) amp,an1
 125  format(2f10.1)


      AMPC=1./ABS(AN1) 
      AMPD=1./AMP
C 
C     CALCULATE THE ERROR VECTOR  
 10   CONTINUE

      DO 20 N=1,ND 
       XTEM=X(n)
       CALL GRAV2D(XTEM,NS,AN)
       CA=(AN-CALCR)*AMPC
       OB=(DATA(N)-TR)*AMPD 
       DIF(N)=OB-CA 
       IF(IREL.EQ.1) DIF(N)=DIF(N)/SIGN(AMAX1(DMIN,ABS(OB)),OB)
 20    CONTINUE 
C 
      if (idebug .eq. 1)  write(13,*) 'leaving DIFERS..'
       RETURN 
       end
C
      SUBROUTINE DIFSET(VAR,DIF,ND) 
      integer idebug
      common /debug/ idebug
      COMMON/REF/XR,TR,AMP,IGEOM
c  
C
      DIMENSION VAR(20),DIF(ND) 
      if (idebug .eq. 1)  write(13,*) 'in DIFSET..'

ccc       XREF=VAR(1)
ccc       DATR=VAR(2)                 
ccc       AMPO=ABS(VAR(3)) 
       IGEOM=ND 
      if (idebug .eq. 1)  write(13,*) 'leaving DIFSET..'
       RETURN 
       END  
c
       SUBROUTINE GRAV2D(X,NS,AN)
C      CALCULATES THE GRAVITY FIELD FROM A
C      2 D  SRUCTURE  
C
      integer idebug
      common /debug/ idebug
c  
       COMMON/MODEL/DELRO,XB(20),ZB(20)
       COMMON/PROF/xp(500),XSTART,DELX,AZIM 
       DIMENSION XS(20),A(20),B(20),D(20),C(20) 
       DATA RAD/0.017453/  
       DATA GRAVC/6.668e-3/  
      if (idebug .eq. 1)  write(13,*) 'in GRAV2D..'
       DO 1 I=1,NS
 1     XS(I)=(XB(I)-X)*SIN(AZIM*RAD)
       XS(NS+1)=XS(1)
       ZB(NS+1)=ZB(1)
       DO 2 I=1,NS
       IF(ZB(I+1).EQ.ZB(I)) GO TO 3
       A(I)=(XS(I+1)-XS(I))/(ZB(I+1)-ZB(I))
       B(I)=(XS(I)*ZB(I+1)-XS(I+1)*ZB(I))/(ZB(I+1)-ZB(I))
       D(I)=B(I)/(1.0+A(I)*A(I))
       C(I)=D(I)*A(I)
       GO TO 2
 3     D(I)=0.0
       C(I)=(XS(I)*ZB(I+1)-XS(I+1)*ZB(I))/(XS(I+1)-XS(I))
 2     CONTINUE  
       AN=0.
       DO 4 I=1,NS 
       F=(XS(I+1)**2+ZB(I+1)**2)/(XS(I)**2+ZB(I)**2)
       AN=AN+0.5*D(I)*ALOG(F)+C(I)*(ATAN2(XS(I+1),ZB(I+1))-ATAN2(XS(I),
     1ZB(I))) 
 4     CONTINUE
       AN=AN*2.0*GRAVC*DELRO      
      if (idebug .eq. 1)  write(13,*) 'leaving GRAV2D..'
       RETURN
       END
c
      SUBROUTINE FORWDR(PAR,VAR,DIF,MAX,ND,NVAR,SUMSQ,IPAR) 
C 
C     ROUTINE TO CALCULATE PARTIAL DERIVATIVES ARRAY AND 
C     ERROR VECTOR FOR S/R NLSQ2 
C
C     ENTRY FORWDR : FILLS ERROR ARRAY DIF FOR ND DATA POINTS
C                   AND RETURNS SUM OF SQUARES OF DIF'S IN
C                   SUMSQ.
C                   IF IPAR.NE.2 THEN CALCULATES PARTIAL
C                   DERIVATIVE VECTORS IN PAR(I,J) WHERE
C                   I=1,2,....,ND FOR J=1,2,..,NVAR WHERE NVAR
C                   IS NUMBER OF VARIABLES IN INVERSION 
C 
C    ENTRY FORSET :  INITIALISES ROUTINE BY SETTING DERIVATIVE 
C                   INCREMENTS DEL(I),I=1,NVAR AND SETS PRINT FLAG
C 
C    NOTE : FORSET MUST BE CALLED  BEFORE USING FORWD 
C             FOR DEL(I)=0. THEN THE CORRESPONDING ENTRIES IN PAR 
C             ARE SET TO 0. I.E. THE PARAMETER DERIVATIVES ARE SET TO 
C             0. AND HENCE PARAMETER WILL CHANGE DURING INVERSION 
C 
C
      integer idebug
      common /debug/ idebug
c  
C
      DIMENSION PAR(MAX,*),VAR(NVAR),DIF(ND)
      COMMON/DELT/DEL(20),IFLAG

      if (idebug .eq. 1)  write(13,*) 'in FORWDR..'
C 
C     RECOVER PARAMETERS TO ORIGINAL FORM 
C 
      if (idebug .eq. 1)  write(13,*) 'calling PARINV',nvar
      CALL PARINV(VAR,NVAR) 
      if (idebug .eq. 1)  write(13,*) 'after PARINV'
C 
C     CALCULATE PARTIAL DERIVATIVE VECTOR 
C 
      IF(IPAR.EQ.2) GO TO 40

      DO 30 N=1,NVAR
C 
C       CHECK FOR CONSTANT PARAMETERS
        IF(DEL(N).NE.0.) GO TO 10 
        CALL SETRA(PAR(1,N),ND,0.)
        GO TO 30
C 
 10     VAR(N)=VAR(N)-DEL(N)
        CALL DIFERS(VAR,PAR(1,N),ND)
        DV=2.*DEL(N)
        VAR(N)=VAR(N)+DV
        CALL DIFERS(VAR,DIF,ND) 
        VAR(N)=VAR(N)-DEL(N)
C 
C       COMPUTE THE PARTIALS DERIVATIVES 
        DKON=1. 
        IF(N.NE.1) DKON=VAR(N)
        DO 20 M=1,ND
        PAR(M,N)=DKON*(DIF(M)-PAR(M,N))/DV 
 20     CONTINUE
C 
 30     CONTINUE
C 
C     CALCULATE THE ERROR VECTOR
 40   CALL DIFERS(VAR,DIF,ND) 
C 
C     CALCULATE THE SUMSQ
      SUMSQ=0.
      DO 42 N=1,ND
 42     SUMSQ=SUMSQ+DIF(N)*DIF(N) 
      RMS=SQRT(SUMSQ/FLOAT(ND))
C 
C     LIST PARAMETERS IF REQUIRED 
C 
      IF(IPAR.NE.2) IFLAG=IFLAG+1 
      IF(IFLAG.GT.0) WRITE(11, 800) IFLAG,RMS,(VAR(I),I=1,NVAR)
 800  FORMAT(I5,10F15.5)
C 
C     RESET PARAMS FOR CONTINUATION OF INVERSION 
      if (idebug .eq. 1)  write(13,*) 'calling PARDIR..',nvar
      CALL PARDIR(VAR,NVAR) 
C 
      if (idebug .eq. 1)  write(13,*) 'leaving FORWDR..'
      RETURN
      end
C 
      SUBROUTINE FORSET(PAR,VAR,DIF,MAX,ND,NVAR,SUMSQ,IPAR) 
C
      integer idebug
      common /debug/ idebug
c  
      DIMENSION PAR(MAX,*),VAR(NVAR),DIF(ND)
      COMMON/DELT/DEL(20),IFLAG

      if (idebug .eq. 1)  write(13,*) 'in FORSET..'
C 

C     SET IFLAG FOR PRINTING
      IFLAG=IPAR
C 
C     TRANSFER PARAM TO DERIVATIVE DELTASP 
      DO 50 N=1,NVAR
 50   DEL(N)=VAR(N) 
C 
      if (idebug .eq. 1)  write(13,*) 'leaving FORSET..'
      RETURN
      END 
c


