hydro/Routing/Noah_distr_routing_subsurface.F

   1 subroutine subsurfaceRouting ( subrt_data, subrt_static, subrt_input, subrt_output)
   2 #ifdef MPP_LAND
   3     use module_mpp_land, only:  mpp_land_com_real, mpp_land_com_integer
   4     use module_subsurface_data
   5     use module_subsurface_static_data
   6     use module_subsurface_input
   7     use module_subsurface_output
   8 #endif
   9     implicit none
  10     type (subsurface_struct), intent(inout) :: subrt_data
  11     type (subsurface_static_interface), intent(inout) :: subrt_static
  12     type (subsurface_input_interface), intent(inout) :: subrt_input
  13     type (subsurface_output_interface), intent(inout) :: subrt_output
  14     !integer, INTENT(IN) :: ixrt, jxrt , nsoil, rt_option
  15     !REAL, INTENT(IN)                          :: DT
  16     !real,INTENT(IN), DIMENSION(NSOIL)      :: SLDPTH
  17     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT) :: sub_resid
  18     !real,INTENT(INOUT), DIMENSION(subrt_static%IXRT,subrt_static%JXRT)::INFXSUBRT
  19     !real,INTENT(INOUT) :: QSUBBDRYTRT
  20     !REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: QSUBBDRYRT, QSUBRT
  21     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,NSOIL) :: SMCRT, SMCWLTRT, SMCMAXRT, SMCREFRT
  22
  23
  24     !INTEGER :: SO8RT_D(IXRT,JXRT,3)
  25     !REAL :: SO8RT(IXRT,JXRT,8)
  26     !REAL, INTENT(IN)                          :: dist(ixrt,jxrt,9)
  27     !  -----local array ----------
  28     !REAL, DIMENSION(IXRT,JXRT)   :: ZWATTABLRT
  29     REAL, DIMENSION(subrt_static%ixrt,subrt_static%jxrt)   :: CWATAVAIL
  30     INTEGER, DIMENSION(subrt_static%ixrt,subrt_static%jxrt) :: SATLYRCHK
  31
  32     CWATAVAIL = 0.
  33     CALL FINDZWAT(subrt_static%ixrt,subrt_static%jxrt,subrt_static%NSOIL, &
  34         subrt_data%grid_transform%smcrt, &
  35         subrt_data%grid_transform%smcmaxrt, &
  36         subrt_data%grid_transform%smcrefrt, &
  37         subrt_data%grid_transform%smcwltrt, subrt_data%properties%zsoil,SATLYRCHK,subrt_data%properties%zwattablrt, &
  38         CWATAVAIL,subrt_data%properties%sldpth)
  39 #ifdef MPP_LAND
  40     call MPP_LAND_COM_REAL(subrt_data%properties%zwattablrt,subrt_static%ixrt,subrt_static%jxrt,99)
  41     call MPP_LAND_COM_REAL(CWATAVAIL,subrt_static%ixrt,subrt_static%jxrt,99)
  42     call MPP_LAND_COM_INTEGER(SATLYRCHK,subrt_static%ixrt,subrt_static%jxrt,99)
  43 #endif
  44
  45
  46     !DJG Second, Call subsurface routing routine...
  47 #ifdef HYDRO_D
  48     print *, "Beginning SUB_routing..."
  49     print *, "Routing method is ",subrt_static%rt_option, " direction."
  50 #endif
  51
  52     !!!! Find saturated layer depth...
  53     ! Loop through domain to determine sat. layers and assign wat tbl depth...
  54     !    and water available for subsfc routing (CWATAVAIL)...
  55     ! This subroutine returns: ZWATTABLRT, CWATAVAIL and SATLYRCHK
  56
  57
  58     CALL SUBSFC_RTNG( subrt_data, subrt_static, subrt_input, subrt_output, CWATAVAIL, SATLYRCHK)
  59
  60 #ifdef HYDRO_D
  61     print *, "SUBROUTE routing called and returned..."
  62 #endif
  63
  64 end subroutine subsurfaceRouting
  65
  66 !DJG ------------------------------------------------
  67 !DJG   SUBROUTINE SUBSFC_RTNG
  68 !DJG ------------------------------------------------
  69
  70 SUBROUTINE SUBSFC_RTNG(subrt_data, subrt_static, subrt_input, subrt_output, CWATAVAIL, SATLYRCHK)
  71
  72     !       use module_mpp_land, only: write_restart_rt_3, write_restart_rt_2, &
  73         !            my_id
  74 #ifdef MPP_LAND
  75     use module_mpp_land, only: MPP_LAND_COM_REAL, sum_real1, &
  76         my_id, io_id, numprocs
  77     use module_subsurface_data
  78     use module_subsurface_static_data
  79     use module_subsurface_input
  80     use module_subsurface_output
  81     use module_hydro_stop, only:HYDRO_stop
  82 #endif
  83     IMPLICIT NONE
  84
  85     !DJG -------- DECLARATIONS ------------------------
  86
  87     type (subsurface_struct), intent(inout) :: subrt_data
  88     type (subsurface_static_interface), intent(inout) :: subrt_static
  89     type (subsurface_input_interface), intent(inout) :: subrt_input
  90     type (subsurface_output_interface), intent(inout) :: subrt_output
  91
  92     !INTEGER, INTENT(IN) :: IXRT,JXRT,NSOIL
  93     !INTEGER, INTENT(IN) :: IXRT,JXRT
  94
  95     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: SOXRT
  96     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: SOYRT
  97     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: LATKSATRT
  98     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: SOLDEPRT
  99
 100     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)   :: ZWATTABLRT
 101     REAL, INTENT(IN), DIMENSION(subrt_static%ixrt,subrt_static%jxrt)   :: CWATAVAIL
 102     INTEGER, INTENT(IN), DIMENSION(subrt_static%ixrt,subrt_static%jxrt) :: SATLYRCHK
 103
 104
 105     !REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: QSUBRT
 106     !REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: QSUBBDRYRT
 107
 108     !REAL, INTENT(IN)                          :: dist(ixrt,jxrt,9)
 109     !REAL, INTENT(IN)                          :: DT
 110     !REAL, INTENT(IN), DIMENSION(subrt_static%nsoil)        :: ZSOIL
 111     !REAL, INTENT(IN), DIMENSION(subrt_static%nsoil)      :: SLDPTH
 112     !REAL, INTENT(INOUT)                       :: QSUBBDRYTRT
 113     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: INFXSUBRT
 114     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,subrt_static%nsoil) :: SMCMAXRT
 115     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,subrt_static%nsoil) :: SMCREFRT
 116     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,subrt_static%nsoil) :: SMCRT
 117     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,subrt_static%nsoil) :: SMCWLTRT
 118
 119     REAL, DIMENSION(subrt_static%IXRT,subrt_static%JXRT)        :: ywtmp
 120     !DJG Local Variables
 121
 122     INTEGER     :: I,J,KK
 123     !djg        INTEGER, DIMENSION(IXRT,JXRT) :: SATLYRCHK
 124
 125     REAL        :: GRDAREA
 126     REAL        :: SUBFLO
 127     REAL        :: WATAVAIL
 128
 129     !INTEGER :: SO8RT_D(IXRT,JXRT,3)
 130     !REAL :: SO8RT(IXRT,JXRT,8)
 131     !INTEGER, INTENT(IN)  ::  rt_option
 132     integer  ::  index
 133
 134     INTEGER :: DT_STEPS             !-- number of timestep in routing
 135     REAL :: SUBDT                !-- subsurface routing timestep
 136     INTEGER :: KRT                  !-- routing counter
 137     REAL, DIMENSION(subrt_static%IXRT,subrt_static%JXRT,subrt_static%nsoil) :: SMCTMP  !--temp store of SMC
 138     REAL, DIMENSION(subrt_static%IXRT,subrt_static%JXRT) :: ZWATTABLRTTMP ! temp store of ZWAT
 139     REAL, DIMENSION(subrt_static%IXRT,subrt_static%JXRT) :: INFXSUBRTTMP ! temp store of infilx
 140     !djg        REAL, DIMENSION(IXRT,JXRT) :: CWATAVAIL ! temp specif. of wat avial
 141
 142
 143
 144     !DJG Debug Variables...
 145     REAL :: qsubchk,qsubbdrytmp
 146     REAL :: junk1,junk2,junk3,junk5,junk6,junk7
 147     INTEGER, PARAMETER :: double=8
 148     REAL (KIND=double) :: smctot1a,smctot2a
 149     INTEGER :: kx,count
 150
 151 #ifdef HYDRO_D
 152     ! ADCHANGE: Water balance variables
 153     real   :: smctot1,smctot2
 154     real   :: suminfxsrt1,suminfxsrt2
 155     real   :: qbdry1,qbdry2
 156     real   :: sumqsubrt1, sumqsubrt2
 157 #endif
 158
 159     !DJG -----------------------------------------------------------------
 160     !DJG  SUBSURFACE ROUTING LOOP
 161     !DJG    - SUBSURFACE ROUTING RUN ON NOAH TIMESTEP
 162     !DJG    - SUBSURFACE ROUITNG ONLY PERFORMED ON SATURATED LAYERS
 163     !DJG -----------------------------------------------------------------
 164
 165     !Z. Cui -----------------------------------------------------------------
 166     !Move SUBDT initialization up to avoid crash when HYDRO_D is defined
 167     SUBDT = subrt_static%dt                !-- initialize the routing timestep to DT
 168     !Z. Cui -----------------------------------------------------------------
 169 #ifdef HYDRO_D
 170     ! ADCHANGE: START Initial water balance variables
 171     ! ALL VARS in MM
 172     suminfxsrt1 = 0.
 173     qbdry1 = 0.
 174     smctot1 = 0.
 175     sumqsubrt1 = 0.
 176     do i=1,subrt_static%IXRT
 177         do j=1,subrt_static%JXRT
 178             suminfxsrt1 = suminfxsrt1 + subrt_input%infiltration_excess(I,J) / float(subrt_static%IXRT*subrt_static%JXRT)
 179             qbdry1 = qbdry1 + subrt_data%state%qsubbdryrt(i,j) / subrt_data%properties%distance_to_neighbor(i,j,9)*SUBDT*1000. / float(subrt_static%IXRT*subrt_static%JXRT)
 180             sumqsubrt1 = sumqsubrt1 + subrt_data%state%qsubrt(i,j) / subrt_data%properties%distance_to_neighbor(i,j,9)*SUBDT*1000. / float(subrt_static%IXRT*subrt_static%JXRT)
 181             do kk=1,subrt_static%nsoil
 182                 smctot1 = smctot1 + subrt_data%grid_transform%smcrt(I,J,KK)*subrt_data%properties%sldpth(KK)*1000. / float(subrt_static%IXRT*subrt_static%JXRT)
 183             end do
 184         end do
 185     end do
 186
 187 #ifdef MPP_LAND
 188     ! not tested
 189     CALL sum_real1(suminfxsrt1)
 190     CALL sum_real1(qbdry1)
 191     CALL sum_real1(sumqsubrt1)
 192     CALL sum_real1(smctot1)
 193     suminfxsrt1 = suminfxsrt1/float(numprocs)
 194     qbdry1 = qbdry1/float(numprocs)
 195     sumqsubrt1 = sumqsubrt1/float(numprocs)
 196     smctot1 = smctot1/float(numprocs)
 197 #endif
 198     ! END Initial water balance variables
 199 #endif
 200
 201
 202     !yw GRDAREA=DXRT*DXRT
 203     ! GRDAREA=subrt_data%properties%distance_to_neighbor(i,j,9)
 204
 205
 206     !DJG debug subsfc...
 207     subflo = 0.0
 208
 209     !DJG Set up mass balance checks...
 210     !         CWATAVAIL = 0.            !-- initialize subsurface watavail
 211
 212
 213     !!!! Find saturated layer depth...
 214     ! Loop through domain to determine sat. layers and assign wat tbl depth...
 215     !    and water available for subsfc routing (CWATAVAIL)...
 216     !
 217     !   This function call appears to have been moved to SubsurfaceRouting()
 218     !
 219     !
 220     !         CALL FINDZWAT(IXRT,JXRT,subrt_static%nsoil,SMCRT,SMCMAXRT,SMCREFRT, &
 221         !                             SMCWLTRT,ZSOIL,SATLYRCHK,subrt_data%properties%zwattablrt, &
 222         !                             CWATAVAIL,SLDPTH)
 223
 224
 225
 226
 227     !DJG debug variable...
 228
 229     !DJG Courant check temp variable setup...
 230     ZWATTABLRTTMP = subrt_data%properties%zwattablrt !-- temporary storage of water table level
 231
 232
 233
 234
 235     !!!! Call subsurface routing subroutine...
 236 #ifdef HYDRO_D
 237     print *, "calling subsurface routing subroutine...Opt. ",subrt_static%rt_option
 238 #endif
 239
 240     ! ADCHANGE: IMPORTANT!
 241     ! 2D subsurface option currently has bug so forcing to option 1 in this routine to
 242     ! allow users to still have option to use 2d overland (both are controlled by same
 243     ! rt_option flag). Remove this hard-coded option when rt_option=2 is fixed for subsurface.
 244     !     if(subrt_static%rt_option .eq. 1) then
 245     CALL ROUTE_SUBSURFACE1(subrt_data%properties%distance_to_neighbor, &
 246                            subrt_data%properties%zwattablrt, &
 247                            subrt_data%properties%lksatrt, subrt_data%properties%nexprt, &
 248                            subrt_data%properties%soldeprt, &
 249                            subrt_static%IXRT, subrt_static%JXRT, &
 250                            subrt_data%properties%surface_slope, &
 251                            subrt_data%properties%max_surface_slope_index, &
 252                            CWATAVAIL, SUBDT, &
 253                            subrt_data%state%qsubbdryrt, subrt_data%state%qsubbdrytrt, &
 254                            subrt_data%state%qsubrt)
 255     !     else
 256     !        CALL ROUTE_SUBSURFACE2(subrt_data%properties%distance_to_neighbor,subrt_data%properties%zwattablrt, &
 257         !               subrt_data%state%qsubrt, subrt_data%properties%surface_slope_x, subrt_data%properties%surface_slope_y, &
 258         !               subrt_data%properties%lksatrt, subrt_data%properties%soldeprt, subrt_static%IXRT, subrt_static%JXRT, subrt_data%state%qsubbdryrt, subrt_data%state%qsubbdrytrt, &
 259         !               CWATAVAIL,SUBDT)
 260     !     end if
 261
 262 #ifdef HYDRO_D
 263     write(6,*) "finish calling ROUTE_SUBSURFACE ", subrt_static%rt_option
 264 #endif
 265
 266
 267     !!!! Update soil moisture fields with subsurface flow...
 268
 269     !!!! Loop through subsurface routing domain...
 270     DO I=1,subrt_static%IXRT
 271         DO J=1,subrt_static%JXRT
 272
 273             !!DJG Check for courant condition violation...put limit on qsub
 274             !!DJG QSUB HAS units of m^3/s SUBFLO has units of m
 275
 276             ! ADCHANGE: Moved this constraint to the ROUTE_SUBSURFACE routines
 277             !IF (CWATAVAIL(i,j).le.ABS(qsubrt(i,j))/subrt_data%properties%distance_to_neighbor(i,j,9)*SUBDT) THEN
 278             !  subrt_data%state%qsubrt(i,j) = -1.0*CWATAVAIL(i,j)
 279             !  SUBFLO = subrt_data%state%qsubrt(i,j)  !Units of qsubrt converted via CWATAVAIL
 280             !ELSE
 281             SUBFLO = subrt_data%state%qsubrt(i,j) / subrt_data%properties%distance_to_neighbor(i,j,9) * SUBDT !Convert qsubrt from m^3/s to m
 282             !END IF
 283
 284             WATAVAIL = 0.  !Initialize to 0. for every cell...
 285
 286
 287             !!DJG Begin loop through soil profile to adjust soil water content
 288             !!DJG based on subsfc flow (SUBFLO)...
 289
 290             IF (SUBFLO.GT.0) THEN ! Increase soil moist for +SUBFLO (Inflow)
 291
 292                 ! Loop through soil layers from bottom to top
 293                 DO KK=subrt_static%nsoil,1,-1
 294
 295
 296                     ! Check for saturated layers
 297                     IF (subrt_data%grid_transform%smcrt(I,J,KK) .GE. subrt_data%grid_transform%smcmaxrt(I,J,KK)) THEN
 298                         IF (subrt_data%grid_transform%smcrt(I,J,KK) .GT. subrt_data%grid_transform%smcmaxrt(I,J,KK)) THEN
 299                             print *, "FATAL ERROR: Subsfc acct. SMCMAX exceeded...", &
 300                                 subrt_data%grid_transform%smcrt(I,J,KK), subrt_data%grid_transform%smcmaxrt(I,J,KK),KK,i,j
 301                             call hydro_stop("In SUBSFC_RTNG() - SMCMAX exceeded")
 302                         ELSE
 303                         END IF
 304                     ELSE
 305                         WATAVAIL = (subrt_data%grid_transform%smcmaxrt(I,J,KK)-subrt_data%grid_transform%smcrt(I,J,KK))*subrt_data%properties%sldpth(KK)
 306                         IF (WATAVAIL.GE.SUBFLO) THEN
 307                             subrt_data%grid_transform%smcrt(I,J,KK) = subrt_data%grid_transform%smcrt(I,J,KK) + SUBFLO/subrt_data%properties%sldpth(KK)
 308                             SUBFLO = 0.
 309                         ELSE
 310                             SUBFLO = SUBFLO - WATAVAIL
 311                             subrt_data%grid_transform%smcrt(I,J,KK) = subrt_data%grid_transform%smcmaxrt(I,J,KK)
 312                         END IF
 313                     END IF
 314
 315                     IF (SUBFLO.EQ.0.) EXIT
 316                     !                IF (SUBFLO.EQ.0.) goto 669
 317
 318                 END DO      ! END DO FOR SOIL LAYERS
 319
 320                 669           continue
 321
 322                 ! If all layers sat. add remaining subflo to infilt. excess...
 323                 IF (KK.eq.0.AND.SUBFLO.gt.0.) then
 324                     subrt_input%infiltration_excess(I,J) = subrt_input%infiltration_excess(I,J) + SUBFLO*1000.    !Units = mm
 325                     SUBFLO=0.
 326                 END IF
 327
 328                 !DJG Error trap...
 329                 if (subflo.ne.0.) then
 330 #ifdef HYDRO_D
 331                     print *, "Subflo (+) not expired...:",subflo, i, j, kk, subrt_data%grid_transform%smcrt(i,j,1), &
 332                         subrt_data%grid_transform%smcrt(i,j,2), subrt_data%grid_transform%smcrt(i,j,3), &
 333                         subrt_data%grid_transform%smcrt(i,j,4), subrt_data%grid_transform%smcrt(i,j,5), &
 334                         subrt_data%grid_transform%smcrt(i,j,6), subrt_data%grid_transform%smcrt(i,j,7), &
 335                         subrt_data%grid_transform%smcrt(i,j,8), "SMCMAX", subrt_data%grid_transform%smcmaxrt(i,j,1)
 336 #endif
 337                 end if
 338
 339
 340             ELSE IF (SUBFLO.LT.0) THEN    ! Decrease soil moist for -SUBFLO (Drainage)
 341
 342
 343                 !DJG loop from satlyr back down and subtract out subflo as necess...
 344                 !    now set to SMCREF, 8/24/07
 345                 !DJG and then using unsat cond as opposed to Ksat...
 346
 347                 DO KK=SATLYRCHK(I,J),subrt_static%nsoil
 348                     WATAVAIL = (subrt_data%grid_transform%smcrt(I,J,KK) - subrt_data%grid_transform%smcrefrt(I,J,KK)) * subrt_data%properties%sldpth(KK)
 349                     IF (WATAVAIL.GE.ABS(SUBFLO)) THEN
 350                         !?yw mod                 IF (WATAVAIL.GE.(ABS(SUBFLO)+0.000001) ) THEN
 351                         subrt_data%grid_transform%smcrt(I,J,KK) = subrt_data%grid_transform%smcrt(I,J,KK) + SUBFLO/subrt_data%properties%sldpth(KK)
 352                         SUBFLO=0.
 353                     ELSE     ! Since subflo is small on a time-step following is unlikely...
 354                         subrt_data%grid_transform%smcrt(I,J,KK) = subrt_data%grid_transform%smcrefrt(I,J,KK)
 355                         SUBFLO=SUBFLO+WATAVAIL
 356                     END IF
 357                     IF (SUBFLO.EQ.0.) EXIT
 358                     !                IF (SUBFLO.EQ.0.) goto 668
 359
 360                 END DO  ! END DO FOR SOIL LAYERS
 361                 668        continue
 362
 363
 364                 !DJG Error trap...
 365                 if(abs(subflo) .le. 1.E-7 )  subflo = 0.0  !truncate residual to 1E-7 prec.
 366
 367                 if (subflo.ne.0.) then
 368 #ifdef HYDRO_D
 369                     print *, "Subflo (-) not expired:",i,j,subflo,CWATAVAIL(i,j)
 370                     print *, "zwatabl = ", subrt_data%properties%zwattablrt(i,j)
 371                     print *, "QSUBRT(I,J)=",subrt_data%state%qsubrt(i,j)
 372                     print *, "WATAVAIL = ",WATAVAIL, "kk=",kk
 373                     print *
 374 #endif
 375                 end if
 376
 377
 378
 379             END IF  ! end if for +/- SUBFLO soil moisture accounting...
 380
 381
 382
 383
 384         END DO        ! END DO X dim
 385     END DO          ! END DO Y dim
 386     !!!! End loop through subsurface routing domain...
 387
 388 #ifdef MPP_LAND
 389     do i = 1, subrt_static%nsoil
 390         call MPP_LAND_COM_REAL(subrt_data%grid_transform%smcrt(:,:,i),subrt_static%IXRT,subrt_static%JXRT,99)
 391     end DO
 392 #endif
 393
 394 #ifdef HYDRO_D
 395     ! ADCHANGE: START Final water balance variables
 396     ! ALL VARS in MM
 397     suminfxsrt2 = 0.
 398     qbdry2 = 0.
 399     smctot2 = 0.
 400     sumqsubrt2 = 0.
 401     do i=1,subrt_static%IXRT
 402         do j=1,subrt_static%JXRT
 403             suminfxsrt2 = suminfxsrt2 + subrt_input%infiltration_excess(I,J) / float(subrt_static%IXRT*subrt_static%JXRT)  !
 404             qbdry2 = qbdry2 + subrt_data%state%qsubbdryrt(i,j)/subrt_data%properties%distance_to_neighbor(i,j,9)*SUBDT*1000. / float(subrt_static%IXRT*subrt_static%JXRT)
 405             sumqsubrt2 = sumqsubrt2 + subrt_data%state%qsubrt(i,j)/subrt_data%properties%distance_to_neighbor(i,j,9)*SUBDT*1000. / float(subrt_static%IXRT*subrt_static%JXRT)
 406             do kk=1,subrt_static%nsoil
 407                 smctot2 = smctot2 + subrt_data%grid_transform%smcrt(I,J,KK)*subrt_data%properties%sldpth(KK)*1000. / float(subrt_static%IXRT*subrt_static%JXRT)
 408             end do
 409         end do
 410     end do
 411
 412 #ifdef MPP_LAND
 413     ! not tested
 414     CALL sum_real1(suminfxsrt2)
 415     CALL sum_real1(qbdry2)
 416     CALL sum_real1(sumqsubrt2)
 417     CALL sum_real1(smctot2)
 418     suminfxsrt2 = suminfxsrt2/float(numprocs)
 419     qbdry2 = qbdry2/float(numprocs)
 420     sumqsubrt2 = sumqsubrt2/float(numprocs)
 421     smctot2 = smctot2/float(numprocs)
 422 #endif
 423
 424 #ifdef MPP_LAND
 425     if (my_id .eq. IO_id) then
 426 #endif
 427         print *, "SUBSFC Routing Mass Bal: "
 428         print *, "WB_SUB!QsubDiff", sumqsubrt2-sumqsubrt1
 429         print *, "WB_SUB!Qsub1", sumqsubrt1
 430         print *, "WB_SUB!Qsub2", sumqsubrt2
 431         print *, "WB_SUB!InfxsDiff", suminfxsrt2-suminfxsrt1
 432         print *, "WB_SUB!Infxs1", suminfxsrt1
 433         print *, "WB_SUB!Infxs2", suminfxsrt2
 434         print *, "WB_SUB!QbdryDiff", qbdry2-qbdry1
 435         print *, "WB_SUB!Qbdry1", qbdry1
 436         print *, "WB_SUB!Qbdry2", qbdry2
 437         print *, "WB_SUB!SMCDiff", smctot2-smctot1
 438         print *, "WB_SUB!SMC1", smctot1
 439         print *, "WB_SUB!SMC2", smctot2
 440         print *, "WB_SUB!Residual", sumqsubrt1 - ( (suminfxsrt2-suminfxsrt1) &
 441             + (smctot2-smctot1) )
 442 #ifdef MPP_LAND
 443     endif
 444 #endif
 445     ! END Final water balance variables
 446 #endif
 447
 448
 449     !DJG ----------------------------------------------------------------
 450 END SUBROUTINE SUBSFC_RTNG
 451 !DJG ----------------------------------------------------------------
 452
 453
 454 !DJG ------------------------------------------------------------------------
 455 !DJG  SUBSURFACE FINDZWAT
 456 !DJG ------------------------------------------------------------------------
 457 SUBROUTINE FINDZWAT(IXRT,JXRT,NSOIL,SMCRT,SMCMAXRT,SMCREFRT, &
 458         SMCWLTRT,ZSOIL,SATLYRCHK,ZWATTABLRT,CWATAVAIL,&
 459         SLDPTH)
 460
 461     IMPLICIT NONE
 462
 463     !DJG -------- DECLARATIONS ------------------------
 464
 465     INTEGER, INTENT(IN) :: IXRT,JXRT,NSOIL
 466     REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCMAXRT
 467     REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCREFRT
 468     REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCRT
 469     REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCWLTRT
 470     REAL, INTENT(IN), DIMENSION(NSOIL)        :: ZSOIL
 471     REAL, INTENT(IN), DIMENSION(NSOIL)        :: SLDPTH
 472     REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: ZWATTABLRT
 473     REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: CWATAVAIL
 474     INTEGER, INTENT(OUT), DIMENSION(IXRT,JXRT) :: SATLYRCHK
 475
 476     !DJG Local Variables
 477     INTEGER :: KK,i,j
 478
 479     !!!! Find saturated layer depth...
 480     ! Loop through domain to determine sat. layers and assign wat tbl depth...
 481
 482
 483     SATLYRCHK = 0  !set flag for sat. layers
 484     CWATAVAIL = 0.  !set wat avail for subsfc rtng = 0.
 485
 486     DO J=1,JXRT
 487         DO I=1,IXRT
 488
 489             ! Loop through soil layers from bottom to top
 490             DO KK=NSOIL,1,-1
 491
 492                 ! Check for saturated layers
 493                 ! Add additional logical check to ensure water is 'available' for routing,
 494                 !  (i.e. not 'frozen' or otherwise immobile)
 495                 !                IF (SMCRT(I,J,KK).GE.SMCMAXRT(I,J,KK).AND.SMCMAXRT(I,J,KK) &
 496                     !                  .GT.SMCWLTRT(I,J,KK)) THEN
 497                 IF ( (SMCRT(I,J,KK).GE.SMCREFRT(I,J,KK)).AND.(SMCREFRT(I,J,KK) &
 498                         .GT.SMCWLTRT(I,J,KK)) ) THEN
 499                     ! Add additional check to ensure saturation from bottom up only...8/8/05
 500                     IF((SATLYRCHK(I,J).EQ.KK+1) .OR. (KK.EQ.NSOIL) ) SATLYRCHK(I,J) = KK
 501                 END IF
 502
 503             END DO
 504
 505
 506             ! Designate ZWATTABLRT based on highest sat. layer and
 507             ! Define amount of water avail for subsfc routing on each gridcell (CWATAVAIL)
 508             !  note: using a 'field capacity' value of SMCREF as lower limit...
 509
 510             IF (SATLYRCHK(I,J).ne.0) then
 511                 IF (SATLYRCHK(I,J).ne.1) then  ! soil column is partially sat.
 512                     ZWATTABLRT(I,J) = -ZSOIL(SATLYRCHK(I,J)-1)
 513                     !DJG 2/16/2016 fix                  DO KK=SATLYRCHK(I,J),NSOIL
 514                     !old                   CWATAVAIL(I,J) = (SMCRT(I,J,SATLYRCHK(I,J))-&
 515                         !old                                    SMCREFRT(I,J,SATLYRCHK(I,J))) * &
 516                         !old                                    (ZSOIL(SATLYRCHK(I,J)-1)-ZSOIL(NSOIL))
 517                     !DJG 2/16/2016 fix                    CWATAVAIL(I,J) = CWATAVAIL(I,J)+(SMCRT(I,J,KK)- &
 518                         !DJG 2/16/2016 fix                                     SMCREFRT(I,J,KK))*SLDPTH(KK)
 519                     !DJG 2/16/2016 fix                  END DO
 520
 521
 522                 ELSE  ! soil column is fully saturated to sfc.
 523                     ZWATTABLRT(I,J) = 0.
 524                     !DJG 2/16/2016 fix                  DO KK=SATLYRCHK(I,J),NSOIL
 525                     !DJG 2/16/2016 fix                    CWATAVAIL(I,J) = (SMCRT(I,J,KK)-SMCREFRT(I,J,KK))*SLDPTH(KK)
 526                     !DJG 2/16/2016 fix                  END DO
 527                 END IF
 528                 !DJG 2/16/2016 fix accumulation of CWATAVAIL...
 529                 DO KK=SATLYRCHK(I,J),NSOIL
 530                     CWATAVAIL(I,J) = CWATAVAIL(I,J)+(SMCRT(I,J,KK)- &
 531                         SMCREFRT(I,J,KK))*SLDPTH(KK)
 532                 END DO
 533             ELSE  ! no saturated layers...
 534                 ZWATTABLRT(I,J) = -ZSOIL(NSOIL)
 535                 SATLYRCHK(I,J) = NSOIL + 1
 536             END IF
 537
 538         END DO
 539     END DO
 540
 541
 542     !DJG ----------------------------------------------------------------
 543 END SUBROUTINE FINDZWAT
 544 !DJG ----------------------------------------------------------------
 545
 546 !===================================================================================================
 547 ! Subroutine Name:
 548 !   subroutine ROUTE_SUBSURFACE1
 549 ! Author(s)/Contact(s):
 550 !   D. Gochis <gochis><ucar><edu>
 551 ! Abstract:
 552 !   Calculates single direction (steepest slope) subsurface flow
 553 ! History Log:
 554 !   3/27/03 -Created, DG.
 555 !   1/05/04 -Modified, DG.
 556 !   1/07/19 -Modified, AD.
 557 ! Usage:
 558 ! Parameters:
 559 ! Input Files:
 560 ! Output Files:
 561 ! Condition codes:
 562 ! User controllable options: None.
 563 ! Notes:
 564 !  - Adapted from Wigmosta, 1994
 565 !  - Returns qsub=DQSUB which in turn becomes SUBFLO in head calc.
 566 !===================================================================================================
 567 SUBROUTINE ROUTE_SUBSURFACE1(dist, z, latksat, nexp, soldep, &
 568                              XX, YY, SO8RT, SO8RT_D, &
 569                              CWATAVAIL, SUBDT, QSUBDRY, QSUBDRYT, qsub)
 570 #ifdef MPP_LAND
 571    use module_mpp_land, only: left_id, down_id, right_id, up_id, &
 572                               mpp_land_com_real8, my_id, mpp_land_com_real
 573 #endif
 574    use module_hydro_stop, only:HYDRO_stop
 575
 576    implicit none
 577
 578    ! Passed variables
 579    integer, intent(in) :: XX,YY
 580    ! Passed arrays
 581    real, intent(in), dimension(XX,YY,9)    :: dist ! distance to neighbor cells (m)
 582    real, intent(in), dimension(XX,YY)      :: z ! depth to water table (m)
 583    real, intent(in), dimension(XX,YY)      :: latksat
 584                                               ! lateral saturated hydraulic conductivity (m/s)
 585    real, intent(in), dimension(XX,YY)      :: nexp ! latksat decay coefficient
 586    real, intent(in), dimension(XX,YY)      :: soldep ! soil depth (m)
 587    real, intent(in), dimension(XX,YY,8)    :: SO8RT ! terrain slope in all directions (m/m)
 588    integer, intent(in), dimension(XX,YY,3) :: SO8RT_D ! steepest terrain slope cell (i, j, index)
 589    real, intent(in), dimension(XX,YY)      :: CWATAVAIL ! water available for routing (m)
 590    real, intent(in)                        :: SUBDT ! subsurface routing timestep (s)
 591    real, intent(inout), dimension(XX,YY)   :: QSUBDRY ! subsurface flow at domain boundary (m3/s)
 592    real, intent(inout)                     :: QSUBDRYT
 593                                               ! total subsurface flow outside of domain (m3/s)
 594    real, intent(out), dimension(XX,YY)     :: qsub ! subsurface flow outside of cell (m3/s)
 595
 596    ! Local variables
 597    integer :: i, j ! indices for local loops
 598    integer :: IXX0, JYY0, index ! i, j, and direction index for steepest slope neighbor
 599    real    :: beta ! total head slope (m/m)
 600    real    :: hh ! interim variable for subsurface solution per DHSVM (dimensionless)
 601    real    :: gamma ! interim variable for subsurface solution per DHSVM (m3/s)
 602    real    :: ksat ! saturated hydraulic conductivity (m/s)
 603    real    :: qqsub ! net subsurface flow out of cell (m3/s?)
 604    real    :: waterToRoute ! total water to route from cell (m)
 605    ! Local arrays
 606    real*8, DIMENSION(XX,YY) :: qsub_tmp, QSUBDRY_tmp ! temp trackers for fluxes (m3/s)
 607    ! Local parameters
 608    real :: tmp_dist(9)
 609
 610    ! Initialize temp variables
 611    qsub_tmp = 0.
 612    QSUBDRY_tmp = 0.
 613
 614 #ifdef HYDRO_D
 615    write(6,*) "call subsurface routing xx= , yy =", yy, xx
 616 #endif
 617
 618    do j=2,YY-1 ! start j loop
 619
 620       do i=2,XX-1 ! start i loop
 621
 622          if (i.ge.2.AND.i.le.XX-1.AND.j.ge.2.AND.j.le.YY-1) then !if grdcl chk
 623             ! check for boundary grid point?
 624
 625             ! Set initial guess values for steepest slope neighbor based on terrain
 626             IXX0 = SO8RT_D(i,j,1)
 627             JYY0 = SO8RT_D(i,j,2)
 628             index = SO8RT_D(i,j,3)
 629             !beta = -1.0
 630
 631             ! ADCHANGE: Always call routine to check for steepest elev+head slope
 632             ! Updates IXX0, JYY0, index, and beta for steepest slope neighbor.
 633             ! Note that beta will be -1 if steepest slope cannot be found.
 634             tmp_dist = dist(i,j,:)
 635             call GETSUB8(i, j, XX, YY, Z, so8rt, tmp_dist, &
 636                          IXX0, JYY0, index, beta)
 637
 638             if (dist(i,j,index) .le. 0) then
 639                 write(6,*) "FATAL ERROR: dist(i,j,index) is <= zero "
 640                 call hydro_stop("In ROUTE_SUBSURFACE1() - dist(i,j,index) is <= zero ")
 641             endif
 642             if (soldep(i,j) .eq. 0) then
 643                 call hydro_stop("In ROUTE_SUBSURFACE1() - soldep is = zero")
 644             endif
 645
 646             if (beta .gt. 0) then            !if-then for flux calc
 647                if (beta .lt. 1E-20 ) then
 648 #ifdef HYDRO_D
 649                   print*, "Message: beta need to be reset to 1E-20. beta = ",beta
 650 #endif
 651                   beta = 1E-20
 652                endif
 653
 654                ! Do the rest if the lowest grid can be found.
 655                hh = ( 1 - ( z(i,j) / soldep(i,j) ) ) ** nexp(i,j)
 656                ksat = latksat(i,j)
 657
 658                if (hh .lt. 0.) then
 659                   print *, "hsub<0 at gridcell...", i,j,hh,z(i+1,j),z(i,j), &
 660                        soldep(i,j)
 661                   call hydro_stop("In ROUTE_SUBSURFACE1() - hsub<0 at gridcell ")
 662                endif
 663
 664                ! Calculate flux from cell
 665                ! AD_NOTE: gamma and qqsub are negative when flow is out of cell
 666                gamma = -1.0 * ( (dist(i,j,index) * ksat * soldep(i,j)) / nexp(i,j) ) * beta
 667                qqsub = gamma * hh
 668
 669                ! AD_NOTE: Moved this water available constraint from outside qsub calc loop
 670                ! to inside to better account for adjustments to adjacent cells.
 671                ! Calculate total water to route (where dist(i,j,9) is cell area):
 672                waterToRoute = ABS(qqsub) / dist(i,j,9) * SUBDT
 673                if ( (qqsub .le. 0.0) .and. (CWATAVAIL(i,j) .lt. waterToRoute) ) THEN
 674                   qqsub = -1.0 * CWATAVAIL(i,j) * dist(i,j,9) / SUBDT
 675                endif
 676
 677                ! Remove from cell qsub to track net fluxes over full i, j loop
 678                ! (remember: qqsub is negative when flow is out of cell!)
 679                qsub(i,j) = qsub(i,j) + qqsub
 680                ! Add to neighbor cell qsub to track net fluxes over full i, j loop
 681                qsub_tmp(ixx0,jyy0) = qsub_tmp(ixx0,jyy0) - qqsub
 682
 683                !DJG Error Checks...
 684                if (qqsub .gt. 0) then
 685                   print*, "FATAL ERROR: qqsub should be negative, qqsub =",qqsub,&
 686                      "gamma=",gamma,"hh=",hh,"beta=",beta,&
 687                      "so8RT=",so8RT(i,j,index),"latksat=",ksat, &
 688                      "tan(beta)=",tan(beta),i,j,z(i,j),z(IXX0,JYY0)
 689                   print*, "ixx0=",ixx0, "jyy0=",jyy0
 690                   print*, "soldep =", soldep(i,j), "nexp=",nexp(i,j)
 691                   call hydro_stop("In ROUTE_SUBSURFACE1() - qqsub should be negative")
 692                endif
 693
 694                ! Make boundary adjustments if cells are on edge of local domain
 695 #ifdef MPP_LAND
 696                if ( ((ixx0.eq.XX).and.(right_id .lt. 0)) .or. &
 697                   ((ixx0.eq.1) .and.(left_id  .lt. 0)) .or. &
 698                   ((jyy0.eq.1) .and.(down_id  .lt. 0)) .or. &
 699                   ((JYY0.eq.YY).and.(up_id    .lt. 0)) ) then
 700 #else
 701                if ((ixx0.eq.1).or.(ixx0.eq.xx).or.(jyy0.eq.1).or.(jyy0.eq.yy)) then
 702 #endif
 703                   ! If on edge, move flux tracking BACK from neighbor cell qsub and into
 704                   ! boundary qsub (note: qsubdry is negative and qqsub is negative, so we
 705                   ! are making it a bigger sink here)
 706                   qsub_tmp(ixx0,jyy0) = qsub_tmp(ixx0,jyy0) + qqsub
 707                   QSUBDRY_tmp(ixx0,jyy0) = QSUBDRY_tmp(ixx0,jyy0) + qqsub
 708                   ! Add to total BOUNDARY qsub
 709                   QSUBDRYT = QSUBDRYT + qqsub
 710                endif
 711
 712             endif  ! end if for flux calc
 713
 714          endif   ! end if for gridcell check
 715
 716       end do  ! end i loop
 717
 718    end do   ! end j loop
 719
 720 #ifdef MPP_LAND
 721    call MPP_LAND_COM_REAL8(qsub_tmp,XX,YY,1)
 722    call MPP_LAND_COM_REAL8(QSUBDRY_tmp,XX,YY,1)
 723 #endif
 724
 725    ! Sum grids to get net of self and neighbor fluxes after looping through all cells
 726    ! (i.e., net of out-flux from one cell and in-flux from its neighbor cell)
 727    qsub = qsub + qsub_tmp
 728    QSUBDRY = QSUBDRY + QSUBDRY_tmp
 729
 730
 731 #ifdef MPP_LAND
 732    call MPP_LAND_COM_REAL(qsub,XX,YY,99)
 733    call MPP_LAND_COM_REAL(QSUBDRY,XX,YY,99)
 734 #endif
 735
 736    return
 737
 738 end subroutine ROUTE_SUBSURFACE1
 739
 740 !===================================================================================================
 741
 742
 743
 744 !DJG ----------------------------------------------------------------
 745 !DJG     SUBROUTINE ROUTE_SUBSURFACE2
 746 !DJG ----------------------------------------------------------------
 747
 748 SUBROUTINE ROUTE_SUBSURFACE2(                                 &
 749         dist,z,qsub,sox,soy,                                   &
 750         latksat,soldep,XX,YY,QSUBDRY,QSUBDRYT,CWATAVAIL,   &
 751         SUBDT)
 752
 753     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 754     !
 755     !  Subroutine to route subsurface flow through the watershed
 756     !DJG ----------------------------------------------------------------
 757     !
 758     !  Called from: main.f (Noah_router_driver)
 759     !
 760     !  Returns: qsub=DQSUB   which in turn becomes SUBFLO in head calc.
 761     !
 762     !  Created:    D. Gochis                           3/27/03
 763     !              Adaptded from Wigmosta, 1994
 764     !
 765     !  Modified:   D. Gochis                           1/05/04
 766     !
 767     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 768 #ifdef MPP_LAND
 769     use module_mpp_land, only: left_id,down_id,right_id,&
 770         up_id,mpp_land_com_real,MPP_LAND_UB_COM, &
 771         MPP_LAND_LR_COM,mpp_land_com_integer
 772 #endif
 773     use module_hydro_stop, only:HYDRO_stop
 774
 775     IMPLICIT NONE
 776
 777
 778     !! Declare Passed variables
 779
 780     INTEGER, INTENT(IN) :: XX,YY
 781
 782     !! Declare passed arrays
 783
 784     REAL, INTENT(IN), DIMENSION(XX,YY) :: z
 785     REAL, INTENT(IN), DIMENSION(XX,YY) :: sox
 786     REAL, INTENT(IN), DIMENSION(XX,YY) :: soy
 787     REAL, INTENT(IN), DIMENSION(XX,YY) :: latksat
 788     REAL, INTENT(IN), DIMENSION(XX,YY) :: CWATAVAIL
 789     REAL, INTENT(IN), DIMENSION(XX,YY) :: soldep
 790     REAL, INTENT(OUT), DIMENSION(XX,YY) :: qsub
 791     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: QSUBDRY
 792     REAL, INTENT(INOUT) :: QSUBDRYT
 793     REAL, INTENT(IN) :: SUBDT
 794     real, intent(in), dimension(xx,yy,9) :: dist
 795
 796     !!! Declare Local Variables
 797
 798     REAL :: dzdx,dzdy,beta,gamma
 799     REAL :: qqsub,hh,ksat, gsize
 800
 801     INTEGER :: i,j
 802     !!! Initialize variables
 803     REAL, PARAMETER :: nexp=1.0      ! local power law exponent
 804     qsub = 0.                        ! initialize flux = 0. !DJG 5 May 2014
 805
 806     !yw        soldep = 2.
 807
 808
 809     ! Begin Subsurface routing
 810
 811     !!! Loop to route water in x-direction
 812     do j=1,YY
 813         do i=1,XX
 814             ! check for boundary grid point?
 815             if (i.eq.XX) GOTO 998
 816             gsize = dist(i,j,3)
 817
 818             dzdx= (z(i,j) - z(i+1,j))/gsize
 819             beta=sox(i,j) + dzdx + 1E-30
 820             if (abs(beta) .lt. 1E-20) beta=1E-20
 821             if (beta.lt.0) then
 822                 !yw            hh=(1-(z(i+1,j)/soldep(i,j)))**nexp
 823                 hh=(1-(z(i+1,j)/soldep(i+1,j)))**nexp
 824                 ! Change later to use mean Ksat of two cells
 825                 ksat=latksat(i+1,j)
 826             else
 827                 hh=(1-(z(i,j)/soldep(i,j)))**nexp
 828                 ksat=latksat(i,j)
 829             end if
 830
 831             if (hh .lt. 0.) then
 832                 print *, "hsub<0 at gridcell...", i,j,hh,z(i+1,j),z(i,j), &
 833                     soldep(i,j),nexp
 834                 call hydro_stop("In ROUTE_SUBSURFACE2() - hsub<0 at gridcell")
 835             end if
 836
 837             !Err. tan slope          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*tan(beta)
 838             !AD_CHANGE: beta is already a slope so no tan (consistent with ROUTE_SUBSURFACE1)
 839             gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*beta
 840             !DJG lacks tan(beta) of original Wigmosta version          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*beta
 841
 842             qqsub = gamma * hh
 843             qsub(i,j) = qsub(i,j) + qqsub
 844             qsub(i+1,j) = qsub(i+1,j) - qqsub
 845
 846             ! Boundary adjustments
 847 #ifdef MPP_LAND
 848             if ((i.eq.1).AND.(beta.lt.0.).and.(left_id.lt.0)) then
 849 #else
 850                 if ((i.eq.1).AND.(beta.lt.0.)) then
 851 #endif
 852                     qsub(i,j) = qsub(i,j) - qqsub
 853                     QSUBDRY(i,j) = QSUBDRY(i,j) - qqsub
 854                     QSUBDRYT = QSUBDRYT - qqsub
 855 #ifdef MPP_LAND
 856                 else if ((i.eq.(xx-1)).AND.(beta.gt.0.) &
 857                         .and.(right_id.lt.0) ) then
 858 #else
 859                 else if ((i.eq.(xx-1)).AND.(beta.gt.0.)) then
 860 #endif
 861                     qsub(i+1,j) = qsub(i+1,j) + qqsub
 862                     QSUBDRY(i+1,j) = QSUBDRY(i+1,j) + qqsub
 863                     QSUBDRYT = QSUBDRYT + qqsub
 864                 end if
 865
 866                 998       continue
 867
 868                 !! End loop to route sfc water in x-direction
 869         end do
 870     end do
 871
 872 #ifdef MPP_LAND
 873     call MPP_LAND_LR_COM(qsub,XX,YY,99)
 874     call MPP_LAND_LR_COM(QSUBDRY,XX,YY,99)
 875 #endif
 876
 877
 878     !!! Loop to route water in y-direction
 879     do j=1,YY
 880         do i=1,XX
 881             ! check for boundary grid point?
 882             if (j.eq.YY) GOTO 999
 883             gsize = dist(i,j,1)
 884
 885             dzdy= (z(i,j) - z(i,j+1))/gsize
 886             beta=soy(i,j) + dzdy + 1E-30
 887             if (abs(beta) .lt. 1E-20) beta=1E-20
 888             if (beta.lt.0) then
 889                 !yw            hh=(1-(z(i,j+1)/soldep(i,j)))**nexp
 890                 hh=(1-(z(i,j+1)/soldep(i,j+1)))**nexp
 891                 ksat=latksat(i,j+1)
 892             else
 893                 hh=(1-(z(i,j)/soldep(i,j)))**nexp
 894                 ksat=latksat(i,j)
 895             end if
 896
 897             if (hh .lt. 0.) GOTO 999
 898
 899             !Err. tan slope          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*tan(beta)
 900             gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*beta
 901
 902             qqsub = gamma * hh
 903             qsub(i,j) = qsub(i,j) + qqsub
 904             qsub(i,j+1) = qsub(i,j+1) - qqsub
 905
 906             ! Boundary adjustments
 907
 908 #ifdef MPP_LAND
 909             if ((j.eq.1).AND.(beta.lt.0.).and.(down_id.lt.0)) then
 910 #else
 911             if ((j.eq.1).AND.(beta.lt.0.)) then
 912 #endif
 913                 qsub(i,j) = qsub(i,j) - qqsub
 914                 QSUBDRY(i,j) = QSUBDRY(i,j) - qqsub
 915                 QSUBDRYT = QSUBDRYT - qqsub
 916 #ifdef MPP_LAND
 917             else if ((j.eq.(yy-1)).AND.(beta.gt.0.)  &
 918                 .and. (up_id.lt.0) ) then
 919 #else
 920             else if ((j.eq.(yy-1)).AND.(beta.gt.0.)) then
 921 #endif
 922                 qsub(i,j+1) = qsub(i,j+1) + qqsub
 923                 QSUBDRY(i,j+1) = QSUBDRY(i,j+1) + qqsub
 924                 QSUBDRYT = QSUBDRYT + qqsub
 925             end if
 926
 927             999       continue
 928
 929                 !! End loop to route sfc water in y-direction
 930         end do
 931     end do
 932
 933 #ifdef MPP_LAND
 934     call MPP_LAND_UB_COM(qsub,XX,YY,99)
 935     call MPP_LAND_UB_COM(QSUBDRY,XX,YY,99)
 936 #endif
 937
 938     return
 939     !DJG------------------------------------------------------------
 940 end subroutine ROUTE_SUBSURFACE2
 941 !DJG------------------------------------------------------------