module module_GW_baseflow
#ifdef MPP_LAND
use module_mpp_land
#endif
implicit none
#include "gw_field_include.inc"
#include "rt_include.inc"
#include "namelist.inc"
contains
!------------------------------------------------------------------------------
!DJG Simple GW Bucket Model
!------------------------------------------------------------------------------
subroutine simp_gw_buck(ix,jx,ixrt,jxrt,numbasns,basns_area,&
gwsubbasmsk, runoff1x, runoff2x, z_gwsubbas, qin_gwsubbas,&
qout_gwsubbas,qinflowbase,gw_strm_msk,gwbas_pix_ct,dist,DT,&
C,ex,z_mx,GWBASESWCRT,OVRTSWCRT)
implicit none
!!!Declarations...
integer, intent(in) :: ix,jx,ixrt,jxrt
integer, intent(in) :: numbasns
integer, intent(in), dimension(ix,jx) :: gwsubbasmsk
real, intent(in), dimension(ix,jx) :: runoff2x
real, intent(in), dimension(ix,jx) :: runoff1x
real, intent(in) :: basns_area(numbasns),dist(ixrt,jxrt,9),DT
real, intent(in),dimension(numbasns) :: C,ex,z_mx
real, intent(out),dimension(numbasns) :: qout_gwsubbas
real, intent(out),dimension(numbasns) :: qin_gwsubbas
real, intent(out),dimension(numbasns) :: z_gwsubbas
real, intent(out),dimension(ixrt,jxrt) :: qinflowbase
integer, intent(in),dimension(ixrt,jxrt) :: gw_strm_msk
integer, intent(in) :: GWBASESWCRT
integer, intent(in) :: OVRTSWCRT
real*8, dimension(numbasns) :: sum_perc8,ct_bas8
real, dimension(numbasns) :: sum_perc
real, dimension(numbasns) :: net_perc
real, dimension(numbasns) :: ct_bas
real, dimension(numbasns) :: gwbas_pix_ct
integer :: i,j,bas
real :: zbastmp
character(len=19) :: header
character(len=1) :: jnk
!!!Initialize variables...
ct_bas8 = 0
sum_perc8 = 0.
net_perc = 0.
qout_gwsubbas = 0.
qin_gwsubbas = 0.
!!!Calculate aggregated percolation from deep runoff into GW basins...
do i=1,ix
do j=1,jx
do bas=1,numbasns
if(gwsubbasmsk(i,j).eq.bas) then
if(OVRTSWCRT.ne.0) then
sum_perc8(bas) = sum_perc8(bas)+runoff2x(i,j) !Add only drainage to bucket...runoff2x in (mm)
else
sum_perc8(bas) = sum_perc8(bas)+runoff1x(i,j)+runoff2x(i,j) !Add sfc water & drainage to bucket...runoff1x and runoff2x in (mm)
end if
ct_bas8(bas) = ct_bas8(bas) + 1
end if
end do
end do
end do
#ifdef MPP_LAND
call sum_real8(sum_perc8,numbasns)
call sum_real8(ct_bas8,numbasns)
#endif
sum_perc = sum_perc8
ct_bas = ct_bas8
!!!Loop through GW basins to adjust for inflow/outflow
DO bas=1,numbasns ! Loop for GW bucket calcs...
#ifdef MPP_LAND
if(ct_bas(bas) .gt. 0) then
#endif
net_perc(bas) = sum_perc(bas) / ct_bas(bas) !units (mm)
qin_gwsubbas(bas) = net_perc(bas)/1000. * ct_bas(bas) * basns_area(bas) !units (m^3)
!Adjust level of GW depth...(conceptual GW bucket units (m))
z_gwsubbas(bas) = z_gwsubbas(bas) + net_perc(bas) / 1000.0 ! (m)
zbastmp = z_gwsubbas(bas)
!Calculate baseflow as a function of GW depth...
if(GWBASESWCRT.eq.1) then !active exponential bucket...
! Assuming and exponential relation between z and Q...
qout_gwsubbas(bas) = C(bas)*EXP(ex(bas)*z_gwsubbas(bas)/z_mx(bas)) !Exp.model. q_out (m^3/s)
!Adjust level of GW depth...
z_gwsubbas(bas) = z_gwsubbas(bas) - qout_gwsubbas(bas)*DT &
/ (ct_bas(bas)*basns_area(bas)) !units(m)
elseif (GWBASESWCRT.eq.2) then !Pass through/steady-state bucket
! Assuming a steady-state (inflow=outflow) model...
qout_gwsubbas(bas) = qin_gwsubbas(bas) !steady-state model...(m^3)
end if
#ifdef MPP_LAND
endif
#endif
END DO ! End loop for GW bucket calcs...
!!!Distribute basin integrated baseflow to stream pixels as stream 'inflow'...
qinflowbase = 0.
do i=1,ixrt
do j=1,jxrt
!!! -simple uniform disaggregation (8.31.06)
if (gw_strm_msk(i,j).gt.0) then
if(GWBASESWCRT.eq.1) then !calc stream inflow from exponential bucket... (m^3/s to mm)
qinflowbase(i,j) = qout_gwsubbas(gw_strm_msk(i,j))*1000.*DT/ &
gwbas_pix_ct(gw_strm_msk(i,j))/dist(i,j,9) !(mm)
elseif (GWBASESWCRT.eq.2) then !calc stream inflow from passthrough/steady-state bucket (m^3 to mm)
qinflowbase(i,j) = qout_gwsubbas(gw_strm_msk(i,j))*1000./ &
gwbas_pix_ct(gw_strm_msk(i,j))/dist(i,j,9) !(mm)
end if
end if
end do
end do
!!! - weighted redistribution...(need to pass accum weights (slope) in...)
! NOT FINISHED just BASIC framework...
! do bas=1,numbasns
! do k=1,gwbas_pix_ct(bas)
! qinflowbase(i,j) = k*slope
! end do
! end do
return
!------------------------------------------------------------------------------
End subroutine simp_gw_buck
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!DJG Wedge-Aquifer Scheme (TBA)
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
!DJG TOPMODEL Scheme (TBA)
!------------------------------------------------------------------------------
#ifdef MPP_LAND
subroutine pix_ct_1(in_gw_strm_msk,ixrt,jxrt,gwbas_pix_ct,numbasns)
USE module_mpp_land
implicit none
integer :: i,j,ixrt,jxrt,numbasns, bas
integer,dimension(ixrt,jxrt) :: in_gw_strm_msk
integer,dimension(global_rt_nx,global_rt_ny) :: gw_strm_msk
real,dimension(numbasns) :: gwbas_pix_ct
gw_strm_msk = 0
call write_IO_rt_int(in_gw_strm_msk, gw_strm_msk)
if(my_id .eq. IO_id) then
gwbas_pix_ct = 0.
do bas = 1,numbasns
do i=1,global_rt_nx
do j=1,global_rt_ny
if(gw_strm_msk(i,j) .eq. bas) then
gwbas_pix_ct(gw_strm_msk(i,j)) = gwbas_pix_ct(gw_strm_msk(i,j)) &
+ 1.0
endif
end do
end do
end do
end if
call mpp_land_bcast_real(numbasns,gwbas_pix_ct)
return
end subroutine pix_ct_1
#endif
!------------------------------------------------------------------------------
! Benjamin Fersch 2d groundwater model
!------------------------------------------------------------------------------
subroutine gw2d_ini(did,dt,dx)
use module_GW_baseflow_data, only: gw2d
implicit none
integer did
real dt,dx
gw2d(did)%dx=dx
gw2d(did)%dt=dt
! bftodo: develop proper landtype mask
gw2d(did)%compres=0. ! currently not implemented
return
end subroutine gw2d_ini
subroutine gw2d_allocate(did, ix, jx, nsoil)
use module_GW_baseflow_data, only: gw2d
implicit none
integer ix, jx, nsoil
integer istatus, did
if(gw2d(did)%allo_status .eq. 1) return
gw2d(did)%allo_status = 1
gw2d(did)%ix = ix
gw2d(did)%jx = jx
allocate(gw2d(did)%ltype (ix,jx))
allocate(gw2d(did)%elev (ix,jx))
allocate(gw2d(did)%bot (ix,jx))
allocate(gw2d(did)%hycond (ix,jx))
allocate(gw2d(did)%poros (ix,jx))
allocate(gw2d(did)%compres(ix,jx))
allocate(gw2d(did)%ho (ix,jx))
allocate(gw2d(did)%h (ix,jx))
allocate(gw2d(did)%convgw (ix,jx))
! allocate(gw2d(did)% (ix,jx))
end subroutine gw2d_allocate
subroutine gwstep(ix, jx, dx, &
ltype, elev, bot, &
hycond, poros, compres, &
ho, h, convgw, &
ebot, eocn, &
dt, istep)
! #else
! dx, istep, dt, & !supplied
! ims,ime,jms,jme,its,ite,jts,jte, & !supplied
! ids,ide,jds,jde,ifs,ife,jfs,jfe) !supplied
! #endif
! New (volug): calling routines use change in head, convgw = d(h-ho)/dt.
! Steps ground-water hydrology (head) through one timestep.
! Modified from Prickett and Lonnquist (1971), basic one-layer aquifer
! simulation program, with mods by Zhongbo Yu(1997).
! Solves S.dh/dt = d/dx(T.dh/dx) + d/dy(T.dh/dy) + "external sources"
! for a single layer, where h is head, S is storage coeff and T is
! transmissivity. 3-D arrays in main program (hycond,poros,h,bot)
! are 2-D here, since only a single (uppermost) layer is solved.
! Uses an iterative time-implicit ADI method.
! use module_hms_constants
integer, intent(in) :: ix, jx
integer, intent(in), dimension(ix,jx) :: ltype ! land-sfc type (supp)
real, intent(in), dimension(ix,jx) :: &
elev, & ! elev/bathymetry of sfc rel to sl (m) (supp)
bot, & ! elev. aquifer bottom rel to sl (m) (supp)
hycond, & ! hydraulic conductivity (m/s per m/m) (supp)
poros, & ! porosity (m3/m3) (supp)
compres, & ! compressibility (1/Pa) (supp)
ho ! head at start of timestep (m) (supp)
real, intent(inout), dimension(ix,jx) :: &
h, & ! head, after ghmcompute (m) (ret)
convgw ! convergence due to gw flow (m/s) (ret)
real, intent(inout) :: ebot, eocn
integer :: istep !, dt
real, intent(in) :: dt, dx
! #endif
! eocn = mean spurious sink for h_ocn = sealev fix (m/s)(ret)
! This equals the total ground-water flow across
! land->ocean boundaries.
! ebot = mean spurious source for "bot" fix (m/s) (returned)
! time = elapsed time from start of run (sec)
! dt = timestep length (sec)
! istep = timestep counter
! Local arrays:
real, dimension(ix,jx) :: sf2 ! storage coefficient (m3 of h2o / bulk m3)
real, dimension(ix,jx,2) :: t ! transmissivity (m2/s)..1 for N-S,..2 for E-W
real, dimension(0:ix+jx) :: b,g ! work arrays
real, parameter :: botinc = 0.01 ! re-wetting increment to fix h < bot
! parameter (botinc = 0. ) ! re-wetting increment to fix h < bot
! (m); else no flow into dry cells
real, parameter :: delskip = 0.005 ! av.|dhead| value for iter.skip out(m)
integer, parameter :: itermax = 10 ! maximum number of iterations
integer, parameter :: itermin = 3 ! minimum number of iterations
real, parameter :: sealev = -1. ! sea-level elevation (m)
! die müssen noch sortiert, geprüft und aufgeräumt werden
integer :: &
iter, &
j, &
i, &
jp, &
ip, &
ii, &
n, &
jj, &
ierr, &
ier
! real :: su, sc, shp, bb, aa, cc, w, zz, tareal, dtoa, dtot
real :: &
dy, &
e, &
su, &
sc, &
shp, &
bb, &
dd, &
aa, &
cc, &
w, &
ha, &
delcur, &
dtot, &
dtoa, &
darea, &
tareal, &
zz
#ifdef MPP_LAND
real mpiDelcur
integer mpiSize
#endif
dy = dx
darea = dx*dy
call scopy (ix*jx, ho, 1, h, 1)
! Top of iterative loop for ADI solution
iter = 0
!~~~~~~~~~~~~~
80 continue
!~~~~~~~~~~~~~
iter = iter+1
#ifdef MPP_LAND
call MPP_LAND_COM_REAL(h, ix, jx, 99)
#endif
e = 0. ! absolute changes in head (for iteration control)
! eocn = 0. ! accumulated fixes for h = 0 over ocean (diag)
! ebot = 0. ! accumulated fixes for h < bot (diagnostic)
! Set storage coefficient (sf2)
! #ifdef HMSWRF
!
tareal = 0.
!
! do j=jfs,jfe
! do i=ifs,ife
!
!
! #else
do j=1,jx
do i=1,ix
if(ltype(i,j) .ge. 1) tareal = tareal + darea
! #endif
! unconfined water table (h < e): V = poros*(h-b)
! dV/dh = poros
! saturated to surface (h >= e) : V = poros*(e-b) + (h-e)
! dV/dh = 1
! (compressibility is ignored)
!
! su = poros(i,j)*(1.-theta(i,j)) ! old (pre-volug)
su = poros(i,j) ! new (volug)
sc = 1.
if (ho(i,j).le.elev(i,j) .and. h(i,j).le.elev(i,j)) then
sf2(i,j) = su
else if (ho(i,j).ge.elev(i,j) .and. h(i,j).ge.elev(i,j)) then
sf2(i,j) = sc
else if (ho(i,j).le.elev(i,j) .and. h(i,j).ge.elev(i,j)) then
shp = sf2(i,j) * (h(i,j) - ho(i,j))
sf2(i,j) = shp * sc / (shp - (su-sc)*(elev(i,j)-ho(i,j)))
else if (ho(i,j).ge.elev(i,j) .and. h(i,j).le.elev(i,j)) then
shp = sf2(i,j) * (ho(i,j) - h(i,j))
sf2(i,j) = shp * su / (shp + (su-sc)*(ho(i,j)-elev(i,j)))
endif
enddo
enddo
#ifdef MPP_LAND
! communicate storage coefficient
call MPP_LAND_COM_REAL(sf2, ix, jx, 99)
#endif
!==========================
! Column calculations
!==========================
! Set transmissivities. Use min(h,elev)-bot instead of h-bot,
! since if h > elev, thickness of groundwater flow is just
! elev-bot.
! #ifdef HMSWRF
!
! do j=jfs,jfe
! jp = min (j+1,jfe)
! do i=ifs,ife
! ip = min (i+1,ife)
!
! #else
do j=1,jx
jp = min (j+1,jx)
do i=1,ix
ip = min (i+1,ix)
! #endif
t(i,j,2) = sqrt( abs( &
hycond(i, j)*(min(h(i ,j),elev(i ,j))-bot(i ,j)) &
*hycond(ip,j)*(min(h(ip,j),elev(ip,j))-bot(ip,j)) &
) ) &
! #ifdef HMSWRF
* (0.5*(dy+dy)) & ! in WRF the dx and dy are usually equal
/ (0.5*(dx+dx))
! #else
! * (0.5*(dy(i,j)+dy(ip,j))) &
! / (0.5*(dx(i,j)+dx(ip,j)))
! #endif
t(i,j,1) = sqrt( abs( &
hycond(i,j )*(min(h(i,j ),elev(i,j ))-bot(i,j )) &
*hycond(i,jp)*(min(h(i,jp),elev(i,jp))-bot(i,jp)) &
) ) &
! #ifdef HMSWRF
* (0.5*(dx+dx)) &
/ (0.5*(dy+dy))
! #else
! * (0.5*(dx(i,j)+dx(i,jp))) &
! / (0.5*(dy(i,j)+dy(i,jp)))
! #endif
enddo
enddo
#ifdef MPP_LAND
! communicate transmissivities in x and y direction
call MPP_LAND_COM_REAL(t(:,:,1), ix, jx, 99)
call MPP_LAND_COM_REAL(t(:,:,2), ix, jx, 99)
#endif
b = 0.
g = 0.
!-------------------
do 190 ii=1,ix
!-------------------
i=ii
if (mod(istep+iter,2).eq.1) i=ix-i+1
! calculate b and g arrays
!>>>>>>>>>>>>>>>>>>>>
do 170 j=1,jx
!>>>>>>>>>>>>>>>>>>>>
! bb = (sf2(i,j)/dt) * darea(i,j)
! dd = ( ho(i,j)*sf2(i,j)/dt ) * darea(i,j)
bb = (sf2(i,j)/dt) * darea
dd = ( ho(i,j)*sf2(i,j)/dt ) * darea
aa = 0.0
cc = 0.0
if (j-1) 90,100,90
90 aa = -t(i,j-1,1)
bb = bb + t(i,j-1,1)
100 if (j-jx) 110,120,110
110 cc = -t(i,j,1)
bb = bb + t(i,j,1)
120 if (i-1) 130,140,130
130 bb = bb + t(i-1,j,2)
dd = dd + h(i-1,j)*t(i-1,j,2)
140 if (i-ix) 150,160,150
150 bb = bb + t(i,j,2)
dd = dd + h(i+1,j)*t(i,j,2)
160 w = bb - aa*b(j-1)
b(j) = cc/w
g(j) = (dd-aa*g(j-1))/w
!>>>>>>>>>>>>>>>
170 continue
!>>>>>>>>>>>>>>>
! re-estimate heads
e = e + abs(h(i,jx)-g(jx))
h(i,jx) = g(jx)
n = jx-1
180 if (n.eq.0) goto 185
ha = g(n) - b(n)*h(i,n+1)
e = e + abs(ha-h(i,n))
h(i,n) = ha
n = n-1
goto 180
185 continue
!-------------
190 continue
!-------------
#ifdef MPP_LAND
call MPP_LAND_COM_REAL(h, ix, jx, 99)
#endif
!=======================
! Row calculations
!=======================
! set transmissivities (same as above)
do j=1,jx
jp = min (j+1,jx)
do i=1,ix
ip = min (i+1,ix)
t(i,j,2) = sqrt( abs( &
hycond(i, j)*(min(h(i ,j),elev(i ,j))-bot(i ,j)) &
*hycond(ip,j)*(min(h(ip,j),elev(ip,j))-bot(ip,j)) &
) ) &
! * (0.5*(dy(i,j)+dy(ip,j))) &
! / (0.5*(dx(i,j)+dx(ip,j)))
* (0.5*(dy+dy)) &
/ (0.5*(dx+dx))
t(i,j,1) = sqrt( abs( &
hycond(i,j )*(min(h(i,j ),elev(i,j ))-bot(i,j )) &
*hycond(i,jp)*(min(h(i,jp),elev(i,jp))-bot(i,jp)) &
) ) &
* (0.5*(dx+dx)) &
/ (0.5*(dy+dy))
enddo
enddo
#ifdef MPP_LAND
! communicate transmissivities in x and y direction
call MPP_LAND_COM_REAL(t(:,:,1), ix, jx, 99)
call MPP_LAND_COM_REAL(t(:,:,2), ix, jx, 99)
#endif
b = 0.
g = 0.
!-------------------
do 300 jj=1,jx
!-------------------
j=jj
if (mod(istep+iter,2).eq.1) j = jx-j+1
! calculate b and g arrays
!>>>>>>>>>>>>>>>>>>>>
do 280 i=1,ix
!>>>>>>>>>>>>>>>>>>>>
! bb = (sf2(i,j)/dt) * darea(i,j)
! dd = ( ho(i,j)*sf2(i,j)/dt ) * darea(i,j)
bb = (sf2(i,j)/dt) * darea
dd = ( ho(i,j)*sf2(i,j)/dt ) * darea
aa = 0.0
cc = 0.0
if (j-1) 200,210,200
200 bb = bb + t(i,j-1,1)
dd = dd + h(i,j-1)*t(i,j-1,1)
210 if (j-jx) 220,230,220
220 dd = dd + h(i,j+1)*t(i,j,1)
bb = bb + t(i,j,1)
230 if (i-1) 240,250,240
240 bb = bb + t(i-1,j,2)
aa = -t(i-1,j,2)
250 if (i-ix) 260,270,260
260 bb = bb + t(i,j,2)
cc = -t(i,j,2)
270 w = bb - aa*b(i-1)
b(i) = cc/w
g(i) = (dd-aa*g(i-1))/w
!>>>>>>>>>>>>>>>
280 continue
!>>>>>>>>>>>>>>>
! re-estimate heads
e = e + abs(h(ix,j)-g(ix))
h(ix,j) = g(ix)
n = ix-1
290 if (n.eq.0) goto 295
ha = g(n)-b(n)*h(n+1,j)
e = e + abs(h(n,j)-ha)
h(n,j) = ha
n = n-1
goto 290
295 continue
!-------------
300 continue
!-------------
! fix head < bottom of aquifer
! #endif
!
! #ifdef HMSWRF
!
! do j=jfs,jfe
! do i=ifs,ife
!
! #else
do j=1,jx
do i=1,ix
! #endif
if (ltype(i,j).eq.1 .and. h(i,j).le.bot(i,j)+botinc) then
! #ifndef HMSWRF
e = e + bot(i,j) + botinc - h(i,j)
! ebot = ebot + (bot(i,j)+botinc-h(i,j))*sf2(i,j)*darea(i,j)
ebot = ebot + (bot(i,j)+botinc-h(i,j))*sf2(i,j)*darea
! #endif
h(i,j) = bot(i,j) + botinc
endif
enddo
enddo
! maintain head = sea level for ocean (only for adjacent ocean,
! rest has hycond=0)
! #ifdef HMSWRF
!
! do j=jfs,jfe
! do i=its,ife
!
! #else
do j=1,jx
do i=1,ix
! #endif
if (ltype(i,j).eq.2) then
! #ifndef HMSWRF
eocn = eocn + (h(i,j)-sealev)*sf2(i,j)*darea
! eocn = eocn + (h(i,j)-sealev)*sf2(i,j)*darea(i,j)
! #endif
h(i,j) = sealev
endif
enddo
enddo
! Loop back for next ADI iteration
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! #ifdef HMSWRF
! delcur = e/(xdim*ydim)
! #else
delcur = e/(ix*jx)
! #endif
#ifdef MPP_LAND
call mpi_reduce(delcur, mpiDelcur, 1, MPI_REAL, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
call MPI_COMM_SIZE( MPI_COMM_WORLD, mpiSize, ierr )
mpiDelcur = mpiDelcur/mpiSize
call mpi_bcast(delcur, 1, mpi_real, 0, MPI_COMM_WORLD, ierr)
#endif
if ( (delcur.gt.delskip*dt/86400. .and. iter.lt.itermax) &
.or. iter.lt.itermin ) then
goto 80
else
endif
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Compute convergence rate due to ground water flow (returned)
! #ifdef HMSWRF
!
! do j=jfs,jfe
! do i=ifs,ife
!
! #else
do j=1,jx
do i=1,ix
! #endif
if (ltype(i,j).eq.1) then
convgw(i,j) = sf2(i,j) * (h(i,j)-ho(i,j)) / dt
else
convgw(i,j) = 0.
endif
enddo
enddo
! Diagnostic water conservation check for this timestep
dtot = 0. ! total change in water storage (m3)
dtoa = 0.
! #ifdef HMSWRF
!
! do j=jts,jte
! do i=its,ite
!
! #else
do j=1,jx
do i=1,ix
! #endif
if (ltype(i,j).eq.1) then
! #ifdef HMSWRF
dtot = dtot + sf2(i,j) *(h(i,j)-ho(i,j)) * darea
dtoa = dtoa + sf2(i,j) * abs(h(i,j)-ho(i,j)) * darea
! #else
! dtot = dtot + sf2(i,j) *(h(i,j)-ho(i,j)) * darea(i,j)
! dtoa = dtoa + sf2(i,j) * abs(h(i,j)-ho(i,j)) * darea(i,j)
! #endif
endif
enddo
enddo
dtot = (dtot/tareal)/dt ! convert to m/s, rel to land area
dtoa = (dtoa/tareal)/dt
eocn = (eocn/tareal)/dt
ebot = (ebot/tareal)/dt
zz = 1.e3 * 86400. ! convert printout to mm/day
#ifdef HYDRO_D
write (*,900) &
dtot*zz, dtoa*zz, -eocn*zz, ebot*zz, &
(dtot-(-eocn+ebot))*zz
#endif
900 format &
(3x,' dh/dt |dh/dt| ocnflx botfix',&
' ',' ghmerror' &
! /3x,4f9.4,2(9x),e14.4)
/3x,5(e14.4))
return
end subroutine gwstep
SUBROUTINE SCOPY (NT, ARR, INCA, BRR, INCB)
!
! Copies array ARR to BRR, incrementing by INCA and INCB
! respectively, up to a total length of NT words of ARR.
! (Same as Cray SCOPY.)
!
real, DIMENSION(*) :: ARR, BRR
integer :: ia, nt, inca, incb, ib
!
IB = 1
DO 10 IA=1,NT,INCA
BRR(IB) = ARR(IA)
IB = IB + INCB
10 CONTINUE
!
RETURN
END SUBROUTINE SCOPY
end module module_GW_baseflow