|
XBeach
|
00001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00002 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00003 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00004 ! Jaap van Thiel de Vries, Robert McCall ! 00005 ! ! 00006 ! d.roelvink@unesco-ihe.org ! 00007 ! UNESCO-IHE Institute for Water Education ! 00008 ! P.O. Box 3015 ! 00009 ! 2601 DA Delft ! 00010 ! The Netherlands ! 00011 ! ! 00012 ! This library is free software; you can redistribute it and/or ! 00013 ! modify it under the terms of the GNU Lesser General Public ! 00014 ! License as published by the Free Software Foundation; either ! 00015 ! version 2.1 of the License, or (at your option) any later version. ! 00016 ! ! 00017 ! This library is distributed in the hope that it will be useful, ! 00018 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00019 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00020 ! Lesser General Public License for more details. ! 00021 ! ! 00022 ! You should have received a copy of the GNU Lesser General Public ! 00023 ! License along with this library; if not, write to the Free Software ! 00024 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00025 ! USA ! 00026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00027 module groundwaterflow 00028 implicit none 00029 save 00030 private ! set private default 00031 ! set these public to users of groundwater module 00032 public gw_init 00033 public gw_bc 00034 public gwflow 00035 00036 contains 00037 00038 00039 00040 subroutine gw_init(s,par) 00041 use params, only: parameters 00042 use xmpi_module 00043 use spaceparams 00044 00045 IMPLICIT NONE 00046 00047 type(parameters) :: par 00048 type(spacepars) :: s 00049 00050 integer :: i,j 00051 00052 if(.not. xmaster) return ! Robert: change in way called from libxbeach 00053 00054 allocate (s%gwhead(s%nx+1,s%ny+1)) 00055 allocate (s%gwheadb(s%nx+1,s%ny+1)) 00056 allocate (s%gwlevel(s%nx+1,s%ny+1)) 00057 allocate (s%gwheight(s%nx+1,s%ny+1)) 00058 allocate (s%gwu(s%nx+1,s%ny+1)) 00059 allocate (s%gwv(s%nx+1,s%ny+1)) 00060 allocate (s%gwqx(s%nx+1,s%ny+1)) 00061 allocate (s%gwqy(s%nx+1,s%ny+1)) 00062 allocate (s%gww(s%nx+1,s%ny+1)) 00063 allocate (s%infil(s%nx+1,s%ny+1)) 00064 allocate (s%gwbottom(s%nx+1,s%ny+1)) 00065 allocate (s%dinfil(s%nx+1,s%ny+1)) 00066 allocate (s%gw0back(2,s%ny+1)) 00067 allocate (s%gwcurv(s%nx+1,s%ny+1)) 00068 allocate (s%Kx(s%nx+1,s%ny+1)) 00069 allocate (s%Ky(s%nx+1,s%ny+1)) 00070 allocate (s%Kz(s%nx+1,s%ny+1)) 00071 allocate (s%Kzinf(s%nx+1,s%ny+1)) 00072 00073 s%gww=0.d0 00074 s%infil = 0.d0 00075 00076 if (par%gwflow==1) then 00077 if (par%aquiferbotfile==' ') then ! Not a filename 00078 s%gwbottom=par%aquiferbot 00079 else 00080 open(31,file=trim(par%aquiferbotfile)) 00081 do j=1,s%ny+1 00082 read(31,*)(s%gwbottom(i,j),i=1,s%nx+1) 00083 end do 00084 close(31) 00085 endif 00086 00087 if (par%gw0file==' ') then ! Not a filename 00088 s%gwhead=par%gw0 00089 else 00090 open(31,file=trim(par%gw0file)) 00091 do j=1,s%ny+1 00092 read(31,*)(s%gwhead(i,j),i=1,s%nx+1) 00093 end do 00094 close(31) 00095 endif 00096 if (xmpi_istop) then 00097 s%gwbottom(1,:) = s%gwbottom(2,:) 00098 endif 00099 if (xmpi_isbot) then 00100 s%gwbottom(s%nx+1,:) = s%gwbottom(s%nx,:) 00101 endif 00102 if (xmpi_isleft .and. s%ny>0) then 00103 s%gwbottom(:,1) = s%gwbottom(:,2) 00104 endif 00105 if (xmpi_isright .and. s%ny>0) then 00106 s%gwbottom(:,s%ny+1) = s%gwbottom(:,s%ny) 00107 endif 00108 00109 s%gwbottom = min(s%gwbottom,s%zb) 00110 00111 s%gwhead=max(s%gwhead,s%gwbottom) 00112 s%gw0back=s%gwhead(s%nx:s%nx+1,:) 00113 00114 s%gwlevel=min(s%zb,s%gwhead) 00115 s%gwlevel=max(s%gwlevel,s%gwbottom) !+par%eps) 00116 00117 s%gwheight=s%gwlevel-s%gwbottom 00118 00119 s%gwu=0.d0 00120 s%gwv=0.d0 00121 s%gww=0.d0 00122 s%gwqx = 0.d0 00123 s%gwqy = 0.d0 00124 s%gwcurv=0.d0 00125 s%infil = 0.d0 00126 s%dinfil=max(par%dwetlayer/3.d0,0.02) ! Centroid of area influenced instantly by free surface level lies at dwetlayer/3 00127 endif 00128 end subroutine gw_init 00129 00130 subroutine gw_bc(s,par) 00131 use params, only: parameters 00132 use xmpi_module 00133 use spaceparams 00134 use paramsconst 00135 00136 IMPLICIT NONE 00137 00138 type(parameters) :: par 00139 type(spacepars) :: s 00140 00141 00142 s%gwbottom=min(s%gwbottom,s%zb) !-par%eps) 00143 00144 if(xmpi_istop) then 00145 ! s%gwhead(1,:)=s%zs0(1,:) 00146 s%gwhead(1,:)=s%gwhead(2,:) 00147 s%gwbottom(1,:) = s%gwbottom(2,:) 00148 endif 00149 if (xmpi_isbot) then 00150 if (par%tideloc==4 .or. (par%tideloc==2 .and. par%paulrevere==PAULREVERE_LAND)) then 00151 ! s%gwhead(s%nx+1,:)=s%zs0(s%nx+1,:) 00152 s%gwhead(s%nx+1,:)=s%gwhead(s%nx,:) 00153 else 00154 s%gwhead(s%nx+1,:)=s%gw0back(2,:) 00155 endif 00156 s%gwbottom(s%nx+1,:) = s%gwbottom(s%nx,:) 00157 endif 00158 if (xmpi_isleft .and. s%ny>0) then 00159 s%gwbottom(:,1) = s%gwbottom(:,2) 00160 s%gwhead(:,1)=s%gwhead(:,2) 00161 s%gwlevel(:,1)=s%gwlevel(:,2) 00162 endif 00163 if (xmpi_isright .and. s%ny>0) then 00164 s%gwbottom(:,s%ny+1) = s%gwbottom(:,s%ny) 00165 s%gwhead(:,s%ny+1)=s%gwhead(:,s%ny) 00166 s%gwlevel(:,s%ny+1)=s%gwlevel(:,s%ny) 00167 endif 00168 00169 #ifdef USEMPI 00170 call xmpi_shift_ee(s%gwhead) 00171 call xmpi_shift_ee(s%gwlevel) 00172 #endif 00173 if(xmpi_istop) then 00174 s%gwlevel(1,:)=min(s%gwhead(1,:),s%zb(1,:)) 00175 endif 00176 if (xmpi_isbot) then 00177 s%gwlevel(s%nx+1,:)=min(s%gwhead(s%nx+1,:),s%zb(s%nx+1,:)) 00178 endif 00179 00180 end subroutine gw_bc 00181 00182 subroutine gwflow(s,par) 00183 00184 use params, only: parameters 00185 use xmpi_module 00186 use spaceparams 00187 use constants, only: pi 00188 use paramsconst 00189 00190 IMPLICIT NONE 00191 00192 type(parameters) :: par 00193 type(spacepars),target :: s 00194 00195 ! internal variables 00196 integer :: i,j 00197 integer :: incl 00198 integer,parameter :: maxiter = 50 00199 real*8 :: dzr,wdt 00200 real*8,dimension(:,:),allocatable,save :: gwhu,gwhv,gwheadtop 00201 real*8,dimension(:,:),allocatable,save :: Kx,Ky,Kz,Kzinf 00202 real*8,dimension(:,:),allocatable,save :: Kxupd,Kyupd,Kzupd 00203 logical,dimension(:,:),allocatable,save :: connected 00204 real*8,dimension(:,:),allocatable,save :: fracdt 00205 real*8,dimension(:,:),allocatable,save :: zsupd,ratio 00206 real*8,dimension(:,:),allocatable,save :: dynpresupd 00207 real*8,dimension(:,:),allocatable,save :: infiluncon 00208 ! wwvv changed allocatable into pointer: 00209 ! wwvv see also later in this file 00210 real*8,dimension(:,:),pointer,save :: infilcon 00211 real*8,dimension(:,:),allocatable,save :: infilhorgw,infilhorsw 00212 real*8,dimension(:,:),allocatable,save :: gwumean 00213 real*8 :: factime 00214 logical :: initial 00215 real*8,save :: connectcrit 00216 00217 00218 initial = .false. 00219 00220 ! allocate variables 00221 if (.not. allocated(gwhu)) then 00222 allocate(gwhu(1:s%nx+1,1:s%ny+1)) 00223 allocate(gwhv(1:s%nx+1,1:s%ny+1)) 00224 allocate(gwheadtop(1:s%nx+1,1:s%ny+1)) 00225 00226 allocate(Kx(1:s%nx+1,1:s%ny+1)) 00227 allocate(Ky(1:s%nx+1,1:s%ny+1)) 00228 allocate(Kz(1:s%nx+1,1:s%ny+1)) 00229 allocate(Kzinf(1:s%nx+1,1:s%ny+1)) 00230 allocate(Kxupd(1:s%nx+1,1:s%ny+1)) 00231 allocate(Kyupd(1:s%nx+1,1:s%ny+1)) 00232 allocate(Kzupd(1:s%nx+1,1:s%ny+1)) 00233 00234 allocate(connected(1:s%nx+1,1:s%ny+1)) 00235 allocate(fracdt(1:s%nx+1,1:s%ny+1)) 00236 allocate(zsupd(1:s%nx+1,1:s%ny+1)) 00237 allocate(ratio(1:s%nx+1,1:s%ny+1)) 00238 00239 allocate(dynpresupd(1:s%nx+1,1:s%ny+1)) 00240 00241 allocate(infiluncon(1:s%nx+1,1:s%ny+1)) 00242 allocate(infilcon(1:s%nx+1,1:s%ny+1)) 00243 allocate(infilhorgw(1:s%nx+1,1:s%ny+1)) 00244 allocate(infilhorsw(1:s%nx+1,1:s%ny+1)) 00245 allocate(gwumean(1:s%nx+1,1:s%ny+1)) 00246 00247 ! initialise variables 00248 Kx = par%kx 00249 Ky = par%ky 00250 Kz = par%kz 00251 Kzinf = par%kz 00252 00253 Kxupd = Kx 00254 Kyupd = Ky 00255 Kzupd = Kz 00256 00257 infiluncon = 0.d0 00258 infilcon = 0.d0 00259 gwumean = 0.d0 00260 00261 dynpresupd = 0.d0 00262 00263 gwhu = 0.d0 00264 gwhv = 0.d0 00265 00266 initial = .true. 00267 if (par%gwnonh==0) then 00268 connectcrit = par%dwetlayer 00269 else 00270 connectcrit = par%eps 00271 endif 00272 00273 ratio = gw_calculate_smoothwetlayer(s,par%dwetlayer,par%px) 00274 endif 00275 ! 00276 ! Definition of horizontal groundwater velocity: 00277 ! Because there is not simple way of measuring flow in pores, Darcy relates to the 00278 ! flow velocity through the ground in terms of surface water (not pore water). 00279 ! Therefore, after the flux has been calculated the groundwater level change is 00280 ! divided by the porosity to get changes in terms of pore water. 00281 ! 00282 ! Definition of vertical groundwater/surface water exchange: 00283 ! infil is defined positive from sea to groundwater, gww is defined positive up 00284 ! (from groundwater to sea!). Both infil and gww are in volumes of surface water 00285 ! (not pore water). gwu and gwv are also in terms of surface water volume. 00286 ! 00287 ! Initialize 00288 s%gww = 0.d0 00289 s%infil = 0.d0 00290 infiluncon = 0.d0 00291 infilcon = 0.d0 00292 infilhorgw = 0.d0 00293 infilhorsw = 0.d0 00294 fracdt = 0.d0 00295 ! 00296 ! If using dynamic pressure (nonh==1), very short wave lengths mess things 00297 ! up because they are not resolved properly in the nonh pressure solver. 00298 ! To surpress this, carry out pressure averaging over time scale ~ 1/4 Trep 00299 if (par%wavemodel==WAVEMODEL_NONH) then 00300 factime = min(4*par%dt/par%Trep,1.d0) 00301 dynpresupd = (1-factime) * dynpresupd + factime * s%pres 00302 endif 00303 ! 00304 ! Infiltration and exfiltration are handled separately for places where the 00305 ! groundwater and surface water are connected (Poisson equation), and where 00306 ! they are disconnected (swash infiltration, seepage). "Connected" is set 00307 ! by the difference between the groundwater level and the bed level, and a 00308 ! numerical constant for stability 00309 where(s%wetz==1 .and. s%zb-s%gwlevel<connectcrit) 00310 connected = .true. 00311 elsewhere 00312 connected = .false. 00313 endwhere 00314 ! 00315 ! Handle infiltration in unconnected cells. Some cells may be 00316 ! unconnected for only part of the timestep, therefore a fraction 00317 ! of timestep is returned which was required to reach connection 00318 if (par%gwnonh==1) then 00319 call gw_unconnected_infil(s,par,Kzinf,connected,fracdt,infiluncon) 00320 ! where ((.not. connected) .and. (gwlevel>zb)) 00321 ! infiluncon = infiluncon - (gwlevel-zb)/par%dt*par%por 00322 ! fracdt = 1.d0 00323 ! endwhere 00324 where (s%gwlevel>s%zb) 00325 infiluncon = infiluncon - (s%gwlevel-s%zb)/par%dt*par%por 00326 endwhere 00327 endif 00328 ! 00329 ! 00330 ! Update groundwater level and water surface level to account for 00331 ! infiltration and exfiltration effects. Note infil is avarage 00332 ! infiltration rate over whole timestep 00333 s%gwlevel = s%gwlevel+par%dt*infiluncon/par%por 00334 zsupd = s%zs-par%dt*infiluncon 00335 ! 00336 ! Recalculate connected property, this is needed further on in the code 00337 where(zsupd-s%zb>=par%eps .and. s%zb-s%gwlevel<connectcrit) 00338 connected = .true. 00339 elsewhere 00340 connected = .false. 00341 endwhere 00342 ! 00343 ! Thickness of groundwater layer 00344 s%gwheight=s%gwlevel-s%gwbottom 00345 ! 00346 ! Determine intermediate aquifer depths (upwind scheme) 00347 call gw_calculate_interfaceheight(s,gwhu,gwhv,initial) 00348 ! 00349 ! Calculate head smoothing function over dwetlayer 00350 ratio = gw_calculate_smoothwetlayer(s,par%dwetlayer,par%px) 00351 ! 00352 ! Calculate top boundary condition for head (pressure on groundwater 00353 ! from surface water, or atmospheric pressure when dry). Takes into 00354 ! account a transition layer, specified by user to smooth transition 00355 ! from surface water head to atmospheric head 00356 if (par%wavemodel==WAVEMODEL_NONH) then 00357 gwheadtop = gwCalculateHeadTop(s%nx,s%ny,zsupd,s%gwlevel,s%wetz,ratio,1,par%g,dynpresupd) 00358 else 00359 gwheadtop = gwCalculateHeadTop(s%nx,s%ny,zsupd,s%gwlevel,s%wetz,ratio,0,par%g,dynpresupd) 00360 endif 00361 ! 00362 ! Compute groundwater head in cell centres 00363 ! In case of hydrostatic this is the same as the surface water head, 00364 ! or groundwater level. In order to ensure small errors do not explode 00365 ! minimum time averaging is done. 00366 ! In case of non-hydrostatic computations, the head is solved by the 00367 ! Poisson solver. In the case of a turbulent/modflow groundwater scheme 00368 ! the solver is forced using the estimate of the hydraulic conductivity 00369 ! from the previous time step. The estimate from the previous time step 00370 ! is also used to compute velocities (else no convervation of mass). The 00371 ! computation of the velocity also returns a new update of the hydraulic 00372 ! conductivity to be used in the next time step 00373 if (par%gwnonh == 0) then 00374 ! hydrostatic pressure assumption, with time delay required for 00375 ! stability 00376 ! factime = min(par%dt/0.2d0,1.d0) 00377 ! gwhead = (1.d0-factime)*gwhead + factime*gwheadtop 00378 s%gwhead = (s%gwhead + gwheadtop)/2 00379 else 00380 ! non-hydrostatic pressure, solve using Poisson solver. Kxupd is the 00381 ! updated hydraulic conductivity from the previous time step 00382 call gw_solver(par,s,gwheadtop,Kxupd,Kyupd,Kzupd,gwhu,gwhv,fracdt) 00383 endif 00384 ! 00385 ! 00386 ! Calculate groundwater velocities 00387 ! We use Kxupd input which is the same value as sent into the Poisson solver in 00388 ! the previous command. In the case of the non-hydrostatic solver, we must force 00389 ! the computation of the velocities to use the hydraulic conductivity estimate 00390 ! put into the Poisson solver. In practice this is only important is the hydraulic 00391 ! conductivity varies due to the turbulent/modflow scheme. 00392 ! 00393 ! Input Kx is the basic initial hydaulic conductivity 00394 ! Input Kxupd is the updated hydraulic conductivity from the previous time step 00395 ! used in the Poisson solver 00396 ! Output Kxupd is the new updated hydraulic conductivity for the next time step 00397 call gw_calculate_velocities(s,par,fracdt,s%gwu,s%gwv,s%gww,Kx,Ky,Kz,Kxupd,Kyupd,Kzupd) 00398 ! 00399 ! 00400 ! For next time step do some smoothing of Kx,Ky,Kx 00401 do j=1,s%ny+1 00402 do i=1,s%nx+1 00403 incl = min(i+3,s%nx+1)-max(1,i-3)+1 00404 Kx(i,j) = sum(Kxupd(max(1,i-3):min(i+3,s%nx+1),j))/incl 00405 Kz(i,j) = sum(Kzupd(max(1,i-3):min(i+3,s%nx+1),j))/incl 00406 Ky(i,j) = sum(Kyupd(max(1,i-3):min(i+3,s%nx+1),j))/incl 00407 enddo 00408 enddo 00409 Kxupd = Kx 00410 Kzupd = Kz 00411 Kyupd = Ky 00412 s%Kx = Kx 00413 s%Ky = Ky 00414 s%Kz = Kz 00415 Kx = par%kx 00416 Ky = par%ky 00417 Kz = par%kz 00418 s%Kzinf = Kzinf 00419 ! 00420 ! 00421 ! Calculate horizontal fluxes 00422 s%gwqx=s%gwu*gwhu 00423 s%gwqy=s%gwv*gwhv 00424 ! 00425 ! 00426 ! Stop cells from drying up 00427 if (s%ny>0) then 00428 do i=2,s%nx 00429 do j=2,s%ny 00430 if (s%gwlevel(i,j)<=s%gwbottom(i,j)+par%eps) then 00431 s%gwqx(i,j) = min(s%gwqx(i,j),0.d0) ! let no water out, only in 00432 s%gwqx(i-1,j) = max(s%gwqx(i-1,j),0.d0) ! let no water out, only in 00433 s%gwqy(i,j) = min(s%gwqy(i,j),0.d0) ! let no water out, only in 00434 s%gwqy(i,j-1) = max(s%gwqy(i,j-1),0.d0) ! let no water out, only in 00435 endif 00436 enddo 00437 enddo 00438 else 00439 do i=2,s%nx 00440 if (s%gwlevel(i,1)<=s%gwbottom(i,1)+par%eps) then 00441 s%gwqx(i,1) = min(s%gwqx(i,1),0.d0) ! let no water out, only in 00442 s%gwqx(i-1,1) = max(s%gwqx(i-1,1),0.d0) ! let no water out, only in 00443 endif 00444 enddo 00445 endif 00446 ! 00447 ! 00448 ! Force continuity, where gww is computed for hydrostatic computation and 00449 ! gww strictly enforced (no roundoff errors) for nonhydrostatic computation 00450 if (s%ny>2) then 00451 do j=2,s%ny 00452 do i=2,s%nx 00453 s%gww(i,j) = (-s%gwqx(i,j)*s%dnu(i,j)+s%gwqx(i-1,j)*s%dnu(i-1,j) & 00454 -s%gwqy(i,j)*s%dsv(i,j)+s%gwqy(i,j-1)*s%dsv(i,j-1))*s%dsdnzi(i,j) 00455 enddo 00456 enddo 00457 else 00458 if (s%ny>0) then 00459 do i=2,s%nx 00460 s%gww(i,2) = (-s%gwqx(i,2)+s%gwqx(i-1,2))/s%dsz(i,2) 00461 enddo 00462 s%gww(:,1) = s%gww(:,2) 00463 s%gww(:,3) = s%gww(:,2) 00464 else 00465 do i=2,s%nx 00466 s%gww(i,1) = (-s%gwqx(i,1)+s%gwqx(i-1,1))/s%dsz(i,1) 00467 enddo 00468 endif 00469 endif 00470 ! 00471 ! 00472 ! updating groundwater level is complicated: 00473 ! - in sections that are unconnected : dz/dt = -w/por == grad(qx,qy)/por 00474 ! - in sections that are connected: 00475 ! if w > 0 (exfiltration) : dz/dt = 0 (rise until gwlevel == zb, remainder is 00476 ! sent to zs) 00477 ! if w < 0 (infiltration) : dz/dt = 0 (unless available zs is less than w*dt, 00478 ! then drop water level on remainder) 00479 ! 00480 if (par%gwnonh==1) then 00481 ! non-hydrostatic computation 00482 do j=min(2,s%ny+1),max(s%ny,1) 00483 do i=2,s%nx 00484 ! exfiltration terms, same for connected and unconnected cells 00485 if (s%gww(i,j)>0.d0) then 00486 if (connected(i,j)) then 00487 dzr = max(s%zb(i,j)-s%gwlevel(i,j),0.d0) 00488 wdt = s%gww(i,j)*par%dt/par%por 00489 if (wdt>=dzr) then 00490 infilcon(i,j) = -s%gww(i,j)*(1.d0-dzr/wdt) 00491 else 00492 ! nothing 00493 endif 00494 endif 00495 ! connected infiltration terms, unconnected cells already computed 00496 ! infiltration at the start of the groundwater flow subroutine 00497 elseif (s%gww(i,j)<0.d0) then 00498 if (connected(i,j)) then 00499 dzr = max(s%zs(i,j)-s%zb(i,j)-par%eps,0.d0) 00500 wdt = -s%gww(i,j)*par%dt 00501 if (wdt>=dzr) then 00502 infilcon(i,j) = -s%gww(i,j)*(dzr/wdt) 00503 else 00504 infilcon(i,j) = -s%gww(i,j) 00505 endif 00506 endif 00507 else 00508 infilcon(i,j) = 0.d0 00509 endif 00510 enddo 00511 enddo 00512 else 00513 ! hydrostatic computation 00514 where(s%wetz==0) 00515 s%dinfil=par%dwetlayer/3.d0 ! Centroid of area influenced instantly by free surface level lies at dwetlayer/3 00516 elsewhere 00517 s%dinfil=min(s%dinfil,s%zb-s%gwlevel) 00518 s%dinfil=max(s%dinfil,par%dwetlayer/3.d0) 00519 endwhere 00520 ! wwvv since gw_calculate_hydrostatic_w returns a pointer: 00521 ! wwvv changed '=' into '=>' 00522 ! wwvv see also the declaration of infilcon 00523 infilcon => gw_calculate_hydrostatic_w(s,par,ratio) 00524 where (infilcon*par%dt>s%hh-par%eps) 00525 infilcon=(s%hh-par%eps)/par%dt 00526 endwhere 00527 s%dinfil=s%dinfil+infilcon*par%dt/par%por 00528 if (par%gwhorinfil==1) then 00529 call gw_horizontal_infil_exfil(s,par,infilhorgw,infilhorsw,Kx,dynpresupd) 00530 endif 00531 endif 00532 ! Compute new water level 00533 s%gwlevel = s%gwlevel + (s%gww+infilcon+infilhorgw)*par%dt/par%por 00534 ! 00535 ! In case of hydrostatic approach, try to fix the connected cells 00536 00537 ! 00538 ! Update total infiltration contribution 00539 s%infil = infiluncon+infilcon+infilhorsw 00540 ! 00541 ! Model boundaries 00542 ! Robert: check if all these are actually needed 00543 #ifdef USEMPI 00544 call xmpi_shift_ee(s%gwlevel) 00545 call xmpi_shift_ee(s%gwhead) 00546 call xmpi_shift_ee(s%gwcurv) 00547 #endif 00548 if (xmpi_istop) then 00549 s%gwlevel(1,:) = s%gwlevel(2,:) 00550 s%gwhead(1,:) = s%gwhead(2,:) 00551 s%gwcurv(1,:) = s%gwcurv(2,:) 00552 s%gwbottom(1,:) = s%gwbottom(2,:) 00553 s%gww(1,:) = s%gww(2,:) 00554 endif 00555 if (xmpi_isbot) then 00556 s%gwlevel(s%nx+1,:) = s%gwlevel(s%nx,:) 00557 s%gwhead(s%nx+1,:) = s%gwhead(s%nx,:) 00558 s%gwcurv(s%nx+1,:) = s%gwcurv(s%nx,:) 00559 s%gwbottom(s%nx+1,:) = s%gwbottom(s%nx,:) 00560 s%gww(s%nx+1,:) = s%gww(s%nx,:) 00561 endif 00562 if (xmpi_isleft .and. s%ny>0) then 00563 s%gwlevel(:,1) = s%gwlevel(:,2) 00564 s%gwhead(:,1) = s%gwhead(:,2) 00565 s%gwcurv(:,1) = s%gwcurv(:,2) 00566 s%gwbottom(:,1) = s%gwbottom(:,2) 00567 s%gww(:,1) = s%gww(:,2) 00568 endif 00569 if (xmpi_isright .and. s%ny>0) then 00570 s%gwlevel(:,s%ny+1) = s%gwlevel(:,s%ny) 00571 s%gwhead(:,s%ny+1) = s%gwhead(:,s%ny) 00572 s%gwcurv(:,s%ny+1) = s%gwcurv(:,s%ny) 00573 s%gwbottom(:,s%ny+1) = s%gwbottom(:,s%ny) 00574 s%gww(:,s%ny+1) = s%gww(:,s%ny) 00575 endif 00576 ! 00577 ! output pressure head at bottom 00578 ! 00579 s%gwheadb = gwCalculateHeadBottom(s%gwcurv,gwheadtop,s%gwheight,s%nx,s%ny,par%gwheadmodel) 00580 s%gwheadb(1,:) = s%gwheadb(2,:) 00581 s%gwheadb(s%nx+1,:) = s%gwheadb(s%nx,:) 00582 end subroutine gwflow 00583 00584 subroutine gw_unconnected_infil(s,par,Kzinf,connected,fracdt,infil) 00585 use params, only: parameters 00586 use xmpi_module 00587 use spaceparamsdef 00588 use paramsconst 00589 00590 IMPLICIT NONE 00591 00592 type(parameters),intent(in) :: par 00593 type(spacepars) :: s 00594 real*8,dimension(s%nx+1,s%ny+1),intent(inout):: Kzinf 00595 logical,dimension(s%nx+1,s%ny+1),intent(in) :: connected 00596 real*8,dimension(s%nx+1,s%ny+1),intent(out) :: fracdt,infil 00597 ! local 00598 integer :: i,j 00599 logical :: turbapprox 00600 real*8,dimension(s%nx+1,s%ny+1) :: Kzupd 00601 00602 ! initialise 00603 fracdt = 0.d0 00604 infil = 0.d0 00605 ! 00606 ! infiltration depends on turbulent or laminar method 00607 select case (par%gwscheme) 00608 case (GWSCHEME_LAMINAR) 00609 turbapprox = .false. 00610 case (GWSCHEME_TURBULENT) 00611 turbapprox = .true. 00612 end select 00613 ! 00614 ! 00615 ! Start Kzupd by seeding with Kz. 00616 Kzupd = Kzinf 00617 do j=1,s%ny+1 00618 do i=1,s%nx+1 00619 if (.not. connected(i,j)) then 00620 if (s%wetz(i,j)==1 .and. s%gwlevel(i,j)<s%zb(i,j)) then 00621 call gw_calc_local_infil(par,s%zb(i,j),s%hh(i,j),s%gwlevel(i,j),Kzinf(i,j),& 00622 s%dinfil(i,j),s%D50top(i,j),turbapprox, & 00623 fracdt(i,j),infil(i,j),Kzupd(i,j)) 00624 endif 00625 endif 00626 enddo 00627 enddo 00628 ! 00629 ! Update infiltration layer thickness 00630 where (infil>0.d0) 00631 s%dinfil = s%dinfil+par%dt*infil/par%por 00632 elsewhere 00633 ! s%dinfil = max(s%dinfil-par%dt*par%kz/par%por,0.d0) 00634 s%dinfil = 0.d0 00635 Kzupd = par%kz 00636 endwhere 00637 ! Ensure nothing went wrong here with fracdt 00638 fracdt = min(max(fracdt,0.d0),1.d0) 00639 ! Update Kz 00640 Kzinf = Kzupd 00641 end subroutine gw_unconnected_infil 00642 00643 subroutine gw_calculate_interfaceheight(s,gwhu,gwhv,initial) 00644 use xmpi_module 00645 use spaceparams 00646 00647 IMPLICIT NONE 00648 00649 type(spacepars),intent(in) :: s 00650 logical,intent(in) :: initial 00651 real*8,dimension(s%nx+1,s%ny+1),intent(out) :: gwhu,gwhv 00652 ! local 00653 integer :: i,j 00654 00655 ! HU 00656 if (initial) then 00657 do j=1,s%ny+1 00658 do i=1,s%nx 00659 if (s%gwhead(i+1,j)>s%gwhead(i,j)) then 00660 gwhu(i,j) = s%gwheight(i+1,j) 00661 elseif (s%gwhead(i+1,j)<s%gwhead(i,j)) then 00662 gwhu(i,j) = s%gwheight(i,j) 00663 else 00664 gwhu(i,j) = 0.5d0*(s%gwheight(i+1,j)+s%gwheight(i,j)) 00665 endif 00666 enddo 00667 enddo 00668 else 00669 do j=1,s%ny+1 00670 do i=1,s%nx 00671 if (s%gwu(i,j)>0.d0) then 00672 gwhu(i,j) = min(s%gwheight(i,j),max(s%zb(i+1,j)-s%gwbottom(i+1,j),0.d0)) 00673 elseif (s%gwu(i,j)<0.d0) then 00674 gwhu(i,j) = min(s%gwheight(i+1,j),max(s%zb(i,j)-s%gwbottom(i,j),0.d0)) 00675 else 00676 gwhu(i,j) = 0.5d0*(s%gwheight(i+1,j)+s%gwheight(i,j)) 00677 endif 00678 enddo 00679 enddo 00680 endif 00681 gwhu = max(gwhu,0.d0) ! in case of drying cells 00682 ! boundaries 00683 if (xmpi_isbot) then 00684 gwhu(s%nx+1,:) = gwhu(s%nx,:) 00685 endif 00686 00687 ! HV 00688 if (s%ny>0) then 00689 if (initial) then 00690 do j=1,s%ny 00691 do i=1,s%nx+1 00692 if (s%gwhead(i,j+1)>s%gwhead(i,j)) then 00693 gwhv(i,j) = s%gwheight(i,j+1) 00694 elseif (s%gwhead(i,j+1)<s%gwhead(i,j)) then 00695 gwhv(i,j) = s%gwheight(i,j) 00696 else 00697 gwhv(i,j) = 0.5d0*(s%gwheight(i,j+1)+s%gwheight(i,j)) 00698 endif 00699 enddo 00700 enddo 00701 else 00702 do j=1,s%ny 00703 do i=1,s%nx+1 00704 if (s%gwv(i,j)>0.d0) then 00705 gwhv(i,j) = min(s%gwheight(i,j),max(s%zb(i,j+1)-s%gwbottom(i,j+1),0.d0)) 00706 elseif (s%gwv(i,j)<0.d0) then 00707 gwhv(i,j) = min(s%gwheight(i,j+1),max(s%zb(i,j)-s%gwbottom(i,j),0.d0)) 00708 else 00709 gwhv(i,j) = 0.5d0*(s%gwheight(i,j+1)+s%gwheight(i,j)) 00710 endif 00711 enddo 00712 enddo 00713 endif 00714 gwhv = max(gwhv,0.d0) ! in case of drying cells 00715 ! Boundaries 00716 if (xmpi_isright) then 00717 gwhv(:,s%ny+1) = gwhv(:,s%ny) 00718 endif 00719 else ! ny==0 00720 gwhv(:,1) = s%gwheight(:,1) 00721 endif 00722 end subroutine gw_calculate_interfaceheight 00723 00724 pure function gw_calculate_smoothwetlayer(s,dwetlayer,px) result(r) 00725 use spaceparams 00726 00727 IMPLICIT NONE 00728 00729 type(spacepars),intent(in) :: s 00730 real*8,intent(in) :: dwetlayer,px 00731 real*8,dimension(s%nx+1,s%ny+1) :: r 00732 00733 ! Free surface head and ratio free surface to groundwater head to be used 00734 r = (s%zb-s%gwlevel)/dwetlayer 00735 ! This ratio should lie between 0 and 1 00736 r = min(max(r,0.d0),1.d0) 00737 ! Convert this to a smoother function at the edges 00738 r = 1.d0-(cos(r*px)+1)/2 00739 end function gw_calculate_smoothwetlayer 00740 00741 pure function gwCalculateHeadTop(nx,ny,zs,gwlevel,wetz,r,nonh,g,pres) result(gwhead) 00742 IMPLICIT NONE 00743 ! input 00744 integer,intent(in) :: nx,ny,nonh 00745 real*8,intent(in) :: g 00746 real*8,dimension(nx+1,ny+1),intent(in) :: zs,gwlevel,r,pres 00747 integer,dimension(nx+1,ny+1),intent(in) :: wetz 00748 ! result 00749 real*8,dimension(nx+1,ny+1) :: gwhead 00750 00751 ! If this point is wet (wetz==1), and 0<r<=1, additional head is added to 00752 ! the groundwater, described as ratio "r" times the difference in head between 00753 ! the surface water and the groundwater. 00754 if(nonh==1) then 00755 gwhead = gwlevel + wetz*(zs-gwlevel+pres/g)*(1.d0-r) 00756 else 00757 gwhead = gwlevel + wetz*(zs-gwlevel)*(1.d0-r) 00758 endif 00759 end function gwCalculateHeadTop 00760 00761 pure subroutine gwCalculateHeadGradient(nx,ny,gwhead,dsu,dnv,dheaddx,dheaddy) 00762 IMPLICIT NONE 00763 ! input/output 00764 integer,intent(in) :: nx,ny 00765 real*8,dimension(nx+1,ny+1),intent(in) :: gwhead,dsu,dnv 00766 real*8,dimension(nx+1,ny+1),intent(out) :: dheaddx,dheaddy 00767 ! internal 00768 integer :: i,j 00769 00770 dheaddx=0.d0 00771 dheaddy=0.d0 00772 00773 do j=1,ny+1 00774 do i=1,nx 00775 dheaddx(i,j)=(gwhead(i+1,j)-gwhead(i,j))/dsu(i,j) 00776 end do 00777 end do 00778 00779 do j=1,ny 00780 do i=1,nx+1 00781 dheaddy(i,j)=(gwhead(i,j+1)-gwhead(i,j))/dnv(i,j) 00782 end do 00783 end do 00784 end subroutine gwCalculateHeadGradient 00785 00786 subroutine gw_solver(par,s,hbc,Kx,Ky,Kz,hu,hv,fracdt) 00787 use params, only: parameters 00788 use xmpi_module 00789 use spaceparamsdef 00790 use solver_module, only: solver_tridiag,solver_sip 00791 use paramsconst 00792 00793 IMPLICIT NONE 00794 00795 type(parameters),intent(in) :: par 00796 type(spacepars) :: s 00797 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: hbc 00798 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: Kx,Ky,Kz 00799 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: hu,hv 00800 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: fracdt 00801 ! internal 00802 real*8,dimension(:,:,:),allocatable,save :: A,work 00803 real*8,dimension(:,:),allocatable,save :: rhs,x,res 00804 integer :: i,j,n,m,k,it 00805 integer :: imin,imax,jmin,jmax 00806 real*8 :: dxum,dxup,dxpc 00807 real*8 :: dyvm,dyvp 00808 real*8 :: dyum,dyup,dxvm,dxvp,dA 00809 real*8 :: hcmx,hcc,hcpx 00810 real*8 :: hcmy,hcpy 00811 real*8 :: twothird 00812 real*8 :: hum,hup,hvm,hvp 00813 00814 twothird = 2.d0/3 00815 00816 i=0 00817 k=0 00818 00819 ! use tridiagonal solver if ny<=2 00820 if (s%ny<3) then 00821 if (par%ny==0) then 00822 j=1 00823 else 00824 j=2 00825 endif 00826 n = s%nx-1 00827 ! allocate Matrix solver coefficients 00828 if(.not.allocated(A)) then 00829 allocate(A(5,n,1)) 00830 allocate(work(1,n,1)) 00831 work = 0.d0 00832 allocate(rhs(n,1)) 00833 allocate(x(n,1)) 00834 endif 00835 ! build Matrix solver coefficients 00836 select case (par%gwheadmodel) 00837 case (GWHEADMODEL_PARABOLIC) 00838 do k=1,n 00839 i = k+1 00840 if (s%zb(i,1)>s%gwbottom(i,1)+par%eps) then 00841 ! abbreviations 00842 dxum = s%dsu(i-1,j) 00843 dxup = s%dsu(i,j) 00844 dxpc = s%dsz(i,j) 00845 hcmx = s%gwheight(i-1,j) 00846 hcc = s%gwheight(i,j) 00847 hcpx = s%gwheight(i+1,j) 00848 hum = hu(i-1,j) 00849 hup = hu(i,j) 00850 if (k==-1) then 00851 ! A(2,k,1) = 0.d0 ! dummy 00852 ! A(1,k,1) = +Kx(i,j)*hup/dxup*twothird*hcc**2 +Kz(i-1,j)*2*dxpc*hcc*(1.d0-fracdt(i,j)) 00853 ! A(3,k,1) = -Kx(i,j)*hup/dxup*twothird*hcpx**2 00854 ! rhs(k,1) = - Kx(i,j)*hbc(i+1,j)*hup/dxup + Kx(i+1,j)*hbc(i,j)*hup/dxup 00855 elseif (k==-n) then 00856 ! A(2,k,1) = -Kx(i,j)*hum/dxum*twothird*hcmx**2 00857 ! A(1,k,1) = Kx(i,j)*hum/dxum*twothird*hcc**2 +Kz(i,j)*2*dxpc*hcc*(1.d0-fracdt(i,j)) 00858 ! A(3,k,1) = 0.d0 ! dummy 00859 ! rhs(k,1) = Kx(i,j)*hbc(i,j)*hum/dxum - Kx(i,j)*hbc(i-1,j)*hum/dxum 00860 else 00861 A(2,k,1) = -Kx(i-1,j)*hum/dxum*twothird*hcmx**2 00862 A(1,k,1) = Kx(i-1,j)*hum/dxum*twothird*hcc**2 +Kx(i,j)*hup/dxup*twothird*hcc**2 + & 00863 Kz(i,j)*2*dxpc*hcc*(1.d0-fracdt(i,j)) 00864 A(3,k,1) = -Kx(i,j)*hup/dxup*twothird*hcpx**2 00865 rhs(k,1) = Kx(i-1,j)*hum*(hbc(i ,j)-hbc(i-1,j))/dxum & 00866 - Kx(i ,j)*hup*(hbc(i+1,j)-hbc(i ,j))/dxup 00867 endif 00868 else 00869 A(:,k,1) = 0.d0 00870 rhs(k,1) = 0.d0 00871 endif 00872 enddo 00873 case (GWHEADMODEL_EXPONENTIAL) 00874 do k=1,n 00875 i = k+1 00876 ! abbreviations 00877 dxum = s%dsu(i-1,j) 00878 dxup = s%dsu(i,j) 00879 dxpc = s%dsz(i,j) 00880 hcmx = s%gwheight(i-1,j) 00881 hcc = s%gwheight(i,j) 00882 hcpx = s%gwheight(i+1,j) 00883 hum = hu(i-1,j) 00884 hup = hu(i,j) 00885 if (k==-1) then 00886 A(2,k,1) = 0.d0 ! dummy 00887 A(1,k,1) = +Kx(i+1,j)*hup/dxup*twothird*hcc**2 +Kz(i,j)*2*dxpc*hcc*(1.d0-fracdt(i,j)) 00888 A(3,k,1) = -Kx(i+1,j)*hup/dxup*twothird*hcpx**2 00889 rhs(k,1) = - Kx(i+1,j)*hbc(i+1,j)*hup/dxup + Kx(i+1,j)*hbc(i,j)*hup/dxup 00890 elseif (k==-n) then 00891 A(2,k,1) = -Kx(i,j)*hum/dxum*twothird*hcmx**2 00892 A(1,k,1) = Kx(i,j)*hum/dxum*twothird*hcc**2 +Kz(i,j)*2*dxpc*hcc*(1.d0-fracdt(i,j)) 00893 A(3,k,1) = 0.d0 ! dummy 00894 rhs(k,1) = Kx(i,j)*hbc(i,j)*hum/dxum - Kx(i,j)*hbc(i-1,j)*hum/dxum 00895 else 00896 A(2,k,1) = -Kx(i-1,j)*hum/dxum*(cosh(hcmx)-1/hcmx*sinh(hcmx)) 00897 A(1,k,1) = Kx(i-1,j)*hum/dxum*(cosh(hcc)-1/hcc*sinh(hcc)) + & 00898 Kx(i,j)*hup/dxup*(cosh(hcc)-1/hcc*sinh(hcc)) + & 00899 Kz(i,j)*dxpc*sinh(hcc)*(1.d0-fracdt(i,j)) 00900 A(3,k,1) = -Kx(i,j)*hup/dxup*(cosh(hcpx)-1/hcpx*sinh(hcpx)) 00901 00902 rhs(k,1) = Kx(i-1,j)*hbc(i,j)*hum/dxum - Kx(i-1,j)*hbc(i-1,j)*hum/dxum & 00903 - Kx(i,j)*hbc(i+1,j)*hup/dxup + Kx(i,j)*hbc(i,j)*hup/dxup 00904 endif 00905 enddo 00906 end select 00907 call solver_tridiag(A,rhs,x,work(1,:,1),n-1,0,fixshallow=.true.) 00908 ! return to global variable 00909 s%gwcurv(2:s%nx,j) = x(:,1) 00910 select case (par%gwheadmodel) 00911 case (GWHEADMODEL_PARABOLIC) 00912 s%gwhead(2:s%nx,j) = hbc(2:s%nx,j) - twothird*x(:,1)*s%gwheight(2:s%nx,j)**2 00913 case (GWHEADMODEL_EXPONENTIAL) 00914 s%gwhead(2:s%nx,j) = hbc(2:s%nx,j)+ & 00915 x(:,1)/s%gwheight(2:s%nx,j)*sinh(s%gwheight(2:s%nx,j)) & 00916 -x(:,1)*cosh(s%gwheight(2:s%nx,j)) 00917 end select 00918 s%gwhead(1,j) = s%gwhead(2,j) 00919 s%gwhead(s%nx+1,j) = s%gwhead(s%nx,j) 00920 ! spread left and right 00921 if (s%ny==2) then 00922 s%gwhead(:,1) = s%gwhead(:,2) 00923 s%gwhead(:,3) = s%gwhead(:,2) 00924 endif 00925 else 00926 if (par%gwfastsolve==1) then 00927 ! Use tridiagonal solver per row with explicit longshore velocities 00928 n = s%nx+1 00929 m = s%ny+1 00930 ! allocate Matrix solver coefficients 00931 if(.not.allocated(A)) then 00932 allocate(A(5,n,1)) 00933 allocate(work(1,n,1)) 00934 work = 0.d0 00935 allocate(rhs(n,1)) 00936 allocate(x(n,1)) 00937 endif 00938 do j=2,m-1 00939 ! build Matrix solver coefficients 00940 select case (par%gwheadmodel) 00941 case (GWHEADMODEL_PARABOLIC) 00942 do i=1,n 00943 imin = max(i-1,1) 00944 imax = min(i+1,s%nx+1) 00945 ! abbreviations 00946 dxum = s%dsu(imin,j) 00947 dxup = s%dsu(i,j) 00948 dxpc = s%dsz(i,j) 00949 dA = s%dsdnzi(i,j) 00950 dyum = s%dnu(imin,j) 00951 dyup = s%dnu(i,j) 00952 dyvm = s%dnv(i,j-1) 00953 dyvp = s%dnv(i,j) 00954 dxvm = s%dsv(i,j-1) 00955 dxvp = s%dsv(i,j) 00956 hcmx = s%gwheight(imin,j) 00957 hcc = s%gwheight(i,j) 00958 hcpx = s%gwheight(imax,j) 00959 hum = hu(imin,j) 00960 hup = hu(i,j) 00961 hvm = hv(i,j-1) 00962 hvp = hv(i,j) 00963 ! wwvv problem what is the value of k here 00964 ! wwvv to avoid warnings about uninitialized value for k 00965 ! wwvv I initialize k at the start of the program 00966 A(2,k,1) = -Kx(imin,j)*hum*dyum/dxum*twothird*hcmx**2 00967 A(1,k,1) = Kx(imin,j)*hum*dyum/dxum*twothird*hcc**2 + & 00968 Kx(i,j)*hup*dyup/dxup*twothird*hcc**2 + & 00969 Kz(i,j)*2*dA*hcc*(1.d0-fracdt(i,j)) 00970 A(3,k,1) = -Kx(i,j)*hup*dyup/dxup*twothird*hcpx**2 00971 rhs(k,1) = Kx(imin,j)*hum*dyum*(hbc(i ,j)-hbc(imin,j))/dxum & 00972 - Kx(i ,j)*hup*dyup*(hbc(imax,j)-hbc(i ,j))/dxup & 00973 + hvm*dxvm*s%gwv(i,j-1) - hvp*dxvp*s%gwv(i,j) 00974 00975 enddo 00976 case (GWHEADMODEL_EXPONENTIAL) 00977 ! wwvv problem detected by compiler: should this not be an i-loop? 00978 ! wwvv compiler suspects uninitialized value for i, I initialize i 00979 ! wwvv to 0 at the start of the subroutine 00980 do k=1,n 00981 imin = max(i-1,1) 00982 imax = min(i+1,s%nx+1) 00983 ! abbreviations 00984 dxum = s%dsu(imin,j) 00985 dxup = s%dsu(i,j) 00986 dxpc = s%dsz(i,j) 00987 dA = s%dsdnzi(i,j) 00988 dyum = s%dnu(imin,j) 00989 dyup = s%dnu(i,j) 00990 dyvm = s%dnv(i,j-1) 00991 dyvp = s%dnv(i,j) 00992 dxvm = s%dsv(i,j-1) 00993 dxvp = s%dsv(i,j) 00994 hcmx = s%gwheight(imin,j) 00995 hcc = s%gwheight(i,j) 00996 hcpx = s%gwheight(imax,j) 00997 hum = hu(imin,j) 00998 hup = hu(i,j) 00999 hvm = hv(i,j-1) 01000 hvp = hv(i,j) 01001 A(2,k,1) = -Kx(imin,j)*hum*dyum/dxum*(cosh(hcmx)-1/hcmx*sinh(hcmx)) 01002 A(1,k,1) = Kx(imin,j)*hum*dyum/dxum*(cosh(hcc)-1/hcc*sinh(hcc)) + & 01003 Kx(i,j)*hup*dyup/dxup*(cosh(hcc)-1/hcc*sinh(hcc)) + & 01004 Kz(i,j)*dA*sinh(hcc)*(1.d0-fracdt(i,j)) 01005 A(3,k,1) = -Kx(i,j)*hup*dyup/dxup*(cosh(hcpx)-1/hcpx*sinh(hcpx)) 01006 01007 rhs(k,1) = Kx(imin,j)*hbc(i,j)*hum*dyum/dxum - Kx(i,j)*hbc(i-1,j)*hum*dyum/dxum & 01008 - Kx(i,j)*hbc(imax,j)*hup*dyup/dxup + Kx(i,j)*hbc(i,j)*hup*dyup/dxup 01009 enddo 01010 end select 01011 call solver_tridiag(A,rhs,x,work(1,:,1),n-1,0,fixshallow=.true.) 01012 s%gwcurv(:,j) = x(:,1) 01013 select case (par%gwheadmodel) 01014 case (GWHEADMODEL_PARABOLIC) 01015 s%gwhead(2:s%nx,j) = hbc(2:s%nx,j) - twothird*x(:,1)*s%gwheight(2:s%nx,j)**2 01016 case (GWHEADMODEL_EXPONENTIAL) 01017 s%gwhead(2:s%nx,j) = hbc(2:s%nx,j)+ & 01018 x(:,1)/s%gwheight(2:s%nx,j)*sinh(s%gwheight(2:s%nx,j)) & 01019 -x(:,1)*cosh(s%gwheight(2:s%nx,j)) 01020 end select 01021 enddo 01022 ! spread left and right 01023 s%gwhead(:,1) = s%gwhead(:,2) 01024 s%gwhead(:,s%ny+1) = s%gwhead(:,s%ny) 01025 else 01026 n = s%nx+1 01027 m = s%ny+1 01028 ! allocate Matrix solver coefficients 01029 if(.not.allocated(A)) then 01030 allocate(A(5,n,m)) 01031 allocate(work(5,n,m)) 01032 work = 0.d0 01033 allocate(rhs(n,m)) 01034 allocate(x(n,m)) 01035 x = 0.d0 01036 allocate(res(n,m)) 01037 res = 0.d0 01038 endif 01039 ! build Matrix solver coefficients 01040 select case (par%gwheadmodel) 01041 case (GWHEADMODEL_PARABOLIC) 01042 do j=1,m 01043 do i=1,n 01044 imin = max(i-1,1) 01045 imax = min(i+1,s%nx+1) 01046 jmin = max(j-1,1) 01047 jmax = min(j+1,s%ny+1) 01048 ! abbreviations (need to include dx and dy for mass conservation 01049 ! on curvilinear grid 01050 dxum = s%dsu(imin,j) 01051 dxup = s%dsu(i,j) 01052 dxpc = s%dsz(i,j) 01053 dA = s%dsdnzi(i,j) 01054 dyum = s%dnu(imin,j) 01055 dyup = s%dnu(i,j) 01056 dyvm = s%dnv(i,jmin) 01057 dyvp = s%dnv(i,j) 01058 dxvm = s%dsv(i,jmin) 01059 dxvp = s%dsv(i,j) 01060 hcmx = s%gwheight(imin,j) 01061 hcpx = s%gwheight(imax,j) 01062 hcc = s%gwheight(i,j) 01063 hcmy = s%gwheight(i,jmin) 01064 hcpy = s%gwheight(i,jmax) 01065 hum = hu(imin,j) 01066 hup = hu(i,j) 01067 hvm = hv(i,jmin) 01068 hvp = hv(i,j) 01069 ! Matrix A 01070 ! main diagonal 01071 A(1,i,j) = Kx(imin,j)*hum*dyum/dxum*twothird*hcc**2+Kx(i,j)*hup*dyup/dxup*twothird*hcc**2 + & 01072 Ky(i,jmin)*hvm*dxvm/dyvm*twothird*hcc**2+Ky(i,j)*hvp*dxvp/dyvp*twothird*hcc**2 + & 01073 Kz(i,j)*2*dA*hcc*(1.d0-fracdt(i,j)) 01074 ! down x 01075 A(2,i,j) = -Kx(imin,j)*hum*dyum/dxum*twothird*hcmx**2 01076 ! up x 01077 A(3,i,j) = -Kx(i,j)*hup*dyup/dxup*twothird*hcpx**2 01078 ! down y 01079 A(4,i,j) = -Ky(i,jmin)*hvm*dxvm/dyvm*twothird*hcmy**2 01080 ! up y 01081 A(5,i,j) = -Ky(i,j)*hvp*dxvp/dyvp*twothird*hcpy**2 01082 ! 01083 ! RHS 01084 rhs(i,j) = Kx(imin,j)*dyum*hum*(hbc(i ,j)-hbc(imin,j))/dxum & 01085 - Kx(i ,j)*dyup*hup*(hbc(imax,j)-hbc(i ,j))/dxup & 01086 + Ky(i,jmin)*dxvm*hvm*(hbc(i,j )-hbc(i,jmin))/dyvm & 01087 - Ky(i,j )*dxvp*hvp*(hbc(i,jmax)-hbc(i,j ))/dyvp 01088 enddo 01089 enddo 01090 res = 0.d0 01091 call solver_sip(A,rhs,x,res,work,it,s%nx,s%ny) 01092 s%gwcurv(:,:) = x 01093 s%gwcurv(1,:) = s%gwcurv(2,:) 01094 s%gwcurv(s%nx+1,:) = s%gwcurv(s%nx,:) 01095 s%gwcurv(:,1) = s%gwcurv(:,2) 01096 s%gwcurv(:,s%ny+1) = s%gwcurv(:,s%ny) 01097 end select 01098 endif ! speedup 2D 01099 01100 select case (par%gwheadmodel) 01101 case (GWHEADMODEL_PARABOLIC) 01102 s%gwhead = hbc - twothird*s%gwcurv*s%gwheight**2 01103 case (GWHEADMODEL_EXPONENTIAL) 01104 s%gwhead = hbc + s%gwcurv/s%gwheight*sinh(s%gwheight) & 01105 - s%gwcurv*cosh(s%gwheight) 01106 end select 01107 endif 01108 end subroutine gw_solver 01109 01110 subroutine gw_horizontal_infil_exfil(s,par,infilhorgw,infilhorsw,Kx,dynpres) 01111 ! compute horizontal part of infiltration/exfiltration 01112 use params, only: parameters 01113 use xmpi_module 01114 use spaceparamsdef 01115 01116 type(parameters),intent(in) :: par 01117 type(spacepars) :: s 01118 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: Kx,dynpres 01119 real*8,dimension(s%nx+1,s%ny+1),intent(out) :: infilhorgw,infilhorsw 01120 ! internal 01121 integer :: i,j 01122 real*8 :: dz,dx,dy,vel,hsurf,hgw,dhead 01123 real*8,parameter :: visc = 1.0d-6 01124 real*8 :: vcr,scaleareab,scaleareas 01125 real*8 :: dum1,dum2,dum3 01126 01127 infilhorgw = 0.d0 01128 infilhorsw = 0.d0 01129 vcr = par%gwReturb/par%D50(1)*par%por*visc 01130 ! 01131 ! compute horizontal flow velocity between bed and surface water VERTICAL interfaces 01132 ! 01133 ! loop over all u-interfaces 01134 do i=1,s%nx 01135 do j=1,s%ny+1 01136 dz = s%zb(i+1,j)-s%zb(i,j) 01137 if (dz>0.d0 .and. s%wetz(i,j)==1) then ! jump up 01138 hsurf = s%zs(i,j)+dynpres(i,j)/par%g 01139 hgw = s%gwhead(i+1,j) 01140 dx = s%dsz(i+1,j) 01141 dhead = hsurf-hgw ! positive gradient if hsurf>hgw 01142 ! negative velocitity if infiltration 01143 call gw_calculate_velocities_local(vel,dum1,dhead/dx,Kx(i+1,j),0.d0,par%gwscheme, & 01144 .false.,.false.,par%gwheadmodel,dum2,dum3,vcr) 01145 vel = -vel ! convert negative down velocity into positive if infiltration, 01146 ! negative if exfiltration 01147 ! horizontal infiltration occurs along small side section, so scale the horizontal 01148 ! infiltration area by the vertical area of the cell used in the mass balance equation 01149 ! and infiltration velocities 01150 ! 01151 ! Scale area for bed cell and surface water cell 01152 scaleareab = dz/s%dsz(i+1,j) 01153 scaleareas = dz/s%dsz(i,j) 01154 ! Limit flow to that available in space in ground and 01155 ! in surface water volume 01156 if (vel>0.d0) then ! infiltration 01157 vel = min(vel, & 01158 s%hh(i,j)*s%dsz(i,j)/par%dt/dz, & 01159 max(s%zb(i+1,j)-s%gwlevel(i+1,j),0.d0)*par%por*s%dsz(i+1,j)/par%dt/dz) 01160 else 01161 vel = max(vel, & 01162 -(s%gwlevel(i+1,j)-s%gwbottom(i+1,j))*s%dsz(i+1,j)*par%por/par%dt/dz) 01163 endif 01164 infilhorgw(i+1,j) = infilhorgw(i+1,j)+vel*scaleareab 01165 infilhorsw(i,j) = infilhorsw(i,j)+vel*scaleareas 01166 elseif (dz<0.d0 .and. s%wetz(i+1,j)==1) then ! jump down 01167 dz = -dz 01168 hsurf = s%zs(i+1,j)+dynpres(i+1,j)/par%g 01169 hgw = s%gwhead(i,j) 01170 dx = s%dsz(i,j) 01171 dhead = hsurf-hgw ! positive gradient if hsurf>hgw 01172 ! negative velocitity if infiltration 01173 call gw_calculate_velocities_local(vel,dum1,dhead/dx,Kx(i,j),0.d0,par%gwscheme, & 01174 .false.,.false.,par%gwheadmodel,dum2,dum3,vcr) 01175 vel = -vel ! convert negative down velocity into positive if infiltration, 01176 ! negative if exfiltration 01177 ! horizontal infiltration occurs along small side section, so scale the horizontal 01178 ! infiltration area by the vertical area of the cell used in the mass balance equation 01179 ! and infiltration velocities 01180 ! 01181 ! Scale area for bed cell and surface water cell 01182 scaleareab = dz/s%dsz(i,j) 01183 scaleareas = dz/s%dsz(i+1,j) 01184 ! Limit flow to that available in space in ground and 01185 ! in surface water volume 01186 if (vel>0.d0) then ! infiltration 01187 vel = min(vel, & 01188 s%hh(i+1,j)*s%dsz(i+1,j)/par%dt/dz, & 01189 max(s%zb(i,j)-s%gwlevel(i,j),0.d0)*par%por*s%dsz(i,j)/par%dt/dz) 01190 else 01191 vel = max(vel, & 01192 -(s%gwlevel(i,j)-s%gwbottom(i,j))*s%dsz(i,j)*par%por/par%dt/dz) 01193 endif 01194 infilhorgw(i,j) = infilhorgw(i,j)+vel*scaleareab 01195 infilhorsw(i+1,j) = infilhorsw(i+1,j)+vel*scaleareas 01196 else ! equal bed height 01197 ! nothing 01198 endif 01199 enddo 01200 enddo 01201 ! 01202 ! 01203 ! loop over all v-interfaces 01204 do i=1,s%nx+1 01205 do j=1,s%ny 01206 dz = s%zb(i,j+1)-s%zb(i,j) 01207 if (dz>0.d0 .and. s%wetz(i,j)==1) then ! jump up 01208 hsurf = s%zs(i,j)+dynpres(i,j)/par%g 01209 hgw = s%gwhead(i,j+1) 01210 dy = s%dnz(i,j+1) 01211 dhead = hsurf-hgw ! positive gradient if hsurf>hgw 01212 ! negative velocitity if infiltration 01213 call gw_calculate_velocities_local(vel,dum1,dhead/dy,Kx(i,j+1),0.d0,par%gwscheme, & 01214 .false.,.false.,par%gwheadmodel,dum2,dum3,vcr) 01215 vel = -vel ! convert negative down velocity into positive if infiltration, 01216 ! negative if exfiltration 01217 ! horizontal infiltration occurs along small side section, so scale the horizontal 01218 ! infiltration area by the vertical area of the cell used in the mass balance equation 01219 ! and infiltration velocities 01220 ! 01221 ! Scale area for bed cell and surface water cell 01222 scaleareab = dz/s%dnz(i,j+1) 01223 scaleareas = dz/s%dnz(i,j) 01224 ! Limit flow to that available in space in ground and 01225 ! in surface water volume 01226 if (vel>0.d0) then ! infiltration 01227 vel = min(vel, & 01228 s%hh(i,j)*s%dsz(i,j)/par%dt/dz, & 01229 max(s%zb(i,j+1)-s%gwlevel(i,j+1),0.d0)*par%por*s%dnz(i,j+1)/par%dt/dz) 01230 else 01231 vel = max(vel, & 01232 -(s%gwlevel(i,j+1)-s%gwbottom(i,j+1))*s%dnz(i,j+1)*par%por/par%dt/dz) 01233 endif 01234 infilhorgw(i,j+1) = infilhorgw(i,j+1)+vel*scaleareab 01235 infilhorsw(i,j) = infilhorsw(i,j)+vel*scaleareas 01236 elseif (dz<0.d0 .and. s%wetz(i,j+1)==1) then ! jump down 01237 dz = -dz 01238 hsurf = s%zs(i,j+1)+dynpres(i,j+1)/par%g 01239 hgw = s%gwhead(i,j) 01240 dy = s%dnz(i,j) 01241 dhead = hsurf-hgw ! positive gradient if hsurf>hgw 01242 ! negative velocitity if infiltration 01243 call gw_calculate_velocities_local(vel,dum1,dhead/dy,Kx(i,j),0.d0,par%gwscheme, & 01244 .false.,.false.,par%gwheadmodel,dum2,dum3,vcr) 01245 vel = -vel ! convert negative down velocity into positive if infiltration, 01246 ! negative if exfiltration 01247 ! horizontal infiltration occurs along small side section, so scale the horizontal 01248 ! infiltration area by the vertical area of the cell used in the mass balance equation 01249 ! and infiltration velocities 01250 ! 01251 ! Scale area for bed cell and surface water cell 01252 scaleareab = dz/s%dnz(i,j) 01253 scaleareas = dz/s%dnz(i,j+1) 01254 ! Limit flow to that available in space in ground and 01255 ! in surface water volume 01256 if (vel>0.d0) then ! infiltration 01257 vel = min(vel, & 01258 s%hh(i,j+1)*s%dnz(i,j+1)/par%dt/dz, & 01259 max(s%zb(i,j)-s%gwlevel(i,j),0.d0)*par%por*s%dnz(i,j)/par%dt/dz) 01260 else 01261 vel = max(vel, & 01262 -(s%gwlevel(i,j)-s%gwbottom(i,j))*s%dnz(i,j)*par%por/par%dt/dz) 01263 endif 01264 infilhorgw(i,j) = infilhorgw(i,j)+vel*scaleareab 01265 infilhorsw(i,j+1) = infilhorsw(i,j+1)+vel*scaleareas 01266 else ! equal bed height 01267 ! nothing 01268 endif 01269 enddo 01270 enddo 01271 01272 end subroutine gw_horizontal_infil_exfil 01273 01274 subroutine gw_calculate_velocities(s,par,fracdt,gwu,gwv,gww,Kx,Ky,Kz,Kxupd,Kyupd,Kzupd) 01275 use params, only: parameters 01276 use xmpi_module 01277 use spaceparamsdef 01278 use paramsconst 01279 01280 IMPLICIT NONE 01281 01282 type(parameters),intent(in) :: par 01283 type(spacepars) :: s 01284 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: fracdt 01285 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: Kx,Ky,Kz 01286 real*8,dimension(s%nx+1,s%ny+1),intent(out) :: gwu,gwv,gww 01287 real*8,dimension(s%nx+1,s%ny+1),intent(out) :: Kxupd,Kyupd,Kzupd 01288 ! internal 01289 real*8,dimension(s%nx+1,s%ny+1) :: dheaddx,dheaddy 01290 real*8 :: vcr ! for turbulent/MODFLOW approximation 01291 real*8,parameter :: visc = 1.0d-6 01292 integer :: i,j 01293 01294 ! Calculate the head gradient 01295 call gwCalculateHeadGradient(s%nx,s%ny,s%gwhead,s%dsu,s%dnv,dheaddx,dheaddy) 01296 if (xmpi_istop) then 01297 dheaddx(1,:) = dheaddx(2,:) 01298 endif 01299 if (xmpi_isleft .and. s%ny>0) then 01300 dheaddy(:,1) = dheaddy(:,2) 01301 endif 01302 ! initialise 01303 Kxupd = Kx 01304 Kyupd = Ky 01305 Kzupd = Kz 01306 vcr = par%gwReturb/par%D50(1)*par%por*visc 01307 ! Compute of all loops 01308 ! u-direction 01309 do j=1,s%ny+1 01310 do i=1,s%nx+1 01311 call gw_calculate_velocities_local(gwu(i,j),Kxupd(i,j),dheaddx(i,j),Kx(i,j),0.d0, & 01312 par%gwscheme,.false.,par%wavemodel==WAVEMODEL_NONH,par%gwheadmodel, & 01313 s%gwcurv(i,j),s%gwheight(i,j),vcr) 01314 enddo 01315 enddo 01316 ! v-direction 01317 do j=1,s%ny+1 01318 do i=1,s%nx+1 01319 call gw_calculate_velocities_local(gwv(i,j),Kyupd(i,j),dheaddy(i,j),Ky(i,j),0.d0, & 01320 par%gwscheme,.false.,par%wavemodel==WAVEMODEL_NONH,par%gwheadmodel, & 01321 s%gwcurv(i,j),s%gwheight(i,j),vcr) 01322 enddo 01323 enddo 01324 ! w-direction 01325 do j=1,s%ny+1 01326 do i=1,s%nx+1 01327 call gw_calculate_velocities_local(gww(i,j),Kzupd(i,j),0.d0,Kz(i,j),fracdt(i,j), & 01328 par%gwscheme,.true.,par%wavemodel==WAVEMODEL_NONH,par%gwheadmodel, & 01329 s%gwcurv(i,j),s%gwheight(i,j),vcr) 01330 enddo 01331 enddo 01332 end subroutine gw_calculate_velocities 01333 01334 pure subroutine gw_calculate_velocities_local(vel,Kupd,headgrad,Kin,fracdt,gwscheme,isvert,isnonh,headmodel,gwc,gwh,vcr) 01335 ! compute local groundwater velocity component 01336 use paramsconst 01337 IMPLICIT NONE 01338 01339 ! input/output 01340 real*8,intent(out) :: vel,Kupd 01341 real*8,intent(in) :: headgrad,Kin,fracdt 01342 integer, intent(in) :: gwscheme 01343 integer, intent(in) :: headmodel 01344 logical,intent(in) :: isvert,isnonh 01345 real*8,intent(in) :: gwc,gwh,vcr 01346 ! internal 01347 real*8 :: vest,err,fac,kest,vupd 01348 01349 select case (gwscheme) 01350 case (GWSCHEME_LAMINAR) 01351 if (isvert .and. isnonh) then ! vertical flow non-hydrostatic 01352 select case(headmodel) 01353 case (GWHEADMODEL_PARABOLIC) 01354 vel = -Kin*2*gwc*gwh*(1.d0-fracdt) 01355 case (GWHEADMODEL_EXPONENTIAL) 01356 vel = -Kin*gwc*sinh(gwh)*(1.d0-fracdt) 01357 end select 01358 elseif (isvert .and. .not. isnonh) then ! vertical flow hydrostatic 01359 vel = 0.d0 01360 else ! horizontal flow 01361 vel = -Kin*headgrad 01362 endif 01363 Kupd = Kin 01364 case (GWSCHEME_TURBULENT) 01365 if (isvert .and. isnonh) then ! vertical flow non-hydrostatic 01366 select case(headmodel) 01367 case (GWHEADMODEL_PARABOLIC) 01368 vest = -Kin*2*gwc*gwh 01369 if (abs(vest)>vcr) then 01370 err = 1.d0 01371 do while (err>0.001d0) 01372 fac = min(sqrt(vcr/max(abs(vest),vcr)),1.d0) 01373 kest = Kin*fac 01374 vupd = -kest*2*gwc*gwh 01375 vest = (vest+vupd)/2 01376 err = abs(vupd-vest)/abs(vest) 01377 enddo 01378 Kupd = kest 01379 else 01380 Kupd = Kin 01381 endif 01382 if (isnonh) then 01383 vel = -Kin*2*gwc*gwh ! result bound by pressure solver estimate, will 01384 ! update next time step with better K approximation 01385 else 01386 vel = -Kupd*2*gwc*gwh 01387 endif 01388 case (GWHEADMODEL_EXPONENTIAL) 01389 vest = -Kin*gwc*sinh(gwh) 01390 if (abs(vest)>vcr) then 01391 err = 1.d0 01392 do while (err>0.001d0) 01393 fac = min(sqrt(vcr/max(abs(vest),vcr)),1.d0) 01394 kest = Kin*fac 01395 vupd = -kest*gwc*sinh(gwh) 01396 vest = (vest+vupd)/2 01397 err = abs(vupd-vest)/abs(vest) 01398 enddo 01399 Kupd = kest 01400 else 01401 Kupd = Kin 01402 endif 01403 if (isnonh) then 01404 vel = -Kin*gwc*sinh(gwh) ! result bound by pressure solver estimate, will 01405 ! update next time step with better K approximation 01406 else 01407 vel = -Kupd*gwc*sinh(gwh) 01408 endif 01409 end select ! head model 01410 vel = vel*(1.d0-fracdt) 01411 elseif (isvert .and. .not. isnonh) then ! vertical flow hydrostatic 01412 vel = 0.d0 01413 Kupd = Kin 01414 else ! horizontal flow 01415 vest = -Kin*headgrad 01416 if (abs(vest)>vcr) then 01417 err = 1.d0 01418 do while (err>0.001d0) 01419 fac = min(sqrt(vcr/max(abs(vest),vcr)),1.d0) 01420 kest = Kin*fac 01421 vupd = -kest*headgrad 01422 vest = (vest+vupd)/2 01423 err = abs(vupd-vest)/abs(vest) 01424 enddo 01425 Kupd = kest 01426 else 01427 Kupd = Kin 01428 endif 01429 if (isnonh) then 01430 vel = -Kin*headgrad ! result bound by pressure solver estimate, will 01431 ! update next time step with better K approximation 01432 else 01433 vel = -Kupd*headgrad 01434 endif 01435 endif 01436 end select ! turbulence model 01437 end subroutine gw_calculate_velocities_local 01438 01439 function gw_calculate_hydrostatic_w(s,par,ratio) result(infil) 01440 use params, only: parameters 01441 use xmpi_module 01442 use spaceparamsdef 01443 use paramsconst 01444 01445 IMPLICIT NONE 01446 01447 type(parameters) :: par 01448 type(spacepars),target :: s 01449 ! wwvv the following line from s.ind : 01450 real*8,dimension(:,:), pointer :: infil 01451 real*8,dimension(s%nx+1,s%ny+1) :: ratio 01452 ! local 01453 integer :: i,j 01454 real*8 :: w1,w2,w3,dummy,dummy2 01455 logical :: turbapprox 01456 01457 01458 ! wwvv the following line is from s.inp : 01459 infil => s%infil 01460 ! Select turbulent approximation of groundwater flow 01461 if (par%gwscheme==GWSCHEME_TURBULENT) then 01462 turbapprox = .true. 01463 else 01464 turbapprox = .false. 01465 endif 01466 ! 01467 ! Velocity in connected cells is based on part instantaneous recovery 01468 ! of volume (w1) and part darcy driven flow due to top pressure (w2). 01469 ! These velocities are weighted according to 'ratio', the closeness 01470 ! of the groundwater surface to the bed. 01471 if (s%ny==0) then 01472 do i=2,s%nx 01473 if (s%gwlevel(i,1)>s%zb(i,1)) then 01474 infil(i,1) = par%por*(s%zb(i,1)-s%gwlevel(i,1))/par%dt 01475 ! infil(i,1) = max(infil(i,1),-par%kz) 01476 elseif (s%wetz(i,1)==1) then 01477 ! subroutine returns infiltration rate 01478 call gw_calc_local_infil(par,s%zb(i,1),s%hh(i,1),s%gwlevel(i,1),par%kz,s%dinfil(i,1), & 01479 s%D50top(i,1),turbapprox,dummy,w1,dummy2) 01480 w1 = min(w1,(1.d0+s%hh(i,1)/(par%dwetlayer/3))*par%kz) 01481 ! part is instant, part is through infiltration 01482 w2 = par%por*(s%zb(i,1)-s%gwlevel(i,1))/par%dt 01483 w2 = min(w2,(s%hh(i,1)-par%eps)/par%dt) 01484 w3 = ratio(i,1)*w1+(1.d0-ratio(i,1))*w2 01485 ! this should not exceed the available space for water 01486 infil(i,1) = min(w3,(s%zb(i,1)-s%gwlevel(i,1))*par%por/par%dt) 01487 else 01488 infil(i,1) = 0.d0 01489 endif 01490 enddo 01491 ! boundaries 01492 if (xmpi_istop) then 01493 infil(1,1) = infil(2,1) 01494 endif 01495 if (xmpi_isbot) then 01496 infil(s%nx+1,1) = infil(s%nx,1) 01497 endif 01498 else ! s%ny>0 01499 do j=2,s%ny 01500 do i=2,s%nx 01501 if (s%gwlevel(i,j)>s%zb(i,j)) then 01502 infil(i,j) = par%por*(s%zb(i,j)-s%gwlevel(i,j))/par%dt 01503 elseif (s%wetz(i,j)==1) then 01504 ! subroutine returns infiltration rate 01505 call gw_calc_local_infil(par,s%zb(i,j),s%hh(i,j),s%gwlevel(i,j),par%kz,s%dinfil(i,j), & 01506 s%D50top(i,j),turbapprox,dummy,w1,dummy2) 01507 w1 = min(w1,(1.d0+s%hh(i,j)/(par%dwetlayer/3))*par%kz) 01508 ! part is instant, part is through infiltration 01509 w2 = par%por*(s%zb(i,j)-s%gwlevel(i,j))/par%dt 01510 w2 = min(w2,(s%hh(i,j)-par%eps)/par%dt) 01511 w3 = ratio(i,j)*w1+(1.d0-ratio(i,j))*w2 01512 ! this should not exceed the available space for water 01513 infil(i,j) = min(w3,(s%zb(i,j)-s%gwlevel(i,j))*par%por/par%dt) 01514 else 01515 infil(i,j) = 0.d0 01516 endif 01517 enddo 01518 enddo 01519 if (xmpi_istop) then 01520 infil(1,:) = infil(2,:) 01521 endif 01522 if (xmpi_isbot) then 01523 infil(s%nx+1,:) = infil(s%nx,:) 01524 endif 01525 if (xmpi_isleft) then 01526 infil(:,1) = infil(:,2) 01527 endif 01528 if (xmpi_isright) then 01529 infil(:,s%ny+1) = infil(:,s%ny) 01530 endif 01531 endif ! s%ny>0 01532 end function gw_calculate_hydrostatic_w 01533 01534 01535 pure subroutine gw_calc_local_infil(par,zb,hh,gwlevel,kzlocal,dinfil,D50top,turb,fracdt,infil,kzupd) 01536 use params, only: parameters 01537 01538 IMPLICIT NONE 01539 01540 type(parameters),intent(in) :: par 01541 real*8,intent(in) :: zb,hh,gwlevel,dinfil,D50top,kzlocal 01542 logical,intent(in) :: turb 01543 real*8,intent(out) :: fracdt,infil,kzupd 01544 ! local 01545 real*8,parameter :: visc = 1.0d-6 01546 real*8 :: dis 01547 real*8 :: vest,vupd,Re,err,kze,newinfil 01548 real*8,parameter :: beta = 3.236067977499790d0 ! golden ratio 01549 real*8,parameter :: alfa = 2.236067977499790d0 ! golden ratio 01550 logical :: firstloop 01551 real*8 :: twodtpdi,fourdtph,dtpsq,twodtp,disq,dtp 01552 ! 01553 ! Determine infiltration velocity in single cell 01554 ! 01555 ! The routine only solves Laminar/Turbulent unsaturated flow, 01556 ! saturated flow handled later by Poisson solver. The fraction of 01557 ! time that unsaturated flow occurs is given by fracdt. 01558 ! 01559 ! Unsaturated flow: 01560 ! w(n+1) = kz * (1+dp/dz) = kz * (1 + (hh)/(dz(n)+w(n+1)*dt/por) ) 01561 ! w -> (-dzold + dt/por k + Sqrt[dzold^2 + 2 dt/por dzold k + 4 dt/por hold k + dt^2 k^2])/(2 dt/por) 01562 ! 01563 ! The routine for turbulent flow is the same as for Darcy flow, 01564 ! but iterate until correct value of vest found 01565 ! 01566 ! 01567 ! Coefficients for implicit solution 01568 dtp = par%dt/par%por 01569 twodtp = 2*dtp 01570 twodtpdi = twodtp*dinfil 01571 fourdtph = 4*dtp*hh 01572 dtpsq = dtp**2 01573 disq = dinfil**2 01574 ! Compute infiltration velocity assuming fully-unsaturated flow 01575 if(turb) then 01576 ! First estimate without turbulence 01577 vest = ( (-dinfil + dtp*par%kz + & 01578 sqrt(disq + & 01579 twodtpdi*par%kz + & 01580 fourdtph*par%kz + & 01581 dtpsq*par%kz**2& 01582 )& 01583 )& 01584 /twodtp & 01585 ) 01586 Re = abs(vest)*D50top/par%por/visc 01587 if (Re>par%gwReturb) then 01588 err = 1.d0 01589 ! itialize with previous (turbulent) Kzinf 01590 kze = kzlocal 01591 firstloop = .true. 01592 do while (err>0.01d0) 01593 vupd = ( (-dinfil + dtp*kze + & 01594 sqrt(disq + & 01595 twodtpdi*kze + & 01596 fourdtph*kze + & 01597 dtpsq*kze**2& 01598 )& 01599 )& 01600 /twodtp & 01601 ) 01602 if (firstloop) then 01603 err = 1.d0 01604 firstloop = .false. 01605 else 01606 err = abs(vest-vupd)/abs(vest) 01607 endif 01608 ! New estimate of the velocity, using golden ratio 01609 vest = (vest+alfa*vupd)/beta 01610 ! update estimate of Kz for next iteration 01611 Re = abs(vest)*D50top/par%por/visc 01612 kze = par%kz*sqrt(par%gwReturb/Re) 01613 enddo 01614 kzupd = kze 01615 else 01616 kzupd = par%kz 01617 endif 01618 else 01619 ! No need to iterate 01620 vest = ( (-dinfil + dtp*par%kz + & 01621 sqrt(disq + & 01622 twodtpdi*par%kz + & 01623 fourdtph*par%kz + & 01624 dtpsq*par%kz**2& 01625 )& 01626 )& 01627 /twodtp & 01628 ) 01629 kzupd = par%kz 01630 endif 01631 ! now compute depth of infiltration at this rate 01632 dis = vest*par%dt/par%por 01633 ! what part of timestep is unsaturated, assuming 01634 ! (incorrectly) that this is linear? 01635 fracdt = max(0.d0,min((zb-gwlevel)/dis,1.d0)) 01636 infil = vest*fracdt ! this is the average over the WHOLE TIMESTEP 01637 ! check to see that the amount of water available is larger than the amount 01638 ! infiltrating 01639 if (infil > (hh-par%eps)/par%dt) then 01640 newinfil = max((hh-par%eps)/par%dt,0.d0) 01641 fracdt = fracdt*newinfil/infil 01642 infil = newinfil 01643 endif 01644 01645 01646 01647 01648 01649 01650 !! Distance from bed to groundwater level 01651 !dis = max(zb-gwlevel,0.d0) 01652 !! how long does it take to fill this distance under free vertical flow? 01653 !dtunsat = (dis*par%por) / (kzlocal*(1+hh/dis)) 01654 !! what part of timestep is unsaturated? 01655 !fracdt = min(dtunsat,par%dt)/par%dt 01656 !subdt = fracdt*par%dt/par%por 01657 !! does this need any unsaturated flow? 01658 !if (subdt>0.d0) then 01659 ! ! estimate vertical velocity due to unsaturated flow 01660 ! vest = ( (-dinfil + subdt*kzlocal + & 01661 ! sqrt(dinfil**2 + & 01662 ! 2*subdt*dinfil*kzlocal + & 01663 ! 4*subdt*hh*kzlocal + & 01664 ! subdt**2*kzlocal**2& 01665 ! )& 01666 ! )& 01667 ! /(2*subdt) & 01668 ! ) 01669 !else 01670 ! vest = 0.d0 01671 !endif 01672 !! Only set error if using turbulent flow approximation AND 01673 !! the Reynolds number is above the critical Re number 01674 !if (turb) then 01675 ! ! Calculate Reynolds number 01676 ! Re = abs(vest)*D50top/par%por/visc 01677 ! if (Re>par%gwReturb) then 01678 ! err = 1.d0 01679 ! endif 01680 !else 01681 ! err=-1.d0 01682 !endif 01683 !! update estimate if needed 01684 !kze = kzlocal ! note that this variable MUST remain separate from Kz(i,j) 01685 ! ! which is used to update the Poisson part of the equation 01686 !do while (err>0.001d0) 01687 ! Re = abs(vest)*D50top/par%por/visc 01688 ! kze = kzlocal*sqrt(par%gwReturb/Re) 01689 ! ! how long does it take to fill distance with new kz 01690 ! dtunsat = (dis*par%por) / (kze*(1+hh/dis)) 01691 ! ! what part of timestep unsaturated and saturated? 01692 ! fracdt = min(dtunsat,par%dt)/par%dt 01693 ! subdt = fracdt*par%dt/par%por 01694 ! ! new vertical velocity estimate 01695 ! vupd = ( (-dinfil + subdt*kze + & 01696 ! sqrt(dinfil**2 + & 01697 ! 2*subdt*dinfil*kze + & 01698 ! 4*subdt*hh*kze + & 01699 ! subdt**2*kze**2& 01700 ! )& 01701 ! )& 01702 ! /(2*subdt) & 01703 ! ) 01704 ! err = abs(vest-vupd)/abs(vest) 01705 ! vest = (vest+vupd)/2 01706 !enddo ! err>0.001 01707 !infil = vest*fracdt ! this is the average over the WHOLE TIMESTEP 01708 !! check to see that the amount of water available is larger than the amount 01709 !! infiltrating 01710 !if (infil > (hh-par%eps)/par%dt) then 01711 ! newinfil = max((hh-par%eps)/par%dt,0.d0) 01712 ! fracdt = fracdt*newinfil/infil 01713 ! infil = newinfil 01714 !endif 01715 !kzupd = kze 01716 end subroutine gw_calc_local_infil 01717 01718 ! TODO: subroutine is not used. Remove this 01719 #if 0 01720 pure subroutine gw_calc_local_connected_infil(par,gwhead,gwlevel,zb,hh,D50top,turb,zs,infil) 01721 use params, only: parameters 01722 01723 IMPLICIT NONE 01724 01725 type(parameters),intent(in) :: par 01726 real*8,intent(in) :: gwhead,gwlevel,zb,hh,D50top,zs 01727 logical,intent(in) :: turb 01728 real*8,intent(out) :: infil 01729 ! local 01730 real*8,parameter :: visc = 1.0d-6 01731 real*8 :: headdif 01732 real*8 :: vest,vupd,Re,err,kze 01733 logical :: infiltrate, exfiltrate 01734 01735 ! Determine infiltration and exfiltration velocity in single cell in case of hydrostatic 01736 ! groundwater computation 01737 ! 01738 ! Since in connected cells hydrostatic mode gwhead == headtop, we use the difference between 01739 ! headtop and gwhead when the cell needs filling, and the difference between gwlevel and zb 01740 ! when the cell needs emptying 01741 ! 01742 ! Saturated flow: 01743 ! infil = kz * (dp/dz) = kz * headdif/dz = kx * headdif / (infil*dt/por) 01744 ! infil = sqrt(kz * headdif * por/dt) 01745 ! 01746 err = 0 01747 01748 if(gwlevel<zb) then 01749 infiltrate = .true. 01750 headdif = min(zs-gwhead,hh) 01751 else 01752 infiltrate = .false. 01753 if(gwlevel>zb) then 01754 exfiltrate = .true. 01755 headdif = gwlevel-zb 01756 else 01757 exfiltrate = .false. 01758 endif 01759 endif 01760 01761 if (infiltrate .or. exfiltrate) then 01762 ! wwvv problem detected by compiler: 01763 ! wwvv it is possible the err does not get initialized 01764 ! wwvv so I initialize err at the beginning at the subroutine to 0 01765 ! compute infiltration rate 01766 kze = par%kx 01767 vest = sqrt(kze*headdif*par%por/par%dt) 01768 if (turb) then 01769 ! Calculate Reynolds number 01770 Re = abs(vest)*D50top/par%por/visc 01771 if (Re>par%gwReturb) then 01772 err = 1.d0 01773 endif 01774 else 01775 err=-1.d0 01776 endif 01777 ! Compute turbulent flow velocity 01778 do while (err>0.001d0) 01779 Re = abs(vest)*D50top/par%por/visc 01780 kze = par%kz*sqrt(par%gwReturb/Re) 01781 ! new vertical velocity estimate 01782 vupd = sqrt(kze*headdif*par%por/par%dt) 01783 err = abs(vest-vupd)/abs(vest) 01784 vest = (vest+vupd)/2 01785 enddo ! err>0.001 01786 else 01787 vest = 0.d0 01788 endif 01789 01790 if (exfiltrate) then 01791 vest = min(vest,(gwlevel-zb)/par%dt*par%por) 01792 infil = -vest 01793 else 01794 infil = min(vest,(hh-par%eps)/par%dt) 01795 infil = min(infil,(zb-gwlevel)/par%dt*par%por) 01796 endif 01797 01798 end subroutine gw_calc_local_connected_infil 01799 #endif 01800 01801 pure function gwCalculateHeadBottom(gwcurv,gwheadtop,gwheight,nx,ny,model) result(gwheadb) 01802 use paramsconst 01803 01804 IMPLICIT NONE 01805 integer,intent(in) :: nx,ny 01806 real*8,dimension(nx+1,ny+1),intent(in) :: gwcurv,gwheadtop,gwheight 01807 integer ,intent(in) :: model 01808 real*8,dimension(nx+1,ny+1) :: gwheadb 01809 ! internal 01810 01811 select case (model) 01812 case (GWHEADMODEL_PARABOLIC) 01813 gwheadb = gwheadtop-gwcurv*gwheight**2 01814 case (GWHEADMODEL_EXPONENTIAL) 01815 gwheadb = gwcurv+gwheadtop-gwcurv*cosh(gwheight) 01816 end select 01817 end function gwCalculateHeadBottom 01818 01819 01820 end module
1.7.6.1
