!  Program Name:
!  Author(s)/Contact(s):
!  Abstract:
!  History Log:
! 
!  Usage:
!  Parameters: <Specify typical arguments passed>
!  Input Files:
!        <list file names and briefly describe the data they include>
!  Output Files:
!        <list file names and briefly describe the information they include>
! 
!  Condition codes:
!        <list exit condition or error codes returned >
!        If appropriate, descriptive troubleshooting instructions or
!        likely causes for failures could be mentioned here with the
!        appropriate error code
! 
!  User controllable options: <if applicable>

!DJG ------------------------------------------------
!DJG   SUBROUTINE RT_PARM
!DJG ------------------------------------------------

	SUBROUTINE RT_PARM(IX,JY,IXRT,JXRT,VEGTYP,RETDP,OVRGH,  &
                      AGGFACTR)
#ifdef MPP_LAND
        use module_mpp_land, only: left_id,down_id,right_id,&
               up_id,mpp_land_com_real,MPP_LAND_UB_COM, &
               MPP_LAND_LR_COM,mpp_land_com_integer 
#endif

	IMPLICIT NONE

!DJG -------- DECLARATIONS -----------------------
 
	INTEGER, INTENT(IN) :: IX,JY,IXRT,JXRT,AGGFACTR

	INTEGER, INTENT(IN), DIMENSION(IX,JY)	:: VEGTYP
	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)	:: RETDP
	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)	:: OVRGH


!DJG Local Variables

	INTEGER	:: I,J,IXXRT,JYYRT
        INTEGER :: AGGFACYRT,AGGFACXRT


!DJG Assign RETDP and OVRGH based on VEGTYP...

	do J=1,JY
          do I=1,IX

             do AGGFACYRT=AGGFACTR-1,0,-1
              do AGGFACXRT=AGGFACTR-1,0,-1

               IXXRT=I*AGGFACTR-AGGFACXRT
               JYYRT=J*AGGFACTR-AGGFACYRT
#ifdef MPP_LAND
       if(left_id.ge.0) IXXRT=IXXRT+1
       if(down_id.ge.0) JYYRT=JYYRT+1
#else
!yw ????
!       IXXRT=IXXRT+1
!       JYYRT=JYYRT+1
#endif

!        if(AGGFACTR .eq. 1) then
!            IXXRT=I
!            JYYRT=J
!        endif



!DJG Urban, rock, playa, snow/ice...
	       IF (VEGTYP(I,J).EQ.1.OR.VEGTYP(I,J).EQ.26.OR.   &
                      VEGTYP(I,J).EQ.26.OR.VEGTYP(I,J).EQ.24) THEN
                 RETDP(IXXRT,JYYRT)=1.3
                 OVRGH(IXXRT,JYYRT)=0.1
!DJG Wetlands and water bodies...
	       ELSE IF (VEGTYP(I,J).EQ.17.OR.VEGTYP(I,J).EQ.18.OR.  &
                      VEGTYP(I,J).EQ.19.OR.VEGTYP(I,J).EQ.16) THEN
                 RETDP(IXXRT,JYYRT)=10.0
                 OVRGH(IXXRT,JYYRT)=0.2
!DJG All other natural covers...
               ELSE 
                 RETDP(IXXRT,JYYRT)=5.0
                 OVRGH(IXXRT,JYYRT)=0.2
               END IF

              end do
             end do

          end do
        end do
#ifdef MPP_LAND
        call MPP_LAND_COM_REAL(RETDP,IXRT,JXRT,99)
        call MPP_LAND_COM_REAL(OVRGH,IXRT,JXRT,99)
#endif

!DJG ----------------------------------------------------------------
  END SUBROUTINE RT_PARM
!DJG ----------------------------------------------------------------





!DJG ------------------------------------------------
!DJG   SUBROUTINE SUBSFC_RTNG
!DJG ------------------------------------------------

	SUBROUTINE SUBSFC_RTNG(dist,ZWATTABLRT,QSUBRT,SOXRT,    &
          SOYRT,LATKSATRT,SOLDEPRT,QSUBBDRYRT,QSUBBDRYTRT,      &    
          NSOIL,SMCRT,INFXSUBRT,SMCMAXRT,SMCREFRT,ZSOIL,IXRT,JXRT,DT,    &
          SMCWLTRT,SO8RT,SO8RT_D, rt_option,SLDPTH,junk4,CWATAVAIL, &
          SATLYRCHK)

!       use module_mpp_land, only: write_restart_rt_3, write_restart_rt_2, &
!            my_id
#ifdef MPP_LAND
        use module_mpp_land, only: MPP_LAND_COM_REAL, sum_real1, &
		my_id, io_id, numprocs
#endif
	IMPLICIT NONE

!DJG -------- DECLARATIONS ------------------------

	INTEGER, INTENT(IN) :: IXRT,JXRT,NSOIL

	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: SOXRT,junk4
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: SOYRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: LATKSATRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)    :: SOLDEPRT

	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)   :: ZWATTABLRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)   :: CWATAVAIL
        INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: SATLYRCHK


	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: QSUBRT
	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: QSUBBDRYRT

	REAL, INTENT(IN)                          :: dist(ixrt,jxrt,9)
	REAL, INTENT(IN)                          :: DT
	REAL, INTENT(IN), DIMENSION(NSOIL)        :: ZSOIL
	REAL, INTENT(IN), DIMENSION(NSOIL) 	  :: SLDPTH
	REAL, INTENT(INOUT)                       :: QSUBBDRYTRT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: INFXSUBRT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,NSOIL) :: SMCMAXRT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,NSOIL) :: SMCREFRT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,NSOIL) :: SMCRT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,NSOIL) :: SMCWLTRT

	REAL, DIMENSION(IXRT,JXRT)	:: ywtmp
!DJG Local Variables

	INTEGER	:: I,J,KK
!djg        INTEGER, DIMENSION(IXRT,JXRT) :: SATLYRCHK

	REAL 	:: GRDAREA
	REAL	:: SUBFLO
	REAL	:: WATAVAIL

        INTEGER :: SO8RT_D(IXRT,JXRT,3)
        REAL :: SO8RT(IXRT,JXRT,8)
        integer ::  rt_option, index

        INTEGER :: DT_STEPS             !-- number of timestep in routing
        REAL :: SUBDT                !-- subsurface routing timestep
        INTEGER :: KRT                  !-- routing counter
        REAL, DIMENSION(IXRT,JXRT,NSOIL) :: SMCTMP  !--temp store of SMC
        REAL, DIMENSION(IXRT,JXRT) :: ZWATTABLRTTMP ! temp store of ZWAT
        REAL, DIMENSION(IXRT,JXRT) :: INFXSUBRTTMP ! temp store of infilx
!djg        REAL, DIMENSION(IXRT,JXRT) :: CWATAVAIL ! temp specif. of wat avial
        
        

!DJG Debug Variables...
        REAL :: qsubchk,qsubbdrytmp
        REAL :: junk1,junk2,junk3,junk5,junk6,junk7
        INTEGER, PARAMETER :: double=8
        REAL (KIND=double) :: smctot1a,smctot2a
	INTEGER :: kx,count

#ifdef HYDRO_D
! ADCHANGE: Water balance variables
       real   :: smctot1,smctot2
       real   :: suminfxsrt1,suminfxsrt2
       real   :: qbdry1,qbdry2
       real   :: sumqsubrt1, sumqsubrt2
#endif
        
!DJG -----------------------------------------------------------------
!DJG  SUBSURFACE ROUTING LOOP
!DJG    - SUBSURFACE ROUTING RUN ON NOAH TIMESTEP
!DJG    - SUBSURFACE ROUITNG ONLY PERFORMED ON SATURATED LAYERS
!DJG -----------------------------------------------------------------

#ifdef HYDRO_D
! ADCHANGE: START Initial water balance variables 
! ALL VARS in MM
       suminfxsrt1 = 0.
       qbdry1 = 0.
       smctot1 = 0.
       sumqsubrt1 = 0.
       do i=1,IXRT
         do j=1,JXRT
           suminfxsrt1 = suminfxsrt1 + INFXSUBRT(I,J) / float(IXRT*JXRT)
           qbdry1 = qbdry1 + QSUBBDRYRT(I,J)/dist(i,j,9)*SUBDT*1000. / float(IXRT*JXRT)
           sumqsubrt1 = sumqsubrt1 + QSUBRT(I,J)/dist(i,j,9)*SUBDT*1000. / float(IXRT*JXRT)
           do kk=1,NSOIL
               smctot1 = smctot1 + SMCRT(I,J,KK)*SLDPTH(KK)*1000. / float(IXRT*JXRT)
           end do
         end do
       end do

#ifdef MPP_LAND
! not tested
       CALL sum_real1(suminfxsrt1)
       CALL sum_real1(qbdry1)
       CALL sum_real1(sumqsubrt1)
       CALL sum_real1(smctot1)
       suminfxsrt1 = suminfxsrt1/float(numprocs)
       qbdry1 = qbdry1/float(numprocs)
       sumqsubrt1 = sumqsubrt1/float(numprocs)
       smctot1 = smctot1/float(numprocs)
#endif
! END Initial water balance variables
#endif


        !yw GRDAREA=DXRT*DXRT
        ! GRDAREA=dist(i,j,9)


!DJG debug subsfc...
         subflo = 0.0

!DJG Set up mass balance checks...
!         CWATAVAIL = 0.            !-- initialize subsurface watavail
         SUBDT = DT                !-- initialize the routing timestep to DT


!!!! Find saturated layer depth...
! Loop through domain to determine sat. layers and assign wat tbl depth...
!    and water available for subsfc routing (CWATAVAIL)...
!
!         CALL FINDZWAT(IXRT,JXRT,NSOIL,SMCRT,SMCMAXRT,SMCREFRT, &
!                             SMCWLTRT,ZSOIL,SATLYRCHK,ZWATTABLRT, &
!                             CWATAVAIL,SLDPTH)
         



!DJG debug variable...

!DJG Courant check temp variable setup...
         ZWATTABLRTTMP = ZWATTABLRT !-- temporary storage of water table level




!!!! Call subsurface routing subroutine...
#ifdef HYDRO_D
     print *, "calling subsurface routing subroutine...Opt. ",rt_option
#endif


     if(rt_option .eq. 1) then
        CALL ROUTE_SUBSURFACE1(dist,ZWATTABLRT,QSUBRT,SOXRT,SOYRT,  &   
               LATKSATRT,SOLDEPRT,IXRT,JXRT,QSUBBDRYRT,QSUBBDRYTRT, &   
               SO8RT,SO8RT_D,CWATAVAIL,SUBDT)
     else 
        CALL ROUTE_SUBSURFACE2(dist,ZWATTABLRT,QSUBRT,SOXRT,SOYRT,      &
               LATKSATRT,SOLDEPRT,IXRT,JXRT,QSUBBDRYRT,QSUBBDRYTRT,     &
               CWATAVAIL,SUBDT)
     end if

#ifdef HYDRO_D
     write(6,*) "finish calling ROUTE_SUBSURFACE ", rt_option
#endif


!!!! Update soil moisture fields with subsurface flow...

!!!! Loop through subsurface routing domain...
	DO I=1,IXRT
          DO J=1,JXRT

!!DJG Check for courant condition violation...put limit on qsub
!!DJG QSUB HAS units of m^3/s SUBFLO has units of m
          
! ADCHANGE: Moved this constraint to the ROUTE_SUBSURFACE routines
           !IF (CWATAVAIL(i,j).le.ABS(qsubrt(i,j))/dist(i,j,9)*SUBDT) THEN
           !  QSUBRT(i,j) = -1.0*CWATAVAIL(i,j)
           !  SUBFLO = QSUBRT(i,j)  !Units of qsubrt converted via CWATAVAIL
           !ELSE
             SUBFLO=QSUBRT(I,J)/dist(i,j,9)*SUBDT !Convert qsubrt from m^3/s to m
           !END IF

           WATAVAIL=0.  !Initialize to 0. for every cell...


!!DJG Begin loop through soil profile to adjust soil water content
!!DJG based on subsfc flow (SUBFLO)...

            IF (SUBFLO.GT.0) THEN ! Increase soil moist for +SUBFLO (Inflow)

! Loop through soil layers from bottom to top
              DO KK=NSOIL,1,-1


! Check for saturated layers
                IF (SMCRT(I,J,KK).GE.SMCMAXRT(I,J,KK)) THEN
                  IF (SMCRT(I,J,KK).GT.SMCMAXRT(I,J,KK)) THEN
                   print *, "FATAL ERROR: Subsfc acct. SMCMAX exceeded...", &
                       SMCRT(I,J,KK), SMCMAXRT(I,J,KK),KK,i,j
                   call hydro_stop("In SUBSFC_RTNG() - SMCMAX exceeded")
                  ELSE
                  END IF
                ELSE
                  WATAVAIL = (SMCMAXRT(I,J,KK)-SMCRT(I,J,KK))*SLDPTH(KK)
                  IF (WATAVAIL.GE.SUBFLO) THEN
                    SMCRT(I,J,KK) = SMCRT(I,J,KK) + SUBFLO/SLDPTH(KK)
                    SUBFLO = 0.
                  ELSE
                    SUBFLO = SUBFLO - WATAVAIL
                    SMCRT(I,J,KK) = SMCMAXRT(I,J,KK)
                  END IF
                END IF

                 IF (SUBFLO.EQ.0.) EXIT
!                IF (SUBFLO.EQ.0.) goto 669

              END DO      ! END DO FOR SOIL LAYERS

669           continue

! If all layers sat. add remaining subflo to infilt. excess...                  
              IF (KK.eq.0.AND.SUBFLO.gt.0.) then
                 INFXSUBRT(I,J) = INFXSUBRT(I,J) + SUBFLO*1000.    !Units = mm
                 SUBFLO=0.
              END IF

!DJG Error trap...
	       if (subflo.ne.0.) then
#ifdef HYDRO_D
                  print *, "Subflo (+) not expired...:",subflo,i,j,kk,SMCRT(i,j,1), &
                           SMCRT(i,j,2),SMCRT(i,j,3),SMCRT(i,j,4),SMCRT(i,j,5),  &
                           SMCRT(i,j,6),SMCRT(i,j,7),SMCRT(i,j,8),"SMCMAX",SMCMAXRT(i,j,1)
#endif
               end if

 
            ELSE IF (SUBFLO.LT.0) THEN    ! Decrease soil moist for -SUBFLO (Drainage)


!DJG loop from satlyr back down and subtract out subflo as necess...
!    now set to SMCREF, 8/24/07
!DJG and then using unsat cond as opposed to Ksat...

	      DO KK=SATLYRCHK(I,J),NSOIL
                 WATAVAIL = (SMCRT(I,J,KK)-SMCREFRT(I,J,KK))*SLDPTH(KK)
                 IF (WATAVAIL.GE.ABS(SUBFLO)) THEN
!?yw mod                 IF (WATAVAIL.GE.(ABS(SUBFLO)+0.000001) ) THEN
                   SMCRT(I,J,KK) = SMCRT(I,J,KK) + SUBFLO/SLDPTH(KK)
                   SUBFLO=0.
                 ELSE     ! Since subflo is small on a time-step following is unlikely...
                   SMCRT(I,J,KK)=SMCREFRT(I,J,KK)
                   SUBFLO=SUBFLO+WATAVAIL
                 END IF
                 IF (SUBFLO.EQ.0.) EXIT
!                IF (SUBFLO.EQ.0.) goto 668

              END DO  ! END DO FOR SOIL LAYERS
668        continue


!DJG Error trap...
              if(abs(subflo) .le. 1.E-7 )  subflo = 0.0  !truncate residual to 1E-7 prec.

	       if (subflo.ne.0.) then
#ifdef HYDRO_D
                  print *, "Subflo (-) not expired:",i,j,subflo,CWATAVAIL(i,j)
                  print *, "zwatabl = ", ZWATTABLRT(I,J)
                  print *, "QSUBRT(I,J)=",QSUBRT(I,J)
                  print *, "WATAVAIL = ",WATAVAIL, "kk=",kk
                  print *
#endif
               end if



            END IF  ! end if for +/- SUBFLO soil moisture accounting...




          END DO        ! END DO X dim
        END DO          ! END DO Y dim
