hydro/Routing/module_gw_gw2d.F

   1 !  Program Name:
   2 !  Author(s)/Contact(s):
   3 !  Abstract:
   4 !  History Log:
   5 !
   6 !  Usage:
   7 !  Parameters: <Specify typical arguments passed>
   8 !  Input Files:
   9 !        <list file names and briefly describe the data they include>
  10 !  Output Files:
  11 !        <list file names and briefly describe the information they include>
  12 !
  13 !  Condition codes:
  14 !        <list exit condition or error codes returned >
  15 !        If appropriate, descriptive troubleshooting instructions or
  16 !        likely causes for failures could be mentioned here with the
  17 !        appropriate error code
  18 !
  19 !  User controllable options: <if applicable>
  20
  21 !------------------------------------------------------------------------------
  22 ! Benjamin Fersch  2d groundwater model
  23 !------------------------------------------------------------------------------
  24
  25
  26 module module_gw_gw2d
  27
  28
  29 #ifdef MPP_LAND
  30    use module_mpp_land
  31 #endif
  32    use module_gw_gw2d_data, only: gw2d
  33    use module_rt_data, only: rt_domain
  34    use config_base, only: nlst
  35
  36    implicit none
  37
  38 #include "gw_field_include.inc"
  39
  40
  41 #ifdef MPP_LAND
  42  integer, private :: ierr
  43  integer, parameter :: rowshift = 0
  44  integer, parameter :: colshift = 1
  45 #endif
  46
  47
  48  contains
  49
  50
  51    subroutine gw2d_ini(did,dt,dx)
  52
  53      use module_HYDRO_io, only: output_gw_spinup
  54
  55      implicit none
  56      integer did
  57      real dt,dx
  58      integer :: jj, ii, iter, itermax
  59
  60
  61
  62
  63       itermax = nlst(did)%GwPreCycles
  64            gw2d(did)%dx=dx
  65            gw2d(did)%dt=dt
  66
  67            gw2d(did)%qgw_chanrt = 0.
  68            gw2d(did)%qsgwrt = 0.
  69            gw2d(did)%qdarcyRT = 0.
  70            gw2d(did)%excess = 0.
  71
  72            gw2d(did)%compres=0. ! currently not implemented
  73            gw2d(did)%istep=0 ! initialize time step
  74            ! reset cells with undefined hydraulic conductivity
  75            where(gw2d(did)%hycond .eq. 100) gw2d(did)%hycond = 5E-4
  76
  77           do iter=1,itermax
  78 #ifdef HYDRO_D
  79 #ifdef MPP_LAND
  80           if(my_id .eq. IO_id) &
  81 #endif
  82           write(6,*) "       GW Pre-cycle", iter, "of", itermax
  83 #endif
  84            call gwstep(gw2d(did)%ix, gw2d(did)%jx, gw2d(did)%dx, &
  85              gw2d(did)%ltype, gw2d(did)%elev, gw2d(did)%bot, &
  86              gw2d(did)%hycond, gw2d(did)%poros, gw2d(did)%compres, &
  87              gw2d(did)%ho, gw2d(did)%h, gw2d(did)%convgw, gw2d(did)%excess, &
  88              gw2d(did)%ebot, gw2d(did)%eocn, gw2d(did)%dt, &
  89              iter)
  90
  91              gw2d(did)%ho = gw2d(did)%h
  92
  93           if((nlst(did)%GwPreDiag .and. iter==1) .or. &
  94               nlst(did)%GwPreDiag .and. (mod(iter, nlst(did)%GwPreDiagInterval) .eq. 0) ) then
  95            call output_gw_spinup(nlst(did)%igrid, 1000000,                &
  96                             RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt,   &
  97                             nlst(did)%startdate, nlst(did)%olddate, &
  98                             gw2d(did)%ho, gw2d(did)%convgw, gw2d(did)%excess,  &
  99                             nlst(did)%geo_finegrid_flnm,nlst(did)%DT,     &
 100                             RT_DOMAIN(did)%LATVAL,        &
 101                             RT_DOMAIN(did)%LONVAL,rt_domain(did)%overland%properties%distance_to_neighbor,          &
 102                             nlst(did)%output_gw)
 103            end if
 104
 105
 106           end do
 107
 108    return
 109    end subroutine gw2d_ini
 110
 111    subroutine gw2d_allocate(did, ix, jx, nsoil)
 112
 113       implicit none
 114       integer ix, jx, nsoil
 115       integer istatus, did
 116
 117       if(gw2d(did)%allo_status .eq. 1) return
 118       gw2d(did)%allo_status = 1
 119
 120       gw2d(did)%ix = ix
 121       gw2d(did)%jx = jx
 122
 123 #ifdef MPP_LAND
 124       if(down_id == -1)  then !  if south border
 125        gw2d(did)%jts = 1
 126       else
 127        gw2d(did)%jts = 2
 128       endif
 129
 130       if(up_id == -1)    then !if north border
 131         gw2d(did)%jte = jx
 132       else
 133         gw2d(did)%jte = jx-1
 134       endif
 135
 136       if(left_id == -1)  then !if west border
 137         gw2d(did)%its = 1
 138       else
 139         gw2d(did)%its = 2
 140       endif
 141
 142       if(right_id == -1) then ! if east border
 143         gw2d(did)%ite = ix
 144       else
 145         gw2d(did)%ite = ix-1
 146       endif
 147
 148 #else
 149       gw2d(did)%its = 1
 150       gw2d(did)%ite = ix
 151       gw2d(did)%jts = 1
 152       gw2d(did)%jte = jx
 153 #endif
 154
 155       allocate(gw2d(did)%ltype  (ix,jx))
 156       allocate(gw2d(did)%elev   (ix,jx))
 157       allocate(gw2d(did)%bot    (ix,jx))
 158       allocate(gw2d(did)%hycond (ix,jx))
 159       allocate(gw2d(did)%poros  (ix,jx))
 160       allocate(gw2d(did)%compres(ix,jx))
 161       allocate(gw2d(did)%ho     (ix,jx))
 162       allocate(gw2d(did)%h      (ix,jx))
 163       allocate(gw2d(did)%convgw (ix,jx))
 164       allocate(gw2d(did)%excess (ix,jx))
 165
 166       allocate(gw2d(did)%qgw_chanrt (ix,jx))
 167
 168
 169       ! TODO allocate only if gwSoilCoupling is active
 170       allocate(gw2d(did)%qsgwrt   (ix,jx))
 171       allocate(gw2d(did)%qsgw     (rt_domain(did)%ix,rt_domain(did)%jx))
 172       allocate(gw2d(did)%qdarcyRT (ix,jx))
 173
 174     end subroutine gw2d_allocate
 175
 176
 177     subroutine gwstep(ix, jx, dx,              &
 178                       ltype, elev, bot,        &
 179                       hycond, poros, compres,  &
 180                       ho, h, convgw, excess,   &
 181                       ebot, eocn,              &
 182                       dt, istep)
 183
 184 ! New (volug): calling routines use change in head, convgw = d(h-ho)/dt.
 185
 186 ! Steps ground-water hydrology (head) through one timestep.
 187 ! Modified from Prickett and Lonnquist (1971), basic one-layer aquifer
 188 ! simulation program, with mods by Zhongbo Yu(1997).
 189 ! Solves S.dh/dt = d/dx(T.dh/dx) + d/dy(T.dh/dy) + "external sources"
 190 ! for a single layer, where h is head, S is storage coeff and T is
 191 ! transmissivity. 3-D arrays in main program (hycond,poros,h,bot)
 192 ! are 2-D here, since only a single (uppermost) layer is solved.
 193 ! Uses an iterative time-implicit ADI method.
 194
 195 ! use module_hms_constants
 196
 197
 198
 199       integer, intent(in) :: ix, jx
 200
 201       integer, intent(in), dimension(ix,jx) ::  ltype     ! land-sfc type  (supp)
 202       real,    intent(in), dimension(ix,jx) ::  &
 203         elev,           &  ! elev/bathymetry of sfc rel to sl (m) (supp)
 204         bot,            &  ! elev. aquifer bottom rel to sl (m)   (supp)
 205         hycond,         &  ! hydraulic conductivity (m/s per m/m) (supp)
 206         poros,          &  ! porosity (m3/m3)                     (supp)
 207         compres,        &  ! compressibility (1/Pa)               (supp)
 208         ho                 ! head at start of timestep (m)        (supp)
 209
 210       real,    intent(inout), dimension(ix,jx) ::  &
 211         h,              &  ! head, after ghmcompute (m)           (ret)
 212         convgw,         &  ! convergence due to gw flow (m/s)     (ret)
 213         excess
 214
 215       real, intent(inout) :: ebot, eocn
 216
 217
 218
 219       integer ::  istep !, dt
 220       real, intent(in) :: dt, dx
 221
 222 ! #endif
 223 !       eocn  = mean spurious sink for h_ocn = sealev fix (m/s)(ret)
 224 !               This equals the total ground-water flow across
 225 !               land->ocean boundaries.
 226 !       ebot  = mean spurious source for "bot" fix (m/s) (returned)
 227 !       time  = elapsed time from start of run (sec)
 228 !       dt = timestep length (sec)
 229 !       istep = timestep counter
 230
 231 ! Local arrays:
 232
 233       real, dimension(ix,jx)   :: sf2    ! storage coefficient (m3 of h2o / bulk m3)
 234       real, dimension(ix,jx,2) ::   t    ! transmissivity (m2/s)..1 for N-S,..2 for E-W
 235
 236 #ifdef MPP_LAND
 237       real, dimension(:,:), allocatable :: aa, &         ! tridiagonal matrix lower diagonal
 238                                            bb, &         ! tridiagonal matrix main diagonal
 239                                            cc, &         ! tridiagonal matrix upper diagonal
 240                                            dd, &         ! right hand side
 241                                            b2, &
 242                                            c2, &
 243                                            rhs, &
 244                                            wk, &
 245                                            hh
 246       real, dimension(:), allocatable ::   xfac, &
 247                                            zfac
 248 #else
 249       real, dimension(:), allocatable :: aa, &         ! tridiagonal matrix lower diagonal
 250                                          bb, &         ! tridiagonal matrix main diagonal
 251                                          cc, &         ! tridiagonal matrix upper diagonal
 252                                          dd, &         ! right hand side
 253                                          hh            ! solution vector
 254 #endif
 255       real, parameter    :: botinc = 0.01  ! re-wetting increment to fix h < bot
 256 !     parameter (botinc = 0.  )  ! re-wetting increment to fix h < bot
 257                                  ! (m); else no flow into dry cells
 258       real, parameter    :: delskip = 0.005 ! av.|dhead| value for iter.skip out(m)
 259       integer, parameter :: itermax = 1    ! maximum number of iterations
 260       integer, parameter :: itermin = 1    ! minimum number of iterations
 261       real, parameter    :: sealev = 1000.     ! sea-level elevation (m)
 262
 263       integer            :: its, ite, jts, jte, ifs, ife, jfs, jfe, &
 264                             xdim, ydim, fxdim, fydim
 265
 266 ! die müssen noch sortiert, geprüft und aufgeräumt werden
 267       integer ::                &
 268         iter,                   &
 269         j,                      &
 270         i,                      &
 271         jp,                     &
 272         ip,                     &
 273         n,                      &
 274         ierr,                   &
 275         ier,                    &
 276         ioffs,                  &
 277         joffs
 278
 279 !       real :: su, sc, shp, bb, aa, cc, w, zz, tareal, dtoa, dtot
 280       real ::                   &
 281         dy,                     &
 282         e,                      &
 283         su,                     &
 284         sc,                     &
 285         shp,                    &
 286         w,                      &
 287         ha,                     &
 288         delcur,                 &
 289         dtot,                   &
 290         dtoa,                   &
 291         darea,                  &
 292         tareal,                 &
 293         zz
 294
 295 #ifdef MPP_LAND
 296       real ::        mpiDelcur, &
 297                      gdtot,     &
 298                      gdtoa,     &
 299                      geocn,     &
 300                      gebot
 301       integer mpiSize
 302 #endif
 303
 304
 305
 306 dy = dx
 307 darea = dx*dy
 308
 309 ! define indexes for parallel execution
 310
 311 #ifdef MPP_LAND
 312 if(down_id == -1)  then !  if south border
 313   jts = 1
 314 else
 315   jts = 2
 316 endif
 317
 318 if(up_id == -1)    then !if north border
 319   jte = jx
 320 else
 321   jte = jx-1
 322 endif
 323
 324 if(left_id == -1)  then !if west border
 325   its = 1
 326 else
 327   its = 2
 328 endif
 329
 330 if(right_id == -1) then ! if east border
 331   ite = ix
 332 else
 333   ite = ix-1
 334 endif
 335
 336 #else
 337 its = 1
 338 ite = ix
 339 jts = 1
 340 jte = jx
 341 #endif
 342
 343 ifs = 1
 344 ife = ix
 345 jfs = 1
 346 jfe = jx
 347
 348
 349 fxdim = ife-ifs+1
 350 fydim = jfe-jfs+1
 351  xdim = ite-its+1
 352  ydim = jte-jts+1
 353
 354
 355       call scopy (fxdim*fydim, ho(ifs:ife,jfs:jfe), 1,    &
 356                   h(ifs:ife,jfs:jfe), 1)
 357
 358
 359 !       Top of iterative loop for (not anymore ADI) solution
 360
 361       iter = 0
 362 !~~~~~~~~~~~~~
 363    80 continue
 364 !~~~~~~~~~~~~~
 365       iter = iter+1
 366
 367
 368 #ifdef MPP_LAND
 369
 370        call MPP_LAND_COM_REAL(h, fxdim, fydim, 99)
 371
 372 #endif
 373       e    = 0.       ! absolute changes in head (for iteration control)
 374 !      eocn = 0.       ! accumulated fixes for h = 0 over ocean (diag)
 375 !      ebot = 0.       ! accumulated fixes for h < bot (diagnostic)
 376
 377 !       Set storage coefficient (sf2)
 378
 379
 380
 381     tareal = 0.
 382       do j=jts,jte
 383         do i=its,ite
 384
 385
 386         if(ltype(i,j) .ge. 1) tareal = tareal + darea
 387
 388 !         unconfined water table (h < e): V = poros*(h-b)
 389 !                                         dV/dh = poros
 390 !         saturated to surface (h >= e) : V = poros*(e-b) + (h-e)
 391 !                                         dV/dh = 1
 392 !         (compressibility is ignored)
 393 !
 394 !         su = poros(i,j)*(1.-theta(i,j))    ! old (pre-volug)
 395           su = poros(i,j)                    ! new (volug)
 396           sc = 1.
 397
 398 !           if      (ho(i,j).le.elev(i,j) .and. h(i,j).le.elev(i,j)) then
 399             sf2(i,j) = su
 400 !           else if (ho(i,j).ge.elev(i,j) .and. h(i,j).ge.elev(i,j)) then
 401 !             sf2(i,j) = sc
 402 !           else if (ho(i,j).le.elev(i,j) .and. h(i,j).ge.elev(i,j)) then
 403 !             shp = sf2(i,j) * (h(i,j) - ho(i,j))
 404 !             sf2(i,j) = shp * sc / (shp - (su-sc)*(elev(i,j)-ho(i,j)))
 405 !           else if (ho(i,j).ge.elev(i,j) .and. h(i,j).le.elev(i,j)) then
 406 !             shp = sf2(i,j) * (ho(i,j) - h(i,j))
 407 !             sf2(i,j) = shp * su / (shp + (su-sc)*(ho(i,j)-elev(i,j)))
 408 !           endif
 409
 410         enddo
 411       enddo
 412
 413 #ifdef MPP_LAND
 414        ! communicate storage coefficient
 415        call MPP_LAND_COM_REAL(sf2, fxdim, fydim, 99)
 416
 417 #endif
 418
 419 !==========================
 420 !       Column calculations
 421 !==========================
 422
 423 !       Set transmissivities. Use min(h,elev)-bot instead of h-bot,
 424 !       since if h > elev, thickness of groundwater flow is just
 425 !       elev-bot. (uses geometric mean)
 426
 427
 428       do j=jts,jte
 429         jp = min (j+1,jfe)
 430         do i=its,ite
 431           ip = min (i+1,ife)
 432
 433           t(i,j,2) = sqrt( abs(                                           &
 434                         hycond(i, j)*(min(h(i ,j),elev(i ,j))-bot(i ,j))  &
 435                        *hycond(ip,j)*(min(h(ip,j),elev(ip,j))-bot(ip,j))  &
 436                          )    )                                           &
 437                    * (0.5*(dy+dy)) & ! in WRF the dx and dy are usually equal
 438                    / (0.5*(dx+dx))
 439
 440           t(i,j,1) = sqrt( abs(                                           &
 441                         hycond(i,j )*(min(h(i,j ),elev(i,j ))-bot(i,j ))  &
 442                        *hycond(i,jp)*(min(h(i,jp),elev(i,jp))-bot(i,jp))  &
 443                          )    )                                           &
 444                    * (0.5*(dx+dx))  &
 445                    / (0.5*(dy+dy))
 446
 447
 448         enddo
 449       enddo
 450
 451
 452
 453
 454
 455 #ifdef MPP_LAND
 456       ! communicate transmissivities in x and y direction
 457        call MPP_LAND_COM_REAL(t(:,:,1), fxdim, fydim, 99)
 458        call MPP_LAND_COM_REAL(t(:,:,2), fxdim, fydim, 99)
 459
 460
 461        allocate(aa(jts:jte,its:ite))
 462        allocate(bb(jts:jte,its:ite))
 463        allocate(cc(jts:jte,its:ite))
 464        allocate(dd(jts:jte,its:ite))
 465        allocate(c2(1:ydim,1:xdim))
 466        allocate(b2(1:ydim,1:xdim))
 467        allocate(wk(1:ydim,1:xdim))
 468        allocate(hh(0:ydim+1,0:xdim+1))
 469        allocate(xfac(1:ydim))
 470        allocate(zfac(1:ydim))
 471 #else
 472   allocate(aa(jfs:jfe))
 473   allocate(bb(jfs:jfe))
 474   allocate(cc(jfs:jfe))
 475   allocate(dd(jfs:jfe))
 476   allocate(hh(jfs:jfe))
 477
 478 !-------------------
 479       do i=ifs,ife
 480 !-------------------
 481
 482 !>>>>>>>>>>>>>>>>>>>>
 483         do j=jfs,jfe
 484 !>>>>>>>>>>>>>>>>>>>>
 485 #endif
 486 #ifndef MPP_LAND
 487           bb(j) = (sf2(i,j)/dt) * darea
 488           dd(j) = ( ho(i,j)*sf2(i,j)/dt ) * darea
 489           aa(j) = 0.0
 490           cc(j) = 0.0
 491
 492           if ((j-jfs) /= 0) then
 493            aa(j) = -t(i,j-1,1)
 494            bb(j) = bb(j) + t(i,j-1,1)
 495           endif
 496
 497           if ((j-jfe) /= 0) then
 498            cc(j) = -t(i,j,1)
 499            bb(j) = bb(j) + t(i,j,1)
 500           endif
 501
 502           if ((i-ifs) /= 0) then
 503            bb(j) = bb(j) + t(i-1,j,2)
 504            dd(j) = dd(j) + h(i-1,j)*t(i-1,j,2)
 505           endif
 506
 507           if ((i-ife) /= 0) then
 508            bb(j) = bb(j) + t(i,j,2)
 509            dd(j) = dd(j) + h(i+1,j)*t(i,j,2)
 510           endif
 511
 512 !>>>>>>>>>>>>>>>
 513         end do
 514 !>>>>>>>>>>>>>>>
 515
 516   call trdiagSolve(aa, bb, cc, dd, hh, fydim)
 517
 518   h(i,:) = hh
 519   end do
 520
 521 deallocate(aa)
 522 deallocate(bb)
 523 deallocate(cc)
 524 deallocate(dd)
 525 deallocate(hh)
 526
 527 #else
 528 !-------------------
 529       do i=its,ite
 530 !-------------------
 531
 532 !>>>>>>>>>>>>>>>>>>>>
 533         do j=jts,jte
 534 !>>>>>>>>>>>>>>>>>>>>
 535           bb(j,i) = (sf2(i,j)/dt) * darea
 536           dd(j,i) = ( ho(i,j)*sf2(i,j)/dt ) * darea
 537           aa(j,i) = 0.0
 538           cc(j,i) = 0.0
 539
 540           if (((j-jfs) /= 0)) then
 541            aa(j,i) = -t(i,j-1,1)
 542            bb(j,i) = bb(j,i) + t(i,j-1,1)
 543           endif
 544
 545           if (((j-jfe) /= 0)) then
 546            cc(j,i) = -t(i,j,1)
 547            bb(j,i) = bb(j,i) + t(i,j,1)
 548           endif
 549
 550           if (((i-ifs) /= 0)) then
 551            bb(j,i) = bb(j,i) + t(i-1,j,2)
 552            dd(j,i) = dd(j,i) + h(i-1,j)*t(i-1,j,2)
 553           endif
 554
 555           if (((i-ife) /= 0)) then
 556            bb(j,i) = bb(j,i) + t(i,j,2)
 557            dd(j,i) = dd(j,i) + h(i+1,j)*t(i,j,2)
 558           endif
 559
 560 !>>>>>>>>>>>>>>>
 561         end do
 562 !>>>>>>>>>>>>>>>
 563
 564 !-------------
 565   end do
 566 !-------------
 567
 568     if(np_up_down .gt. 1) then
 569         call sub_n_form(xdim, ydim, aa, &
 570                         bb, cc, &
 571                         dd, &
 572                         c2, b2, hh, wk, xfac, zfac, &
 573                         p_up_down+1, np_up_down, 2)
 574
 575
 576         call parysolv1(c2, b2, hh, 1., my_id+1, p_up_down+1, &
 577                         xdim, ydim, np_left_right, np_up_down)
 578
 579     else
 580         call sub_tri_solv(xdim,ydim,aa(jts:jte,its:ite), &
 581                           bb(jts:jte,its:ite), cc(jts:jte,its:ite), &
 582                           dd(jts:jte,its:ite), &
 583                           hh, wk,xfac,zfac,2)
 584     endif
 585
 586 ioffs = its-1
 587 joffs = jts-1
 588 !-------------------
 589       do i=its,ite
 590 !-------------------
 591
 592 !>>>>>>>>>>>>>>>>>>>>
 593         do j=jts,jte
 594 !>>>>>>>>>>>>>>>>>>>>
 595
 596               h(i,j) = hh(j-joffs,i-ioffs)
 597
 598          end do
 599      end do
 600
 601 #endif
 602
 603 #ifdef MPP_LAND
 604
 605        call MPP_LAND_COM_REAL(h, fxdim, fydim, 99)
 606
 607 #endif
 608
 609
 610 !=======================
 611 !       Row calculations
 612 !=======================
 613
 614 !       set transmissivities (same as above)
 615
 616
 617       do j=jts,jte
 618         jp = min (j+1,jfe)
 619         do i=its,ite
 620           ip = min (i+1,ife)
 621           t(i,j,2) = sqrt( abs(                                            &
 622                         hycond(i, j)*(min(h(i ,j),elev(i ,j))-bot(i ,j))   &
 623                        *hycond(ip,j)*(min(h(ip,j),elev(ip,j))-bot(ip,j))   &
 624                          )    )                                            &
 625                    * (0.5*(dy+dy))                                         &
 626                    / (0.5*(dx+dx))
 627
 628           t(i,j,1) = sqrt( abs(                                            &
 629                         hycond(i,j )*(min(h(i,j ),elev(i,j ))-bot(i,j ))   &
 630                        *hycond(i,jp)*(min(h(i,jp),elev(i,jp))-bot(i,jp))   &
 631                          )    )                                            &
 632                    * (0.5*(dx+dx))                                         &
 633                    / (0.5*(dy+dy))
 634
 635
 636         enddo
 637       enddo
 638
 639 #ifdef MPP_LAND
 640       ! communicate transmissivities in x and y direction
 641        call MPP_LAND_COM_REAL(t(:,:,1), fxdim, fydim, 99)
 642        call MPP_LAND_COM_REAL(t(:,:,2), fxdim, fydim, 99)
 643 #endif
 644
 645 #ifndef MPP_LAND
 646 allocate(aa(ifs:ife))
 647 allocate(bb(ifs:ife))
 648 allocate(cc(ifs:ife))
 649 allocate(dd(ifs:ife))
 650 allocate(hh(ifs:ife))
 651
 652
 653 !-------------------
 654       do j=jfs,jfe
 655 !-------------------
 656
 657
 658 !>>>>>>>>>>>>>>>>>>>>
 659         do i=ifs,ife
 660 !>>>>>>>>>>>>>>>>>>>>
 661           bb(i) = (sf2(i,j)/dt) * darea
 662           dd(i) = ( ho(i,j)*sf2(i,j)/dt ) * darea
 663           aa(i) = 0.0
 664           cc(i) = 0.0
 665
 666           if ((j-jfs) /= 0) then
 667            bb(i) = bb(i) + t(i,j-1,1)
 668            dd(i) = dd(i) + h(i,j-1)*t(i,j-1,1)
 669           endif
 670
 671           if ((j-jfe) /= 0) then
 672            dd(i) = dd(i) + h(i,j+1)*t(i,j,1)
 673            bb(i) = bb(i) + t(i,j,1)
 674           endif
 675
 676           if ((i-ifs) /= 0) then
 677            bb(i) = bb(i) + t(i-1,j,2)
 678            aa(i) = -t(i-1,j,2)
 679           endif
 680
 681           if ((i-ife) /= 0) then
 682            bb(i) = bb(i) + t(i,j,2)
 683            cc(i) = -t(i,j,2)
 684           endif
 685
 686 !>>>>>>>>>>>>>>>
 687         end do
 688 !>>>>>>>>>>>>>>>
 689
 690   call trdiagSolve(aa, bb, cc, dd, hh, fxdim)
 691
 692   h(:,j) = hh
 693   end do
 694
 695 #else
 696 !-------------------
 697       do i=its,ite
 698 !-------------------
 699
 700 !>>>>>>>>>>>>>>>>>>>>
 701         do j=jts,jte
 702 !>>>>>>>>>>>>>>>>>>>>
 703           bb(j,i) = (sf2(i,j)/dt) * darea
 704           dd(j,i) = ( ho(i,j)*sf2(i,j)/dt ) * darea
 705           aa(j,i) = 0.0
 706           cc(j,i) = 0.0
 707
 708           if (((j-jfs) /= 0)) then
 709            bb(j,i) = bb(j,i) + t(i,j-1,1)
 710            dd(j,i) = dd(j,i) + h(i,j-1)*t(i,j-1,1)
 711           endif
 712
 713           if (((j-jfe) /= 0)) then
 714            dd(j,i) = dd(j,i) + h(i,j+1)*t(i,j,1)
 715            bb(j,i) = bb(j,i) + t(i,j,1)
 716           endif
 717
 718           if (((i-ifs) /= 0)) then
 719            bb(j,i) = bb(j,i) + t(i-1,j,2)
 720            aa(j,i) = -t(i-1,j,2)
 721           endif
 722
 723           if (((i-ife) /= 0)) then
 724            bb(j,i) = bb(j,i) + t(i,j,2)
 725            cc(j,i) = -t(i,j,2)
 726           endif
 727
 728 !>>>>>>>>>>>>>>>
 729         end do
 730 !>>>>>>>>>>>>>>>
 731
 732 !-------------
 733 end do
 734 !-------------
 735
 736     if(np_left_right .gt. 1) then
 737
 738 ! 3 c(,)  -- subdiagonal elements of tridiagonal systems
 739 ! 4 a(,)  -- diagonal elements of tridiagonal systems
 740 ! 5 b(,)  -- superdiagonal elements of tridiagonal systems
 741 ! 6 r(,)  -- right-hand side elements of tridiagonal systems
 742 ! 7 c2(,) -- front-leg elements of N-systems
 743 ! 8 b2(,) -- back-leg elements of N-systems
 744 ! 9 r2(,) -- right-hand side elements of N-systems (0:ydim+1,0:xdim+1)
 745 ! 10 wk(,) -- work array with same dimensions as a, b, c, etc.
 746
 747         call sub_n_form(xdim, ydim, aa, &
 748                         bb, cc, &
 749                         dd, &
 750                         c2, b2, hh, wk, xfac, zfac, &
 751                         p_left_right+1, np_left_right, 1)
 752
 753         call parxsolv1(c2, b2, hh, 1., my_id+1, p_left_right+1, &
 754                         xdim, ydim, np_left_right, np_up_down)
 755
 756     else
 757         call sub_tri_solv(xdim,ydim,aa, &
 758                           bb, cc, &
 759                           dd, &
 760                           hh, wk,xfac,zfac,1)
 761     endif
 762 ioffs = its-1
 763 joffs = jts-1
 764 !-------------------
 765       do i=its,ite
 766 !-------------------
 767
 768 !>>>>>>>>>>>>>>>>>>>>
 769         do j=jts,jte
 770 !>>>>>>>>>>>>>>>>>>>>
 771
 772                h(i,j) = hh(j-joffs,i-ioffs)
 773
 774       end do
 775      end do
 776
 777 deallocate(b2)
 778 deallocate(c2)
 779 deallocate(wk)
 780 deallocate(xfac)
 781 deallocate(zfac)
 782 #endif
 783 deallocate(aa)
 784 deallocate(bb)
 785 deallocate(cc)
 786 deallocate(dd)
 787 deallocate(hh)
 788
 789 ! fix head < bottom of aquifer
 790
 791       do j=jts,jte
 792         do i=its,ite
 793           if (ltype(i,j).eq.1 .and. h(i,j).le.bot(i,j)+botinc) then
 794
 795             e = e +  bot(i,j) + botinc - h(i,j)
 796 !             ebot = ebot + (bot(i,j)+botinc-h(i,j))*sf2(i,j)*darea(i,j)
 797             ebot = ebot + (bot(i,j)+botinc-h(i,j))*sf2(i,j)*darea
 798
 799             h(i,j) = bot(i,j) + botinc
 800           endif
 801         enddo
 802       enddo
 803 !        maintain head = sea level for ocean (only for adjacent ocean,
 804 !        rest has hycond=0)
 805
 806       do j=jts,jte
 807         do i=its,ite
 808           if (ltype(i,j).eq.2) then
 809
 810             eocn = eocn + (h(i,j)-sealev)*sf2(i,j)*darea
 811 !             eocn = eocn + (h(i,j)-sealev)*sf2(i,j)*darea(i,j)
 812
 813 !             h(i,j) = sealev (no update of outer boundary cells)
 814           endif
 815         enddo
 816       enddo
 817
 818 !        Loop back for next ADI iteration
 819
 820 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 821       delcur = e/(xdim*ydim)
 822
 823 !       print*, 'delcur before mpi:', delcur
 824
 825 #ifdef MPP_LAND
 826
 827 call mpi_reduce(delcur, mpiDelcur, 1, MPI_REAL, MPI_SUM, 0, HYDRO_COMM_WORLD, ierr)
 828 call MPI_COMM_SIZE( HYDRO_COMM_WORLD, mpiSize, ierr )
 829
 830 if(my_id .eq. IO_id) delcur = mpiDelcur/mpiSize
 831
 832 call mpi_bcast(delcur, 1, mpi_real, 0, HYDRO_COMM_WORLD, ierr)
 833
 834 #endif
 835
 836 !       if ( (delcur.gt.delskip*dt/86400. .and. iter.lt.itermax)      &
 837       if ( (delcur.gt.delskip .and. iter.lt.itermax)      &
 838            .or. iter.lt.itermin ) then
 839
 840 #ifdef HYDRO_D
 841
 842 #ifdef MPP_LAND
 843 if(my_id .eq. IO_id)  write(6,*) "Iteration", iter, "of", itermax, "error:", delcur
 844 #else
 845                       write(6,*) "Iteration", iter, "of", itermax, "error:", delcur
 846 #endif
 847
 848 #endif
 849
 850       goto 80
 851       endif
 852
 853 #ifdef MPP_LAND
 854
 855        call MPP_LAND_COM_REAL(h, fxdim, fydim, 99)
 856
 857 #endif
 858
 859
 860
 861 !     Compute exfiltration amount and
 862 !     convergence rate due to ground water
 863 !     flow
 864
 865       do j=jts,jte
 866         do i=its,ite
 867
 868           if((elev(i,j) - h(i,j)) .lt. 0.) then
 869             excess(i,j) = sf2(i,j)*(h(i,j) - elev(i,j))
 870                  h(i,j) = elev(i,j)
 871           else
 872             excess(i,j) = 0.
 873           end if
 874
 875           if(ltype(i,j).eq.1) then
 876             convgw(i,j) = sf2(i,j) * (h(i,j)-ho(i,j)) / dt
 877           else
 878             convgw(i,j) = 0.
 879           endif
 880         enddo
 881       enddo
 882
 883 !      call MPP_LAND_COM_REAL(convgw, fxdim, fydim, 99)
 884
 885 !        Diagnostic water conservation check for this timestep
 886
 887       dtot = 0.     ! total change in water storage (m3)
 888       dtoa = 0.
 889
 890       do j=jts,jte
 891         do i=its,ite
 892           if (ltype(i,j).eq.1) then
 893
 894             dtot = dtot + sf2(i,j) *(h(i,j)-ho(i,j)) * darea
 895             dtoa = dtoa + sf2(i,j) * abs(h(i,j)-ho(i,j)) * darea
 896
 897 !             dtot = dtot + sf2(i,j) *(h(i,j)-ho(i,j)) * darea(i,j)
 898 !             dtoa = dtoa + sf2(i,j) * abs(h(i,j)-ho(i,j)) * darea(i,j)
 899           endif
 900         enddo
 901       enddo
 902
 903       dtot = (dtot/tareal)/dt   ! convert to m/s, rel to land area
 904       dtoa = (dtoa/tareal)/dt
 905       eocn = (eocn/tareal)/dt
 906       ebot = (ebot/tareal)/dt
 907
 908       zz = 1.e3 * 86400.                    ! convert printout to mm/day
 909 #ifdef HYDRO_D
 910 #ifdef MPP_LAND
 911
 912       call MPI_REDUCE(dtot,gdtot,1, MPI_REAL, MPI_SUM, IO_id, HYDRO_COMM_WORLD, ierr)
 913       call MPI_REDUCE(dtoa,gdtoa,1, MPI_REAL, MPI_SUM, IO_id, HYDRO_COMM_WORLD, ierr)
 914       call MPI_REDUCE(eocn,geocn,1, MPI_REAL, MPI_SUM, IO_id, HYDRO_COMM_WORLD, ierr)
 915       call MPI_REDUCE(ebot,gebot,1, MPI_REAL, MPI_SUM, IO_id, HYDRO_COMM_WORLD, ierr)
 916
 917       if(my_id .eq. IO_id) then
 918         write (*,900)                         &
 919           gdtot*zz, gdtoa*zz, -geocn*zz, gebot*zz,     &
 920           (gdtot-(-geocn+gebot))*zz
 921        endif
 922
 923 #else
 924
 925         write (*,900)                         &
 926           dtot*zz, dtoa*zz, -eocn*zz, ebot*zz,     &
 927           (dtot-(-eocn+ebot))*zz
 928 #endif
 929 #endif
 930   900 format                                       &
 931         (3x,'    dh/dt       |dh/dt|        ocnflx        botfix',&
 932             '        ghmerror'  &
 933 !         /3x,4f9.4,2(9x),e14.4)
 934         /3x,5(e14.4))
 935
 936       return
 937       end subroutine gwstep
 938
 939
 940       SUBROUTINE SCOPY (NT, ARR, INCA, BRR, INCB)
 941 !
 942 !        Copies array ARR to BRR, incrementing by INCA and INCB
 943 !        respectively, up to a total length of NT words of ARR.
 944 !        (Same as Cray SCOPY.)
 945 !
 946       real, DIMENSION(*) :: ARR, BRR
 947       integer :: ia, nt, inca, incb, ib
 948 !
 949       IB = 1
 950       DO 10 IA=1,NT,INCA
 951          BRR(IB) = ARR(IA)
 952          IB = IB + INCB
 953    10 CONTINUE
 954 !
 955       RETURN
 956       END SUBROUTINE SCOPY
 957
 958
 959 subroutine trdiagSolve(a,b,c,rhs,x,n)
 960
 961       implicit none
 962
 963       integer,intent(in) :: n
 964       real,dimension(n),intent(in) :: a, b, c, rhs
 965       real,dimension(n),intent(out) :: x
 966       real,dimension(n) :: cp, dp
 967       real :: m
 968       integer i
 969
 970 ! initialize c-prime and d-prime
 971         cp(1) = c(1)/b(1)
 972         dp(1) = rhs(1)/b(1)
 973 ! solve for vectors c-prime and d-prime
 974          do i = 2,n
 975            m = b(i)-cp(i-1)*a(i)
 976            cp(i) = c(i)/m
 977            dp(i) = (rhs(i)-dp(i-1)*a(i))/m
 978          enddo
 979 ! initialize x
 980          x(n) = dp(n)
 981 ! solve for x from the vectors c-prime and d-prime
 982         do i = n-1, 1, -1
 983           x(i) = dp(i)-cp(i)*x(i+1)
 984         end do
 985
 986
 987 end subroutine trdiagSolve
 988
 989
 990 subroutine gwSoilFlux(did)
 991
 992
 993   implicit none
 994
 995   integer, intent(in)   :: did
 996
 997
 998   real, dimension(rt_domain(did)%ixrt,rt_domain(did)%jxrt) :: smcrel, ztrans, headChange
 999   real :: frac, zres
