c-----twodgw: two-dimensional unsteady groundwater model.
c-----Developed by Dr. V. M. Ponce
c-----Units in meters and seconds (hours,days)
c-----Version 1.81, 08oct01
program twodgw
parameter (mz=120)
parameter (mz1= mz-1)
parameter (mz2= mz1**2)
parameter (ny=9)
character*9 date
parameter (vn=1.8,date='26Sep01')
parameter (nsm=60,nmh=60,nhd=24,xndy=365.25)
parameter (nsh=nsm*nmh,nsd=nsm*nmh*nhd,nmd=nmh*nhd)
parameter (nsy=nsd*xndy)
dimension head(0:mz,0:mz),heada(0:mz,0:mz) !head old, head new
dimension rain(mz1,mz1) ! rain
dimension xirr(mz1,mz1) ! irrigation
dimension pump(mz1,mz1) ! pumping
dimension bott(0:mz,0:mz),depth(0:mz,0:mz)
dimension mapbin(0:mz,0:mz)
dimension kleft(0:mz),kright(0:mz)
dimension kc(0:mz)
character*4 gw01,gw02,gw03,gw10,gw20,gw30,gw40,gw50,gw60,gw70,gw80
character*70 me10,me20,me30,me40,me50,me60,me70,me80
character*70 label1,label2,label3
data ts,itd /0.,0/
data rain,xirr,pump /mz2*0.,mz2*0.,mz2*0./
data sumpumpvolume,cumvolumechange /2*0./
data iflag,kflag /2*0/
c-----open statement------------------------------------------------------------------
open(07,file='twodgw.data',status='unknown')
open(08,file='twodgw.results',status='unknown')
open(09,file='twodgw.headdata',status='unknown')
open(10,file='twodgw.pumpdata',status='unknown')
open(11,file='twodgw.gridconfig',status='unknown')
open(12,file='twodgw.irrigation',status='unknown')
open(21,file='twodgw.gridecho',status='unknown')
open(23,file='twodgw.hcenter',status='unknown')
c-----reading input-------------------------------------------------------------------
c-----input of record gw01-------------------------------------------------------------
read(7,*) gw01,label1 ! label1= project
c-----input of record gw02-------------------------------------------------------------
read(7,*) gw02,label2 ! label2= date, run number
c-----input of record gw03-------------------------------------------------------------
read(7,*) gw03,label3 ! label3= run characteristics
c-----input of record gw10-------------------------------------------------------------
read(7,*) gw10,me10,tr,st,ds,dthr,nz
c-----tr= transmissivity, in m2/s
c-----st= storativity, in dimensionless units
c-----ds= grid space interval, in m
c-----dthr= time step, in hr
c-----nz= actual square grid size (0:nz)
c-----input of record gw20-------------------------------------------------------------
read(7,*) gw20,me20,indm
c-----indm= mapping indicator:
c----- indm=0 for regular grid configuration
c----- indm=1 for irregular grid configuration
c-----input of record gw30-------------------------------------------------------------
read(7,*) gw30,me30,indb
c-----indb= input indicator:
c----- indb=0 for boundary inflow (linear extrapolation from interior points)
c----- indb=1 for boundary noflow (repetition of adjacent interior point)
c-----input of record gw40-------------------------------------------------------------
read(7,*) gw40,me40,indi,href,bottref,hdref,locleft,locright
c-----indi= input indicator:
c----- indi=0 for initial conditions input from auxiliary input file
c----- indi=1 for initial conditions equal to reference head value
c----- indi=2 for initial condition including reference depletion
c-----href= reference piezometric head, in m
c-----bottref= reference rock bottom elevation, in m
c-----hdref= reference depletion head, in m
c-----locleft= location of start of reference depletion
c-----locright= location of end of reference depletion
c-----input of record gw40-------------------------------------------------------------
read(7,*) gw50,me50,indr,rainref
c-----indr= input indicator:
c----- indr=0 rain percolation feature off
c----- indr=1 rain percolation feature on
c-----rainref= reference percolation from rain, in mm/yr
c-----input of record gw50-------------------------------------------------------------
read(7,*) gw60,me60,indx,xirrref,irrleft,irrright
c-----indx= input indicator:
c----- indx=0 irrigation percolation feature off
c----- indx=1 irrigation percolation feature on
c-----xirrref= reference percolation from irrigation, in mm/yr
c-----irrleft= left location for irrigation
c-----irrright= right location for irrigation
c-----input of record gw60-------------------------------------------------------------
read(7,*) gw70,me70,indp,pumpref
c-----indp= input indicator:
c----- indp=0 standard (regular) pumping feature off
c----- indp=1 standard (regular) pumping feature on
c----- indp=2 pumping from actual pumps, read from twodgw.pumpdata
c-----pumpref= reference pumping rate (yearly equivalent), in L/s
c-----input of record gw70-------------------------------------------------------------
read(7,*) gw80,me80,tty,tpd
c-----tty= total simulation time, in yr
c-----tpd= time interval for printing report, in days
c-----definition of nz values----------------------------------------------------------
nz10=nz/10
nzh=nz/2
nzp=nz/nz10
nz1=nz-1
nzs=nz*nz
nz2=nz1**2
c-----input and echo of mapped binary codes--------------------------------------------
if(indm.eq.1) then
do j= 0,nz
read(11,'(121i1)') (mapbin(j,k),k=0,nz)
write(21,'(121i1)') (mapbin(j,k),k=0,nz) !echo to twodgw.gridecho
enddo
c-----calculation of kleft and kright limits of integration----------------------------
do j= 0,nz
do k= 0,nz1
if(k.eq.0.and.mapbin(j,k).eq.1) then
kleft(j)= 0
elseif(mapbin(j,k).eq.0.and.mapbin(j,k+1).eq.1) then
kleft(j)= k+1
endif
enddo
write(21,'(2i4)') j,kleft(j)
enddo
do j= 0,nz
do k= 0,nz1
if(mapbin(j,k).eq.1.and.mapbin(j,k+1).eq.0) then
kright(j)= k
elseif(k.eq.nz1.and.mapbin(j,k+1).eq.1) then
kright(j)= nz
endif
enddo
write(21,'(2i4)') j,kright(j)
enddo
endif
c-----default values------------------------------------------------------------------
if(tr.eq.0.) then
tr= 0.0091 ! transmissivity, in m2/sec
endif
if(st.eq.0.) then
st= 0.1 ! storativity, dimensionless
endif
if(ds.eq.0.) then
ds= 100. ! space step, in meters
endif
if(nz.gt.mz) then
iflag= 1
nz= mz
nz10=nz/10
nzh=nz/2
nzp=nz/nz10
nz1=nz-1
nzs=nz*nz
nz2=nz1**2
elseif(nz.lt.60) then
kflag= 1
endif
if(href.eq.0.) then
href= 700. ! reference head, meters
endif
if(int(rainref-0.001).eq.-1) then
rainref= 0. ! reference percolation from rain, in mm/yr
endif
if(int(xirrref-0.001).eq.-1) then
xirrref= 0. ! reference percolation from irrigation, in mm/yr
endif
if(int(pumpref-0.001).eq.-1) then
pumpref= 100. ! 100 Liters/sec
endif
if(locleft.eq.-1) then
locleft= 0.25*nz + 0.1
endif
if(locright.eq.-1) then
locright= 0.75*nz + 0.1
endif
if(irrleft.eq.-1) then
irrleft= 0.25*nz + 0.1
endif
if(irrright.eq.-1) then
irrright= 0.75*nz + 0.1
endif
if(tty.eq.0.) then
tty= 20. ! total simulation time, in year(s)
endif
if(tpd.eq.0.) then
tpd= 365. ! time to print report, in day(s)
endif
c-----printing header-----------------------------------------------------------------
100 FORMAT(1X,162(1H-))
write(8,100)
write(6,100)
write(23,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
c write(23,'((1h-))')
write(8,'(a)') 'SDSU Two-Dimensional Unsteady Groundwater Model'
write(6,'(a)') 'SDSU Two-Dimensional Unsteady Groundwater Model'
write(23,'(a)') 'SDSU Two-Dimensional Unsteady Groundwater Model'
write(8,'(a,f3.1,2a)') 'Version ',vn,', Date: ',date
write(6,'(a,f3.1,2a)') 'Version ',vn,', Date: ',date
write(23,'(a,f3.1,2a)') 'Version ',vn,', Date: ',date
write(8,100)
write(6,100)
write(23,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
c write(23,'((1h-))')
write(8,'(a)') label1
write(6,'(a)') label1
write(23,'(a)') label1
write(8,'(a)') label2
write(6,'(a)') label2
write(23,'(a)') label2
write(8,'(a)') label3
write(6,'(a)') label3
write(23,'(a)') label3
write(8,100)
write(6,100)
write(23,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
c write(23,'((1h-))')
write(23,'(a)') 'HEAD AT CENTER FIELD'
write(23,100)
c write(23,'((1h-))')
if(iflag.eq.1) then
write(6,'(a,i3)')
1'***Warning: nz value has been reset to maximum array size mz= ',
2mz
write(8,'(a,i3)')
1'***Warning: nz value has been reset to maximum array size mz= ',
2mz
endif
if(kflag.eq.1) then
write(6,'(a,i3)')
1'***Fatal error: nz cannot be less than 60. '
write(8,'(a,i3)')
1'***Fatal error: nz cannot be less than 60. '
stop 'nz too low'
endif
write(8,'(a,i3)') 'Indicator for mapped input geometry= ',indm
write(6,'(a,i3)') 'Indicator for mapped input geometry= ',indm
write(8,'(a,i3)') 'Indicator for type of boundary inflow= ',indb
write(6,'(a,i3)') 'Indicator for type of boundary inflow= ',indb
write(8,'(a,i3)') 'Indicator for initial conditions= ',indi
write(6,'(a,i3)') 'Indicator for initial conditions= ',indi
write(8,'(a,i3)') 'Indicator for rain percolation presence= ',indr
write(6,'(a,i3)') 'Indicator for rain percolation presence= ',indr
write(8,'(a,i3)') 'Indicator for irrigation perc. presence= ',indx
write(6,'(a,i3)') 'Indicator for irrigation perc. presence= ',indx
write(8,'(a,i3)') 'Indicator for pumping presence= ',indp
write(6,'(a,i3)') 'Indicator for pumping presence= ',indp
write(8,'(a,f10.4,a)') 'Transmissivity= ',tr,' m2/s'
write(6,'(a,f10.4,a)') 'Transmissivity= ',tr,' m2/s'
write(8,'(a,f10.2,a)') 'Storativity= ',st,' '
write(6,'(a,f10.2,a)') 'Storativity= ',st,' '
write(8,'(a,i10,a)') 'Grid space interval= ',int(ds),' m'
write(6,'(a,i10,a)') 'Grid space interval= ',int(ds),' m'
write(8,'(a,i10,a)') 'Actual array size= ',nz,' units'
write(6,'(a,i10,a)') 'Actual array size= ',nz,' units'
write(8,'(a,i10,a)') 'Reference head= ',int(href),' m'
write(6,'(a,i10,a)') 'Reference head= ',int(href),' m'
write(8,'(a,f10.3,a)') 'Simulation time= ',tty,' yr'
write(6,'(a,f10.3,a)') 'Simulation time= ',tty,' yr'
write(8,'(a,i10,a)')
1'Reference percolation from rainfall= ',int(rainref),' mm/yr'
write(6,'(a,i8,a)')
2'Reference percolation from rainfall= ',int(rainref),' mm/yr'
write(8,'(a,i8,a)')
1'Reference percolation from irrigation= ',int(xirrref),' mm/yr'
write(6,'(a,i8,a)')
2'Reference percolation from irrigation= ',int(xirrref),' mm/yr'
write(8,'(a,i8,a)')
1'Reference pumping rate= ',int(pumpref),' L/s'
write(6,'(a,i8,a)')
2'Reference pumping rate= ',int(pumpref),' L/s'
c-----calculation of input aquifer diffusivity----------------------------------------
diff= tr/st
write(8,'(a,f10.4,a)') 'Input diffusivity= ',diff,' m2/s'
write(6,'(a,f10.4,a)') 'Input diffusivity= ',diff,' m2/s'
c-----calculation of time interval----------------------------------------------------
dtsec= dthr*3600.
dt= dtsec
dssq= ds**2
drey= 4.*diff*dt/dssq
drey1= 1. - drey
write(8,'(a,f10.4)')
1'Calculated cell Reynolds number= ',drey
write(6,'(a,f10.4)')
1'Calculated cell Reynolds number= ',drey
tts= nsy*tty
nt= tts/dt + 0.001
np= tpd*nsd/dt + 0.001
write(8,'(a,i8,a)') 'Calculated time interval= ',int(dt),' s'
write(6,'(a,i8,a)') 'Calculated time interval= ',int(dt),' s'
dth= dt/nsh + 0.001
write(8,'(a,i8,a)') 'Calculated time interval= ',int(dth),' hr'
write(6,'(a,i8,a)') 'Calculated time interval= ',int(dth),' hr'
write(6,'(a,i9)') 'Number of time steps = ',nt
write(8,'(a,i9)') 'Number of time steps = ',nt
c-----calculation of source/sink constant---------------------------------------------
diffc= 0.25*drey*dssq/dt
trc= diffc*st
ssco= 0.25/trc
c-----conversion of sources and sinks to m3/s-----------------------------------------
rainco= ssco*0.001*dssq/nsy
xirrico= ssco*0.001*dssq/nsy
pumpco= ssco*0.001
rainref= rainco*rainref ! conversion to m3/sec, and then to L units
xirrref= xirrico*xirrref ! conversion to m3/sec, and then to L units
pumpref= pumpco*pumpref ! conversion to m3/sec, and then to L units
c-----referencing rock bottom elevation-----------------------------------------------
if(indm.eq.0) then !300
if(indi.ge.1) then !302
do j= 0,nz
do k= 0,nz
bott(j,k)= bottref
enddo
enddo
endif !302
endif !300
c-----initial conditions (first level)------------------------------------------------
if(indm.le.1) then !310
if(indi.eq.0) then !320
do j= 0,nz
read(9,*) (head(j,k),k=0,nz)
enddo
elseif(indi.eq.1) then !320
do j= 0,nz
do k= 0,nz
head(j,k)= href
enddo
enddo
elseif(indi.eq.2) then !320
do j= 0,nz
do k= 0,nz
head(j,k)= href
enddo
enddo
do j= locleft,locright
do k= locleft,locright
head(j,k)= hdref
enddo
enddo
endif !320
endif !310
c-----calculating initial volume----------------------------------------------------
if(indm.eq.0) then !350
do j= 0,nz
do k= 0,nz
depth(j,k)= head(j,k) - bott(j,k)
enddo
enddo
volinit= volall(depth,dssq,mz,nz,nz1)
volcalc0= volinit
endif !350
c-----specifying percolation from rain----------------------------------------------
if(indr.eq.1) then !360
if(indm.le.1) then !370
do j= 1,nz1
do k= 1,nz1
rain(j,k)= rainref
enddo
enddo
endif !370
endif !360
c-----specifying percolation from irrigation----------------------------------------
if(indx.eq.1) then !380
if(indm.eq.0) then !390
do k= irrleft,irrright
do j= irrleft,irrright
xirr(j,k)= xirrref
enddo
enddo
elseif(indm.eq.1) then !390
do j= 0,nz
read(12,*) (xirr(j,k),k=0,nz)
enddo
do j= 0,nz
do k= 0,nz
xirr(j,k)= irrico*xirr(j,k)
enddo
enddo
endif !390
endif !380
c-----specifying pumping locations--------------------------------------------------
if(indm.eq.0) then !420
if(indp.eq.0) then !410
continue
elseif(indp.eq.1) then !410
c-----reference pumping locations---------------------------------------------------
do j= nzp,nz-nzp,nzp
pump(j,j)= pumpref
enddo
do j= nzp,nz-nzp,nzp
k= nz - j
pump(j,k)= pumpref
enddo
c-----pumping locations read from input file----------------------------------------
elseif(indp.eq.2) then !410
do j= 0,nz
read(10,*) (pump(j,k),k=0,nz)
enddo
do j= 0,nz
do k= 0,nz
pump(j,k)= pumpco*pump(j,k)
enddo
enddo
endif !410
else !420
if(indp.eq.2) then !425
do j= 0,nz
read(10,*) (pump(j,k),k=0,nz)
enddo
do j= 0,nz
do k= 0,nz
pump(j,k)= pumpco*pump(j,k)
enddo
enddo
endif !425
endif !420
c-----calculation of pumped volume----------------------------------------------------
do j= 1,nz1
do k= 1,nz1
if(pump(j,k).gt.0.) then !520
sumpumpvolume= sumpumpvolume + dt*pump(j,k)/(ssco*1.E06)
endif !520
enddo
enddo
sumpumpvolumeaqui= sumpumpvolume/st
c-----printing of initial conditions--------------------------------------------------
if(indm.eq.0) then !430
kc(0)= 0
do j= 1,nz
kc(j)= kc(j-1) + 1
enddo
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
write(8,'(a,i5,a)') 'Time= ',itd,' days'
write(6,'(a,i5,a)') 'Time= ',itd,' days'
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
write(8,'(3x,11i9)') (kc(j),j=0,nz,nzp)
write(6,'(3x,11i9)') (kc(j),j=0,nz,nzp)
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
do k= 0,nz,nzp
write(8,'(i3,11f9.3)') k,(head(j,k),j=0,nz,nzp)
write(6,'(i3,11f9.3)') k,(head(j,k),j=0,nz,nzp)
enddo
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
write(8,'(a,f23.3,a)')'Initial aquifer volume= ',volinit,' hm3'
write(6,'(a,f23.3,a)')'Initial aquifer volume= ',volinit,' hm3'
write(8,'(a,f23.3,a)')'Pumped volume (per Dt)= ',sumpumpvolume,
1' hm3'
write(6,'(a,f23.3,a)')'Pumped volume (per Dt)= ',sumpumpvolume,
1' hm3'
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
else
stop 'indm=1'
endif !430
c-----main time loop------------------------------------------------------------------
n= 0
write(23,*)n,head(nz/2,nz/2)
do 50 n= 1,nt
ts= ts + dt
td= ts/nsd
c-----calculation of head at advanced time level--------------------------------------
if(indm.eq.0) then !490
do j= 1,nz1
do k= 1,nz1
heada(j,k)=
1drey*(0.25*(head(j-1,k)+head(j,k+1)+head(j+1,k)+head(j,k-1)))
2+ drey1*head(j,k)
enddo
enddo
if(indr.eq.0.and.indx.eq.0) then !500
continue
else !500
do j= 1,nz1
do k= 1,nz1
heada(j,k)= heada(j,k) + rain(j,k) + xirr(j,k)
enddo
enddo
endif !500
if(indp.eq.0) then !510
continue
else !510
do k= 1,nz1
do j= 1,nz1
heada(j,k)= heada(j,k) - pump(j,k)
enddo
enddo
endif !510
else !490
do j= 1,nz1
do k= kleft(j)+1,kright(j)-1
heada(j,k)= 0.25*(head(j-1,k)+head(j,k+1)+head(j+1,k)+head(j,k-1))
enddo
enddo
if(indr.eq.0.and.indx.eq.0) then !515
continue
else !515
do j= 1,nz1
do k= kleft(j)+1,kright(j)-1
heada(j,k)= heada(j,k) + rain(j,k) + xirr(j,k)
enddo
enddo
endif !515
if(indp.eq.0) then !518
continue
else !518
do j= 1,nz1
do k= kleft(j)+1,kright(j)-1
heada(j,k)= heada(j,k) - pump(j,k)
enddo
enddo
endif !518
endif !490
c-----calculation of boundary conditions----------------------------------------------
c-----inflow: linear extrapolation from interior points------------------------------
if(indm.eq.1) then !600
if(indb.eq.0) then !610
do k= kleft(0)+1,kright(0)-1
heada(0,k)= 2*heada(1,k) - heada(2,k)
enddo
do k= kleft(nz)+1,kright(nz)-1
heada(nz,k)= 2*heada(nz1,k) - heada(nz-2,k)
enddo
do j= 1,nz1
heada(j,kleft(j))= 2*heada(j,kleft(j+1)) - heada(j,kleft(j+2))
enddo
do j= 1,nz1
heada(j,kright(nz))= 2*heada(j,kright(nz-1))
1 - heada(j,kright(nz-2))
enddo
c-----noflow: repetition from adjacent interior point--------------------------------
else !610
do k= kleft(0)+1,kright(0)-1
heada(0,k)= heada(1,k)
enddo
do k= kleft(nz)+1,kright(nz)-1
heada(nz,k)= heada(nz1,k)
enddo
do j= 1,nz1
heada(j,kleft(j))= heada(j,kleft(j)+1)
enddo
do j= 1,nz1
heada(j,kright(nz))= heada(j,kright(nz-1))
enddo
endif !610
elseif(indm.eq.0) then !600
c-----inflow: linear extrapolation from interior points------------------------------
if(indb.eq.0) then !611
do k= 1,nz1
heada(0,k)= 2*heada(1,k) - heada(2,k)
enddo
do k= 1,nz1
heada(nz,k)= 2*heada(nz1,k) - heada(nz-2,k)
enddo
do j= 1,nz1
heada(j,0)= 2*heada(j,1) - heada(j,2)
enddo
do j= 1,nz1
heada(j,nz)= 2*heada(j,nz1) - heada(j,nz-2)
enddo
sumlatinflow= funvolatinf(head,heada,mz,nz,nz1,trc,dt)
sumlatinflowaqui= sumlatinflow/st
c-----noflow: repetition from adjacent interior point--------------------------------
else !611
do k= 1,nz1
heada(0,k)= heada(1,k)
enddo
do k= 1,nz1
heada(nz,k)= heada(nz1,k)
enddo
do j= 1,nz1
heada(j,0)= heada(j,1)
enddo
do j= 1,nz1
heada(j,nz)= heada(j,nz1)
enddo
endif !611
endif !600
c-----calculation of corner points----------------------------------------------------
if(indm.eq.0) then !690
if(indb.eq.0) then !700
heada(0,0)= 0.5*(heada(0,1) + heada(1,0))
heada(nz,0)= 0.5*(heada(nz1,0) + heada(nz,1))
heada(0,nz)= 0.5*(heada(0,nz1) + heada(1,nz))
heada(nz,nz)= 0.5*(heada(nz1,nz) + heada(nz,nz1))
else !700
heada(0,0)= heada(1,1)
heada(0,nz)= heada(1,nz1)
heada(nz,0)= heada(nz1,1)
heada(nz,nz)= heada(nz1,nz1)
endif !700
else !690
continue ! no corner point calculation when indm=1
endif !690
c-----matrix relocation (ha into h)---------------------------------------------------
if(indm.eq.0) then !800
do j= 0,nz
do k= 0,nz
head(j,k)= heada(j,k)
enddo
enddo
else !800
do j= 0,nz
do k= kleft(j),kright(j)
head(j,k)= heada(j,k)
enddo
enddo
endif !800
c-----calculating aquifer volume----------------------------------------------------
do j= 0,nz
do k= 0,nz
depth(j,k)= head(j,k) - bott(j,k)
enddo
enddo
volcalc= volall(depth,dssq,mz,nz,nz1)
volumechange= volcalc0 - volcalc - sumpumpvolumeaqui
1 + sumlatinflowaqui
cumvolumechange= cumvolumechange + volumechange
percvolcalc= 100.*volcalc/volinit
percvolume= 100.*(volinit + cumvolumechange)/volinit
volcalc0= volcalc
c-----printing every np intervals-----------------------------------------------------
if(indm.eq.0) then !810
if(mod(n,np).eq.0) then !880
tyear= td/xndy + 1
write(8,'(a,i6,a,i2,a)')
1'Time= ',int(td),' days (Year ',int(tyear),')'
write(6,'(a,i6,a,i2,a)')
1'Time= ',int(td),' days (Year ',int(tyear),')'
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
write(8,'(3x,11i9)') (kc(j),j=0,nz,nzp)
write(6,'(3x,11i9)') (kc(j),j=0,nz,nzp)
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
do k= 0,nz,nzp
write(8,'(i3,11f9.3)') k,(head(j,k),j=0,nz,nzp)
write(6,'(i3,11f9.3)') k,(head(j,k),j=0,nz,nzp)
enddo
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
c-----printing volume data-----------------------------------------------------------
write(8,'(a,f12.3,a)')
1'Calculated aquifer volume= ',volcalc,' hm3'
write(6,'(a,f12.3,a)')
1'Calculated aquifer volume= ',volcalc,' hm3'
write(8,'(a,f12.3,a)')
1'Calculated to initial aquifer volume=',percvolcalc,' %'
write(6,'(a,f12.3,a)')
1'Calculated to initial aquifer volume=',percvolcalc,' %'
if(indp.eq.1) then !790
write(8,'(a,f12.3,a)')
1'Cumulative volume change= ',cumvolumechange,' hm3'
write(6,'(a,f12.3,a)')
1'Cumulative volume change= ',cumvolumechange,' hm3'
write(8,'(a,f12.3,a)')
1'Volume conservation= ',percvolume,' %'
write(6,'(a,f12.3,a)')
1'Volume conservation= ',percvolume,' %'
endif !790
write(23,*)n,head(nz/2,nz/2)
write(8,100)
write(6,100)
c write(8,'((1h-))')
c write(6,'((1h-))')
endif !880
endif !810
50 enddo
end
c-----end of program twodgw-----------------------------------------------------------
function volall(depth,dssq,mz,nz,nz1)
dimension depth(0:mz,0:mz)
dssq2= dssq/2
dssq4= dssq/4
volall= 0.
sumdepth= 0.
do k= 1,nz1
do j= 1,nz1
sumdepth= sumdepth + depth(j,k)
enddo
enddo
volall= volall + sumdepth*dssq
sumdepth= 0.
do j= 1,nz1
sumdepth= sumdepth + depth(j,0)
enddo
do j= 1,nz1
sumdepth= sumdepth + depth(j,nz)
enddo
do k= 1,nz1
sumdepth= sumdepth + depth(0,k)
enddo
do k= 1,nz1
sumdepth= sumdepth + depth(nz,k)
enddo
volall= volall + sumdepth*dssq2
sumdepth= depth(0,0) + depth(nz,0) + depth(0,nz) + depth(nz,nz)
volall= volall + sumdepth*dssq4
volall= volall/1.E+6 ! conversion to hm3
end
c-----end of function volall----------------------------------------------------------
function funvolatinf(head,heada,mz,nz,nz1,trc,dt)
dimension head(0:mz,0:mz)
dimension heada(0:mz,0:mz)
fun= 0.
do j= 1,nz
fun= fun + 0.5*trc*(head(j,0)-head(j,1)+heada(j,0)-heada(j,1))
enddo
do j= 1,nz
fun= fun
1 + 0.5*trc*(head(j,nz)-head(j,nz1)+heada(j,nz)-heada(j,nz1))
enddo
do k= 1,nz
fun= fun + 0.5*trc*(head(0,k)-head(1,k)+heada(0,k)-heada(1,k))
enddo
do k= 1,nz
fun= fun
1 + 0.5*trc*(head(nz,k)-head(nz1,k)+heada(nz,k)-heada(nz1,k))
enddo
funvolatinf= fun*dt/1.E+6
end
c-----end of function funvolatinf-----------------------------------------------------