!!!! End loop through subsurface routing domain...

#ifdef MPP_LAND
     do i = 1, NSOIL
        call MPP_LAND_COM_REAL(SMCRT(:,:,i),IXRT,JXRT,99)
     end DO
#endif

#ifdef HYDRO_D
! ADCHANGE: START Final water balance variables
! ALL VARS in MM
        suminfxsrt2 = 0.
        qbdry2 = 0.
        smctot2 = 0.
        sumqsubrt2 = 0.
        do i=1,IXRT
         do j=1,JXRT
            suminfxsrt2 = suminfxsrt2 + INFXSUBRT(I,J) / float(IXRT*JXRT)
            qbdry2 = qbdry2 + QSUBBDRYRT(I,J)/dist(i,j,9)*SUBDT*1000. / float(IXRT*JXRT)
            sumqsubrt2 = sumqsubrt2 + QSUBRT(I,J)/dist(i,j,9)*SUBDT*1000. / float(IXRT*JXRT)
            do kk=1,NSOIL
                smctot2 = smctot2 + SMCRT(I,J,KK)*SLDPTH(KK)*1000. / float(IXRT*JXRT)
            end do
         end do
        end do

#ifdef MPP_LAND
! not tested
        CALL sum_real1(suminfxsrt2)
        CALL sum_real1(qbdry2)
        CALL sum_real1(sumqsubrt2)
        CALL sum_real1(smctot2)
        suminfxsrt2 = suminfxsrt2/float(numprocs)
        qbdry2 = qbdry2/float(numprocs)
        sumqsubrt2 = sumqsubrt2/float(numprocs)
        smctot2 = smctot2/float(numprocs)
#endif

#ifdef MPP_LAND   
       if (my_id .eq. IO_id) then
#endif
       print *, "SUBSFC Routing Mass Bal: "
       print *, "WB_SUB!QsubDiff", sumqsubrt2-sumqsubrt1
       print *, "WB_SUB!Qsub1", sumqsubrt1
       print *, "WB_SUB!Qsub2", sumqsubrt2
       print *, "WB_SUB!InfxsDiff", suminfxsrt2-suminfxsrt1
       print *, "WB_SUB!Infxs1", suminfxsrt1
       print *, "WB_SUB!Infxs2", suminfxsrt2
       print *, "WB_SUB!QbdryDiff", qbdry2-qbdry1
       print *, "WB_SUB!Qbdry1", qbdry1
       print *, "WB_SUB!Qbdry2", qbdry2
       print *, "WB_SUB!SMCDiff", smctot2-smctot1
       print *, "WB_SUB!SMC1", smctot1
       print *, "WB_SUB!SMC2", smctot2
       print *, "WB_SUB!Residual", sumqsubrt1 - ( (suminfxsrt2-suminfxsrt1) &
                       + (smctot2-smctot1) )
#ifdef MPP_LAND
       endif
#endif
! END Final water balance variables
#endif


!DJG ----------------------------------------------------------------
  END SUBROUTINE SUBSFC_RTNG 
!DJG ----------------------------------------------------------------


!DJG ------------------------------------------------------------------------
!DJG  SUBSURFACE FINDZWAT
!DJG ------------------------------------------------------------------------
         SUBROUTINE FINDZWAT(IXRT,JXRT,NSOIL,SMCRT,SMCMAXRT,SMCREFRT, &
                             SMCWLTRT,ZSOIL,SATLYRCHK,ZWATTABLRT,CWATAVAIL,&
                             SLDPTH)

	IMPLICIT NONE

!DJG -------- DECLARATIONS ------------------------

	INTEGER, INTENT(IN) :: IXRT,JXRT,NSOIL
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCMAXRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCREFRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT,NSOIL) :: SMCWLTRT
	REAL, INTENT(IN), DIMENSION(NSOIL)        :: ZSOIL
	REAL, INTENT(IN), DIMENSION(NSOIL)        :: SLDPTH
	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: ZWATTABLRT
	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: CWATAVAIL
        INTEGER, INTENT(OUT), DIMENSION(IXRT,JXRT) :: SATLYRCHK
       
!DJG Local Variables
        INTEGER :: KK,i,j


!!!! Find saturated layer depth...
! Loop through domain to determine sat. layers and assign wat tbl depth...


        SATLYRCHK = 0  !set flag for sat. layers
        CWATAVAIL = 0.  !set wat avail for subsfc rtng = 0.

        DO J=1,JXRT
          DO I=1,IXRT

! Loop through soil layers from bottom to top
              DO KK=NSOIL,1,-1

! Check for saturated layers
! Add additional logical check to ensure water is 'available' for routing,
!  (i.e. not 'frozen' or otherwise immobile)
!                IF (SMCRT(I,J,KK).GE.SMCMAXRT(I,J,KK).AND.SMCMAXRT(I,J,KK) &
!                  .GT.SMCWLTRT(I,J,KK)) THEN
                IF ( (SMCRT(I,J,KK).GE.SMCREFRT(I,J,KK)).AND.(SMCREFRT(I,J,KK) &
                  .GT.SMCWLTRT(I,J,KK)) ) THEN
! Add additional check to ensure saturation from bottom up only...8/8/05
                  IF((SATLYRCHK(I,J).EQ.KK+1) .OR. (KK.EQ.NSOIL) ) SATLYRCHK(I,J) = KK
                END IF

              END DO


! Designate ZWATTABLRT based on highest sat. layer and
! Define amount of water avail for subsfc routing on each gridcell (CWATAVAIL)
!  note: using a 'field capacity' value of SMCREF as lower limit...

              IF (SATLYRCHK(I,J).ne.0) then
                IF (SATLYRCHK(I,J).ne.1) then  ! soil column is partially sat.
                  ZWATTABLRT(I,J) = -ZSOIL(SATLYRCHK(I,J)-1)
!DJG 2/16/2016 fix                  DO KK=SATLYRCHK(I,J),NSOIL
!old                   CWATAVAIL(I,J) = (SMCRT(I,J,SATLYRCHK(I,J))-&
!old                                    SMCREFRT(I,J,SATLYRCHK(I,J))) * &
!old                                    (ZSOIL(SATLYRCHK(I,J)-1)-ZSOIL(NSOIL))
!DJG 2/16/2016 fix                    CWATAVAIL(I,J) = CWATAVAIL(I,J)+(SMCRT(I,J,KK)- &
!DJG 2/16/2016 fix                                     SMCREFRT(I,J,KK))*SLDPTH(KK)
!DJG 2/16/2016 fix                  END DO


                ELSE  ! soil column is fully saturated to sfc.
                  ZWATTABLRT(I,J) = 0.
!DJG 2/16/2016 fix                  DO KK=SATLYRCHK(I,J),NSOIL
!DJG 2/16/2016 fix                    CWATAVAIL(I,J) = (SMCRT(I,J,KK)-SMCREFRT(I,J,KK))*SLDPTH(KK)
!DJG 2/16/2016 fix                  END DO
                END IF
!DJG 2/16/2016 fix accumulation of CWATAVAIL...
                  DO KK=SATLYRCHK(I,J),NSOIL
                    CWATAVAIL(I,J) = CWATAVAIL(I,J)+(SMCRT(I,J,KK)- &
                                     SMCREFRT(I,J,KK))*SLDPTH(KK)
                  END DO
              ELSE  ! no saturated layers...
                ZWATTABLRT(I,J) = -ZSOIL(NSOIL)
                SATLYRCHK(I,J) = NSOIL + 1
              END IF


	   END DO
         END DO


!DJG ----------------------------------------------------------------
  END SUBROUTINE FINDZWAT 
!DJG ----------------------------------------------------------------


!DJG ----------------------------------------------------------------
!DJG ----------------------------------------------------------------
!DJG     SUBROUTINE ROUTE_SUBSURFACE2
!DJG ----------------------------------------------------------------

          SUBROUTINE ROUTE_SUBSURFACE2(                                 &
                dist,z,qsub,sox,soy,                                   &
                latksat,soldep,XX,YY,QSUBDRY,QSUBDRYT,CWATAVAIL,   &
                SUBDT)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Subroutine to route subsurface flow through the watershed
!DJG ----------------------------------------------------------------
!
!  Called from: main.f (Noah_router_driver)
!
!  Returns: qsub=DQSUB   which in turn becomes SUBFLO in head calc.
!
!  Created:    D. Gochis                           3/27/03
!              Adaptded from Wigmosta, 1994
!
!  Modified:   D. Gochis                           1/05/04
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef MPP_LAND
        use module_mpp_land, only: left_id,down_id,right_id,&
               up_id,mpp_land_com_real,MPP_LAND_UB_COM, &
               MPP_LAND_LR_COM,mpp_land_com_integer
#endif

        IMPLICIT NONE


!! Declare Passed variables

        INTEGER, INTENT(IN) :: XX,YY

!! Declare passed arrays

        REAL, INTENT(IN), DIMENSION(XX,YY) :: z
        REAL, INTENT(IN), DIMENSION(XX,YY) :: sox
        REAL, INTENT(IN), DIMENSION(XX,YY) :: soy
        REAL, INTENT(IN), DIMENSION(XX,YY) :: latksat
        REAL, INTENT(IN), DIMENSION(XX,YY) :: CWATAVAIL
        REAL, INTENT(IN), DIMENSION(XX,YY) :: soldep
        REAL, INTENT(OUT), DIMENSION(XX,YY) :: qsub
        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: QSUBDRY
        REAL, INTENT(INOUT) :: QSUBDRYT
        REAL, INTENT(IN) :: SUBDT
        real, intent(in), dimension(xx,yy,9) :: dist 

!!! Declare Local Variables

        REAL :: dzdx,dzdy,beta,gamma
        REAL :: qqsub,hh,ksat, gsize

        INTEGER :: i,j
!!! Initialize variables
        REAL, PARAMETER :: nexp=1.0      ! local power law exponent
        qsub = 0.                        ! initialize flux = 0. !DJG 5 May 2014

!yw        soldep = 2.
        

! Begin Subsurface routing

!!! Loop to route water in x-direction
        do j=1,YY
          do i=1,XX
! check for boundary grid point?
          if (i.eq.XX) GOTO 998
          gsize = dist(i,j,3)

          dzdx= (z(i,j) - z(i+1,j))/gsize
          beta=sox(i,j) + dzdx + 1E-30
          if (abs(beta) .lt. 1E-20) beta=1E-20
          if (beta.lt.0) then
!yw            hh=(1-(z(i+1,j)/soldep(i,j)))**nexp
            hh=(1-(z(i+1,j)/soldep(i+1,j)))**nexp
! Change later to use mean Ksat of two cells
            ksat=latksat(i+1,j)
          else
            hh=(1-(z(i,j)/soldep(i,j)))**nexp
            ksat=latksat(i,j)
          end if

          if (hh .lt. 0.) then
            print *, "hsub<0 at gridcell...", i,j,hh,z(i+1,j),z(i,j), &
                      soldep(i,j),nexp
            call hydro_stop("In ROUTE_SUBSURFACE2() - hsub<0 at gridcell")
          end if

!Err. tan slope          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*tan(beta)
!AD_CHANGE: beta is already a slope so no tan (consistent with ROUTE_SUBSURFACE1) 
          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*beta
!DJG lacks tan(beta) of original Wigmosta version          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*beta

          qqsub = gamma * hh
          qsub(i,j) = qsub(i,j) + qqsub
          qsub(i+1,j) = qsub(i+1,j) - qqsub

! Boundary adjustments
#ifdef MPP_LAND
          if ((i.eq.1).AND.(beta.lt.0.).and.(left_id.lt.0)) then
#else
          if ((i.eq.1).AND.(beta.lt.0.)) then
#endif
            qsub(i,j) = qsub(i,j) - qqsub
            QSUBDRY(i,j) = QSUBDRY(i,j) - qqsub
            QSUBDRYT = QSUBDRYT - qqsub
#ifdef MPP_LAND
          else if ((i.eq.(xx-1)).AND.(beta.gt.0.) &
              .and.(right_id.lt.0) ) then
#else
          else if ((i.eq.(xx-1)).AND.(beta.gt.0.)) then
#endif
            qsub(i+1,j) = qsub(i+1,j) + qqsub
            QSUBDRY(i+1,j) = QSUBDRY(i+1,j) + qqsub
            QSUBDRYT = QSUBDRYT + qqsub
          end if

998       continue

!! End loop to route sfc water in x-direction
          end do
        end do

#ifdef MPP_LAND
       call MPP_LAND_LR_COM(qsub,XX,YY,99)
       call MPP_LAND_LR_COM(QSUBDRY,XX,YY,99)
#endif


!!! Loop to route water in y-direction
        do j=1,YY
          do i=1,XX
! check for boundary grid point?
          if (j.eq.YY) GOTO 999
          gsize = dist(i,j,1)

          dzdy= (z(i,j) - z(i,j+1))/gsize
          beta=soy(i,j) + dzdy + 1E-30
          if (abs(beta) .lt. 1E-20) beta=1E-20
          if (beta.lt.0) then
!yw            hh=(1-(z(i,j+1)/soldep(i,j)))**nexp
            hh=(1-(z(i,j+1)/soldep(i,j+1)))**nexp
            ksat=latksat(i,j+1)
          else
            hh=(1-(z(i,j)/soldep(i,j)))**nexp
            ksat=latksat(i,j)
          end if

          if (hh .lt. 0.) GOTO 999

!Err. tan slope          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*tan(beta)
          gamma=-1.*((gsize*ksat*soldep(i,j))/nexp)*beta

          qqsub = gamma * hh
          qsub(i,j) = qsub(i,j) + qqsub
          qsub(i,j+1) = qsub(i,j+1) - qqsub

! Boundary adjustments

#ifdef MPP_LAND
          if ((j.eq.1).AND.(beta.lt.0.).and.(down_id.lt.0)) then
#else
          if ((j.eq.1).AND.(beta.lt.0.)) then
#endif
            qsub(i,j) = qsub(i,j) - qqsub
            QSUBDRY(i,j) = QSUBDRY(i,j) - qqsub
            QSUBDRYT = QSUBDRYT - qqsub
#ifdef MPP_LAND
          else if ((j.eq.(yy-1)).AND.(beta.gt.0.)  &
                .and. (up_id.lt.0) ) then
#else
          else if ((j.eq.(yy-1)).AND.(beta.gt.0.)) then
#endif
            qsub(i,j+1) = qsub(i,j+1) + qqsub
            QSUBDRY(i,j+1) = QSUBDRY(i,j+1) + qqsub
            QSUBDRYT = QSUBDRYT + qqsub
          end if

999       continue

!! End loop to route sfc water in y-direction
          end do
        end do

#ifdef MPP_LAND
       call MPP_LAND_UB_COM(qsub,XX,YY,99)
       call MPP_LAND_UB_COM(QSUBDRY,XX,YY,99)
#endif

        return
!DJG------------------------------------------------------------
        end subroutine ROUTE_SUBSURFACE2
!DJG------------------------------------------------------------



!DJG ------------------------------------------------
!DJG   SUBROUTINE OV_RTNG
!DJG ------------------------------------------------

	SUBROUTINE OV_RTNG(DT,DTRT_TER,IXRT,JXRT,INFXSUBRT,      &
          SFCHEADSUBRT,DHRT,CH_NETRT,RETDEPRT,OVROUGHRT,      &
          QSTRMVOLRT,QBDRYRT,QSTRMVOLTRT,QBDRYTRT,SOXRT,     &
          SOYRT,dist,LAKE_MSKRT,LAKE_INFLORT,LAKE_INFLOTRT,  &
          SO8RT,SO8RT_D,rt_option,q_sfcflx_x,q_sfcflx_y)

!yyww 
#ifdef MPP_LAND
        use module_mpp_land, only: left_id,down_id,right_id, &
              up_id,mpp_land_com_real, my_id, &
             mpp_land_sync
#endif

	IMPLICIT NONE

