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