hydro/Routing/module_gw_gw2d.F90

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