!DJG --------DECLARATIONS----------------------------

	INTEGER, INTENT(IN)			:: IXRT,JXRT
	REAL, INTENT(IN)			:: DT,DTRT_TER

	INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETRT
	INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: LAKE_MSKRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)	:: INFXSUBRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)	:: SOXRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)	:: SOYRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT,9):: dist 
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)	:: RETDEPRT
	REAL, INTENT(IN), DIMENSION(IXRT,JXRT)	:: OVROUGHRT

	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)	:: SFCHEADSUBRT
	REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)	:: DHRT

	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: QSTRMVOLRT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: LAKE_INFLORT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: QBDRYRT
	REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: q_sfcflx_x,q_sfcflx_y
	REAL, INTENT(INOUT)     :: QSTRMVOLTRT,QBDRYTRT,LAKE_INFLOTRT
        REAL, INTENT(IN), DIMENSION(IXRT,JXRT,8)  :: SO8RT

!DJG Local Variables

	INTEGER :: KRT,I,J,ct

	REAL, DIMENSION(IXRT,JXRT)	:: INFXS_FRAC
	REAL	:: DT_FRAC,SUM_INFXS,sum_head
        INTEGER SO8RT_D(IXRT,JXRT,3), rt_option
	
	


!DJG ----------------------------------------------------------------------
! DJG BEGIN 1-D or 2-D OVERLAND FLOW ROUTING LOOP
!DJG ---------------------------------------------------------------------
!DJG  Loop over 'routing time step'
!DJG  Compute the number of time steps based on NOAH DT and routing DTRT_TER

       DT_FRAC=INT(DT/DTRT_TER)

#ifdef HYDRO_D
       write(6,*) "OV_RTNG  DT_FRAC, DT, DTRT_TER",DT_FRAC, DT, DTRT_TER
       write(6,*) "IXRT, JXRT = ",ixrt,jxrt
#endif

!DJG NOTE: Applying all infiltration excess water at once then routing
!DJG       Pre-existing SFHEAD gets combined with Precip. in the
!DJG       calculation of INFXS1 during subroutine SRT.f.
!DJG debug


!DJG Assign all infiltration excess to surface head...
            SFCHEADSUBRT=INFXSUBRT

!DJG Divide infiltration excess over all routing time-steps
!	     INFXS_FRAC=INFXSUBRT/(DT/DTRT_TER)

!DJG Set flux accumulation fields to 0. before each loop...
      q_sfcflx_x = 0.
      q_sfcflx_y = 0.
      ct =0


!DJG Execute routing time-step loop...


      DO KRT=1,DT_FRAC

        DO J=1,JXRT
          DO I=1,IXRT

!DJG Removed 4_29_05, sfhead now updated in route_overland subroutine...
!           SFCHEADSUBRT(I,J)=SFCHEADSUBRT(I,J)+DHRT(I,J)
!!           SFCHEADSUBRT(I,J)=SFCHEADSUBRT(I,J)+DHRT(I,J)+INFXS_FRAC(I,J)
!           DHRT(I,J)=0.

!DJG ERROR Check...

	   IF (SFCHEADSUBRT(I,J).lt.0.) THEN 
#ifdef HYDRO_D
		print *, "ywcheck 2 ERROR!!!: Neg. Surface Head Value at (i,j):",    &
                    i,j,SFCHEADSUBRT(I,J)
                print *, "RETDEPRT(I,J) = ",RETDEPRT(I,J), "KRT=",KRT
                print *, "INFXSUBRT(i,j)=",INFXSUBRT(i,j)
                print *, "jxrt=",jxrt," ixrt=",ixrt
#endif
           END IF

!DJG Remove surface water from channel cells
!DJG Channel inflo cells specified as nonzeros from CH_NET
!DJG 9/16/04  Channel Extractions Removed until stream model implemented...



!yw            IF (CH_NETRT(I,J).ne.-9999) THEN
           IF (CH_NETRT(I,J).ge.0) THEN
             ct = ct +1

!DJG Temporary test to up the retention depth of channel grid cells to 'soak' 
!more water into valleys....set retdep = retdep*100 (=5 mm)

!	     RETDEPRT(I,J) = RETDEPRT(I,J) * 100.0    !DJG TEMP HARDWIRE!!!!
!	     RETDEPRT(I,J) = 10.0    !DJG TEMP HARDWIRE!!!!

             IF (SFCHEADSUBRT(I,J).GT.RETDEPRT(I,J)) THEN
!!               QINFLO(CH_NET(I,J)=QINFLO(CH_NET(I,J)+SFCHEAD(I,J) - RETDEPRT(I,J)
               QSTRMVOLTRT = QSTRMVOLTRT + (SFCHEADSUBRT(I,J) - RETDEPRT(I,J))
               QSTRMVOLRT(I,J) = QSTRMVOLRT(I,J)+SFCHEADSUBRT(I,J)-RETDEPRT(I,J)

             ! if(QSTRMVOLRT(I,J) .gt. 0) then 
             !     print *, "QSTRVOL GT 0", QSTRMVOLRT(I,J),I,J 
             !  endif

               SFCHEADSUBRT(I,J) = RETDEPRT(I,J)
             END IF
           END IF

!DJG Lake inflow withdrawl from surface head...(4/29/05)
           

           IF (LAKE_MSKRT(I,J).gt.0) THEN
             IF (SFCHEADSUBRT(I,J).GT.RETDEPRT(I,J)) THEN
               LAKE_INFLOTRT = LAKE_INFLOTRT + (SFCHEADSUBRT(I,J) - RETDEPRT(I,J))
               LAKE_INFLORT(I,J) = LAKE_INFLORT(I,J)+SFCHEADSUBRT(I,J)-RETDEPRT(I,J)
               SFCHEADSUBRT(I,J) = RETDEPRT(I,J)
              
             END IF
           END IF



         END DO
        END DO

!yw check         call MPP_LAND_COM_REAL(QSTRMVOLRT,IXRT,JXRT,99)
!DJG----------------------------------------------------------------------
!DJG CALL OVERLAND FLOW ROUTING SUBROUTINE
!DJG----------------------------------------------------------------------

!DJG Debug...


           if(rt_option .eq. 1) then
              CALL ROUTE_OVERLAND1(DTRT_TER,dist,SFCHEADSUBRT,DHRT,SOXRT,   &
		SOYRT,RETDEPRT,OVROUGHRT,IXRT,JXRT,QBDRYRT,QBDRYTRT,    & 
                SO8RT,SO8RT_D,q_sfcflx_x,q_sfcflx_y)
            else
              CALL ROUTE_OVERLAND2(DTRT_TER,dist,SFCHEADSUBRT,DHRT,SOXRT,   &
                  SOYRT,RETDEPRT,OVROUGHRT,IXRT,JXRT,QBDRYRT,QBDRYTRT,  &
                  q_sfcflx_x,q_sfcflx_y)    
            end if
             
        END DO          ! END routing time steps

#ifdef HYDRO_D
 	print *, "End of OV_routing call..."
#endif

!----------------------------------------------------------------------
! END OVERLAND FLOW ROUTING LOOP
!     CHANNEL ROUTING TO FOLLOW 
!----------------------------------------------------------------------

!DJG ----------------------------------------------------------------
  END SUBROUTINE OV_RTNG 
!DJG ----------------------------------------------------------------

!DJG     SUBROUTINE ROUTE_OVERLAND1
!DJG ----------------------------------------------------------------

          SUBROUTINE ROUTE_OVERLAND1(dt,                                &
     &          gsize,h,qsfc,sox,soy,                                   &
     &     retent_dep,dist_rough,XX,YY,QBDRY,QBDRYT,SO8RT,SO8RT_D,      &
     &     q_sfcflx_x,q_sfcflx_y)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Subroutine to route excess rainfall over the watershed
!     using a 2d diffusion routing scheme.
!
!  Called from: main.f
!
!      Will try to formulate this to be called from NOAH
!
!  Returns: qsfc=DQOV   which in turn becomes DH in head calc.
!
!  Created:  Adaptded from CASC2D source code
!  NOTE: dh from original code has been replaced by qsfc
!        dhh replaced by qqsfc
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#ifdef MPP_LAND
        use module_mpp_land, only: left_id,down_id,right_id, &
              up_id,mpp_land_com_real, my_id, mpp_land_com_real8,&
             mpp_land_sync
#endif

        IMPLICIT NONE


!! Declare Passed variables

        INTEGER, INTENT(IN) :: XX,YY
        REAL, INTENT(IN) :: dt, gsize(xx,yy,9)

!! Declare passed arrays

        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: h
        REAL, INTENT(IN), DIMENSION(XX,YY) :: qsfc
        REAL, INTENT(IN), DIMENSION(XX,YY) :: sox
        REAL, INTENT(IN), DIMENSION(XX,YY) :: soy
        REAL, INTENT(IN), DIMENSION(XX,YY) :: retent_dep
        REAL, INTENT(IN), DIMENSION(XX,YY) :: dist_rough
        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: QBDRY
        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: q_sfcflx_x, q_sfcflx_y
        REAL, INTENT(INOUT) :: QBDRYT
        REAL, INTENT(IN), DIMENSION(XX,YY,8) :: SO8RT
        REAL*8, DIMENSION(XX,YY) :: QBDRY_tmp, DH
        REAL*8, DIMENSION(XX,YY) :: DH_tmp

!!! Declare Local Variables

        REAL :: dhdx,dhdy,alfax,alfay
        REAL :: hh53,qqsfc,hh,dt_new,hmax
        REAL :: sfx,sfy
        REAL :: tmp_adjust

        INTEGER :: i,j
        REAL IXX8,IYY8
        INTEGER  IXX0,JYY0,index, SO8RT_D(XX,YY,3)
        REAL  tmp_gsize,hsum

!!! Initialize variables



!!! Begin Routing of Excess Rainfall over the Watershed

        DH=0.
        DH_tmp=0.
        QBDRY_tmp =0.

!!! Loop to route water
        do j=2,YY-1
          do i=2,XX-1
          if (h(I,J).GT.retent_dep(I,J)) then 
             IXX0 = SO8RT_D(i,j,1)
             JYY0 = SO8RT_D(i,j,2)
             index = SO8RT_D(i,j,3)
             tmp_gsize = 1.0/gsize(i,j,index)
             sfx = so8RT(i,j,index)-(h(IXX0,JYY0)-h(i,j))*0.001*tmp_gsize
             hmax = h(i,j)*0.001  !Specify max head for mass flux limit...
             if(sfx .lt. 1E-20) then
               call GETMAX8DIR(IXX0,JYY0,I,J,H,RETENT_DEP,so8rt,gsize(i,j,:),sfx,XX,YY)
             end if
             if(IXX0 > 0) then  ! do the rest if the lowest grid can be found.
                 if(sfx .lt. 1E-20) then
#ifdef HYDRO_D
                      print*, "Message: sfx reset to 1E-20. sfx =",sfx
                      print*, "i,j,index,IXX0,JYY0",i,j,index,IXX0,JYY0
                      print*, "so8RT(i,j,index), h(IXX0,JYY0), h(i,j), gsize(i,j,index) ", &
                         so8RT(i,j,index), h(IXX0,JYY0), h(i,j), gsize(i,j,index)
#endif
                      sfx = 1E-20
                 end if
                 alfax = sqrt(sfx) / dist_rough(i,j) 
                 hh=(h(i,j)-retent_dep(i,j)) * 0.001
                 hh53=hh**(5./3.)

! Calculate q-flux...
                 qqsfc = alfax*hh53*dt * tmp_gsize

!Courant check (simple mass limit on overland flow)...
                 if (qqsfc.ge.(hmax*dt*tmp_gsize)) qqsfc = hmax*dt*tmp_gsize

! Accumulate directional fluxes on routing subgrid...
                 if (IXX0.gt.i) then
                   q_sfcflx_x(I,J) = q_sfcflx_x(I,J) + qqsfc * &
                         (1.0 - 0.5 * (ABS(j-JYY0)))
                 else if (IXX0.lt.i) then
                   q_sfcflx_x(I,J) = q_sfcflx_x(I,J) - 1.0 * &
                         qqsfc * (1.0 - 0.5 * (ABS(j-JYY0)))
                 else
                   q_sfcflx_x(I,J) = q_sfcflx_x(I,J) + 0.
                 end if
                 if (JYY0.gt.j) then
                   q_sfcflx_y(I,J) = q_sfcflx_y(I,J) + qqsfc * &
                          (1.0 - 0.5 * (ABS(i-IXX0)))
                 elseif (JYY0.lt.j) then
                   q_sfcflx_y(I,J) = q_sfcflx_y(I,J) - 1.0 * &
                          qqsfc * (1.0 - 0.5 * (ABS(i-IXX0)))
                 else
                   q_sfcflx_y(I,J) = q_sfcflx_y(I,J) + 0.
                 end if


!DJG put adjustment in for (h) due to qqsfc

!yw changed as following:
                 tmp_adjust=qqsfc*1000
                 if((h(i,j) - tmp_adjust) <0 )  then
#ifdef HYDRO_D
                   print*, "Error Warning: surface head is negative:  ",i,j,ixx0,jyy0, &
                       h(i,j) - tmp_adjust
#endif
                     tmp_adjust = h(i,j)
                 end if
 	         DH(i,j) = DH(i,j)-tmp_adjust
                 DH_tmp(ixx0,jyy0) = DH_tmp(ixx0,jyy0) + tmp_adjust
      !yw end change
                  
      !DG Boundary adjustments here
            !DG Constant Flux Condition
#ifdef MPP_LAND
      if( ((ixx0.eq.XX).and.(right_id .lt. 0)) .or. &
          ((ixx0.eq.1) .and.(left_id  .lt. 0)) .or. &
          ((jyy0.eq.1) .and.(down_id  .lt. 0)) .or. &
          ((JYY0.eq.YY).and.(up_id    .lt. 0)) ) then 
!              QBDRY_tmp(IXX0,JYY0)=QBDRY_tmp(IXX0,JYY0) - qqsfc*1000.
#else
                if ((ixx0.eq.XX).or.(ixx0.eq.1).or.(jyy0.eq.1)   &
                     .or.(JYY0.eq.YY )) then
!                     QBDRY(IXX0,JYY0)=QBDRY(IXX0,JYY0) - qqsfc*1000.
#endif
                     QBDRY_tmp(IXX0,JYY0)=QBDRY_tmp(IXX0,JYY0) - qqsfc*1000.
                     QBDRYT=QBDRYT - qqsfc
                     DH_tmp(IXX0,JYY0)= DH_tmp(IXX0,JYY0)-tmp_adjust
                end if
             end if
!! End loop to route sfc water 
          end if
          end do
        end do

#ifdef MPP_LAND
! use double precision to solve the underflow problem.
       call MPP_LAND_COM_REAL8(DH_tmp,XX,YY,1)
       call MPP_LAND_COM_REAL8(QBDRY_tmp,XX,YY,1)
#endif
       QBDRY = QBDRY + QBDRY_tmp
       DH = DH+DH_tmp 

#ifdef MPP_LAND
       call MPP_LAND_COM_REAL8(DH,XX,YY,99)
       call MPP_LAND_COM_REAL(QBDRY,XX,YY,99)
#endif

        H = H + DH

        return

!DJG ----------------------------------------------------------------------
        end subroutine ROUTE_OVERLAND1


!DJG ----------------------------------------------------------------
        SUBROUTINE GETMAX8DIR(IXX0,JYY0,I,J,H,RETENT_DEP,sox,tmp_gsize,max,XX,YY)
          implicit none
          INTEGER:: IXX0,JYY0,IXX8,JYY8, XX, YY
          INTEGER, INTENT(IN) :: I,J

          REAL,INTENT(IN) :: H(XX,YY),RETENT_DEP(XX,YY),sox(XX,YY,8),tmp_gsize(9)
          REAL  max
          IXX0 = -1
          max = 0
          if (h(I,J).LE.retent_dep(I,J)) return

          IXX8 = I
          JYY8 = J+1
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,1),IXX0,JYY0,max,tmp_gsize(1),XX,YY)

          IXX8 = I+1
          JYY8 = J+1
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,2),IXX0,JYY0,max,tmp_gsize(2),XX,YY)

          IXX8 = I+1
          JYY8 = J
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,3),IXX0,JYY0,max,tmp_gsize(3),XX,YY)

          IXX8 = I+1
          JYY8 = J-1
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,4),IXX0,JYY0,max,tmp_gsize(4),XX,YY)

          IXX8 = I
          JYY8 = J-1
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,5),IXX0,JYY0,max,tmp_gsize(5),XX,YY)

          IXX8 = I-1
          JYY8 = J-1
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,6),IXX0,JYY0,max,tmp_gsize(6),XX,YY)

          IXX8 = I-1
          JYY8 = J
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,7),IXX0,JYY0,max,tmp_gsize(7),XX,YY)

          IXX8 = I-1
          JYY8 = J+1
          call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,8),IXX0,JYY0,max,tmp_gsize(8),XX,YY)
        RETURN
        END SUBROUTINE GETMAX8DIR

        SUBROUTINE GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox   &
            ,IXX0,JYY0,max,tmp_gsize,XX,YY)
        implicit none
        integer,INTENT(INOUT) ::IXX0,JYY0
        INTEGER, INTENT(IN) :: I,J,IXX8,JYY8,XX,YY
        REAL,INTENT(IN) :: H(XX,YY),RETENT_DEP(XX,YY),sox(XX,YY)
        REAL, INTENT(INOUT) ::max
        real, INTENT(IN) :: tmp_gsize
        real :: sfx

             sfx = sox(i,j)-(h(IXX8,JYY8)-h(i,j))*0.001/tmp_gsize
             if(sfx .le. 0 ) return
             if(max < sfx ) then
                   IXX0 = IXX8
                   JYY0 = JYY8
                   max = sfx
             end if

        END SUBROUTINE GET8DIR
