MAIN
PROGRAM
! Inversion of Apparent Resistivity Data
Using Evolutionary
! Programming Technique
! By: Ayeni Gboyega O. (GLY/99/017),
! [B.Sc. Research project]
! Department of Geology,
!
!
! DATE: FEBRUARY 2004
! PROGRAMMING LANGUAGES: FORTRAN 77, FORTRAN 90 & VISUAL BASIC
! COMPILER: MICROSOFT FORTRAN POWERSTATION
! Declaration of Global Variables
integer MaxIteration,Generation
real ParamCheck,Fittest(100),Message,ModelRes(50),ncount,xval(50)
real FieldRes(50),ElecSep(50),Low_Lmt(50),Upp_Lmt(50),n
real Rraf(1200),Xxval(50),Nncount
real RMS_error,Fittestl(100)
integer No_Layer,No_Param,i,j,No_Point,seed,printgbo
character
start*1,iterate*1,CEXIT*1,Pprint,restart*1,Invert1*1
character infile*20,Invert*1,location*40,Saved*20,StnNum*10
character Instr*20,VESDate*20
! Parameters common to
common/val/bak_parent(50,100),bak_point(50,100)
common/val1/parent(50,100),chi(100),Rchi(100),Generation1
common/val3/infile
common/val4/Fittest1(100)
common/val5/iterate
common/val6/location
common/val7/Saved
common/val8/StnNum
common/val9/Rchi1(100),Rchil(100)
common/val10/Instr,VESDate
common/val11/Nnopoint,Inopoint,Nelecsep(50),Nrapf(50)
common/val12/XNrapf(50),XNelecsep(50)
common/val13/DataRes(50),DataElecSep(50),k
common/val21/RMS,raf(1200)
common/val22/Comp(50),x_axis(50)
! Title output
write(*,*)'*****************************************************'
write(*,*)'*
WELCOME TO GOA_RES 1.0
*'
write(*,*)'*****************************************************'
write(*,*)'* AN
INVERSION SOFTWARE FOR INTERPRETATION OF
*'
write(*,*)'*
SCHLUBERGER RESISTIVITY SOUNDING DATA
*'
write(*,*)'*****************************************************'
write(*,*)'*****************************************************'
write(*,*)'*
*'
write(*,*)'*
B.Sc Research Thesis *'
write(*,*)'* Ayeni Gboyega O. (GLY/99/017) *'
write(*,*)'*
write(*,*)'* 2004 *'
write(*,*)'*
*'
write(*,*)'*****************************************************'
write(*,*)'*****
SUPERVISOR: PROFESSOR M.O. OLORUNFEMI
*****'
write(*,*)'*****************************************************'
write(*,*)'*****************************************************'
write(*,*)'*
SCHLUMBERGER ARRAY
*'
write(*,*)'*****************************************************'
! For Help in using the program,
! open the Read-me file (readme.txt) in
program folder.
write(*,*)
write(*,*)'For HELP, REFER TO THE READ-ME FILE'
write(*,*)
! Set ParamCheck to minimum chi-square
value expected.
2 ParamCheck=5
! Open all input and output files.
13 call OpenFiles
! Write the array type to Data output file.
write(9,*)'ARRAY TYPE: SCHLUMBERGER'
write(9,*)
write(*,*)'______________________________________________________
write(*,*)
write(*,*)'ARRAY TYPE: SCHLUMBERGER'
write(*,*)
!*********************************************************************
! Acquire all user defined Input Data
11 call DataInput(No_Layer,No_Param,Low_Lmt,Upp_Lmt,No_Point,ElecSep, &
FieldRes,seed,MaxIteration,location)
!*********************************************************************
No_Point1=No_Point
! Initialize control for initial graph output.
999 Generation1=1.
DataElecSep = ElecSep
DataRes = FieldRes
Inopoint=No_Point
No_Point=No_Point
! Initialize all variables used in the 'forward' subroutine
raf=0.
ncount=0.
xval=0.
Rchil=0.
call CloseFiles
call OpenFiles
go to 35
!*********************************************************************
! Refresh all files
17 call CloseFiles
call OpenFiles
raf=0.
ncount=0.
xval=0.
Rchil=Rchi
Fittestl=Fittest
! Compute the lower and upper limits of new parameters for next
iteration.
do i=1,No_Param
Low_Lmt(i)=Fittest(i)-(.05*Fittest(i))
Upp_Lmt(i)=Fittest(i)+(.05*Fittest(i))
end do
!*************** GENERATION OF POPULATION ************************
! Generation of 50 models using the lower and upper limits of the
! layer parameters
35 call PopulationGen(No_Param,Low_Lmt,Upp_Lmt,seed)
!*************** 1st STAGE INVERSION
OF DATA ***********************
! Run the forward program to generate values
of Apparent Resistivity and AB/2 for each model generated
by PopulationGen() above.
call ForwardRes(No_Layer,No_Point,1)
!*************** COMPUTATION OF FITNESS ***********************
! Calculate r.m.s.
and chi-square errors of each model.
call cubspl(XNelecsep,XNrapf,k,raf,ncount,xval)
call Error(raf,k,1)
!*************** OUTPUT GRAPH OF FIELD DATA *****************
if(Generation1.eq.1.)then
Generation1=Generation1+1
! Close and then open all input and output files
call CloseFiles
call OpenFiles
! Open files needed for graph output
open(unit=8,file=infile,status='old')
! Call graph() to output graph of Field
Data
! call graph1(No_Point)
iterVB=2
call graphicplot(2.)
open(unit=25,file='iter.gbo',status='unknown')
! write(*,*)'b4',iterVB
read (25,26)iterVB
! write(*,*)iterVB
write(*,*)
write(*,*)'------------------------------------------'
write(*,*)
close (25)
if(iterVB.eq.2)go to 83
if(iterVB.ne.0)iterVB=2
if(iterVB.ne.2)iterVB=0
write(*,*)
write(*,*)'I=Inversion of Data'
write(*,*)'C=Change Model Entries'
write(*,*)
82 read(*,*)Invert1
if(Invert1.eq.'I'.or.Invert1.eq.'i')go to 83
if(Invert1.eq.'C'.or.Invert1.eq.'c') go to 11
if(Invert1.ne.'I'.and.Invert1.ne.'i')then
if(Invert1.ne.'C'.and.Invert1.ne.'c')then
write(*,*)'PLEASE CHOOSE:'
write(*,*)' I: Invert Data'
write(*,*)' C: Change
Model Parameters'
go to 82
end if
end if
end if
!*********************************************************************
! Start the iteration process from 1 to MaxIteration.
83 do
Generation=1,MaxIteration
!*********************************************************************
! Set the first model as the fittest model
do i=1,No_Param
Fittest(i)=parent(i,1)
end do
!********************* Mutation PROCESS **************************
! Call subroutine Mutation() to create 50
sets of mutated models
! from the 50 sets of models
generated by PopulationGen()
call Mutation(No_Param,seed)
!*********************************************************************
! Add the newly generated 50 models (mutants) to previous 50 models
! to create 100 models.
do j=1,50
do i=1,No_Param
parent(i,j+50)=bak_parent(i,j)
end do
end do
!*************** 2nd STAGE INVERSION
OF DATA ***********************
! Run the forward program to generate values of Apparent Resistivity and AB/2 for each new model (mutants) generated
by Mutation() above.
call ForwardRes(No_Layer,No_Point,51)
!*************** COMPUTATION OF FITNESS ***********************
! Calculate r.m.s.
and chi-square errors of each new model (mutants).
call cubspl(XNelecsep,XNrapf,k,raf,ncount,xval)
! call Error(raf,No_Point,51)
! do i=1,No_Point
! write(*,*)i,k,raf(i),DataRes(i)
! end do
call Error(raf,k,51)
!*********************************************************************
! Sort the models (new and old) according to their fitness.
call Sort_Fittness(No_Param)
!*********************************************************************
! Set the model with least r.m.s.
error as the fittest model.
do i=1,No_Param
Fittest(i)=parent(i,1)
end do
! Print iteration and least r.m.s and
chi-square errors
! write(*,*)Generation,Rchi(1),chi(1)
!*********************************************************************
! Compute parameter (Message) used to check for percentage done.
Message=float(Generation)/float(MaxIteration)
!************* OUTPUT GRAPH OF INITIAL MODEL *****************
! For the initial model (based on user input),
if (Generation1.eq.2.)then
Generation1=Generation1+1
! call the forward program to generate
values of App. Res. & AB/2,
call cubspl(XNelecsep,XNrapf,k,raf,ncount,xval)
call forward(ModelRes,No_Layer,24,Fittest1,1)
! Calculate r.m.s and chi-square errors
between Initial Model and Field Data
call Error(raf,k,1)
Rchi1=Rchi
Rchil=Rchi1
Fittestl=Fittest1
! Check for maximum allowable error between field and model
if(Rchi(1).ge.50.)then
write(*,*)'__________________________________________________'
write(*,*)
write(*,*)' R.M.S Error is large. Please Check Model Params '
write(*,*)'__________________________________________________'
write(*,*)
end if
write(*,*)
write(*,*)'Results of Modelling
with Initial Parameters'
write(*,*)'__________________________________________________'
write(*,*)
! Call result() to output 1st Generation
Model Parameters
call result(Rchi,Fittest1,No_Param,No_Layer)
! Close all files
call CloseFiles
! Open files needed for graph output
open(unit=14,file='TPrint.out',status='unknown')
open(unit=12,file='Generate.out',status='unknown')
open(unit=18,file='InterPlot.out',status='unknown')
open(unit=8,file=infile,status='old')
! Call graph() to output graph of 1st
Generation (Initial) Model
! call graph(24,No_Point1)
!*********************************************************************
! Close and then open all files
call CloseFiles
call OpenFiles
!*********************************************************************
! Continue the Inversion process or change model parameters
iterVB=2
39 call graphicplot(1.)
open(unit=25,file='iter.gbo',status='unknown')
write(*,*)'------------------------------------------'
read (25,26)iterVB
write(*,*)
close (25)
if(iterVB.eq.2)go to 80
if(iterVB.ne.2)iterVB=0
! if(iterVB.ne.0)iterVB=2
printgbo=0.
open(unit=26,file='print.gbo',status='unknown')
read(26,26)printgbo
close(26)
if(printgbo.eq.1.)then
! call graphicplot(3.)
if(Rchi(1).gt.Rchil(1))then
call Rresult1(Rchil,No_Point,Fittestl,No_Param,No_Layer)
else
call Rresult1(Rchi,No_Point,Fittest,No_Param,No_Layer)
end if
call graphicplot(3.)
go to 39
end if
! write(*,*)iterVB
write(*,*)
write(*,*)'I=Continue Inversion'
write(*,*)'C=Change Model Parameters'
write(*,*)
81 read(*,*)Invert
if(Invert.eq.'I'.or.Invert.eq.'i')go
to 80
if(Invert.eq.'C'.or.Invert.eq.'c')
go to 13
if(Invert.ne.'I'.and.Invert.ne.'i')then
if(Invert.ne.'C'.and.Invert.ne.'c')then
write(*,*)'PLEASE CHOOSE:'
write(*,*)' I: Invert Data'
write(*,*)' C: Change Model Parameters'
go to 81
end if
end if
end if
!********************************************************************
! Output Message indicating stage of
inversion
80 if (Generation.eq.10)then
write(*,*)'Please wait....'
end if
if (Message.eq.0.25)then
write(*,*)'25% Completed. Please Wait....'
end if
if (Message.eq.0.5)then
write(*,*)'50% Completed. Please Wait......'
end if
if (Message.eq.0.75)then
write(*,*)'75% Completed. Please Wait.........'
write(*,*)
end if
!*********************************************************************
! Check the termination criteria
if(Generation.eq.MaxIteration)then
!*********************************************************************
! Refresh file containing the Field Data
close(8)
open(unit=8,file=infile,status='old')
! Call the forward program to generate values of App. Res. &
AB/2
call forward(ModelRes,No_Layer,24,Fittest,1)
call cubspl(XNelecsep,XNrapf,k,raf,ncount,xval)
call Rrms(k,No_Param)
! Rchi(1)=RMS
! call ForwardRes(No_Layer,No_Point,1)
! Ensure present error is less than previous. If not, print
previous model
if(Rchi(1).gt.Rchil(1))then
call result(Rchil(1),Fittestl,No_Param,No_Layer)
else
call result(Rchi(1),Fittest,No_Param,No_Layer)
end if
! Close and open files containing results of inversion
close(12)
close(13)
open(unit=18,file='InterPlot.out',status='unknown')
open(unit=12,file='Generate.out',status='unknown')
open(unit=13,file='InterFieldData.out',status='unknown')
write(*,*)
! Check for maximum allowable error between field and model
if(chi(1).ge.5.or.Rchi(1).ge.10.)then
write(*,*)
write(*,*)'R.M.S-Error too large. Check Model Params'
write(*,*)
end if
! Check for maximum allowable error between field and model
if(Rchi(1).le.ParamCheck)then
write(*,*)'You have attained an optimum model'
end if
! Close all files
call CloseFiles
!************* OUTPUT GRAPH OF FINAL MODEL ******************
! Open files needed for final graph output
call OpenFiles
open(unit=8,file=infile,status='old')
!22 call graph(24,No_Point1)
27 call graphicplot(1.)
open(unit=25,file='iter.gbo',status='unknown')
! write(*,*)'b4',iterVB
read (25,26)iterVB
! write(*,*)iterVB
write(*,*)
close (25)
printgbo=0.
open(unit=26,file='print.gbo',status='unknown')
read(26,26)printgbo
close(26)
if(printgbo.eq.1.)then
! call graphicplot(3.)
if(Rchi(1).gt.Rchil(1))then
call Rresult(Rchil,No_Point,Fittestl,No_Param,No_Layer)
else
call Rresult(Rchi,No_Point,Fittest,No_Param,No_Layer)
end if
call graphicplot(3.)
go to 27
end if
26 format(I1)
close(26)
if(iterVB.eq.2)go to 21
if(iterate.eq.'P'.or.iterate.eq.'p')then
18 if(Rchi(1).gt.Rchil(1))then
call Rresult(Rchil,No_Point,Fittestl,No_Param,No_Layer)
else
call Rresult(Rchi,No_Point,Fittest,No_Param,No_Layer)
call graph(24,No_Point1)
end if
end if
if(iterate.eq.'S'.or.iterate.eq.'s')then
19 if(Rchi(1).gt.Rchil(1))then
call Sresult(Rchil(1),No_Point,Fittestl,No_Param,No_Layer)
else
call Sresult(Rchi(1),No_Point,Fittest,No_Param,No_Layer)
end if
call graph(24,No_Point1)
end if
! Refresh all files
21 call CloseFiles
call OpenFiles
if(Rchi(1).lt.Rchil(1))then
Rchil=Rchi
Fittestl=Fittest
end if
!********************************************************************
if(Rchi(1).gt.Rchil(1))then
write(*,*)
write(*,*)'* You have reached attained an Optimum Model *'
end if
write(*,*)'-----------------------------------------------'
write(*,*)
if(iterVB.eq.2)go to 17
write(*,*)
write(*,*)'I=ITERATE'
write(*,*)'C=CONTINUE'
write(*,*)'P=PRINT RESULTS'
write(*,*)'S=SAVE RESULTS'
write(*,*)
iterate=''
15 read(*,*)iterate
write(*,*)
open(unit=8,file=infile,status='old')
if(iterate.eq.'I'.or.iterate.eq.'i')go
to 17
if(iterate.eq.'C'.or.iterate.eq.'c') go to 16
if(iterate.eq.'P'.or.iterate.eq.'p')
go to 18
if(iterate.eq.'S'.or.iterate.eq.'s')
go to 19
if(iterate.ne.'I'.and.iterate.ne.'i')then
if(iterate.ne.'C'.and.iterate.ne.'c')then
if(iterate.ne.'P'.and.iterate.ne.'p')then
if(iterate.ne.'S'.and.iterate.ne.'s')then
write(*,*)'PLEASE CHOOSE:'
write(*,*)' I:
ITERATE'
write(*,*)' C:
CONTINUE'
write(*,*)' P:
PRINT RESULTS'
write(*,*)' S:
SAVE RESULTS'
write(*,*)
go to 15
end if
end if
end if
end if
! call Presult(chi(1),Fittest,No_Param,No_Layer)
!********************************************************************
16 go
to 99
1000 end if
end do
!*********************************************************************
Call
subroutine to close all files before termination
99 call CloseFiles
!********************************************************************
write(*,*)'______________________________________________________'
write(*,*)
write(*,*)'Do you want to Model another set of Data? '
write(*,*)' Y= YES'
write(*,*)' N= NO'
write(*,*)
12 read(*,*)restart
if(restart.eq.'Y'.or.restart.eq.'y')
go to 13
if(restart.eq.'N'.or.restart.eq.'n')
go to 14
if(start.ne.'Y'.and.start.ne.'y')then
if(start.ne.'N'.and.start.ne.'n')then
write(*,*)'PLEASE CHOOSE:'
write(*,*)' Y: MODEL ANOTHER SET OF DATA'
write(*,*)' N: CONTINUE'
go to 12
end if
end if
!********************************************************************
14 write(*,*)
write(*,*)'Do you want to exit the program? Yes=Y; No=N '
6 read(*,*)start
if(start.eq.'Y'.or.start.eq.'y')
go to 4
if(start.eq.'N'.or.start.eq.'n')
go to 5
if(start.ne.'Y'.and.start.ne.'y')then
if(start.ne.'N'.and.start.ne.'n')then
write(*,*)'PLEASE CHOOSE:'
write(*,*)' Y: EXIT PROGRAM'
write(*,*)' N: CONTINUE'
go to 6
end if
end if
write(*,*)
write(*,*)'______________________________________________________'
5 write(*,*)
write(*,*)'DO YOU
WANT TO PRINT THE RESULTS? P=PRINT, C=CONTINUE'
write(*,*)
32 read(*,*)Pprint
if(Pprint.eq.'P'.or.Pprint.eq.'p')
go to 30
if(Pprint.eq.'C'.or.Pprint.eq.'c')
go to 31
if(Pprint.ne.'P'.and.Pprint.ne.'p')then
if(Pprint.ne.'C'.and.Pprint.ne.'c')then
write(*,*)'PLEASE CHOOSE P(PRINT); OR C(CONTINUE)'
go to 32
end if
end if
30 call Presult(Rchi(1),Fittest,No_Param,No_Layer)
write(*,*)
write(*,*)'______________________________________________________'
write(*,*)
write(*,*)'TO PRINT RESULTS OPEN THE FILE (DataOutput) and Print'
write(*,*)'______________________________________________________'
4 write(*,*)
write(*,*)
write(*,*)
write(*,*)
31 write(*,*)'########### Thanks for using this program
###########'
write(*,*)
write(*,*)'----------------------------------------------------'
write(*,*)' TO MODEL
ANOTHER DATA PLEASE RESTART THE PROGRAM
'
write(*,*)'----------------------------------------------------'
write(*,*)
write(*,*)'Press C to close'
close(14)
8 read(*,*)CEXIT
if(CEXIT.eq.'C'.or.CEXIT.eq.'c')
go to 9
if(CEXIT.ne.'C'.and.CEXIT.ne.'y')then
write(*,*)'PLEASE CHOOSE C TO CLOSE THE PROGRAM'
go to 8
end if
9 stop
end