1000   integer :: nsoil, i, j, k
1001
1002   gw2d(did)%qsgwrt = 0.
1003   gw2d(did)%qdarcyRT = 0.
1004
1005 ! Step 1, collect data
1006
1007 ! relative soil moisture content of lowest soil layer (1 = saturated)
1008   nsoil = nlst(did)%nsoil
1009   smcrel = rt_domain(did)%subsurface%grid_transform%smcrt(:,:,nsoil) / RT_DOMAIN(did)%subsurface%grid_transform%smcmaxrt(:,:,nsoil)
1010
1011 ! depth of transition zone from lowest soil layer to groundwater head (in cm)
1012 ! postivie ztrans -> head below LSM soil layer
1013 ! negative ztrans -> head within LSM soil layers
1014   ztrans = (rt_domain(did)%elrt + nlst(did)%zsoil8(nsoil)) - gw2d(did)%ho
1015   ztrans = ztrans * 100
1016
1017   ! darcyGwSoil not defined for ztran = 0
1018   where(ztrans == 0) ztrans = -5
1019
1020 ! Step 2, compute flux either up or down
1021
1022   do j=gw2d(did)%jts, gw2d(did)%jte
1023     do i=gw2d(did)%its, gw2d(did)%ite
1024
1025         if((ztrans(i,j) > 0) .and. (rt_domain(did)%soiltypRT(i,j) < 13)) then
1026         ! if groundwater head < soil layers
1027           call  darcyGwSoil(ztrans(i,j), smcrel(i,j), rt_domain(did)%soiltypRT(i,j), gw2d(did)%qdarcyRT(i,j))
1028
1029           gw2d(did)%qsgwrt(i,j) = gw2d(did)%qdarcyRT(i,j)
1030
1031           ! check and correct for mass balance
1032           if(((gw2d(did)%ho(i,j)-gw2d(did)%bot(i,j)) &
1033             *gw2d(did)%poros(i,j)) < (gw2d(did)%qsgwrt(i,j)*gw2d(did)%dt)) then
1034
1035                 gw2d(did)%qdarcyRT(i,j) = 0.
1036                 gw2d(did)%qsgwrt(i,j) = 0.
1037
1038            end if
1039
1040         else if(ztrans(i,j) < 0 .and. (rt_domain(did)%soiltypRT(i,j) < 13)) then
1041         ! if groundwater head > soil layers
1042           zres = -ztrans(i,j)
1043           do k=nsoil,1,-1
1044
1045              if(zres >= rt_domain(did)%subsurface%properties%sldpth(k)*100.) then
1046              ! complete filling of a LSM soil layer if groundwater head > layer top
1047
1048 !              gw2d(did)%qsgwrt(i,j) = (rt_domain(did)%subsurface%properties%sldpth(k) &
1049 !                                      * (rt_domain(did)%subsurface%grid_transform%smcmaxrt(i,j,k) - RT_DOMAIN(did)%subsurface%grid_transform%smcrt(i,j,k)) &
1050 !                                      + gw2d(did)%qsgwrt(i,j)) / gw2d(did)%dt
1051
1052                rt_domain(did)%subsurface%grid_transform%smcrt(i,j,k) = RT_DOMAIN(did)%subsurface%grid_transform%smcmaxrt(i,j,k)
1053
1054                zres = zres - rt_domain(did)%subsurface%properties%sldpth(k)*100.
1055
1056              else
1057              ! partial filling of a soil layer if not completely below groundwater head
1058
1059                if(zres > (0.5 * rt_domain(did)%subsurface%properties%sldpth(k)*100.)) then
1060
1061                  frac = zres / (rt_domain(did)%subsurface%properties%sldpth(k) * 100.)
1062
1063
1064 !                 gw2d(did)%qsgwrt(i,j) = (rt_domain(did)%subsurface%properties%sldpth(k) &
1065 !                                       * (rt_domain(did)%subsurface%grid_transform%smcmaxrt(i,j,k) - RT_DOMAIN(did)%subsurface%grid_transform%smcrt(i,j,k)) &
1066 !                                       * frac + gw2d(did)%qsgwrt(i,j)) / gw2d(did)%dt
1067
1068                   rt_domain(did)%subsurface%grid_transform%smcrt(i,j,k) = RT_DOMAIN(did)%subsurface%grid_transform%smcmaxrt(i,j,k) * frac
1069
1070                end if
1071
1072              end if
1073           end do
1074         end if
1075     end do
1076   end do
1077
1078           ! sign convention
1079           ! qsgwrt < 0 -> downward flux
1080           ! qsgwrt > 0 -> upward flux
1081
1082 ! TOcheck Step 3, adapt groundwater head (assuming not time lag for percolation / capillary rise flow)
1083
1084 !          modify gw-head before gwstep call with respect to specific yield of the
1085 !          aquifer and the computed flux (qsgwrt)
1086
1087
1088  headChange = (-gw2d(did)%qdarcyRT) * gw2d(did)%dt / gw2d(did)%poros
1089  gw2d(did)%ho = gw2d(did)%ho + headChange
1090
1091 end subroutine gwSoilFlux
1092
1093 subroutine darcyGwSoil(Z, s, soil, q_darcy)
1094
1095 implicit none
1096
1097 INTEGER, INTENT (IN)  :: soil ! soiltype
1098
1099 REAL :: sig_a, sig_b, sig_c
1100
1101 REAL, DIMENSION(9)    :: k_para
1102 REAL, INTENT (IN)     :: Z, s
1103 REAL, INTENT (OUT)    :: q_darcy
1104 real                  :: beta,alpha,q_cap,b,ks,aep,c,q_grav,y,fac
1105
1106 real, dimension(9,12) :: &
1107       k_soil = reshape((/&
1108 0.0778, 3.9939, 0.2913, 4.0801, 0.1386, 4.0500, -12.10, 0.3950, 1.0560,&
1109 0.0924, 4.8822, 0.2674, 3.8915, 0.1365, 4.3800, -09.00, 0.4100, 0.9380,&
1110 0.0367, 4.5259, 0.2446, 4.2849, 0.1208, 4.9000, -21.80, 0.4350, 0.2080,&
1111 0.0101, 3.6896, 0.2153, 4.2765, 0.0887, 5.3000, -78.60, 0.4850, 0.0432,&
1112 0.0101, 3.6896, 0.2153, 4.2765, 0.0887, 5.3000, -78.60, 0.4850, 0.0432,&
1113 0.0169, 2.9936, 0.2858, 4.3738, 0.1026, 5.3900, -47.80, 0.4510, 0.0417,&
1114 0.0271, 4.4743, 0.2587, 3.9055, 0.0920, 7.1200, -29.90, 0.4200, 0.0378,&
1115 0.0227, 4.3768, 0.2658, 3.8234, 0.0843, 7.7500, -35.60, 0.4770, 0.0102,&
1116 0.0127, 6.6836, 0.1725, 3.7512, 0.0703, 8.5200, -63.00, 0.4760, 0.0147,&
1117 0.0530, 9.2423, 0.1859, 3.3688, 0.0728, 10.400, -15.30, 0.4260, 0.0130,&
1118 0.0165, 5.3972, 0.2479, 3.5549, 0.0641, 10.400, -49.00, 0.4920, 0.0062,&
1119 0.0200, 6.0106, 0.2474, 3.4788, 0.0622, 11.400, -40.50, 0.4820, 0.0077/),(/9,12/))
1120
1121
1122
1123  k_para  = k_soil(:,soil)
1124  sig_a   = 1 - exp( -1 * k_para(1) * Z)
1125  sig_b   = k_para(2) * Z**k_para(3)
1126  sig_c   = k_para(4) * exp( -1 * Z**k_para(5))
1127  y       = sig_a/(1  + exp(sig_b * (s - sig_c))) !solving equation (20) in Boogart et al.
1128
1129  b   =   k_para(6)
1130  ks  =   k_para(9)
1131  aep =  -k_para(7)
1132
1133  c      =  2 * b  + 3
1134  q_grav = -1 * ks * s**c
1135
1136 ! alp is constant from equation (13) of paper
1137 beta  = 2 + 3 / b
1138 alpha = 1 + 1.5 /  (beta - 1)
1139 q_cap = ks * alpha * (aep / Z)**beta
1140
1141
1142 q_darcy = y * q_cap + q_grav ![cm/min]
1143
1144 ! limit for exteme gradients with q >> saturated hydraulic conductivity
1145 ! if(q_cap > ks) q_cap = ks
1146 ! if(q_grav < -ks) q_grav = -ks
1147
1148 ! if(q_darcy > ks) q_darcy = ks
1149 ! if(q_darcy < ks) q_darcy = -ks
1150
1151
1152 fac     = 1./6000.
1153 q_darcy = q_darcy * fac
1154 q_cap   = q_cap   * fac
1155 q_grav  = q_grav  * fac
1156
1157 !returns q_darcy in [m/s]
1158
1159 end subroutine darcyGwSoil
1160
1161
1162
1163 subroutine aggregateQsgw(did)
1164
1165
1166
1167   implicit none
1168
1169    integer, intent(in) :: did
1170    integer :: j,i, ixxRT, jyyRT, m,n
1171    real :: agg
1172
1173
1174     do j=1,rt_domain(did)%jx
1175      do i=1,rt_domain(did)%ix
1176
1177        agg= 0.
1178
1179        do m=nlst(did)%aggfactRT-1,0,-1
1180          do n=nlst(did)%aggfactRT-1,0,-1
1181
1182
1183             ixxRT = i * nlst(did)%aggfactRT-n
1184             jyyRT = j * nlst(did)%aggfactRT-m
1185
1186
1187 #ifdef MPP_LAND
1188             if(left_id.ge.0) ixxRT=ixxRT+1
1189             if(down_id.ge.0) jyyRT=jyyRT+1
1190 #endif
1191              agg = agg + gw2d(did)%qdarcyRT(ixxRT, jyyRT)
1192            end do
1193          end do
1194
1195             gw2d(did)%qsgw(i,j) = agg/(nlst(did)%aggfactRT**2)
1196        end do
1197      end do
1198
1199
1200
1201 end subroutine aggregateQsgw
1202
1203 ! Parallel tridiagonal solver useful for domain decomposed ADI
1204 ! Author(s): Mike Lambert
1205 ! Year: 1996
1206 ! Institution: Lawrence Livermore National Laboratory
1207 ! Publication: Lambert, Rodrigue, and Hewett, "A parallel DSDADI method
1208 !                      for solution of the steady state diffusion equation",
1209 !                      Parallel Computing 23 (1997) 2041-2065
1210 ! Ported to MPI by Benjamin Fersch, Karlsruhe Institute of Technology (2013)
1211
1212 #ifdef MPP_LAND
1213       subroutine parysolv1(c,b,r,ct,pid,z_pid, &
1214                             xsps, zsps, xdns, zdns)
1215
1216       implicit none
1217
1218       integer, intent(in) :: XSPS, &
1219                              ZSPS, &
1220                              XDNS, &
1221                              ZDNS
1222
1223       real, dimension(ZSPS, XSPS), intent(inout) ::  c, &
1224                                                      b
1225       real      CLK_PER
1226       parameter (CLK_PER = 6.66666667e-9)
1227
1228       real, dimension(0:ZSPS+1, 0:XSPS+1), intent(inout) ::  r
1229
1230       real, dimension(XSPS,2) :: zn, zntmp
1231
1232       real, dimension(XSPS)   :: t1, t2, fac
1233
1234       real :: clockdt, click
1235       real :: ct, ti, tf, dt
1236
1237       integer :: pid, z_pid
1238       integer :: i, j, sndr_pid, msg_type, cnt, ackn
1239       integer :: sendReq, recvReq
1240
1241       integer   ZN_REC
1242       parameter (ZN_REC = 46)
1243
1244       integer :: source, dest
1245 #ifdef TIMING
1246       dt = clockdt()
1247 #endif
1248
1249       cnt = 2*XSPS
1250
1251       if (z_pid .eq. 1) then
1252
1253 ! Load (ZSPS,j)th equations into passing arrays.
1254         do 10 j = 1, XSPS
1255           zntmp(j,1) = b(ZSPS,j)
1256           zntmp(j,2) = r(ZSPS,j)
1257    10   continue
1258
1259
1260 #ifdef TIMING
1261         ti = click()
1262 #endif
1263
1264 ! ! Send (ZSPS,j)th equations.
1265 ! ! Receive (ZSPS+1,j)th equations.
1266
1267  call mpi_cart_shift(cartGridComm, rowshift, 1, source, dest, ierr)
1268  call MPI_ISEND(zntmp, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, sendReq, ierr)
1269  call MPI_IRECV(   zn, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, recvReq, ierr)
1270  call mpi_wait(sendReq, mpp_status, ierr)
1271  call mpi_wait(recvReq, mpp_status, ierr)
1272
1273 #ifdef TIMING
1274         tf = click()
1275         call add_dt(ct,tf,ti,dt)
1276 #endif
1277
1278         do 20 j = 1, XSPS
1279 ! Backward elimination in (ZSPS,j)th equations to get
1280 ! r(ZSPS,j).
1281         fac(j) = 1./(1. - b(ZSPS,j)*zn(j,1))
1282         r(ZSPS,j) = (r(ZSPS,j)-b(ZSPS,j)*zn(j,2))*fac(j)
1283 ! Forward elimination in (ZSPS+1,j)th equations to get
1284 ! r(ZSPS+1,j).
1285         r(ZSPS+1,j) = zn(j,2) - zn(j,1)*r(ZSPS,j)
1286 ! Completion of backward elimination to get remaining unknowns.
1287         do 30 i = 1, ZSPS-1
1288           r(i,j) = r(i,j) - b(i,j)*r(ZSPS,j)
1289    30   continue
1290    20   continue
1291
1292       else if (z_pid .le. ZDNS/2) then
1293
1294 #ifdef TIMING
1295         ti = click()
1296 #endif
1297 ! ! Receive (0,j)th equations.
1298
1299  call mpi_cart_shift(cartGridComm, rowshift, -1, source, dest, ierr)
1300  call MPI_IRECV(   zn, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, recvReq, ierr)
1301  call mpi_wait(recvReq, mpp_status, ierr)
1302
1303 #ifdef TIMING
1304         tf = click()
1305         call add_dt(ct,tf,ti,dt)
1306 #endif
1307
1308 ! Forward elimination in (j,1)th equations.
1309         do 40 j = 1, XSPS
1310           fac(j) = 1./(1. - c(1,j)*zn(j,1))
1311 ! Check for singular matrix (debugging only)
1312           b(1,j) = b(1,j)*fac(j)
1313           r(1,j) = (r(1,j) - c(1,j)*zn(j,2))*fac(j)
1314 ! Forward elimination in (ZSPS,j)th equations.
1315           fac(j) = 1./(1. - c(ZSPS,j)*b(1,j))
1316 ! Check for singular matrix (debugging only)
1317           b(ZSPS,j) = b(ZSPS,j)*fac(j)
1318           r(ZSPS,j) = (r(ZSPS,j)-c(ZSPS,j)*r(1,j))*fac(j)
1319 ! Store (0,j)th equations for later recovery of r(0,j).
1320           t1(j) = zn(j,1)
1321           t2(j) = zn(j,2)
1322 ! Load (ZSPS,j)th equations into passing arrays.
1323           zntmp(j,1) = b(ZSPS,j)
1324           zntmp(j,2) = r(ZSPS,j)
1325    40   continue
1326
1327 #ifdef TIMING
1328         ti = click()
1329 #endif
1330 ! ! Send (ZSPS,j)th equations.
1331 ! ! Receive (ZSPS+1,j)th equations.
1332
1333  call mpi_cart_shift(cartGridComm, rowshift, 1, source, dest, ierr)
1334  call MPI_ISEND(zntmp, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, sendReq, ierr)
1335  call MPI_IRECV(   zn, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, recvReq, ierr)
1336  call mpi_wait(sendReq, mpp_status, ierr)
1337  call mpi_wait(recvReq, mpp_status, ierr)
1338 #ifdef TIMING
1339         tf = click()
1340         call add_dt(ct,tf,ti,dt)
1341 #endif
1342
1343         do 50 j = 1, XSPS
1344 ! Backward elimination in (ZSPS,j)th equations.
1345           fac(j) = 1./(1. - b(ZSPS,j)*zn(j,1))
1346 ! Check for singular matrix (debugging only)
1347           r(ZSPS,j) = (r(ZSPS,j) - b(ZSPS,j)*zn(j,2))*fac(j)
1348 ! Backward elimination in (ZSPS+1,j)th equations.
1349           r(ZSPS+1,j) = zn(j,2) - zn(j,1)*r(ZSPS,j)
1350 ! Backward elimination in (ZSPS,j)th equations.
1351           r(1,j) = r(1,j) - b(1,j)*r(ZSPS,j)
1352 ! Load (1,j)th equations into passing arrays.
1353           zntmp(j,1) = 0.
1354           zntmp(j,2) = r(1,j)
1355    50   continue
1356
1357 #ifdef TIMING
1358         ti = click()
1359 #endif
1360 ! ! Send (1,j)th equations.
1361
1362 #ifdef TIMING
1363         tf = click()
1364         call add_dt(ct,tf,ti,dt)
1365 #endif
1366
1367  call mpi_cart_shift(cartGridComm, rowshift, -1, source, dest, ierr)
1368  call MPI_ISEND(zntmp, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, sendReq, ierr)
1369
1370         do 60 j = 1, XSPS
1371 ! Backward elimination in (0,j)th equations.
1372         r(0,j) = t2(j) - t1(j)*r(1,j)
1373         do 70 i = 2, ZSPS-1
1374 ! Completion of forward and backward elimination to get remaining
1375 ! unknowns.
1376           r(i,j) = r(i,j) - b(i,j)*r(ZSPS,j) - c(i,j)*r(1,j)
1377    70   continue
1378    60   continue
1379
1380  call mpi_wait(sendReq, mpp_status, ierr)
1381
1382
1383       else if (z_pid .lt. ZDNS) then
1384
1385 #ifdef TIMING
1386         ti = click()
1387 #endif
1388 ! ! Receive (ZSPS+1,j)th equations.
1389
1390  call mpi_cart_shift(cartGridComm, rowshift, 1, source, dest, ierr)
1391  call MPI_IRECV(   zn, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, recvReq, ierr)
1392  call mpi_wait(recvReq, mpp_status, ierr)
1393
1394 #ifdef TIMING
1395         tf = click()
1396         call add_dt(ct,tf,ti,dt)
1397 #endif
1398
1399         do 80 j = 1, XSPS
1400 ! Backward elimination in (ZSPS,j)th equations.
1401           fac(j) = 1./(1. - b(ZSPS,j)*zn(j,1))
1402 ! Check for singular matrix (debugging only)
1403           c(ZSPS,j) = c(ZSPS,j)*fac(j)
1404           r(ZSPS,j) = (r(ZSPS,j)-b(ZSPS,j)*zn(j,2))*fac(j)
1405 ! Backward elimination in (1,j)th equations.
1406           fac(j) = 1./(1. - b(1,j)*c(ZSPS,j))
1407 ! Check for singular matrix (debugging only)
1408           c(1,j) = c(1,j)*fac(j)
1409           r(1,j) = (r(1,j) - b(1,j)*r(ZSPS,j))*fac(j)
1410 ! Store (ZSPS+1,j)th equations for later recovery of
1411 ! r(ZSPS+1,j).
1412           t1(j) = zn(j,1)
1413           t2(j) = zn(j,2)
1414 ! Load passing arrays with (1,j)th equations.
1415           zntmp(j,1) = c(1,j)
1416           zntmp(j,2) = r(1,j)
1417    80   continue
1418
1419 #ifdef TIMING
1420         ti = click()
1421 #endif
1422 ! ! Send (1,j)th equations.
1423 ! ! Receive (0,j)th equations.
1424
1425  call mpi_cart_shift(cartGridComm, rowshift, -1, source, dest, ierr)
1426  call MPI_ISEND(zntmp, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, sendReq, ierr)
1427  call MPI_IRECV(   zn, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, recvReq, ierr)
1428  call mpi_wait(sendReq, mpp_status, ierr)
1429  call mpi_wait(recvReq, mpp_status, ierr)
1430
1431 #ifdef TIMING
1432         tf = click()
1433         call add_dt(ct,tf,ti,dt)
1434 #endif
1435
1436         do 90 j = 1, XSPS
1437 ! Forward elimination in (1,j)th equations
1438           fac(j) = 1./(1. - c(1,j)*zn(j,1))
1439 ! Check for singular matrix (debugging only)
1440           r(1,j) = (r(1,j) - c(1,j)*zn(j,2))*fac(j)
1441 ! Backward elimination in (0,j)th equations.
1442           r(0,j) = zn(j,2) - zn(j,1)*r(1,j)
1443 ! Forward elimination in (ZSPS,j)th equations.
1444           r(ZSPS,j) = r(ZSPS,j) - c(ZSPS,j)*r(1,j)
1445 ! Load (ZSPS,j)th equations into passing arrays.
1446           zntmp(j,1) = 0.
1447           zntmp(j,2) = r(ZSPS,j)
1448    90   continue
1449
1450 #ifdef TIMING
1451         ti = click()
1452 #endif
1453 ! ! Send (ZSPS,j)th equations.
1454
1455  call mpi_cart_shift(cartGridComm, rowshift, 1, source, dest, ierr)
1456  call MPI_ISEND(zntmp, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, sendReq, ierr)
1457
1458 #ifdef TIMING
1459         tf = click()
1460         call add_dt(ct,tf,ti,dt)
1461 #endif
1462
1463         do 100 j = 1, XSPS
1464 ! Forward elimination in (ZSPS+1,j)th equations to get
1465 ! r(ZSPS+1,j).
1466         r(ZSPS+1,j) = t2(j) - t1(j)*r(ZSPS,j)
1467         do 110 i = 2, ZSPS-1
1468 ! Completion of forward and backward elimination to get remaining unknowns.
1469           r(i,j) = r(i,j) - c(i,j)*r(1,j) - b(i,j)*r(ZSPS,j)
1470   110   continue
1471   100   continue
1472
1473  call mpi_wait(sendReq, mpp_status, ierr)
1474
1475       else
1476
1477 ! Load (1,j)th equations into passing arrays.
1478         do 120 j = 1, XSPS
1479           zntmp(j,1) = c(1,j)
1480           zntmp(j,2) = r(1,j)
1481   120   continue
1482
1483 #ifdef TIMING
1484         ti = click()
1485 #endif
1486 ! ! Send (1,j)th equations.
1487 ! ! Receive (0,j)th equations.
1488
1489  call mpi_cart_shift(cartGridComm, rowshift, -1, source, dest, ierr)
1490  call MPI_ISEND(zntmp, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, sendReq, ierr)
1491  call MPI_IRECV(   zn, cnt, MPI_REAL, dest, ZN_REC, cartGridComm, recvReq, ierr)
1492  call mpi_wait(sendReq, mpp_status, ierr)
1493  call mpi_wait(recvReq, mpp_status, ierr)
1494
1495 #ifdef TIMING
1496         tf = click()
1497         call add_dt(ct,tf,ti,dt)
1498 #endif
1499
1500         do 130 j = 1, XSPS
1501 ! Forward elimination in (1,j)th equations to get r(1,j).
1502         fac(j) = 1./(1. - c(1,j)*zn(j,1))
1503 ! Check for singular matrix (debugging only)
1504         r(1,j) = (r(1,j) - c(1,j)*zn(j,2))*fac(j)
1505 ! Backward elimination in (0,j)th equations to get remaining unknowns.
1506         r(0,j) = zn(j,2) - zn(j,1)*r(1,j)
1507         do 140 i = 2, ZSPS
1508 ! Completion of forward elimination to get remaining unknowns.
1509           r(i,j) = r(i,j) - c(i,j)*r(1,j)
1510   140   continue
1511   130   continue
1512
1513       endif
1514
1515       return
1516       end subroutine
1517
1518
1519 ! Parallel tridiagonal solver useful for domain decomposed ADI
1520 ! Author(s): Mike Lambert
1521 ! Year: 1996
1522 ! Institution: Lawrence Livermore National Laboratory
1523 ! Publication: Lambert, Rodrigue, and Hewett, "A parallel DSDADI method
1524 !                      for solution of the steady state diffusion equation",
1525 !                      Parallel Computing 23 (1997) 2041-2065
1526 ! Ported to MPI by Benjamin Fersch, Karlsruhe Institute of Technology (2013)
1527
1528       subroutine parxsolv1(c,b,r,ct,pid,x_pid, &
1529                             xsps, zsps, xdns, zdns)
1530
1531       implicit none
1532
1533        integer, intent(in) :: XSPS, &
1534                               ZSPS, &
1535                               XDNS, &
1536                               ZDNS
1537
1538       real, dimension(ZSPS, XSPS), intent(inout) ::  c, &
1539                                                      b
1540
1541
1542       real, dimension(0:ZSPS+1, 0:XSPS+1), intent(inout) ::  r
1543
1544       real, dimension(ZSPS,2) :: xn, xntmp
1545
1546       integer   XN_REC
1547       parameter (XN_REC = 45)
1548
1549       real, dimension(ZSPS)   :: t1, t2, fac
1550       real :: clockdt, click
1551       real :: ct, ti, tf, dt
1552
1553       integer :: pid, x_pid
1554       integer :: i, j, sndr_pid, msg_type, cnt, ackn
1555       integer :: sendReq, recvReq
1556
1557       integer :: source, dest
1558
1559
1560 #ifdef TIMING
1561       dt = clockdt()
1562 #endif
1563
1564       if (x_pid .eq. 1) then
1565
1566 ! Load passing (i,XSPS)th equations into passing arrays.
1567         do 10 i = 1, ZSPS
1568           xntmp(i,1) = b(i,XSPS)
1569           xntmp(i,2) = r(i,XSPS)
1570    10   continue
1571
1572         cnt = 2*ZSPS
1573 #ifdef TIMING
1574         ti = click()
1575 #endif
1576 ! ! Send (i,XSPS)th equations.
1577 ! ! Receive (i,(XSPS + 1))th equations.
1578
1579  call mpi_cart_shift(cartGridComm, colshift, 1, source, dest, ierr)
1580  call MPI_ISEND(xntmp, cnt, MPI_REAL, dest, XN_REC, cartGridComm, sendReq, ierr)
1581  call MPI_IRECV(   xn, cnt, MPI_REAL, dest, XN_REC, cartGridComm, recvReq, ierr)
1582  call mpi_wait(sendReq, mpp_status, ierr)
1583  call mpi_wait(recvReq, mpp_status, ierr)
1584
1585 #ifdef TIMING
1586         tf = click()
1587         call add_dt(ct,tf,ti,dt)
1588 #endif
1589
1590         do 20 i = 1, ZSPS
1591 ! Backward elimination in (i,XSPS)th equations to get
1592 ! r(i,XSPS)
1593           fac(i) = 1./(1. - b(i,XSPS)*xn(i,1))
1594           r(i,XSPS) = (r(i,XSPS)-b(i,XSPS)*xn(i,2))*fac(i)
1595 ! Forward elimination in (i,XSPS+1)th equations to get
1596 ! r(i,XSPS+1)
1597           r(i,XSPS+1) = xn(i,2) - xn(i,1)*r(i,XSPS)
1598    20   continue
1599
1600 ! Completion of backward elimination to get remaining unknowns.
1601         do 30 j = 1, XSPS-1
1602         do 30 i = 1, ZSPS
1603           r(i,j) = r(i,j) - b(i,j)*r(i,XSPS)
1604    30   continue
1605
1606       else if (x_pid .le. XDNS/2) then
1607
1608         cnt = 2*ZSPS
1609 #ifdef TIMING
1610         ti = click()
1611 #endif
1612 ! ! Receive (i,0)th equations.
1613
1614  call mpi_cart_shift(cartGridComm, colshift, -1, source, dest, ierr)
1615  call MPI_IRECV(   xn, cnt, MPI_REAL, dest, XN_REC, cartGridComm, recvReq, ierr)
1616  call mpi_wait(recvReq, mpp_status, ierr)
1617
1618 #ifdef TIMING
1619         tf = click()
1620         call add_dt(ct,tf,ti,dt)
1621 #endif
1622
1623 ! Forward elimination in (i,1)th equations of subdomain.
1624         do 40 i = 1, ZSPS
1625           fac(i) = 1./(1. - c(i,1)*xn(i,1))
1626           b(i,1) = b(i,1)*fac(i)
1627           r(i,1) = (r(i,1) - c(i,1)*xn(i,2))*fac(i)
1628 ! Forward elimination in (i,XSPS)th equations of subdomain.
1629           fac(i) = 1./(1. - c(i,XSPS)*b(i,1))
1630           b(i,XSPS) = b(i,XSPS)*fac(i)
1631           r(i,XSPS)=(r(i,XSPS)-c(i,XSPS)*r(i,1))*fac(i)
1632 ! Store (i,0)th equations for later recovery of r(i,0).
1633           t1(i) = xn(i,1)
1634           t2(i) = xn(i,2)
1635 ! Load (i,XSPS)th equations into passing arrays.
1636           xntmp(i,1) = b(i,XSPS)
1637           xntmp(i,2) = r(i,XSPS)
1638    40   continue
1639
1640         cnt = 2*ZSPS
1641 #ifdef TIMING
1642         ti = click()
1643 #endif
1644 ! ! Send (i,XSPS)th equations.
1645 ! ! Receive (i,(XSPS + 1))th equations.
1646
1647  call mpi_cart_shift(cartGridComm, colshift, 1, source, dest, ierr)
1648  call MPI_ISEND(xntmp, cnt, MPI_REAL, dest, XN_REC, cartGridComm, sendReq, ierr)
1649  call MPI_IRECV(   xn, cnt, MPI_REAL, dest, XN_REC, cartGridComm, recvReq, ierr)
1650  call mpi_wait(sendReq, mpp_status, ierr)
1651  call mpi_wait(recvReq, mpp_status, ierr)
1652 #ifdef TIMING
1653         tf = click()
1654         call add_dt(ct,tf,ti,dt)
1655 #endif
1656
1657         do 50 i = 1, ZSPS
1658 ! Backward elimination in (i,XSPS)th equations.
1659           fac(i) = 1./(1. - b(i,XSPS)*xn(i,1))
1660           r(i,XSPS) = (r(i,XSPS) - b(i,XSPS)*xn(i,2))*fac(i)
1661 ! Backward elimination in (i,XSPS+1)th equations.
1662           r(i,XSPS+1) = xn(i,2) - xn(i,1)*r(i,XSPS)
1663 ! Backward elimination in (i,1)th equations to get r(i,1).
1664           r(i,1) = r(i,1) - b(i,1)*r(i,XSPS)
1665 ! Load (i,1)th equations into passing array.
1666           xntmp(i,1) = 0.
1667           xntmp(i,2) = r(i,1)
1668    50   continue
1669
1670         cnt = 2*ZSPS
1671 #ifdef TIMING
1672         ti = click()
1673 #endif
1674 ! ! Send (i,1)th equations.
1675
1676 #ifdef TIMING
1677         tf = click()
1678         call add_dt(ct,tf,ti,dt)
1679 #endif
1680  call mpi_cart_shift(cartGridComm, colshift, -1, source, dest, ierr)
1681  call MPI_ISEND(xntmp, cnt, MPI_REAL, dest, XN_REC, cartGridComm, sendReq, ierr)
1682
1683         do 60 i = 1, ZSPS
1684 ! Backward elimination in (i,0)th equations.
1685           r(i,0) = t2(i) - t1(i)*r(i,1)
1686    60   continue
1687
1688 ! Completion of forward and backward elimination for solution of
1689 ! unknowns.
1690         do 70 j = 2, XSPS-1
1691         do 70 i = 1, ZSPS
1692           r(i,j) = r(i,j) - b(i,j)*r(i,XSPS) - c(i,j)*r(i,1)
1693    70   continue
1694
1695  call mpi_wait(sendReq, mpp_status, ierr)
1696
1697       else if (x_pid .lt. XDNS) then
1698
1699         cnt = 2*ZSPS
1700 #ifdef TIMING
1701         ti = click()
1702 #endif
1703 ! ! Receive (i,XSPS+1)th equations.
1704
1705  call mpi_cart_shift(cartGridComm, colshift, 1, source, dest, ierr)
1706  call MPI_IRECV(   xn, cnt, MPI_REAL, dest, XN_REC, cartGridComm, recvReq, ierr)
1707  call mpi_wait(recvReq, mpp_status, ierr)
1708
1709 #ifdef TIMING
1710         tf = click()
1711         call add_dt(ct,tf,ti,dt)
1712 #endif
1713
1714         do 80 i = 1, ZSPS
1715 ! Backward elimination in (i,XSPS)th equations.
1716           fac(i) = 1./(1. - b(i,XSPS)*xn(i,1))
1717           c(i,XSPS) = c(i,XSPS)*fac(i)
1718           r(i,XSPS) = (r(i,XSPS) - b(i,XSPS)*xn(i,2))*fac(i)
1719 ! Backward elimination in (i,1)th equations.
1720           fac(i) = 1./(1. - b(i,1)*c(i,XSPS))
1721           c(i,1) = c(i,1)*fac(i)
1722           r(i,1) = (r(i,1) - b(i,1)*r(i,XSPS))*fac(i)
1723 ! Store (i,XSPS+1)th equations for later recovery of r(i,XSPS+1).
1724           t1(i) = xn(i,1)
1725           t2(i) = xn(i,2)
1726 ! Load passing arrays with (i,1)th equations.
1727           xntmp(i,1) = c(i,1)
1728           xntmp(i,2) = r(i,1)
1729    80   continue
1730
1731         cnt = 2*ZSPS
1732 #ifdef TIMING
1733         ti = click()
1734 #endif
1735 ! ! Send (i,1)th equations.
1736 ! ! Receive (i,0)th equations.
1737  call mpi_cart_shift(cartGridComm, colshift, -1, source, dest, ierr)
1738  call MPI_ISEND(xntmp, cnt, MPI_REAL, dest, XN_REC, cartGridComm, sendReq, ierr)
1739  call MPI_IRECV(   xn, cnt, MPI_REAL, dest, XN_REC, cartGridComm, recvReq, ierr)
1740  call mpi_wait(sendReq, mpp_status, ierr)
1741  call mpi_wait(recvReq, mpp_status, ierr)
1742
1743 #ifdef TIMING
1744         tf = click()
1745         call add_dt(ct,tf,ti,dt)
1746 #endif
1747
1748         do 90 i = 1, ZSPS
1749 ! Forward elimination in (i,1)th equations
1750           fac(i) = 1./(1. - c(i,1)*xn(i,1))
1751           r(i,1) = (r(i,1) - c(i,1)*xn(i,2))*fac(i)
1752 ! Backward elimination in (i,0)th equations.
1753           r(i,0) = xn(i,2) - xn(i,1)*r(i,1)
1754 ! Forward elimination in (i,XSPS)th equations.
1755           r(i,XSPS) = r(i,XSPS) - c(i,XSPS)*r(i,1)
1756 ! Load (i,XSPS)th equations into passing arrays.
1757           xntmp(i,1) = 0.
1758           xntmp(i,2) = r(i,XSPS)
1759    90   continue
1760
1761         cnt = 2*ZSPS
1762 #ifdef TIMING
1763         ti = click()
1764 #endif
1765 ! ! Send (i,XSPS)th equations.
1766
1767  call mpi_cart_shift(cartGridComm, colshift, 1, source, dest, ierr)
1768  call MPI_ISEND(xntmp, cnt, MPI_REAL, dest, XN_REC, cartGridComm, sendReq, ierr)
1769 #ifdef TIMING
1770         tf = click()
1771         call add_dt(ct,tf,ti,dt)
1772 #endif
1773
1774 ! Forward elimination in (i,XSPS)th equations to get
1775 ! r(i,XSPS+1).
1776         do 100 i = 1, ZSPS
1777           r(i,XSPS+1) = t2(i) - t1(i)*r(i,XSPS)
1778   100   continue
1779
1780 ! Completion of forward and backward elimination to get remaining unknowns.
1781         do 110 j = 2, XSPS-1
1782         do 110 i = 1, ZSPS
1783           r(i,j) = r(i,j) - c(i,j)*r(i,1) - b(i,j)*r(i,XSPS)
1784   110   continue
1785
1786  call mpi_wait(sendReq, mpp_status, ierr)
1787
1788       else
1789
1790 ! Load (i,1)th equations into passing arrays.
1791         do 120 i = 1, ZSPS
1792           xntmp(i,1) = c(i,1)
1793           xntmp(i,2) = r(i,1)
1794   120   continue
1795
1796         cnt = 2*ZSPS
1797 #ifdef TIMING
1798         ti = click()
1799 #endif
1800 ! ! Send (i,1)th equations.
1801 ! ! Receive (i,0)th equations.
1802
1803  call mpi_cart_shift(cartGridComm, colshift, -1, source, dest, ierr)
1804  call MPI_ISEND(xntmp, cnt, MPI_REAL, dest, XN_REC, cartGridComm, sendReq, ierr)
1805  call MPI_IRECV(   xn, cnt, MPI_REAL, dest, XN_REC, cartGridComm, recvReq, ierr)
1806  call mpi_wait(sendReq, mpp_status, ierr)
1807  call mpi_wait(recvReq, mpp_status, ierr)
1808
1809 #ifdef TIMING
1810         tf = click()
1811         call add_dt(ct,tf,ti,dt)
1812 #endif
1813
1814         do 130 i = 1, ZSPS
1815 ! Forward elimination in (i,1)th equations to get r(i,1).
1816           fac(i) = 1./(1. - c(i,1)*xn(i,1))
1817           r(i,1) = (r(i,1) - c(i,1)*xn(i,2))*fac(i)
1818 ! Backward elimination in (i,0)th equations to get r(i,0).
1819           r(i,0) = xn(i,2) - xn(i,1)*r(i,1)
1820   130   continue
1821
1822 ! Completion of forward elimination to get remaining unknowns.
1823         do 140 j = 2, XSPS
1824         do 140 i = 1, ZSPS
1825           r(i,j) = r(i,j) - c(i,j)*r(i,1)
1826   140   continue
1827
1828       endif
1829
1830       return
1831       end subroutine
1832
1833
1834 ! Parallel tridiagonal solver useful for domain decomposed ADI
1835 ! Author(s): Mike Lambert
1836 ! Year: 1996
1837 ! Institution: Lawrence Livermore National Laboratory
1838 ! Publication: Lambert, Rodrigue, and Hewett, "A parallel DSDADI method
1839 !                      for solution of the steady state diffusion equation",
1840 !                      Parallel Computing 23 (1997) 2041-2065
1841
1842       subroutine sub_n_form(n_xs,n_zs,c,a,b,r,c2,b2,r2,wk,xfac,zfac, &
1843                             dpid,dn_subs,dir)
1844
1845       implicit none
1846
1847       integer n_xs,n_zs
1848
1849 !  c(,)  -- subdiagonal elements of tridiagonal systems
1850 !  a(,)  -- diagonal elements of tridiagonal systems
1851 !  b(,)  -- superdiagonal elements of tridiagonal systems
1852 !  r(,)  -- right-hand side elements of tridiagonal systems
1853 !  c2(,) -- front-leg elements of N-systems
1854 !  b2(,) -- back-leg elements of N-systems
1855 !  r2(,) -- right-hand side elements of N-systems
1856 !  wk(,) -- work array with same dimensions as a, b, c, etc.
1857       real c(n_zs,n_xs)
1858       real a(n_zs,n_xs)
1859       real b(n_zs,n_xs)
1860       real r(n_zs,n_xs)
1861       real c2(n_zs,n_xs)
1862       real b2(n_zs,n_xs)
1863       real r2(0:n_zs+1,0:n_xs+1)
1864       real wk(n_zs,n_xs)
1865       real fac
1866       real xfac(n_zs)
1867       real zfac(n_xs)
1868
1869       integer dpid,dn_subs,dir
1870       integer i, j, XDIR, ZDIR
1871       parameter (XDIR = 1, ZDIR = 2)
1872
1873       if (dir .eq. XDIR) then
1874
1875 ! Forward elimination of subdiagonal elements
1876         if (dpid .eq. 1) then
1877
1878           do 10 i = 1, n_zs
1879             xfac(i) = 1./a(i,1)
1880             c2(i,1) = 0.
1881             r2(i,1) = r(i,1)*xfac(i)
1882    10     continue
1883
1884           do 20 j = 2, n_xs
1885           do 20 i = 1, n_zs
1886             wk(i,j-1) = b(i,j-1)*xfac(i)
1887             xfac(i) = 1./(a(i,j) - c(i,j)*wk(i,j-1))
1888             c2(i,j) = 0.
1889             r2(i,j) = (r(i,j) - c(i,j)*r2(i,j-1))*xfac(i)
1890    20     continue
1891
1892           do 40 i = 1, n_zs
1893             b2(i,n_xs) = b(i,n_xs)*xfac(i)
1894    40     continue
1895
1896         else
1897
1898           do 50 i = 1, n_zs
1899             xfac(i) = 1./a(i,1)
1900             c2(i,1) = c(i,1)*xfac(i)
1901             wk(i,1) = b(i,1)*xfac(i)
1902             r2(i,1) = r(i,1)*xfac(i)
1903             xfac(i) = 1./a(i,2)
1904             c2(i,2) = c(i,2)*xfac(i)
1905             r2(i,2) = r(i,2)*xfac(i)
1906    50     continue
1907
1908           do 60 j = 3, n_xs
1909           do 60 i = 1, n_zs
1910             wk(i,j-1) = b(i,j-1)*xfac(i)
1911             xfac(i) = 1./(a(i,j) - c(i,j)*wk(i,j-1))
1912             c2(i,j) = -c(i,j)*c2(i,j-1)*xfac(i)
1913             r2(i,j) = (r(i,j) - c(i,j)*r2(i,j-1))*xfac(i)
1914    60     continue
1915
1916           do 80 i = 1, n_zs
1917             b2(i,n_xs) = b(i,n_xs)*xfac(i)
1918    80     continue
1919
1920         endif
1921
1922 ! Backward elimination of superdiagonal elements
1923         if (dpid .eq. dn_subs) then
1924
1925           do 90 j = n_xs-1, 2, -1
1926           do 90 i = 1, n_zs
1927             c2(i,j) = c2(i,j) - wk(i,j)*c2(i,j+1)
1928             b2(i,j) = 0.
1929             r2(i,j) = r2(i,j) - wk(i,j)*r2(i,j+1)
1930    90     continue
1931
1932           do 100 i = 1, n_zs
1933             fac = 1./(1. - wk(i,1)*c2(i,2))
1934             c2(i,1) = c2(i,1)*fac
1935             b2(i,1) = 0.
1936             r2(i,1) = (r2(i,1) - wk(i,1)*r2(i,2))*fac
1937   100     continue
1938
1939         else
1940
1941           do 110 i = 1, n_zs
1942             b2(i,n_xs-1) = wk(i,n_xs-1)
1943   110     continue
1944
1945           do 120 j = n_xs-2, 2, -1
1946           do 120 i = 1, n_zs
1947             c2(i,j) = c2(i,j) - wk(i,j)*c2(i,j+1)
1948             b2(i,j) = -wk(i,j)*b2(i,j+1)
1949             r2(i,j) = r2(i,j) - wk(i,j)*r2(i,j+1)
1950   120     continue
1951
1952 ! If only 2 points in X-direction, do not execute these statements.
1953           if (n_xs .gt. 2) then
1954             do 130 i = 1, n_zs
1955               fac = 1./(1. - wk(i,1)*c2(i,2))
1956               c2(i,1) = c2(i,1)*fac
1957               r2(i,1) = (r2(i,1) - wk(i,1)*r2(i,2))*fac
1958               b2(i,1) = -wk(i,1)*b2(i,2)*fac
1959   130       continue
1960           endif
1961
1962         endif
1963
1964       else if (dir .eq. ZDIR) then
1965
1966 ! Forward elimination of subdiagonal elements
1967         if (dpid .eq. 1) then
1968
1969           do 140 j = 1, n_xs
1970             zfac(j) = 1./a(1,j)
1971             c2(1,j) = 0.
1972             r2(1,j) = r(1,j)*zfac(j)
1973   140     continue
1974
1975           do 150 i = 2, n_zs
1976           do 150 j = 1, n_xs
1977             wk(i-1,j) = b(i-1,j)*zfac(j)
1978             zfac(j) = 1./(a(i,j) - c(i,j)*wk(i-1,j))
1979             c2(i,j) = 0.
1980             r2(i,j) = (r(i,j) - c(i,j)*r2(i-1,j))*zfac(j)
1981   150     continue
1982
1983           do 170 j = 1, n_xs
1984             b2(n_zs,j) = b(n_zs,j)*zfac(j)
1985   170     continue
1986
1987         else
1988
1989           do 180 j = 1, n_xs
1990             zfac(j) = 1./a(1,j)
1991             c2(1,j) = c(1,j)*zfac(j)
1992             wk(1,j) = b(1,j)*zfac(j)
1993             r2(1,j) = r(1,j)*zfac(j)
1994             zfac(j) = 1./a(2,j)
1995             c2(2,j) = c(2,j)*zfac(j)
1996             r2(2,j) = r(2,j)*zfac(j)
1997   180     continue
1998
1999           do 190 i = 3, n_zs
2000           do 190 j = 1, n_xs
2001             wk(i-1,j) = b(i-1,j)*zfac(j)
2002             zfac(j) = 1./(a(i,j) - c(i,j)*wk(i-1,j))
2003             c2(i,j) = -c(i,j)*c2(i-1,j)*zfac(j)
2004             r2(i,j) = (r(i,j) - c(i,j)*r2(i-1,j))*zfac(j)
2005   190     continue
2006
2007           do 210 j = 1, n_xs
2008             b2(n_zs,j) = b(n_zs,j)*zfac(j)
2009   210     continue
2010
2011         endif
2012
2013 ! Backward elimination of superdiagonal elements
2014         if (dpid .eq. dn_subs) then
2015
2016           do 220 j = 1, n_xs
2017           do 220 i = n_zs - 1, 2, -1
2018             c2(i,j) = c2(i,j) - wk(i,j)*c2(i+1,j)
2019             b2(i,j) = 0.
2020             r2(i,j) = r2(i,j) - wk(i,j)*r2(i+1,j)
2021   220     continue
2022
2023           do 230 j = 1, n_xs
2024             fac = 1./(1. - wk(1,j)*c2(2,j))
2025             c2(1,j) = c2(1,j)*fac
2026             b2(1,j) = 0.
2027             r2(1,j) = (r2(1,j) - wk(1,j)*r2(2,j))*fac
2028   230     continue
2029
2030         else
2031
2032           do 240 j = 1, n_xs
2033             b2(n_zs-1,j) = wk(n_zs-1,j)
2034   240     continue
2035
2036           do 250 j = 1, n_xs
2037           do 250 i = n_zs - 2, 2, -1
2038             c2(i,j) = c2(i,j) - wk(i,j)*c2(i+1,j)
2039             b2(i,j) = -wk(i,j)*b2(i+1,j)
2040             r2(i,j)  = r2(i,j) - wk(i,j)*r2(i+1,j)
2041   250     continue
2042
2043 ! If only 2 points in Z-direction, do not execute these statements.
2044           if (n_zs .gt. 2) then
2045             do 260 j = 1, n_xs
2046               fac = 1./(1. - wk(1,j)*c2(2,j))
2047               c2(1,j) = c2(1,j)*fac
2048               r2(1,j) = (r2(1,j) - wk(1,j)*r2(2,j))*fac
2049               b2(1,j) = -wk(1,j)*b2(2,j)*fac
2050   260       continue
2051           endif
2052
2053         endif
2054
2055 ! Announce bad direction specifier (debugging only)
2056 !     else
2057 !       write(*,*) 'sub_n_form:  What direction?'
2058 !       stop
2059       endif
2060
2061       return
2062       end subroutine
2063 #endif
2064
2065 ! Tridiagonal solver useful for domain decomposed ADI
2066 ! Author(s): Mike Lambert
2067 ! Year: 1996
2068 ! Institution: Lawrence Livermore National Laboratory
2069 ! Publication: Lambert, Rodrigue, and Hewett, "A parallel DSDADI method
2070 !                      for solution of the steady state diffusion equation",
2071 !                      Parallel Computing 23 (1997) 2041-2065
2072
2073       subroutine sub_tri_solv(n_xs,n_zs,c,a,b,r,x,wk,xfac,zfac,dir)
2074
2075       implicit none
2076
2077       integer n_xs,n_zs
2078
2079 !  c(,)  -- subdiagonal elements of tridiagonal systems
2080 !  a(,)  -- diagonal elements of tridiagonal systems
2081 !  b(,)  -- superdiagonal elements of tridiagonal systems
2082 !  r(,)  -- right-hand side elements of tridiagonal systems
2083 !  x(,)  -- solutions
2084 !  wk(,) -- work array w/ same dimensions as c, a, b, etc.
2085
2086       real c(n_zs,n_xs)
2087       real a(n_zs,n_xs)
2088       real b(n_zs,n_xs)
2089       real r(n_zs,n_xs)
2090       real x(0:n_zs+1,0:n_xs+1)
2091       real wk(n_zs,n_xs)
2092       real xfac(n_zs)
2093       real zfac(n_xs)
2094
2095       integer dir
2096       integer i,j,XDIR,ZDIR
2097
2098       parameter (XDIR = 1, ZDIR = 2)
2099
2100       if (dir .eq. XDIR) then
2101
2102         do 10 i = 1, n_zs
2103 ! Check for need to pivot (debugging only)
2104         xfac(i) = 1./a(i,1)
2105         x(i,1)  = r(i,1)*xfac(i)
2106    10   continue
2107
2108 ! Forward subdiagonal elimination
2109         do 20 j = 2, n_xs
2110         do 20 i = 1, n_zs
2111         wk(i,j-1) = b(i,j-1)*xfac(i)
2112         xfac(i) = 1./(a(i,j) - c(i,j)*wk(i,j-1))
2113 ! Check for need to pivot (debugging only)
2114         x(i,j) = (r(i,j) - c(i,j)*x(i,j-1))*xfac(i)
2115    20   continue
2116
2117 ! Backsubstitution
2118         do 30 j = n_xs - 1, 1, -1
2119         do 30 i = 1, n_zs
2120         x(i,j)  = x(i,j) - wk(i,j)*x(i,j+1)
2121    30   continue
2122
2123
2124       else if (dir .eq. ZDIR) then
2125
2126        do j = 1, n_xs
2127 ! Check for need to pivot (debugging only)
2128         zfac(j) = 1./a(1,j)
2129         x(1,j)  = r(1,j)*zfac(j)
2130        end do
2131
2132 ! Forward subdiagonal elimination
2133       do j = 1, n_xs
2134        do i = 2, n_zs
2135         wk(i-1,j) = b(i-1,j)*zfac(j)
2136         zfac(j) = 1./(a(i,j) - c(i,j)*wk(i-1,j))
2137 ! Check for need to pivot (debugging only)
2138         x(i,j) = (r(i,j) - c(i,j)*x(i-1,j))*zfac(j)
2139        end do
2140       end do
2141
2142 ! Backsubstitution
2143       do j = 1, n_xs
2144        do i = n_zs - 1, 1, -1
2145         x(i,j)  =  x(i,j) - wk(i,j)*x(i+1,j)
2146        end do
2147       end do
2148
2149 ! Announce bad direction specifier (debugging only)
2150 !     else
2151 !       write(*,*) 'sub_tri_solv:  What direction?'
2152 !       stop
2153       endif
2154
2155       return
2156       end  subroutine
2157
2158
2159 end module module_gw_gw2d