!DJG ----------------------------------------------------------------
!DJG     SUBROUTINE ROUTE_SUBSURFACE1
!DJG ----------------------------------------------------------------

          SUBROUTINE ROUTE_SUBSURFACE1(                                 &
                dist,z,qsub,sox,soy,                                   &
                latksat,soldep,XX,YY,QSUBDRY,QSUBDRYT,SO8RT,SO8RT_D,    &
                CWATAVAIL,SUBDT)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Subroutine to route subsurface flow through the watershed
!
!  Called from: main.f (Noah_router_driver)
!
!  Returns: qsub=DQSUB   which in turn becomes SUBFLO in head calc.
!
!  Created:    D. Gochis                           3/27/03
!              Adaptded from Wigmosta, 1994
!
!  Modified:   D. Gochis                           1/05/04
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#ifdef MPP_LAND
        use module_mpp_land, only: left_id,down_id,right_id,&
           up_id,mpp_land_com_real8,my_id,mpp_land_com_real
#endif

        IMPLICIT NONE


!! Declare Passed variables

        INTEGER, INTENT(IN) :: XX,YY

!! Declare passed arrays

        REAL, INTENT(IN), DIMENSION(XX,YY) :: z
        REAL, INTENT(IN), DIMENSION(XX,YY) :: sox
        REAL, INTENT(IN), DIMENSION(XX,YY) :: soy
        REAL, INTENT(IN), DIMENSION(XX,YY) :: latksat
        REAL, INTENT(IN), DIMENSION(XX,YY) :: CWATAVAIL
        REAL, INTENT(IN), DIMENSION(XX,YY) :: soldep
        REAL, INTENT(OUT), DIMENSION(XX,YY) :: qsub
        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: QSUBDRY
        REAL, INTENT(INOUT) :: QSUBDRYT
        REAL*8, DIMENSION(XX,YY) :: qsub_tmp,QSUBDRY_tmp
!yw        INTEGER, INTENT(OUT) :: flag
        REAL, INTENT(IN) :: dist(xx,yy,9),SUBDT

!!! Declare Local Variables

        REAL :: dzdx,dzdy,beta,gamma
        REAL :: qqsub,hh,ksat

        REAL :: SO8RT(XX,YY,8)
        INTEGER :: SO8RT_D(XX,YY,3), rt_option
        

!!! Initialize variables

        REAL, PARAMETER :: nexp=1.0      ! local power law exponent
        integer IXX0,JYY0,index,i,j
        real tmp_gsize

!     temporary set it to be 2. Should be passed in.
!yw        soldep = 2.
! Begin Subsurface routing



!!! Loop to route water in x-direction
        qsub_tmp = 0.
        QSUBDRY_tmp = 0.

#ifdef HYDRO_D
        write(6,*) "call subsurface routing xx= , yy =", yy, xx
#endif

        do j=2,YY-1
          do i=2,XX-1


        if(i.ge.2.AND.i.le.XX-1.AND.j.ge.2.AND.j.le.YY-1) then !if grdcl chk
! check for boundary grid point?
          IXX0 = SO8RT_D(i,j,1)
          JYY0 = SO8RT_D(i,j,2)

          index = SO8RT_D(i,j,3)

            if(dist(i,j,index) .le. 0) then
               write(6,*) "FATAL ERROR: dist(i,j,index) is <= zero "   
               call hydro_stop("In ROUTE_SUBSURFACE1() - dist(i,j,index) is <= zero ")
            endif
            if(soldep(i,j) .eq. 0) then
               call hydro_stop("In ROUTE_SUBSURFACE1() - soldep is = zero")
            endif

          tmp_gsize = 1.0/dist(i,j,index)

       
          dzdx= (z(i,j) - z(IXX0,JYY0) )* tmp_gsize
          beta=so8RT(i,j,index) + dzdx 

          if(beta .lt. 1E-20 ) then   !if-then for direction...
            call GETSUB8(IXX0,JYY0,I,J,Z,so8rt,dist(i,j,:),beta,XX,YY)
          end if
          if(beta .gt. 0) then            !if-then for flux calc 
              if(beta .lt. 1E-20 ) then
#ifdef HYDRO_D
                   print*, "Message: beta need to be reset to 1E-20. beta = ",beta
#endif
                   beta = 1E-20
              end if

! do the rest if the lowest grid can be found.
              hh=(1-(z(i,j)/soldep(i,j)))**nexp
              ksat=latksat(i,j)

              if (hh .lt. 0.) then
                 print *, "hsub<0 at gridcell...", i,j,hh,z(i+1,j),z(i,j), &
                      soldep(i,j)
                 call hydro_stop("In ROUTE_SUBSURFACE1() - hsub<0 at gridcell ") 
              end if

!err. tan slope     gamma=-1.0*((gsize*ksat*soldep(i,j))/nexp)*tan(beta)
              gamma=-1.0*((dist(i,j,index)*ksat*soldep(i,j))/nexp)*beta
              qqsub = gamma * hh

! ADCHANGE: Moved this water available constraint from outside qsub calc loop to inside
!           to better account for adjustments to adjacent cells
              if( qqsub .le. 0 .and. CWATAVAIL(i,j).lt.ABS(qqsub)/dist(i,j,9)*SUBDT) THEN
                  qqsub = -1.0*CWATAVAIL(i,j)*dist(i,j,9)/SUBDT
              end if

              qsub(i,j) = qsub(i,j) + qqsub
              qsub_tmp(ixx0,jyy0) = qsub_tmp(ixx0,jyy0) - qqsub

!!DJG Error Checks...
              if(qqsub .gt. 0) then
                    print*, "FATAL ERROR: qqsub should be negative, qqsub =",qqsub,&
                       "gamma=",gamma,"hh=",hh,"beta=",beta,"dzdx=",dzdx,&
                       "so8RT=",so8RT(i,j,index),"latksat=",ksat, &
                       "tan(beta)=",tan(beta),i,j,z(i,j),z(IXX0,JYY0)
                    print*, "ixx0=",ixx0, "jyy0=",jyy0
                    print*, "soldep =", soldep(i,j), "nexp=",nexp
                 call hydro_stop("In ROUTE_SUBSURFACE1() - qqsub should be negative") 
              end if




! Boundary adjustments
#ifdef MPP_LAND
      if( ((ixx0.eq.XX).and.(right_id .lt. 0)) .or. &
          ((ixx0.eq.1) .and.(left_id  .lt. 0)) .or. &
          ((jyy0.eq.1) .and.(down_id  .lt. 0)) .or. &
          ((JYY0.eq.YY).and.(up_id    .lt. 0)) ) then 
#else
              if ((ixx0.eq.1).or.(ixx0.eq.xx).or.(jyy0.eq.1).or.(jyy0.eq.yy)) then
#endif
                qsub_tmp(ixx0,jyy0) = qsub_tmp(ixx0,jyy0) + qqsub
                QSUBDRY_tmp(ixx0,jyy0) = QSUBDRY_tmp(ixx0,jyy0) + qqsub

                QSUBDRYT = QSUBDRYT + qqsub
              end if

998           continue

!! End loop to route sfc water in x-direction
      end if  !endif for flux calc

          endif   !! Endif for gridcell check...


          end do  !endif for i-dim
!CRNT debug          if(flag.eq.-99) exit !exit loop for courant violation...
        end do   !endif for j-dim

#ifdef MPP_LAND

       call MPP_LAND_COM_REAL8(qsub_tmp,XX,YY,1)
       call MPP_LAND_COM_REAL8(QSUBDRY_tmp,XX,YY,1)
#endif
       qsub = qsub + qsub_tmp
       QSUBDRY= QSUBDRY + QSUBDRY_tmp 

!ADNOTE: Moved this check to inside qsub calc loop, so no need for additional loop
!        do j=2,YY-1
!          do i=2,XX-1
!            if(dist(i,j,9) .le. 0) then
!               call hydro_stop("In ROUTE_SUBSURFACE1() - dist(i,j,9) is <= zero")
!            endif
!!DJG Feb 16, 2016...comment out to debug...line is identical to line 255
!!            if(CWATAVAIL(i,j).lt.ABS(qsub(i,j))/dist(i,j,9)*SUBDT) THEN
!!              qsub(i,j) = -1.0*CWATAVAIL(i,j)
!!            end if
!          end do
!        end do

#ifdef MPP_LAND
       call MPP_LAND_COM_REAL(qsub,XX,YY,99)
       call MPP_LAND_COM_REAL(QSUBDRY,XX,YY,99)
#endif


        return
!DJG------------------------------------------------------------
        end subroutine ROUTE_SUBSURFACE1
!DJG------------------------------------------------------------

!DJG------------------------------------------------------------


      SUBROUTINE GETSUB8(IXX0,JYY0,I,J,Z,sox,tmp_gsize,max,XX,YY)
          implicit none
          INTEGER:: IXX0,JYY0,IXX8,JYY8, XX, YY
          INTEGER, INTENT(IN) :: I,J

          REAL,INTENT(IN) :: Z(XX,YY),sox(XX,YY,8),tmp_gsize(9)
          REAL  max
          max = -1

          IXX8 = I
          JYY8 = J+1
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,1),IXX0,JYY0,max,tmp_gsize(1),XX,YY)

          IXX8 = I+1
          JYY8 = J+1
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,2),IXX0,JYY0,max,tmp_gsize(2),XX,YY)

          IXX8 = I+1
          JYY8 = J
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,3),IXX0,JYY0,max,tmp_gsize(3),XX,YY)

          IXX8 = I+1
          JYY8 = J-1
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,4),IXX0,JYY0,max,tmp_gsize(4),XX,YY)

          IXX8 = I
          JYY8 = J-1
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,5),IXX0,JYY0,max,tmp_gsize(5),XX,YY)

          IXX8 = I-1
          JYY8 = J-1
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,6),IXX0,JYY0,max,tmp_gsize(6),XX,YY)

          IXX8 = I-1
          JYY8 = J
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,7),IXX0,JYY0,max,tmp_gsize(7),XX,YY)

          IXX8 = I-1
          JYY8 = J+1
          call GETSUB8DIR(IXX8,JYY8,I,J,Z,sox(:,:,8),IXX0,JYY0,max,tmp_gsize(8),XX,YY)
        RETURN
        END SUBROUTINE GETSUB8

        SUBROUTINE GETSUB8DIR(IXX8,JYY8,I,J,Z,sox,IXX0,JYY0,max,tmp_gsize,XX,YY)
        implicit none
        integer,INTENT(INOUT) ::IXX0,JYY0
        INTEGER, INTENT(IN) :: I,J,IXX8,JYY8,XX,YY
        REAL,INTENT(IN) :: Z(XX,YY),sox(XX,YY)
        REAL, INTENT(INOUT) ::max
        real, INTENT(IN) :: tmp_gsize
        real :: beta , dzdx

          dzdx= (z(i,j) - z(IXX0,JYY0) )/tmp_gsize
          beta=sox(i,j) + dzdx 
          if(max < beta ) then
                   IXX0 = IXX8
                   JYY0 = JYY8
                   max = beta 
          end if

        END SUBROUTINE GETSUB8DIR
!DJG ----------------------------------------------------------------------

!DJG     SUBROUTINE ROUTE_OVERLAND2
!DJG ----------------------------------------------------------------

          SUBROUTINE ROUTE_OVERLAND2 (dt,                               &
     &          dist,h,qsfc,sox,soy,                                   &
     &          retent_dep,dist_rough,XX,YY,QBDRY,QBDRYT,               &
     &          q_sfcflx_x,q_sfcflx_y)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Subroutine to route excess rainfall over the watershed
!     using a 2d diffusion routing scheme.
!
!  Called from: main.f
!
!      Will try to formulate this to be called from NOAH
!
!  Returns: qsfc=DQOV   which in turn becomes DH in head calc.
!
!  Created:  Adaptded from CASC2D source code
!  NOTE: dh from original code has been replaced by qsfc
!        dhh replaced by qqsfc
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef MPP_LAND
        use module_mpp_land, only: left_id,down_id,right_id,&
               up_id,mpp_land_com_real,MPP_LAND_UB_COM, &
               MPP_LAND_LR_COM,mpp_land_com_integer
#endif

        IMPLICIT NONE


!! Declare Passed variables

        real :: gsize
        INTEGER, INTENT(IN) :: XX,YY
        REAL, INTENT(IN) :: dt , dist(XX,YY,9)

!! Declare passed arrays

        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: h
        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: qsfc
        REAL, INTENT(IN), DIMENSION(XX,YY) :: sox
        REAL, INTENT(IN), DIMENSION(XX,YY) :: soy
        REAL, INTENT(IN), DIMENSION(XX,YY) :: retent_dep
        REAL, INTENT(IN), DIMENSION(XX,YY) :: dist_rough
        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: QBDRY
        REAL, INTENT(INOUT), DIMENSION(XX,YY) :: q_sfcflx_x,q_sfcflx_y
        REAL, INTENT(INOUT) :: QBDRYT
        REAL  :: DH(XX,YY)

!!! Declare Local Variables

        REAL :: dhdx,dhdy,alfax,alfay
        REAL :: hh53,qqsfc,hh,dt_new
        REAL :: sfx,sfy
        REAL :: tmp_adjust

        INTEGER :: i,j

!!! Initialize variables




!!! Begin Routing of Excess Rainfall over the Watershed


        DH = 0
!!! Loop to route water in x-direction
        do j=1,YY
          do i=1,XX


! check for boundary gridpoint?
          if (i.eq.XX) GOTO 998
           gsize = dist(i,j,3)


! check for detention storage?
          if (h(i,j).lt.retent_dep(i,j).AND.     &
              h(i+1,j).lt.retent_dep(i+1,j)) GOTO 998

          dhdx = (h(i+1,j)/1000. - h(i,j)/1000.) / gsize  ! gisze-(m),h-(mm)

          sfx = (sox(i,j)-dhdx+1E-30)
          if (abs(sfx).lt.1E-20) sfx=1E-20
          alfax = ((abs(sfx))**0.5)/dist_rough(i,j)
          if (sfx.lt.0.) then
              hh=(h(i+1,j)-retent_dep(i+1,j))/1000.
          else
              hh=(h(i,j)-retent_dep(i,j))/1000.
          end if

          if ((retent_dep(i,j).gt.0.).AND.(hh.le.0.)) GOTO 998
          if (hh.lt.0.) then
          GOTO 998
          end if

          hh53=hh**(5./3.)


! Calculate q-flux... (units (m))
          qqsfc = (sfx/abs(sfx))*alfax*hh53*dt/gsize
          q_sfcflx_x(I,J) = q_sfcflx_x(I,J) + qqsfc

!DJG put adjustment in for (h) due to qqsfc

