real head1(200,300,61),head0(200,300),haveobs(61),havegmd3(61)
real rtt(200,300,2),fsg(200,300,2),fhd(200,300,2),ffr(200,300,2),ffl(200,300,2)
real fff(200,300,2),fwl(200,300,2),fet(200,300,2),frc(200,300,2),fst(200,300,2)
integer location,iarea(200,300),wellarea(200,300)
character text*16,flname*80,flnm*80
allocatable location(:,:)

! Read area cells
iarea=0
open(16,file='..\modflow\data\newgmd3_cells.txt',status='old')
read(16,*) na
do i=1,na
  read(16,*) k1,k2
  iarea(k1,k2)=1
enddo
close(16)

open(16,file='..\modflow\data\2023h.txt')
do i=1,200
do j=1,300
  read(16,*) head0(i,j)
enddo
enddo
close(16)

open(16,file='..\input\GMD3GUIpumpin.txt',status='old')
read(16,*) nwell
wellarea=0
allocate(location(nwell,2))
do i=1,nwell
 read(16,*) location(i,1:2)
 wellarea(location(i,1),location(i,2))=1
enddo
close(16)

flname='..\modflow\1948W\1948W.HDS'
do k=1948,2008
 if(k.gt.1945)then
   write(flname(12:15),'(i4.4)') k
   flname(16:16)='W'
   write(flname(18:21),'(i4.4)') k
   flname(22:22)='W'
 endif
 open(16,file=flname,FORM='BINARY',status='old')
 read(16) NT1,NT2,RT1,NT3,TEXT,NT4,NT5,NT6
 READ(16) ((head1(i,j,k-1947), J=1,300),I=1,200)
 close(16)
 if(k.eq.1948) then
  do i=1,200
  do j=1,300
    if(head1(i,j,1).gt.100) head1(i,j,1)=head0(i,j)
  enddo
  enddo
 endif

 if(k.eq.1953) then
  open(17,file='simH2028.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 elseif(k.eq.1958) then
  open(17,file='simH2033.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 elseif(k.eq.1963) then
  open(17,file='simH2038.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 elseif(k.eq.1968) then
  open(17,file='simH2043.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 elseif(k.eq.1978) then
  open(17,file='simH2053.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 elseif(k.eq.1988) then
  open(17,file='simH2063.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 elseif(k.eq.1998) then
  open(17,file='simH2073.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 elseif(k.eq.2008) then
  open(17,file='simH2083.out')
  write(17,*) 'Row, Column, Head_ft'
  do i=1,200
  do j=1,300
   write(17,'(2i5,f12.2)') i,j,head1(i,j,k-1947)
  enddo
  enddo
  close(17)
 endif
enddo

open(17,file='hydrographobswells.out')
write(17,*)'Row Col Year (from 2023 to 2083)'
do i=1,nwell
  i1=location(i,1)
  i2=location(i,2)
  write(17,'(2i5,61(1x,f11.3))') i1,i2,head1(i1,i2,1:61)
enddo
write(17,*) 'Below are the average for all wells and the district:'
haveobs=0.0
do k=1,61
 nsumw=0
 rtw=0.
 do i=1,200
 do j=1,300
  if(head1(i,j,k).gt.100.and.iarea(i,j).eq.1) then
   nsumw=nsumw+1
   rtw=rtw+head1(i,j,k)/10000
  endif
  havegmd3(k)=rtw/nsumw*10000
 enddo
 enddo

 nsumw=0
 rtw=0.
 do i=1,nwell
  i1=location(i,1)
  i2=location(i,2)
  if(head1(i1,i2,k).gt.100) then
   nsumw=nsumw+1
   rtw=rtw+head1(i1,i2,k)/10000
  endif
  if(nsumw.gt.0) haveobs(k)=rtw/nsumw*10000
 enddo
enddo
write(17,'(2i5,61(1x,f11.3))') 0,0,haveobs(1:61)
write(17,'(2i5,61(1x,f11.3))') 0,0,havegmd3(1:61)
close(17)

!Now working on budgets

open(17,file='budgets_gmd3.out')
write(17,'(A89)') 'Year  Specified_head  Lateral_X Lateral_Y Vertical_Z Well ET Recharge Stream Storage Area'
fsg=0.
fhd=0.  
ffr=0.  
fff=0.  
ffl=0.  
fet=0.  
frc=0.  
fst=0.  
!Start of transient
flnm='..\modflow\1948W\1948W.HDS'
do k=2,122
 if(mod(k,2).eq.0) then
  write(flnm(12:15),'(i4.4)') k/2+1947
  flnm(16:16)='W'
  write(flnm(18:21),'(i4.4)') k/2+1947
  flnm(22:22)='W'
  flnm(24:26)='HDS'
 else
  flnm(16:16)='S'
  flnm(22:22)='S'
  flnm(24:26)='HDS'
 endif
 
 open(16,file=flnm,form='binary',status='old')
 read (16) it1,it2,rt1,it3,text,it4,it5,it6
 read(16) ((head0(i,j),j=1,300),i=1,200)
 close(16)

 flnm(24:26)='BUD'
  
 open(16,file=flnm,form='binary',status='old')
 do
  read(16,IOSTAT=iend) it1,it2,text,it3,it4,it5
  if(iend.lt.0) exit
  read(16) ((rtt(i,j,1),j=1,300),i=1,200)
  read(16) ((rtt(i,j,2),j=1,300),i=1,200)

  if(text.eq.'         STORAGE') then
   fsg=rtt  !         storage in cfd
  elseif(text.eq.'   CONSTANT HEAD') then
   fhd=rtt  !   constant head in cfd
  elseif(text.eq.'FLOW RIGHT FACE ')then
   ffr=rtt  ! FLOW RIGHT FACE in cfd
  elseif(text.eq.'FLOW FRONT FACE ')then
   fff=rtt  ! FLOW FRONT FACE in cfd
  elseif(text.eq.'FLOW LOWER FACE ')then
   ffl=rtt  ! FLOW LOWER FACE in cfd
  elseif(text.eq.'           WELLS')then
   fwl=rtt  !           WELLS in cfd
  elseif(text.eq.'              ET')then
   fet=rtt  !              ET in cfd
  elseif(text.eq.'        RECHARGE')then
   frc=rtt  !        RECHARGE in cfd
  elseif(text.eq.'  STREAM LEAKAGE')then
   fst=rtt  !  STREAM LEAKAGE in cfd
  endif
 enddo
 close(16)

!Compute budget for the area - hpa layer only
 cfsg=0.
 cfhd=0.
 cffr=0.
 cfff=0.
 cffl=0.
 cfwl=0.
 cfet=0.
 cfrc=0.
 cfst=0.
 areaT=0.
 do i=1,200
 do j=1,300
  if(head0(i,j).gt.0.and.iarea(i,j).eq.1) then
   areaT=areaT+160

   cfsg=cfsg+fsg(i,j,1)

   cfhd=cfhd+fhd(i,j,1)

   if(j.eq.1) then 
    rt1=0.
   else 
    rt1=ffr(i,j-1,1)
   endif 
   cffr=cffr+rt1-ffr(i,j,1)

   if(i.eq.1) then 
    rt2=0.
   else 
    rt2=fff(i-1,j,1)
   endif 
   cfff=cfff+rt2-fff(i,j,1)

! in case certain model cells have unreasonably high or low fluxed out of the low face
   if(abs(fsg(i,j,1)+fhd(i,j,1)+rt1-ffr(i,j,1)+rt2-fff(i,j,1)+ffl(i,j,1)+fwl(i,j,1)+ &
   fet(i,j,1)+frc(i,j,1)+fst(i,j,1)).lt.abs(fsg(i,j,1))) cffl=cffl-ffl(i,j,1)

   cfwl=cfwl+fwl(i,j,1)
  
   cfet=cfet+fet(i,j,1)

   cfrc=cfrc+frc(i,j,1)

   cfst=cfst+fst(i,j,1)
  endif
 enddo
 enddo

 write(17,'(f6.1,10es15.7)') k/2.+2022,cfhd,cffr,cfff,cffl,cfwl,cfet,cfrc,cfst,-cfsg,areaT
 fsg=0.
 fhd=0.  
 ffr=0.  
 fff=0.  
 ffl=0.  
 fet=0.  
 frc=0.  
 fst=0.  
enddo
close(17)

open(17,file='budgets_obswells.out')
write(17,'(A89)') 'Year  Specified_head  Lateral_X Lateral_Y Vertical_Z Well ET Recharge Stream Storage Area'
fsg=0.
fhd=0.  
ffr=0.  
fff=0.  
ffl=0.  
fet=0.  
frc=0.  
fst=0.  
!Start of transient
flnm='..\modflow\1948W\1948W.HDS'
do k=2,122
 if(mod(k,2).eq.0) then
  write(flnm(12:15),'(i4.4)') k/2+1947
  flnm(16:16)='W'
  write(flnm(18:21),'(i4.4)') k/2+1947
  flnm(22:22)='W'
  flnm(24:26)='HDS'
 else
  flnm(16:16)='S'
  flnm(22:22)='S'
  flnm(24:26)='HDS'
 endif
 
 open(16,file=flnm,form='binary',status='old')
 read (16) it1,it2,rt1,it3,text,it4,it5,it6
 read(16) ((head0(i,j),j=1,300),i=1,200)
 close(16)

 flnm(24:26)='BUD'
  
 open(16,file=flnm,form='binary',status='old')
 do
  read(16,IOSTAT=iend) it1,it2,text,it3,it4,it5
  if(iend.lt.0) exit
  read(16) ((rtt(i,j,1),j=1,300),i=1,200)
  read(16) ((rtt(i,j,2),j=1,300),i=1,200)

  if(text.eq.'         STORAGE') then
   fsg=rtt  !         storage in cfd
  elseif(text.eq.'   CONSTANT HEAD') then
   fhd=rtt  !   constant head in cfd
  elseif(text.eq.'FLOW RIGHT FACE ')then
   ffr=rtt  ! FLOW RIGHT FACE in cfd
  elseif(text.eq.'FLOW FRONT FACE ')then
   fff=rtt  ! FLOW FRONT FACE in cfd
  elseif(text.eq.'FLOW LOWER FACE ')then
   ffl=rtt  ! FLOW LOWER FACE in cfd
  elseif(text.eq.'           WELLS')then
   fwl=rtt  !           WELLS in cfd
  elseif(text.eq.'              ET')then
   fet=rtt  !              ET in cfd
  elseif(text.eq.'        RECHARGE')then
   frc=rtt  !        RECHARGE in cfd
  elseif(text.eq.'  STREAM LEAKAGE')then
   fst=rtt  !  STREAM LEAKAGE in cfd
  endif
 enddo
 close(16)

!Compute budget for the area - hpa layer only
 cfsg=0.
 cfhd=0.
 cffr=0.
 cfff=0.
 cffl=0.
 cfwl=0.
 cfet=0.
 cfrc=0.
 cfst=0.
 areaT=0.
 do i=1,200
 do j=1,300
  if(head0(i,j).gt.0.and.wellarea(i,j).eq.1) then
   areaT=areaT+160

   cfsg=cfsg+fsg(i,j,1)

   cfhd=cfhd+fhd(i,j,1)

   if(j.eq.1) then 
    rt1=0.
   else 
    rt1=ffr(i,j-1,1)
   endif 
   cffr=cffr+rt1-ffr(i,j,1)

   if(i.eq.1) then 
    rt2=0.
   else 
    rt2=fff(i-1,j,1)
   endif 
   cfff=cfff+rt2-fff(i,j,1)

! in case certain model cells have unreasonably high or low fluxed out of the low face
   if(abs(fsg(i,j,1)+fhd(i,j,1)+rt1-ffr(i,j,1)+rt2-fff(i,j,1)+ffl(i,j,1)+fwl(i,j,1)+ &
   fet(i,j,1)+frc(i,j,1)+fst(i,j,1)).lt.abs(fsg(i,j,1))) cffl=cffl-ffl(i,j,1)

   cfwl=cfwl+fwl(i,j,1)
  
   cfet=cfet+fet(i,j,1)

   cfrc=cfrc+frc(i,j,1)

   cfst=cfst+fst(i,j,1)
  endif
 enddo
 enddo

 write(17,'(f6.1,10es15.7)') k/2.+2022,cfhd,cffr,cfff,cffl,cfwl,cfet,cfrc,cfst,-cfsg,areaT
 fsg=0.
 fhd=0.  
 ffr=0.  
 fff=0.  
 ffl=0.  
 fet=0.  
 frc=0.  
 fst=0.  
enddo
close(17)

end