1
2      SUBROUTINE feval(DIM,W,welfile,mfnfile,hedfile,FV,FLAG)
3      IMPLICIT NONE
4C      include "mpif.h"   (commented out by GG 08/04/03)
5
6C*************************************************************************
7!     SUBROUTINE TO EVALUATE THE OBJECTIVE FUNCTION FOR OPTIMIZATION
8!
9!     THIS OBJECTIVE FUNCTION WILL:
10!       1. SET THE VALUES FOR THE SIMULATION
11!       2. CALL FLOW AND TRANSPORT SIMULATORS
12!       3. EXCTRACT NECESSARY INFORMATION FROM OUTPUT FILES
13!          BY CALLING THE SUBROUTINES WRITEWEL AND READHEAD
14!       4. COMPUTE THE COST BY CALLING THE SUBROUTINE COSTA(C,O)
15C************************************************************************
16!     Program created by Katie Kavanagh 9/11/02
17!     Last modified: 8/04/03 by Genetha Gray
18C***********************************************************************
19!     INPUT:
20!             DIM: DIMENSION OF THE PROBLEM
21!               V: POINT TO EVALUATE THE FUNCTION AT
22!         welfile: name of the .wel file (added by GG 08/04/03)
23!         mfnfile: name of the .mfn file (added by GG 08/04/03)
24!         hedfile: name of the .hed file (added by GG 08/04/03)
25!
26!     OUTPUT:
27!             FV: FUNCTION VALUE
28!           FLAG: 0 FOR SUCCESS, >0 FOR FAILURE
29C************************************************************************
30C Description of variables:
31!
32! NROW, NLAY, NCOL: number of rows, layers, columns
33! WELRATE: array of pumping rates
34! COSTC, COSTO capital and operation costs
35! TOTALCOST = COSTC + COSTO
36! HEAD array containing values of head, output from modflow
37C**********************************************************************
38      CHARACTER*80 modflw96
39C*** next line added by GG 08/04/03
40      CHARACTER*80 mfnfile,welfile,hedfile
41
42      INTEGER DIM,I,J,K,FLAG,istat, mype,npes
43      INTEGER NROW,NLAY,NCOL
44      integer modflag
45C!!!! NROW, NCOL, NLAY must match the values in otpmain.f
46      PARAMETER(NROW = 50)
47      PARAMETER(NCOL = 50)
48      PARAMETER(NLAY = 10)
49
50      REAL QTMIN, QT
51      DOUBLE PRECISION MINRATE
52      DOUBLE PRECISION V(DIM),FV,WELRATE(DIM), COSTC,COSTO,TOTALCOST
53      DOUBLE PRECISION W(100)
54      REAL*8 HEAD(NROW,NCOL,NLAY), DRWDOWN(DIM), STARTHEAD(DIM)
55      INTEGER X(DIM), Y(DIM), Z(DIM)
56      common modflag
57C***********************************************************
58C  Modification by GG 08/04/03 to get stand alone executable
59C    Comment out MPI rank calls
60C
61C      call MPI_COMM_rank(MPI_COMM_WORLD,mype,istat);
62C***********************************************************
63
64C Set iterate V  (modified by GG 08/04/03)
65      DO I = 1,DIM
66        V(I) = W(I)
67      ENDDO
68
69C Set values (MINRATE is the shut off pumping rate value)
70      MINRATE = 0.0001
71      QTMIN = -0.031998
72      QT = 0.0
73      COSTO = 0.0
74      COSTC = 0.0
75C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
76C Set your well locations
77       DO 90 , I=1,DIM
78         Z(I) = 10
79 90      CONTINUE
80
81
82
83       X(1) = NINT((1000-V(2))/20.0)
84       Y(1) = NINT(V(1)/20.0)
85
86       X(2) = NINT((1000-V(4))/20.0)
87       Y(2) = NINT(V(3)/20.0)
88
89       X(3) = NINT((1000-V(6))/20.0)
90       Y(3) = NINT(V(5)/20.0)
91
92       X(4) = NINT((1000-V(8))/20.0)
93       Y(4) = NINT(V(7)/20.0)
94
95       X(5) = NINT((1000-V(10))/20.0)
96       Y(5) = NINT(V(9)/20.0)
97
98C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
99C Set the pumping rates at the 5 wells
100       DO 666 , I=1,DIM/2
101       WELRATE(I) = -0.0064
102 666   CONTINUE
103C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
104C Set FLAG = 0 for success
105      FLAG = 0
106C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
107C Check that 2 wells do not occupy the same location
108      call location_constraint(DIM/2,X,Y,FLAG)
109
110C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
111C Check the constraints on the net pumping rate
112      call linear_constraint(DIM/2,X,Y,Z,WELRATE,NLAY,
113     & NROW, NCOL,FLAG)
114C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
115      IF(FLAG == 1)THEN
116      modflag = 1
117      ENDIF
118
119C Compute the cost:
120      IF(FLAG .EQ. 0)THEN
121
122C***********************************************************
123C  Modification by GG 08/04/03 to get stand alone executable
124C    mfn file will be generated by C script that calls stand alone
125C    executable; mype no longer used
126C
127C Set the .mfn file for modflow
128!      call writemfn(mype)
129!      if(mype<10)then
130!      write(mfnfile,'(a2,i1,a4)') 'A1',mype,'.mfn'
131!      endif
132!      if(mype>9)then
133!      write(mfnfile,'(a2,i2,a4)') 'A1',mype,'.mfn'
134!      endif
135C***********************************************************
136
137C Change the .wel file before calling MODFLOW
138C***********************************************************
139C  Modification by GG 08/04/03 to get stand alone executable
140C    name of the .wel file will be predetermined and passed into the
141C    the write wel subroutine instead of passing in mype
142C
143!      call writewel(DIM/2,X,Y,Z,WELRATE,mype)
144C***********************************************************
145      call writewel(DIM/2,X,Y,Z,WELRATE,welfile)
146      CLOSE(13)
147
148C Call MODFLOW
149      call system("echo "//mfnfile//" | modflw96")
150
151C Read in the heads from the .hed file that MODFLOW outputs
152C***********************************************************
153C  Modification by GG 08/04/03 to get stand alone executable
154C    name of .hed file passed in instead of using mype
155C
156!      call readhead(NROW,NCOL,NLAY,HEAD,mype)
157C***********************************************************
158      call readhead(NROW,NCOL,NLAY,HEAD,hedfile)
159      DO 10 , K=1,NLAY
160C      WRITE(*,*)  ((HEAD(I,J,K),J=1,NCOL),I=1,NROW)
161 10   CONTINUE
162
163
164C Compute the installation cost
165      call costac(DIM/2,X,Y,Z,WELRATE,NLAY,NROW,NCOL,COSTC,FLAG)
166C      WRITE(*,*) 'costc', COSTC
167
168C Compute the operational cost
169      call costao(DIM/2,X,Y,Z,WELRATE,HEAD,NLAY,NROW,NCOL,COSTO,FLAG)
170C      WRITE(*,*) 'costo', COSTO
171
172      FV = COSTC + COSTO
173
174C Check to see if the constraints are violated
175      call head_constraint(DIM/2,X,Y,Z,WELRATE,HEAD,NLAY,NROW,
176     & NCOL,COSTO,FLAG)
177
178C***********************************************************
179C  Modification by GG 08/04/03 to get stand alone executable
180C    write the function value and flag to a file in EvalFn.f
181C    instead of writing it to the screen here
182C
183!      WRITE(*,*) '********************************'
184!      WRITE(*,100) ' Function value =', FV
185!      WRITE(*,*) 'FLAG=', FLAG
186!      WRITE(*,*) '********************************'
187!100   FORMAT(a17,E13.5)
188C***********************************************************
189
190      ENDIF
191
192      END
193
194
195
196