!yw changed as following:
           tmp_adjust=qqsfc*1000
          if(tmp_adjust .le. 0 ) GOTO 998
           if((h(i,j) - tmp_adjust) <0 )  then
#ifdef HYDRO_D
               print*, "WARNING: surface head is negative:  ",i,j
#endif
               tmp_adjust = h(i,j)
           end if
           if((h(i+1,j) + tmp_adjust) <0) then 
#ifdef HYDRO_D
               print*, "WARNING: surface head is negative: ",i+1,j
#endif
               tmp_adjust = -1*h(i+1,j)
           end if
 	   Dh(i,j) = Dh(i,j)-tmp_adjust
           Dh(i+1,j) = Dh(i+1,j) + tmp_adjust
!yw end change



!DG Boundary adjustments here
!DG Constant Flux Condition
#ifdef MPP_LAND
          if ((i.eq.1).AND.(sfx.lt.0).and. & 
                (left_id .lt. 0) ) then
#else
          if ((i.eq.1).AND.(sfx.lt.0)) then
#endif
             Dh(i,j) = Dh(i,j) + qqsfc*1000.
            QBDRY(I,J)=QBDRY(I,J) + qqsfc*1000.
            QBDRYT=QBDRYT + qqsfc*1000.
#ifdef MPP_LAND
          else if ( (i.eq.(XX-1)).AND.(sfx.gt.0) &
             .and. (right_id .lt. 0) ) then
#else
          else if ((i.eq.(XX-1)).AND.(sfx.gt.0)) then
#endif
             tmp_adjust = qqsfc*1000.
             if(h(i+1,j).lt.tmp_adjust) tmp_adjust = h(i+1,j)
             Dh(i+1,j) = Dh(i+1,j) - tmp_adjust
!DJG Re-assign h(i+1) = 0.0 when <0.0 (from rounding/truncation error)
            QBDRY(I+1,J)=QBDRY(I+1,J) - tmp_adjust
            QBDRYT=QBDRYT - tmp_adjust
          end if


998     continue

!! End loop to route sfc water in x-direction
          end do
        end do

        H = H + DH
#ifdef MPP_LAND
       call MPP_LAND_LR_COM(H,XX,YY,99)
       call MPP_LAND_LR_COM(QBDRY,XX,YY,99)
#endif


        DH = 0
!!!! Loop to route water in y-direction
        do j=1,YY
          do i=1,XX

!! check for boundary grid point?
          if (j.eq.YY) GOTO 999
           gsize = dist(i,j,1)


!! check for detention storage?
          if (h(i,j).lt.retent_dep(i,j).AND.     & 
              h(i,j+1).lt.retent_dep(i,j+1)) GOTO 999

          dhdy = (h(i,j+1)/1000. - h(i,j)/1000.) / gsize

          sfy = (soy(i,j)-dhdy+1E-30)
          if (abs(sfy).lt.1E-20) sfy=1E-20
          alfay = ((abs(sfy))**0.5)/dist_rough(i,j)
          if (sfy.lt.0.) then
              hh=(h(i,j+1)-retent_dep(i,j+1))/1000.
          else
              hh=(h(i,j)-retent_dep(i,j))/1000.
          end if

          if ((retent_dep(i,j).gt.0.).AND.(hh.le.0.)) GOTO 999
          if (hh.lt.0.) then
            GOTO 999
          end if

         hh53=hh**(5./3.)

! Calculate q-flux...
          qqsfc = (sfy/abs(sfy))*alfay*hh53*dt / gsize
          q_sfcflx_y(I,J) = q_sfcflx_y(I,J) + qqsfc


!DJG put adjustment in for (h) due to qqsfc
!yw	  h(i,j) = h(i,j)-qqsfc*1000.
!yw          h(i,j+1) = h(i,j+1) + qqsfc*1000.
!yw changed as following:
           tmp_adjust=qqsfc*1000
          if(tmp_adjust .le. 0 ) GOTO 999

           if((h(i,j) - tmp_adjust) <0 )  then
#ifdef HYDRO_D
               print *, "WARNING: surface head is negative:  ",i,j
#endif
               tmp_adjust = h(i,j)
           end if
           if((h(i,j+1) + tmp_adjust) <0) then
#ifdef HYDRO_D
               print *, "WARNING: surface head is negative: ",i,j+1
#endif
               tmp_adjust = -1*h(i,j+1)
           end if
	  Dh(i,j) = Dh(i,j)-tmp_adjust
          Dh(i,j+1) = Dh(i,j+1) + tmp_adjust
!yw end change

!          qsfc(i,j) = qsfc(i,j)-qqsfc
!          qsfc(i,j+1) = qsfc(i,j+1) + qqsfc
!!DG Boundary adjustments here
!!DG Constant Flux Condition
#ifdef MPP_LAND
          if ((j.eq.1).AND.(sfy.lt.0)   &
             .and. (down_id .lt. 0) ) then
#else
          if ((j.eq.1).AND.(sfy.lt.0)) then
#endif
            Dh(i,j) = Dh(i,j) + qqsfc*1000.
            QBDRY(I,J)=QBDRY(I,J) + qqsfc*1000.
            QBDRYT=QBDRYT + qqsfc*1000.
#ifdef MPP_LAND
          else if ((j.eq.(YY-1)).AND.(sfy.gt.0) &
             .and. (up_id .lt. 0) ) then
#else
          else if ((j.eq.(YY-1)).AND.(sfy.gt.0)) then
#endif
             tmp_adjust = qqsfc*1000.
             if(h(i,j+1).lt.tmp_adjust) tmp_adjust = h(i,j+1)
             Dh(i,j+1) = Dh(i,j+1) - tmp_adjust
!DJG Re-assign h(j+1) = 0.0 when <0.0 (from rounding/truncation error)
            QBDRY(I,J+1)=QBDRY(I,J+1) - tmp_adjust
            QBDRYT=QBDRYT - tmp_adjust
          end if

999     continue

!!!! End loop to route sfc water in y-direction
          end do
        end do

        H = H +DH
#ifdef MPP_LAND
       call MPP_LAND_UB_COM(H,XX,YY,99)
       call MPP_LAND_UB_COM(QBDRY,XX,YY,99)
#endif
        return

!DJG ----------------------------------------------------------------------
        end subroutine ROUTE_OVERLAND2




!DJG-----------------------------------------------------------------------
!DJG SUBROUTINE TER_ADJ_SOL    - Terrain adjustment of incoming solar radiation
!DJG-----------------------------------------------------------------------
	SUBROUTINE TER_ADJ_SOL(IX,JX,SO8LD_D,TSLP,SHORT,XLAT,XLONG,olddate,DT)

#ifdef MPP_LAND
        use module_mpp_land, only:  my_id, io_id, &
             mpp_land_bcast_int1 
#endif
          implicit none
          integer,INTENT(IN)     :: IX,JX
          INTEGER,INTENT(in), DIMENSION(IX,JX,3)   :: SO8LD_D
          real,INTENT(IN), DIMENSION(IX,JX)  :: XLAT,XLONG
 	  real,INTENT(IN) :: DT
          real,INTENT(INOUT), DIMENSION(IX,JX)  :: SHORT
          character(len=19) :: olddate

! Local Variables...
          real, dimension(IX,JX) ::TSLP,TAZI
          real, dimension(IX,JX) ::SOLDN
	  real :: SOLDEC,DGRD,ITIME2,HRANGLE
	  real :: BINSH,SOLZANG,SOLAZI,INCADJ
	  real :: TAZIR,TSLPR,LATR,LONR,SOLDNADJ
          integer :: JULDAY0,HHTIME0,MMTIME0,YYYY0,MM0,DD0
          integer :: JULDAY,HHTIME,MMTIME,YYYY,MM,DD
	  integer :: I,J
          

!----------------------------------------------------------------------
!  SPECIFY PARAMETERS and VARIABLES
!----------------------------------------------------------------------

       JULDAY = 0	
       SOLDN = SHORT
       DGRD = 3.14159/180.
       
! Set up time variables...
#ifdef MPP_LAND   
       if(my_id .eq. IO_id) then
#endif
          read(olddate(1:4),"(I4)") YYYY0 ! real-time year (GMT)
          read(olddate(6:7),"(I2.2)") MM0 ! real-time month (GMT)
          read(olddate(9:10),"(I2.2)") DD0 ! real-time day (GMT)
          read(olddate(12:13),"(I2.2)") HHTIME0 ! real-time hour (GMT)
          read(olddate(15:16),"(I2.2)") MMTIME0 ! real-time minutes (GMT)
#ifdef MPP_LAND   
       endif
       call mpp_land_bcast_int1(YYYY0) 
       call mpp_land_bcast_int1(MM0) 
       call mpp_land_bcast_int1(DD0) 
       call mpp_land_bcast_int1(HHTIME0) 
       call mpp_land_bcast_int1(MMTIME0) 
#endif


! Set up terrain variables...(returns TSLP&TAZI in radians) 
	call SLOPE_ASPECT(IX,JX,SO8LD_D,TAZI)

!----------------------------------------------------------------------
!  BEGIN LOOP THROUGH GRID
!----------------------------------------------------------------------
        DO J=1,JX
          DO I=1,IX
             YYYY = YYYY0
             MM  = MM0
             DD  = DD0
             HHTIME = HHTIME0
             MMTIME = MMTIME0
      	     call GMT2LOCAL(1,1,XLONG(i,j),YYYY,MM,DD,HHTIME,MMTIME,DT)
             call JULDAY_CALC(YYYY,MM,DD,JULDAY)

! Convert to radians...
           LATR = XLAT(I,J)   !send solsub local lat in deg
           LONR = XLONG(I,J)   !send solsub local lon in deg
           TSLPR = TSLP(I,J)/DGRD !send solsub local slp in deg
           TAZIR = TAZI(I,J)/DGRD !send solsub local azim in deg

!Call SOLSUB to return terrain adjusted incoming solar radiation...
! SOLSUB taken from Whiteman and Allwine, 1986, Environ. Software.

          call SOLSUB(LONR,LATR,TAZIR,TSLPR,SOLDN(I,J),YYYY,MM,         &
               DD,HHTIME,MMTIME,SOLDNADJ,SOLZANG,SOLAZI,INCADJ)

         SOLDN(I,J)=SOLDNADJ

          ENDDO
        ENDDO

	SHORT = SOLDN

        return
	end SUBROUTINE TER_ADJ_SOL  
!DJG-----------------------------------------------------------------------
!DJG END SUBROUTINE TER_ADJ_SOL
!DJG-----------------------------------------------------------------------


!DJG-----------------------------------------------------------------------
!DJG SUBROUTINE GMT2LOCAL
!DJG-----------------------------------------------------------------------
	subroutine GMT2LOCAL(IX,JX,XLONG,YY,MM,DD,HH,MIN,DT)

       implicit none

!!! Declare Passed Args.

        INTEGER, INTENT(INOUT) :: yy,mm,dd,hh,min
        INTEGER, INTENT(IN) :: IX,JX
        REAL,INTENT(IN), DIMENSION(IX,JX)  :: XLONG
        REAL,INTENT(IN) :: DT

!!! Declare local variables

        integer :: i,j,minflag,hhflag,ddflag,mmflag,yyflag
        integer :: adj_min,lst_adj_min,lst_adj_hh,adj_hh
        real, dimension(IX,JX) :: TDIFF
        real :: tmp
        integer :: yyinit,mminit,ddinit,hhinit,mininit

!!! Initialize flags
        hhflag=0
        ddflag=0
        mmflag=0
        yyflag=0

!!! Set up constants...
        yyinit = yy
   	mminit = mm
        ddinit = dd
        hhinit = hh
        mininit = min


! Loop through data...
     do j=1,JX
      do i=1,IX

! Reset yy,mm,dd...
        yy = yyinit
	mm = mminit
        dd = ddinit
        hh = hhinit
	min = mininit

!!! Set up adjustments...
!   - assumes +E , -W  longitude and 0.06667 hr/deg (=24/360)
       TDIFF(I,J) = XLONG(I,J)*0.06667   ! time offset in hr
       tmp = TDIFF(I,J)
       lst_adj_hh = INT(tmp)
       lst_adj_min = NINT(MOD(int(tmp),1)*60.) + int(DT/2./60.)  ! w/ 1/2 timestep adjustment...

!!! Process Minutes...
        adj_min = min+lst_adj_min
        if (adj_min.lt.0) then
          min=60+adj_min
          lst_adj_hh = lst_adj_hh - 1
        else if (adj_min.ge.0.AND.adj_min.lt.60) then
          min=adj_min
        else if (adj_min.ge.60) then
          min=adj_min-60
          lst_adj_hh = lst_adj_hh + 1
        end if

!!! Process Hours
        adj_hh = hh+lst_adj_hh
	if (adj_hh.lt.0) then
          hh = 24+adj_hh
          ddflag=1
        else if (adj_hh.ge.0.AND.adj_hh.lt.24) then
          hh=adj_hh
        else if (adj_hh.ge.24) then
          hh=adj_hh-24
          ddflag = 2
        end if



!!! Process Days, Months, Years
! Subtract a day
        if (ddflag.eq.1) then
          if (dd.gt.1) then
            dd=dd-1
          else
            if (mm.eq.1) then
              mm=12
              yy=yy-1
            else
              mm=mm-1
            end if
            if ((mm.eq.1).or.(mm.eq.3).or.(mm.eq.5).or. &
                (mm.eq.7).or.(mm.eq.8).or.(mm.eq.10).or. &
                 (mm.eq.12)) then
              dd=31
            else

!!! Adjustment for leap years!!!
                if(mm.eq.2) then
                  if(MOD(yy,4).eq.0) then
                    dd=29
                  else
                    dd=28
                  end if
                end if
                if(mm.ne.2) dd=30
            end if
          end if
        end if

! Add a day
        if (ddflag.eq.2) then
          if ((mm.eq.1).or.(mm.eq.3).or.(mm.eq.5).or. &
                (mm.eq.7).or.(mm.eq.8).or.(mm.eq.10).or. &
                 (mm.eq.12)) then
            if (dd.eq.31) then
              dd=1
              if (mm.eq.12) then
                mm=1
                yy=yy+1
              else
                mm=mm+1
              end if
            else
              dd=dd+1
            end if

!!! Adjustment for leap years!!!
          else if (mm.eq.2) then
            if(MOD(yy,4).eq.0) then
              if (dd.eq.29) then
                dd=1
                mm=3
              else
                dd=dd+1
              end if
            else
              if (dd.eq.28) then
                dd=1
                mm=3
              else
                dd=dd+1
              end if
            end if
          else
            if (dd.eq.30) then
              dd=1
              mm=mm+1
            else
              dd=dd+1
            end if
          end if

        end if

       end do   !i-loop
      end do   !j-loop

        return
        end subroutine

!DJG-----------------------------------------------------------------------
!DJG END SUBROUTINE GMT2LOCAL
!DJG-----------------------------------------------------------------------



!DJG-----------------------------------------------------------------------
!DJG SUBROUTINE JULDAY_CALC
!DJG-----------------------------------------------------------------------
      subroutine JULDAY_CALC(YYYY,MM,DD,JULDAY)

	implicit none
	integer,intent(in) :: YYYY,MM,DD
        integer,intent(out) :: JULDAY

        integer :: resid
        integer julm(13)
        DATA JULM/0, 31, 59, 90, 120, 151, 181, 212, 243, 273, &
           304, 334, 365 /

        integer LPjulm(13)
        DATA LPJULM/0, 31, 60, 91, 121, 152, 182, 213, 244, 274, &
           305, 335, 366 /

      resid = MOD(YYYY,4) !Set up leap year check...

      if (resid.ne.0) then    !If not a leap year....
        JULDAY = JULM(MM) + DD
      else                    !If a leap year...
        JULDAY = LPJULM(MM) + DD
      end if

      RETURN
      END subroutine JULDAY_CALC
!DJG-----------------------------------------------------------------------
!DJG END SUBROUTINE JULDAY
!DJG-----------------------------------------------------------------------

!DJG-----------------------------------------------------------------------
!DJG SUBROUTINE SLOPE_ASPECT
!DJG-----------------------------------------------------------------------
	subroutine SLOPE_ASPECT(IX,JX,SO8LD_D,TAZI)

	implicit none
        integer, INTENT(IN)		   :: IX,JX
!	real,INTENT(in),DIMENSION(IX,JX)   :: TSLP  !terrain slope (m/m)
	real,INTENT(OUT),DIMENSION(IX,JX)   :: TAZI  !terrain aspect (deg)

        INTEGER, DIMENSION(IX,JX,3)   :: SO8LD_D
	real :: DGRD
	integer :: i,j

!	TSLP = 0.  !Initialize as flat
	TAZI = 0.  !Initialize as north facing

! Find steepest descent slope and direction...
        do j=1,JX
          do i=1,IX
!	TSLP(I,J) = TANH(Vmax(i,j)) ! calculate slope in radians...

! Convert steepest slope and aspect to radians...
        IF (SO8LD_D(i,j,3).eq.1) then
          TAZI(I,J) = 0.0
        ELSEIF (SO8LD_D(i,j,3).eq.2) then
          TAZI(I,J) = 45.0
        ELSEIF (SO8LD_D(i,j,3).eq.3) then
          TAZI(I,J) = 90.0
        ELSEIF (SO8LD_D(i,j,3).eq.4) then
          TAZI(I,J) = 135.0
        ELSEIF (SO8LD_D(i,j,3).eq.5) then
          TAZI(I,J) = 180.0
        ELSEIF (SO8LD_D(i,j,3).eq.6) then
          TAZI(I,J) = 225.0
        ELSEIF (SO8LD_D(i,j,3).eq.7) then
          TAZI(I,J) = 270.0
        ELSEIF (SO8LD_D(i,j,3).eq.8) then
          TAZI(I,J) = 315.0
	END IF

        DGRD = 3.141593/180.
	TAZI(I,J) = TAZI(I,J)*DGRD ! convert azimuth to radians...

        END DO
      END DO

      RETURN
      END  subroutine SLOPE_ASPECT
!DJG-----------------------------------------------------------------------
!DJG END SUBROUTINE SLOPE_ASPECT
!DJG-----------------------------------------------------------------------

!DJG----------------------------------------------------------------
!DJG    SUBROUTINE SOLSUB
!DJG----------------------------------------------------------------
        SUBROUTINE SOLSUB(LONG,LAT,AZ,IN,SC,YY,MO,IDA,IHR,MM,OUT1, &
                          OUT2,OUT3,INCADJ)


! Notes....

        implicit none
          logical               :: daily, first
          integer               :: yy,mo,ida,ihr,mm,d
          integer,dimension(12) :: nday
          real                  :: lat,long,longcor,longsun,in,inslo
          real :: az,sc,out1,out2,out3,cosbeta,dzero,eccent,pi,calint
          real :: rtod,decmax,omega,onehr,omd,omdzero,rdvecsq,sdec
          real :: declin,cdec,arg,declon,sr,stdmrdn,b,em,timnoon,azslo
          real :: slat,clat,caz,saz,sinc,cinc,hinc,h,cosz,extra,extslo
          real :: t1,z,cosa,a,cosbeta_flat,INCADJ
          integer :: HHTIME,MMTIME,i,ik
          real, dimension(4) :: ACOF,BCOF

! Constants
       daily=.FALSE.
       ACOF(1) = 0.00839
       ACOF(2) = -0.05391
       ACOF(3) = -0.00154
       ACOF(4) = -0.0022
       BCOF(1) = -0.12193
       BCOF(2) = -0.15699
       BCOF(3) = -0.00657
       BCOF(4) = -0.00370
       DZERO = 80.
       ECCENT = 0.0167
       PI = 3.14159
       CALINT = 1.
       RTOD = PI / 180.
       DECMAX=(23.+26./60.)*RTOD
       OMEGA=2*PI/365.
       ONEHR=15.*RTOD

! Calculate Julian Day...
       D = 0
       call JULDAY_CALC(YY,MO,IDA,D)

! Ratio of radius vectors squared...
       OMD=OMEGA*D
       OMDZERO=OMEGA*DZERO
!       RDVECSQ=1./(1.-ECCENT*COS(OMD))**2
       RDVECSQ = 1.    ! no adjustment for orbital changes when coupled to HRLDAS...

! Declination of sun...
       LONGSUN=OMEGA*(D-Dzero)+2.*ECCENT*(SIN(OMD)-SIN(OMDZERO))
       DECLIN=ASIN(SIN(DECMAX)*SIN(LONGSUN))
       SDEC=SIN(DECLIN)
       CDEC=COS(DECLIN)

! Check for Polar Day/night...
       ARG=((PI/2.)-ABS(DECLIN))/RTOD
       IF(ABS(LAT).GT.ARG) THEN
         IF((LAT.GT.0..AND.DECLIN.LT.0) .OR.       &
             (LAT.LT.0..AND.DECLON.GT.0.)) THEN
               OUT1 = 0.
               OUT2 = 0.
               OUT3 = 0.
               RETURN
         ENDIF
         SR=-1.*PI
       ELSE

! Calculate sunrise hour angle...
         SR=-1.*ABS(ACOS(-1.*TAN(LAT*RTOD)*TAN(DECLIN)))
       END IF

! Find standard meridian for site
       STDMRDN=NINT(LONG/15.)*15.
       LONGCOR=(LONG-STDMRDN)/15.

! Compute time correction from equation of time...
       B=2.*PI*(D-.4)/365
       EM=0.
       DO I=1,4
         EM=EM+(BCOF(I)*SIN(I*B)+ACOF(I)*COS(I*B))
       END DO

! Compute time of solar noon...
       TIMNOON=12.-EM-LONGCOR

! Set up a few more terms...
       AZSLO=AZ*RTOD
       INSLO=IN*RTOD
       SLAT=SIN(LAT*RTOD)
       CLAT=COS(LAT*RTOD)
       CAZ=COS(AZSLO)
       SAZ=SIN(AZSLO)
       SINC=SIN(INSLO)
       CINC=COS(INSLO)

! Begin solar radiation calculations...daily first, else instantaneous...
       IF (DAILY) THEN   ! compute daily integrated values...(Not used in HRLDAS!)
         IHR=0
         MM=0
         HINC=CALINT*ONEHR/60.
         IK=(2.*ABS(SR)/HINC)+2.
         FIRST=.TRUE.
         OUT1=0.
         DO I=1,IK
           H=SR+HINC*FLOAT(I-1)
           COSZ=SLAT*SDEC+CLAT*CDEC*COS(H)
           COSBETA=CDEC*((SLAT*COS(H))*(-1.*CAZ*SINC)- &
                SIN(H)*(SAZ*SINC)+(CLAT*COS(H))*CINC)+ &
                SDEC*(CLAT*(CAZ*SINC)+SLAT*CINC)
           EXTRA=SC*RDVECSQ*COSZ
           IF(EXTRA.LE.0.) EXTRA=0.
           EXTSLO=SC*RDVECSQ*COSBETA
           IF(EXTRA.LE.0. .OR. EXTSLO.LT.0.) EXTSLO=0.
           IF(FIRST .AND. EXTSLO.GT.0.) THEN
             OUT2=(H-HINC)/ONEHR+TIMNOON
             FIRST = .FALSE.
           END IF
           IF(.NOT.FIRST .AND. EXTSLO.LE.0.) OUT3=H/ONEHR+TIMNOON
           OUT1=EXTSLO+OUT1
         END DO
         OUT1=OUT1*CALINT*60./1000000.

       ELSE   ! Compute instantaneous values...(Is used in HRLDAS!)

         T1=FLOAT(IHR)+FLOAT(MM)/60.
         H=ONEHR*(T1-TIMNOON)
         COSZ=SLAT*SDEC+CLAT*CDEC*COS(H)

! Assuming HRLDAS forcing already accounts for season, time of day etc,
! subtract out the component of adjustment that would occur for
! a flat surface, this should leave only the sloped component remaining

         COSBETA=CDEC*((SLAT*COS(H))*(-1.*CAZ*SINC)-  &
              SIN(H)*(SAZ*SINC)+(CLAT*COS(H))*CINC)+ &
              SDEC*(CLAT*(CAZ*SINC)+SLAT*CINC)

         COSBETA_FLAT=CDEC*CLAT*COS(H)+SDEC*SLAT

         INCADJ = COSBETA+(1-COSBETA_FLAT)

         EXTRA=SC*RDVECSQ*COSZ
         IF(EXTRA.LE.0.) EXTRA=0.
         EXTSLO=SC*RDVECSQ*INCADJ
!         IF(EXTRA.LE.0. .OR. EXTSLO.LT.0.) EXTSLO=0.  !remove check for HRLDAS.
         OUT1=EXTSLO
         Z=ACOS(COSZ)
         COSA=(SLAT*COSZ-SDEC)/(CLAT*SIN(Z))
         IF(COSA.LT.-1.) COSA=-1.
         IF(COSA.GT.1.) COSA=1.
         A=ABS(ACOS(COSA))
         IF(H.LT.0.) A=-A
         OUT2=Z/RTOD
         OUT3=A/RTOD+180

       END IF    ! End if for daily vs instantaneous values...

!DJG-----------------------------------------------------------------------
       RETURN
       END SUBROUTINE SOLSUB
!DJG-----------------------------------------------------------------------
       
       subroutine seq_land_SO8(SO8LD_D,Vmax,TERR,dx,ix,jx)
         implicit none
         integer :: ix,jx,i,j
         REAL, DIMENSION(IX,JX,8)      :: SO8LD
         INTEGER, DIMENSION(IX,JX,3)   :: SO8LD_D
         real,DIMENSION(IX,JX)      :: TERR
         real                       :: dx(ix,jx,9),Vmax(ix,jx)
         SO8LD_D = -1
         do j = 2, jx -1
            do i = 2, ix -1
               SO8LD(i,j,1) = (TERR(i,j)-TERR(i,j+1))/dx(i,j,1)
               SO8LD_D(i,j,1) = i
               SO8LD_D(i,j,2) = j + 1
               SO8LD_D(i,j,3) = 1
               Vmax(i,j) = SO8LD(i,j,1)

               SO8LD(i,j,2) = (TERR(i,j)-TERR(i+1,j+1))/DX(i,j,2)
               if(SO8LD(i,j,2) .gt. Vmax(i,j) ) then
                 SO8LD_D(i,j,1) = i + 1
                 SO8LD_D(i,j,2) = j + 1
                 SO8LD_D(i,j,3) = 2
                 Vmax(i,j) = SO8LD(i,j,2)
               end if
               SO8LD(i,j,3) = (TERR(i,j)-TERR(i+1,j))/DX(i,j,3)
               if(SO8LD(i,j,3) .gt. Vmax(i,j) ) then
                 SO8LD_D(i,j,1) = i + 1
                 SO8LD_D(i,j,2) = j
                 SO8LD_D(i,j,3) = 3
                 Vmax(i,j) = SO8LD(i,j,3)
               end if
               SO8LD(i,j,4) = (TERR(i,j)-TERR(i+1,j-1))/DX(i,j,4)
               if(SO8LD(i,j,4) .gt. Vmax(i,j) ) then
                 SO8LD_D(i,j,1) = i + 1
                 SO8LD_D(i,j,2) = j - 1
                 SO8LD_D(i,j,3) = 4
                 Vmax(i,j) = SO8LD(i,j,4)
               end if
               SO8LD(i,j,5) = (TERR(i,j)-TERR(i,j-1))/DX(i,j,5)
               if(SO8LD(i,j,5) .gt. Vmax(i,j) ) then
                 SO8LD_D(i,j,1) = i
                 SO8LD_D(i,j,2) = j - 1
                 SO8LD_D(i,j,3) = 5
                 Vmax(i,j) = SO8LD(i,j,5)
               end if
               SO8LD(i,j,6) = (TERR(i,j)-TERR(i-1,j-1))/DX(i,j,6)
               if(SO8LD(i,j,6) .gt. Vmax(i,j) ) then
                 SO8LD_D(i,j,1) = i - 1
                 SO8LD_D(i,j,2) = j - 1
                 SO8LD_D(i,j,3) = 6
                 Vmax(i,j) = SO8LD(i,j,6)
               end if
               SO8LD(i,j,7) = (TERR(i,j)-TERR(i-1,j))/DX(i,j,7)
               if(SO8LD(i,j,7) .gt. Vmax(i,j) ) then
                 SO8LD_D(i,j,1) = i - 1
                 SO8LD_D(i,j,2) = j
                 SO8LD_D(i,j,3) = 7
                 Vmax(i,j) = SO8LD(i,j,7)
               end if
               SO8LD(i,j,8) = (TERR(i,j)-TERR(i-1,j+1))/DX(i,j,8)
               if(SO8LD(i,j,8) .gt. Vmax(i,j) ) then
                 SO8LD_D(i,j,1) = i - 1
                 SO8LD_D(i,j,2) = j + 1
                 SO8LD_D(i,j,3) = 8
                 Vmax(i,j) = SO8LD(i,j,8)
               end if
             enddo
          enddo
          Vmax = TANH(Vmax)    
          return
          end  subroutine seq_land_SO8

#ifdef MPP_LAND
       subroutine MPP_seq_land_SO8(SO8LD_D,Vmax,TERRAIN,dx,ix,jx,&
         global_nx,global_ny)

         use module_mpp_land, only:  my_id, io_id, &
              write_io_real,decompose_data_int,decompose_data_real

         implicit none
         integer,intent(in) :: ix,jx,global_nx,global_ny
         INTEGER, intent(inout),DIMENSION(IX,JX,3)   :: SO8LD_D
!         real,intent(in), DIMENSION(IX,JX)   :: TERRAIN
         real,DIMENSION(IX,JX)   :: TERRAIN
         real,intent(out),dimension(ix,jx) ::  Vmax
         real,intent(in)                     :: dx(ix,jx,9)
         real                     :: g_dx(ix,jx,9)

         real,DIMENSION(global_nx,global_ny)      :: g_TERRAIN
         real,DIMENSION(global_nx,global_ny)      :: g_Vmax
         integer,DIMENSION(global_nx,global_ny,3)      :: g_SO8LD_D
         integer :: k

         g_SO8LD_D = 0
         g_Vmax    = 0
       
         do k = 1, 9 
            call write_IO_real(dx(:,:,k),g_dx(:,:,k)) 
         end do

         call write_IO_real(TERRAIN,g_TERRAIN)
         if(my_id .eq. IO_id) then
            call seq_land_SO8(g_SO8LD_D,g_Vmax,g_TERRAIN,g_dx,global_nx,global_ny)
         endif
          call decompose_data_int(g_SO8LD_D(:,:,3),SO8LD_D(:,:,3))
          call decompose_data_real(g_Vmax,Vmax)
         return
         end subroutine MPP_seq_land_SO8

#endif



      subroutine disaggregateDomain_drv(did)
           use module_RT_data, only: rt_domain
           use module_namelist, only: nlst_rt
           integer :: did
           call disaggregateDomain( RT_DOMAIN(did)%IX,RT_DOMAIN(did)%JX,nlst_rt(did)%NSOIL,&
             RT_DOMAIN(did)%IXRT,RT_DOMAIN(did)%JXRT,nlst_rt(did)%AGGFACTRT,RT_DOMAIN(did)%SICE, &
             RT_DOMAIN(did)%SMC,RT_DOMAIN(did)%SH2OX,RT_DOMAIN(did)%INFXSRT, &
             rt_domain(did)%dist_lsm(:,:,9),RT_DOMAIN(did)%SMCMAX1,RT_DOMAIN(did)%SMCREF1, &
             RT_DOMAIN(did)%SMCWLT1,RT_DOMAIN(did)%VEGTYP,RT_DOMAIN(did)%LKSAT,RT_DOMAIN(did)%dist, &
             RT_DOMAIN(did)%INFXSWGT, RT_DOMAIN(did)%OVROUGHRTFAC,RT_DOMAIN(did)%LKSATFAC, &
             RT_DOMAIN(did)%CH_NETRT,RT_DOMAIN(did)%SH2OWGT,RT_DOMAIN(did)%SMCREFRT,       &
             RT_DOMAIN(did)%INFXSUBRT,RT_DOMAIN(did)%SMCMAXRT, RT_DOMAIN(did)%SMCWLTRT,    &
             RT_DOMAIN(did)%SMCRT, &
             RT_DOMAIN(did)%OVROUGHRT, RT_DOMAIN(did)%LAKE_MSKRT, &
             RT_DOMAIN(did)%LKSATRT, RT_DOMAIN(did)%OV_ROUGH,RT_DOMAIN(did)%SLDPTH, &
	     RT_DOMAIN(did)%soiltypRT, RT_DOMAIN(did)%soiltyp)

      end subroutine disaggregateDomain_drv

      subroutine disaggregateDomain(IX,JX,NSOIL,IXRT,JXRT,AGGFACTRT, &
                     SICE, SMC,SH2OX, INFXSRT, area_lsm, SMCMAX1,SMCREF1, &
               SMCWLT1, VEGTYP, LKSAT, dist,INFXSWGT,OVROUGHRTFAC, &
               LKSATFAC, CH_NETRT,SH2OWGT,SMCREFRT, INFXSUBRT,SMCMAXRT, &
               SMCWLTRT,SMCRT, OVROUGHRT, LAKE_MSKRT, LKSATRT,OV_ROUGH,  &
               SLDPTH, soiltypRT, soiltyp                               &
            )
#ifdef MPP_LAND
        use module_mpp_land, only: left_id,down_id,right_id, &
              up_id,mpp_land_com_real, my_id, io_id, numprocs, &
             mpp_land_sync,mpp_land_com_integer,mpp_land_max_int1, &
             sum_real1
#endif
     implicit none
        integer,INTENT(IN) :: IX,JX,NSOIL,IXRT,JXRT,AGGFACTRT
        real, INTENT(OUT), DIMENSION(IX,JX,NSOIL) :: SICE
        real, INTENT(IN),  DIMENSION(IX,JX,NSOIL) :: SMC,SH2OX
        real, INTENT(IN),  DIMENSION(IX,JX) :: INFXSRT, area_lsm, SMCMAX1,SMCREF1, &
                                               SMCWLT1,  LKSAT
        integer, INTENT(IN), DIMENSION(IX,JX)      :: VEGTYP, soiltyp

        real,INTENT(IN),DIMENSION(IXRT,JXRT,9)::dist
        real,INTENT(IN),DIMENSION(IXRT,JXRT)::INFXSWGT,OVROUGHRTFAC, &
               LKSATFAC
        integer,INTENT(INOUT), DIMENSION(IXRT,JXRT)     ::CH_NETRT, soiltypRT
        real,INTENT(IN),DIMENSION(IXRT,JXRT,NSOIL)::SH2OWGT
        real,INTENT(OUT),DIMENSION(IXRT,JXRT,NSOIL)::SMCREFRT, SMCMAXRT, &
               SMCWLTRT,SMCRT
        real,INTENT(OUT),DIMENSION(IXRT,JXRT)::INFXSUBRT
        real,INTENT(INOUT),DIMENSION(IXRT,JXRT)::OVROUGHRT, LKSATRT
        integer,INTENT(INOUT), DIMENSION(IXRT,JXRT)  ::LAKE_MSKRT
                  

        real,INTENT(IN), DIMENSION(NSOIL)      :: SLDPTH
        REAL    OV_ROUGH(*)



        integer :: i, j, AGGFACYRT, AGGFACXRT, IXXRT, JYYRT,KRT, KF
        REAL  ::  LSMVOL,SMCEXCS, WATHOLDCAP

#ifdef HYDRO_D
! ADCHANGE: Water balance variables
       integer :: kk
       real    :: smctot1,smcrttot2
       real    :: sicetot1
       real    :: suminfxs1,suminfxsrt2
#endif

!-------------------------------------



	SICE=SMC-SH2OX
        SMCREFRT = 0

!DJG First, Disaggregate a few key fields for routing...
!DJG Debug...
#ifdef HYDRO_D
	print *, "Beginning Disaggregation..."
#endif
	
!DJG Mass balance check for disagg...

#ifdef HYDRO_D
! ADCHANGE: START Initial water balance variables
! ALL VARS in MM
        suminfxs1 = 0.
        smctot1 = 0.
        sicetot1 = 0.
        do i=1,IX
         do j=1,JX
            suminfxs1 = suminfxs1 + INFXSRT(I,J) / float(IX*JX)
            do kk=1,NSOIL
                smctot1 = smctot1 + SMC(I,J,KK)*SLDPTH(KK)*1000. / float(IX*JX)
                sicetot1 = sicetot1 + SICE(I,J,KK)*SLDPTH(KK)*1000. / float(IX*JX)
            end do
         end do
        end do

#ifdef MPP_LAND
! not tested
        CALL sum_real1(suminfxs1)
        CALL sum_real1(smctot1)
        CALL sum_real1(sicetot1)
        suminfxs1 = suminfxs1/float(numprocs)
        smctot1 = smctot1/float(numprocs)
        sicetot1 = sicetot1/float(numprocs)
#endif
! END Initial water balance variables
#endif

!DJG Weighting alg. alteration...(prescribe wghts if time = 1)


        do J=1,JX
          do I=1,IX

!DJG Weighting alg. alteration...
              LSMVOL=INFXSRT(I,J)*area_lsm(I,J)


             do AGGFACYRT=AGGFACTRT-1,0,-1
              do AGGFACXRT=AGGFACTRT-1,0,-1

               IXXRT=I*AGGFACTRT-AGGFACXRT
               JYYRT=J*AGGFACTRT-AGGFACYRT
#ifdef MPP_LAND
       if(left_id.ge.0) IXXRT=IXXRT+1
       if(down_id.ge.0) JYYRT=JYYRT+1
#else
!yw ????
!       IXXRT=IXXRT+1
!       JYYRT=JYYRT+1
#endif
!        if(AGGFACTRT .eq. 1) then
!            IXXRT=I
!            JYYRT=J
!        endif


!DJG Implement subgrid weighting routine...
               INFXSUBRT(IXXRT,JYYRT)=LSMVOL*     &
                   INFXSWGT(IXXRT,JYYRT)/dist(IXXRT,JYYRT,9)
  

            do KRT=1,NSOIL  !Do for soil profile loop
               IF(SICE(I,J,KRT).gt.0) then  !...adjust for soil ice
!DJG Adjust SMCMAX for SICE when subsfc routing...make 3d variable
                 SMCMAXRT(IXXRT,JYYRT,KRT)=SMCMAX1(I,J)-SICE(I,J,KRT)
                 SMCREFRT(IXXRT,JYYRT,KRT)=SMCREF1(I,J)-SICE(I,J,KRT)
                 WATHOLDCAP = SMCMAX1(I,J) - SMCWLT1(I,J)
                 IF (SICE(I,J,KRT).le.WATHOLDCAP)    then
                        SMCWLTRT(IXXRT,JYYRT,KRT) = SMCWLT1(I,J)      
                 else
                    if(SICE(I,J,KRT).lt.SMCMAX1(I,J)) &
                          SMCWLTRT(IXXRT,JYYRT,KRT) = SMCWLT1(I,J) - &
                          (SICE(I,J,KRT)-WATHOLDCAP)
                    if(SICE(I,J,KRT).ge.SMCMAX1(I,J)) SMCWLTRT(IXXRT,JYYRT,KRT) = 0.
                 end if
               ELSE
                 SMCMAXRT(IXXRT,JYYRT,KRT)=SMCMAX1(I,J)
                 SMCREFRT(IXXRT,JYYRT,KRT)=SMCREF1(I,J)
                 WATHOLDCAP = SMCMAX1(I,J) - SMCWLT1(I,J)
                 SMCWLTRT(IXXRT,JYYRT,KRT) = SMCWLT1(I,J) 
               END IF   !endif adjust for soil ice...


!Now Adjust soil moisture
!DJG Use SH2O instead of SMC for 'liquid' water...
                 IF(SMCMAXRT(IXXRT,JYYRT,KRT).GT.0) THEN !Check for smcmax data (=0 over water)
                   SMCRT(IXXRT,JYYRT,KRT)=SH2OX(I,J,KRT)*SH2OWGT(IXXRT,JYYRT,KRT)
!old                   SMCRT(IXXRT,JYYRT,KRT)=SMC(I,J,KRT)
                 ELSE
                   SMCRT(IXXRT,JYYRT,KRT) = 0.001  !will be skipped w/ landmask
                   SMCMAXRT(IXXRT,JYYRT,KRT) = 0.001
                 END IF
!DJG Check/Adjust so that subgrid cells do not exceed saturation...
                 IF (SMCRT(IXXRT,JYYRT,KRT).GT.SMCMAXRT(IXXRT,JYYRT,KRT)) THEN
                   SMCEXCS = (SMCRT(IXXRT,JYYRT,KRT) - SMCMAXRT(IXXRT,JYYRT,KRT)) &
                             * SLDPTH(KRT)*1000.  !Excess soil water in units of (mm)
                   SMCRT(IXXRT,JYYRT,KRT) = SMCMAXRT(IXXRT,JYYRT,KRT)
                   DO KF = KRT-1,1, -1  !loop back upward to redistribute excess water from disagg.
                     SMCRT(IXXRT,JYYRT,KF) = SMCRT(IXXRT,JYYRT,KF) + SMCEXCS/(SLDPTH(KF)*1000.) 
                     IF (SMCRT(IXXRT,JYYRT,KF).GT.SMCMAXRT(IXXRT,JYYRT,KF)) THEN  !Recheck new lyr sat.
                       SMCEXCS = (SMCRT(IXXRT,JYYRT,KF) - SMCMAXRT(IXXRT,JYYRT,KF)) &
                           * SLDPTH(KF)*1000.  !Excess soil water in units of (mm)
                       SMCRT(IXXRT,JYYRT,KF) = SMCMAXRT(IXXRT,JYYRT,KF)
                     ELSE  ! Excess soil water expired
                       SMCEXCS = 0.
                       EXIT
                     END IF
                   END DO
                   IF (SMCEXCS.GT.0) THEN  !If not expired by sfc then add to Infil. Excess
                     INFXSUBRT(IXXRT,JYYRT) = INFXSUBRT(IXXRT,JYYRT) + SMCEXCS
                     SMCEXCS = 0.
                   END IF
                 END IF  !End if for soil moisture saturation excess


             end do !End do for soil profile loop



             do KRT=1,NSOIL  !debug loop
               IF (SMCRT(IXXRT,JYYRT,KRT).GT.SMCMAXRT(IXXRT,JYYRT,KRT)) THEN
                      print *, "FATAL ERROR: SMCMAX exceeded upon disaggregation3...", ixxrt,jyyrt,krt,&
                       SMCRT(IXXRT,JYYRT,KRT),SMCMAXRT(IXXRT,JYYRT,KRT)
                      call hydro_stop("In disaggregateDomain() - SMCMAX exceeded upon disaggregation3")
               ELSE IF (SMCRT(IXXRT,JYYRT,KRT).LE.0.) THEN
                       print *, "FATAL ERROR: SMCRT fully depleted upon disaggregation...", ixxrt,jyyrt,krt,&
                       "SMCRT=",SMCRT(IXXRT,JYYRT,KRT),"SH2OWGT=",SH2OWGT(IXXRT,JYYRT,KRT),&
                       "SH2O=",SH2OX(I,J,KRT)
                       print*, "SMC=", SMC(i,j,KRT), "SICE =", sice(i,j,KRT)
                       print *, "VEGTYP = ", VEGTYP(I,J)
                       print *, "i,j,krt, nsoil",i,j,krt,nsoil
! ADCHANGE: If values are close but not exact, end up with a crash. Force values to match.
                       !IF (SMC(i,j,KRT).EQ.sice(i,j,KRT)) THEN
                       IF (ABS(SMC(i,j,KRT) - sice(i,j,KRT)) .LE. 0.00001) THEN
                               print *, "SMC = SICE, soil layer totally frozen, proceeding..."
			       SMCRT(IXXRT,JYYRT,KRT) = 0.001
			       sice(i,j,KRT) = SMC(i,j,KRT)
                       ELSE
                               call hydro_stop("In disaggregateDomain() - SMCRT depleted")
                       END IF
               END IF
             end do !debug loop



!DJG map ov roughness as function of land use provided in VEGPARM.TBL...
! --- added extra check for VEGTYP for 'masked-out' locations...
! --- out of basin locations (VEGTYP=0) have OVROUGH hardwired to 0.1
            IF (VEGTYP(I,J).LE.0) then
              OVROUGHRT(IXXRT,JYYRT) = 0.1     !COWS mask test
            ELSE
               OVROUGHRT(IXXRT,JYYRT)=OV_ROUGH(VEGTYP(I,J))*OVROUGHRTFAC(IXXRT,JYYRT)  ! Distributed calibration...1/17/2012
            END IF



!DJG 6.12.08 Map lateral hydraulic conductivity and apply distributed scaling
! ---        factor that will be read in from hires terrain file
!              LKSATRT(IXXRT,JYYRT) = LKSAT(I,J) 
!              LKSATRT(IXXRT,JYYRT) = LKSAT(I,J) * LKSATFAC(IXXRT,JYYRT) * &  !Apply scaling factor...
! ...and scale from max to 0 when SMC decreases from SMCMAX to SMCREF...
!!DJG error found from KIT,improper scaling       ((SMCMAXRT(IXXRT,JYYRT,NSOIL) - SMCRT(IXXRT,JYYRT,NSOIL)) / &
!                                    (max(0.,(SMCMAXRT(IXXRT,JYYRT,NSOIL) - SMCRT(IXXRT,JYYRT,NSOIL))) / &
!                                    (SMCMAXRT(IXXRT,JYYRT,NSOIL)-SMCREFRT(IXXRT,JYYRT,NSOIL)) )

!AD_CHANGE: 
!New model corrected to scale from 0 at SMCREF to full LKSAT*LKSATFAC at SMCMAX:
		LKSATRT(IXXRT,JYYRT) = LKSAT(I,J) * LKSATFAC(IXXRT,JYYRT) * &
				min (1., &     !just in case, make sure scale factor doesn't go over 1
					( max(0.,(SMCRT(IXXRT,JYYRT,NSOIL) - SMCREFRT(IXXRT,JYYRT,NSOIL))) / &     !becomes 0 if less than SMCREF
				(SMCMAXRT(IXXRT,JYYRT,NSOIL)-SMCREFRT(IXXRT,JYYRT,NSOIL)) ) )

!DJG set up lake mask...
!--- modify to make lake mask large here, but not one of the routed lakes!!!
!--            IF (VEGTYP(I,J).eq.16) then
               IF (VEGTYP(I,J).eq.16 .and. &
                        CH_NETRT(IXXRT,JYYRT).le.0) then
                 !--LAKE_MSKRT(IXXRT,JYYRT) = 1
!yw                 LAKE_MSKRT(IXXRT,JYYRT) = 9999
                 LAKE_MSKRT(IXXRT,JYYRT) = -9999
               end if
               ! BF disaggregate soiltype information for gw-soil-coupling
	       ! TODO: move this disaggregation code line to lsm_init section because soiltype is time-invariant
               soiltypRT(ixxrt,jyyrt) = soiltyp(i,j)


              end do
             end do

          end do
        end do

#ifdef HYDRO_D
! ADCHANGE: START Final water balance variables
! ALL VARS in MM
        suminfxsrt2 = 0.
        smcrttot2 = 0.
        do i=1,IXRT
          do j=1,JXRT
            suminfxsrt2 = suminfxsrt2 + INFXSUBRT(I,J) / float(IXRT*JXRT)
            do kk=1,NSOIL
                smcrttot2 = smcrttot2 + SMCRT(I,J,KK)*SLDPTH(KK)*1000. / float(IXRT*JXRT)
            end do
          end do
        end do

#ifdef MPP_LAND
! not tested
        CALL sum_real1(suminfxsrt2)
        CALL sum_real1(smcrttot2)
       suminfxsrt2 = suminfxsrt2/float(numprocs)
       smcrttot2 = smcrttot2/float(numprocs)
#endif
#ifdef MPP_LAND   
       if(my_id .eq. IO_id) then
#endif
        print *, "Disagg Mass Bal: "
        print *, "WB_DISAG!InfxsDiff", suminfxsrt2-suminfxs1
        print *, "WB_DISAG!Infxs1", suminfxs1
        print *, "WB_DISAG!Infxs2", suminfxsrt2
        print *, "WB_DISAG!SMCDIff", smcrttot2-(smctot1-sicetot1)
        print *, "WB_DISAG!SMC1", smctot1
        print *, "WB_DISAG!SICE1", sicetot1
        print *, "WB_DISAG!SMC2", smcrttot2
        print *, "WB_DISAG!Residual", (suminfxsrt2-suminfxs1) + &
                         (smcrttot2-(smctot1-sicetot1))
#ifdef MPP_LAND
       endif
#endif
! END Final water balance variables
#endif

#ifdef HYDRO_D
	print *, "After Disaggregation..."
#endif

#ifdef MPP_LAND
        call MPP_LAND_COM_REAL(INFXSUBRT,IXRT,JXRT,99)
        call MPP_LAND_COM_REAL(LKSATRT,IXRT,JXRT,99)
        call MPP_LAND_COM_REAL(OVROUGHRT,IXRT,JXRT,99)
        call MPP_LAND_COM_INTEGER(LAKE_MSKRT,IXRT,JXRT,99)
     do i = 1, NSOIL
        call MPP_LAND_COM_REAL(SMCMAXRT(:,:,i),IXRT,JXRT,99)
        call MPP_LAND_COM_REAL(SMCRT(:,:,i),IXRT,JXRT,99)
        call MPP_LAND_COM_REAL(SMCWLTRT(:,:,i),IXRT,JXRT,99)
     end DO
#endif

     end subroutine disaggregateDomain

         subroutine SubsurfaceRouting_drv(did)

             use module_RT_data, only: rt_domain
             use module_namelist, only: nlst_rt
             implicit none
             integer :: did
             IF (nlst_rt(did)%SUBRTSWCRT.EQ.1) THEN
                call subsurfaceRouting (RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt , nlst_rt(did)%nsoil, &
                     RT_DOMAIN(did)%SMCRT,RT_DOMAIN(did)%SMCMAXRT,RT_DOMAIN(did)%SMCREFRT,&
                     RT_DOMAIN(did)%SMCWLTRT, RT_DOMAIN(did)%ZSOIL, RT_DOMAIN(did)%SLDPTH, & 
                     nlst_rt(did)%DT,RT_DOMAIN(did)%ZWATTABLRT,RT_DOMAIN(did)%SOXRT,       &
                     RT_DOMAIN(did)%SOYRT,RT_DOMAIN(did)%LKSATRT, RT_DOMAIN(did)%SOLDEPRT, &
                     RT_DOMAIN(did)%INFXSUBRT,RT_DOMAIN(did)%QSUBBDRYTRT, RT_DOMAIN(did)%QSUBBDRYRT,&
                     RT_DOMAIN(did)%QSUBRT ,nlst_rt(did)%rt_option, RT_DOMAIN(did)%dist, &
                     RT_DOMAIN(did)%sub_resid,RT_DOMAIN(did)%SO8RT_D, RT_DOMAIN(did)%SO8RT)
             endif

         end subroutine SubsurfaceRouting_drv
     
         subroutine subsurfaceRouting (ixrt, jxrt , nsoil, &
                  SMCRT,SMCMAXRT,SMCREFRT,SMCWLTRT, &
                  ZSOIL, SLDPTH, &
                  DT,ZWATTABLRT,SOXRT,SOYRT,LKSATRT,&
                  SOLDEPRT,INFXSUBRT,QSUBBDRYTRT, QSUBBDRYRT,&
                  QSUBRT ,rt_option, dist,sub_resid,SO8RT_D, SO8RT)
#ifdef MPP_LAND
        use module_mpp_land, only:  mpp_land_com_real, mpp_land_com_integer
#endif
         implicit none
         integer, INTENT(IN) :: ixrt, jxrt , nsoil, rt_option
         REAL, INTENT(IN)                          :: DT
         real,INTENT(IN), DIMENSION(NSOIL)      :: ZSOIL, SLDPTH
         REAL, INTENT(IN), DIMENSION(IXRT,JXRT)   :: SOXRT,SOYRT,LKSATRT, SOLDEPRT , sub_resid
         real,INTENT(INOUT), DIMENSION(IXRT,JXRT)::INFXSUBRT
         real,INTENT(INOUT) :: QSUBBDRYTRT
         REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)   :: QSUBBDRYRT, QSUBRT
         REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT,NSOIL) :: SMCRT, SMCWLTRT, SMCMAXRT,SMCREFRT


         INTEGER :: SO8RT_D(IXRT,JXRT,3)
         REAL :: SO8RT(IXRT,JXRT,8)
         REAL, INTENT(IN)                          :: dist(ixrt,jxrt,9)
!  -----local array ----------
         REAL, DIMENSION(IXRT,JXRT)   :: ZWATTABLRT
         REAL, DIMENSION(IXRT,JXRT)   :: CWATAVAIL
         INTEGER, DIMENSION(IXRT,JXRT) :: SATLYRCHK
 



         CWATAVAIL = 0.
         CALL FINDZWAT(IXRT,JXRT,NSOIL,SMCRT,SMCMAXRT,SMCREFRT, &
                             SMCWLTRT,ZSOIL,SATLYRCHK,ZWATTABLRT, &
                             CWATAVAIL,SLDPTH)
#ifdef MPP_LAND
        call MPP_LAND_COM_REAL(ZWATTABLRT,IXRT,JXRT,99)
        call MPP_LAND_COM_REAL(CWATAVAIL,IXRT,JXRT,99)
        call MPP_LAND_COM_INTEGER(SATLYRCHK,IXRT,JXRT,99)
#endif


!DJG Second, Call subsurface routing routine...
#ifdef HYDRO_D
	print *, "Beginning SUB_routing..."
        print *, "Routing method is ",rt_option, " direction."
#endif

!!!! Find saturated layer depth...
! Loop through domain to determine sat. layers and assign wat tbl depth...
!    and water available for subsfc routing (CWATAVAIL)...
! This subroutine returns: ZWATTABLRT, CWATAVAIL and SATLYRCHK


    CALL SUBSFC_RTNG(dist,ZWATTABLRT,QSUBRT,SOXRT,SOYRT,  &
          LKSATRT,SOLDEPRT,QSUBBDRYRT,QSUBBDRYTRT,NSOIL,SMCRT,     &
          INFXSUBRT,SMCMAXRT,SMCREFRT,ZSOIL,IXRT,JXRT,DT,SMCWLTRT,SO8RT,    &
          SO8RT_D, rt_option,SLDPTH,SUB_RESID,CWATAVAIL,SATLYRCHK)

#ifdef HYDRO_D
    print *, "SUBROUTE routing called and returned..."
#endif

    end subroutine subsurfaceRouting 

   
       subroutine OverlandRouting_drv(did)
             use module_RT_data, only: rt_domain
             use module_namelist, only: nlst_rt
             implicit none
             integer :: did
             if(nlst_rt(did)%OVRTSWCRT .eq. 1) then
                 call OverlandRouting (nlst_rt(did)%DT, nlst_rt(did)%DTRT_TER, nlst_rt(did)%rt_option, &
                          rt_domain(did)%ixrt, rt_domain(did)%jxrt,rt_domain(did)%LAKE_MSKRT, &
                          rt_domain(did)%INFXSUBRT, rt_domain(did)%RETDEPRT,rt_domain(did)%OVROUGHRT, &
                          rt_domain(did)%SOXRT, rt_domain(did)%SOYRT, rt_domain(did)%SFCHEADSUBRT,  &
                          rt_domain(did)%DHRT, rt_domain(did)%CH_NETRT, rt_domain(did)%QSTRMVOLRT, &
                          rt_domain(did)%LAKE_INFLORT,rt_domain(did)%QBDRYRT, &
                          rt_domain(did)%QSTRMVOLTRT,rt_domain(did)%QBDRYTRT, rt_domain(did)%LAKE_INFLOTRT,&
                          rt_domain(did)%q_sfcflx_x,rt_domain(did)%q_sfcflx_y, &
                          rt_domain(did)%dist, rt_domain(did)%SO8RT, rt_domain(did)%SO8RT_D , &
                          rt_domain(did)%SMCTOT2,rt_domain(did)%suminfxs1,rt_domain(did)%suminfxsrt, &
                          rt_domain(did)%smctot1,rt_domain(did)%dsmctot )
! ADCHANGE: If overland routing is called, INFXSUBRT is moved to SFCHEADSUBRT, so 
!           zeroing out just in case
             rt_domain(did)%INFXSUBRT = 0.0
             endif
       end subroutine OverlandRouting_drv



       subroutine OverlandRouting (DT, DTRT_TER, rt_option, ixrt, jxrt,LAKE_MSKRT, &
                  INFXSUBRT, RETDEPRT,OVROUGHRT,SOXRT, SOYRT, SFCHEADSUBRT,DHRT, &
                  CH_NETRT, QSTRMVOLRT,LAKE_INFLORT,QBDRYRT, &
                  QSTRMVOLTRT,QBDRYTRT, LAKE_INFLOTRT, q_sfcflx_x,q_sfcflx_y, &
                  dist, SO8RT, SO8RT_D, &
                  SMCTOT2,suminfxs1,suminfxsrt,smctot1,dsmctot )
#ifdef MPP_LAND
        use module_mpp_land, only:  mpp_land_max_int1,  sum_real1, my_id, io_id, numprocs
#endif
       implicit none

       REAL, INTENT(IN) :: DT, DTRT_TER
       integer, INTENT(IN) :: ixrt, jxrt, rt_option
       INTEGER, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: LAKE_MSKRT

       REAL, INTENT(IN), DIMENSION(IXRT,JXRT)   :: INFXSUBRT,  &
                 RETDEPRT,OVROUGHRT,SOXRT, SOYRT
       REAL, INTENT(OUT), DIMENSION(IXRT,JXRT) :: SFCHEADSUBRT,DHRT
       INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETRT
       REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: QSTRMVOLRT,LAKE_INFLORT,QBDRYRT, &
                QSTRMVOLTRT,QBDRYTRT, LAKE_INFLOTRT, q_sfcflx_x,q_sfcflx_y

       REAL, INTENT(IN), DIMENSION(IXRT,JXRT,9):: dist
       REAL, INTENT(IN), DIMENSION(IXRT,JXRT,8)  :: SO8RT
       INTEGER SO8RT_D(IXRT,JXRT,3)

       integer  :: i,j
      

       real            :: smctot2,smctot1,dsmctot
       real            :: suminfxsrt,suminfxs1
! local variable
       real            :: chan_in1,chan_in2
       real            :: lake_in1,lake_in2
       real            :: qbdry1,qbdry2
       integer :: sfcrt_flag



!DJG Third, Call Overland Flow Routing Routine...
#ifdef HYDRO_D
	print *, "Beginning OV_routing..."
        print *, "Routing method is ",rt_option, " direction."
#endif

!DJG debug...OV Routing...
	suminfxs1=0.
        chan_in1=0.
        lake_in1=0.
        qbdry1=0.
        do i=1,IXRT
         do j=1,JXRT
            suminfxs1=suminfxs1+INFXSUBRT(I,J)/float(IXRT*JXRT)
            chan_in1=chan_in1+QSTRMVOLRT(I,J)/float(IXRT*JXRT)
            lake_in1=lake_in1+LAKE_INFLORT(I,J)/float(IXRT*JXRT)
            qbdry1=qbdry1+QBDRYRT(I,J)/float(IXRT*JXRT)
         end do
        end do

#ifdef MPP_LAND
! not tested
        CALL sum_real1(suminfxs1)
        CALL sum_real1(chan_in1)
        CALL sum_real1(lake_in1)
        CALL sum_real1(qbdry1)
        suminfxs1 = suminfxs1/float(numprocs)
        chan_in1 = chan_in1/float(numprocs)
        lake_in1 = lake_in1/float(numprocs)
        qbdry1 = qbdry1/float(numprocs)
#endif


!DJG.7.20.2007 - Global check for infxs>retdep & skip if not...(set sfcrt_flag)
!DJG.7.20.2007 - this check will skip ov rtng when no flow is present...
        
        sfcrt_flag = 0
        
        do j=1,jxrt
          do i=1,ixrt
            if(INFXSUBRT(i,j).gt.RETDEPRT(i,j)) then
              sfcrt_flag = 1
              exit
            end if
          end do
          if(sfcrt_flag.eq.1) exit
        end do   

#ifdef MPP_LAND
       call mpp_land_max_int1(sfcrt_flag)            
#endif
!DJG.7.20.2007 - Global check for infxs>retdep & skip if not...(IF)

    if (sfcrt_flag.eq.1) then  !If/then for sfc_rt check...
#ifdef HYDRO_D
      write(6,*) "calling OV_RTNG "
#endif
      CALL OV_RTNG(DT,DTRT_TER,IXRT,JXRT,INFXSUBRT,SFCHEADSUBRT,DHRT,      &
        CH_NETRT,RETDEPRT,OVROUGHRT,QSTRMVOLRT,QBDRYRT,              &
        QSTRMVOLTRT,QBDRYTRT,SOXRT,SOYRT,dist,                       &
        LAKE_MSKRT,LAKE_INFLORT,LAKE_INFLOTRT,SO8RT,SO8RT_D,rt_option,&
        q_sfcflx_x,q_sfcflx_y) 
    else
      SFCHEADSUBRT = INFXSUBRT
#ifdef HYDRO_D
      print *, "No water to route overland..."
#endif
    end if  !Endif for sfc_rt check...

!DJG.7.20.2007 - Global check for infxs>retdep & skip if not...(ENDIF)

#ifdef HYDRO_D
    print *, "OV routing called and returned..."
#endif

!DJG Debug...OV Routing...
	suminfxsrt=0.
        chan_in2=0.
        lake_in2=0.
        qbdry2=0.
        do i=1,IXRT
         do j=1,JXRT
            suminfxsrt=suminfxsrt+SFCHEADSUBRT(I,J)/float(IXRT*JXRT)
            chan_in2=chan_in2+QSTRMVOLRT(I,J)/float(IXRT*JXRT)
            lake_in2=lake_in2+LAKE_INFLORT(I,J)/float(IXRT*JXRT)
            qbdry2=qbdry2+QBDRYRT(I,J)/float(IXRT*JXRT)
         end do
        end do
#ifdef MPP_LAND
! not tested
        CALL sum_real1(suminfxsrt)
        CALL sum_real1(chan_in2)
        CALL sum_real1(lake_in2)
        CALL sum_real1(qbdry2)
        suminfxsrt = suminfxsrt/float(numprocs)
        chan_in2 = chan_in2/float(numprocs)
        lake_in2 = lake_in2/float(numprocs)
        qbdry2 = qbdry2/float(numprocs)
#endif

#ifdef HYDRO_D
#ifdef MPP_LAND   
       if(my_id .eq. IO_id) then
#endif
	print *, "OV Routing Mass Bal: "
        print *, "WB_OV!InfxsDiff", suminfxsrt-suminfxs1
        print *, "WB_OV!Infxs1", suminfxs1
        print *, "WB_OV!Infxs2", suminfxsrt
        print *, "WB_OV!ChaninDiff", chan_in2-chan_in1
        print *, "WB_OV!Chanin1", chan_in1
        print *, "WB_OV!Chanin2", chan_in2
        print *, "WB_OV!LakeinDiff", lake_in2-lake_in1
        print *, "WB_OV!Lakein1", lake_in1
        print *, "WB_OV!Lakein2", lake_in2
        print *, "WB_OV!QbdryDiff", qbdry2-qbdry1
        print *, "WB_OV!Qbdry1", qbdry1
        print *, "WB_OV!Qbdry2", qbdry2
        print *, "WB_OV!Residual", (suminfxs1-suminfxsrt)-(chan_in2-chan_in1) &
                      -(lake_in2-lake_in1)-(qbdry2-qbdry1)
#ifdef MPP_LAND
       endif
#endif
#endif


       end subroutine OverlandRouting


      subroutine time_seconds(i3)
          integer time_array(8)
          real*8 i3
          call date_and_time(values=time_array)
          i3 = time_array(4)*24*3600+time_array(5) * 3600 + time_array(6) * 60 + &
                time_array(7) + 0.001 * time_array(8)
          return
      end subroutine time_seconds