hydro/Routing/Noah_distr_routing_overland.F

   1 !DJG ------------------------------------------------
   2 !DJG   SUBROUTINE OverlandRouting
   3 !DJG ------------------------------------------------
   4
   5 !subroutine OverlandRouting (DT, DTRT_TER, rt_option, ixrt, jxrt,LAKE_MSKRT, &
   6     !           INFXSUBRT, RETDEPRT,OVROUGHRT,SOXRT, SOYRT, SFCHEADSUBRT,DHRT, &
   7     !           CH_NETRT, QSTRMVOLRT,LAKE_INFLORT,QBDRYRT, &
   8     !           QSTRMVOLTRT,QBDRYTRT, LAKE_INFLOTRT, q_sfcflx_x,q_sfcflx_y, &
   9     !           dist, SO8RT, SO8RT_D, &
  10     !           SMCTOT2,suminfxs1,suminfxsrt,smctot1,dsmctot )
  11
  12
  13 subroutine OverlandRouting( &
  14     ovrt_data, & ! overland data structure
  15     DT, &
  16     DTRT_TER, &
  17     rt_option, &
  18     ixrt, &       ! routing grid x size
  19     jxrt, &      ! routing grid y size
  20     q_sfcflx_x, &! accumulated x flux
  21     q_sfcflx_y  &! accumulated y flux
  22     )
  23 #ifdef MPP_LAND
  24     use module_mpp_land, only:  mpp_land_max_int1,  sum_real1, my_id, io_id, numprocs
  25     use overland_data
  26 #endif
  27     implicit none
  28
  29     type (overland_struct), intent(inout) :: ovrt_data
  30     REAL, INTENT(IN) :: DT, DTRT_TER
  31     integer, INTENT(IN) :: ixrt, jxrt, rt_option
  32     !REAL, INTENT(out), DIMENSION(IXRT,JXRT)   :: SFCHEADSUBRT  ! moved into overland_data%control
  33     real, dimension(IXRT,JXRT), intent(inout) :: q_sfcflx_x,q_sfcflx_y
  34
  35     ! REMOVED VARIABLES
  36     !INTEGER, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: LAKE_MSKRT
  37
  38     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)   :: INFXSUBRT,  &
  39         !          RETDEPRT,OVROUGHRT,SOXRT, SOYRT
  40     !REAL, INTENT(OUT), DIMENSION(IXRT,JXRT) :: SFCHEADSUBRT,DHRT
  41     !INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETRT
  42     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: QSTRMVOLRT,LAKE_INFLORT,QBDRYRT, &
  43         !         QSTRMVOLTRT,QBDRYTRT, LAKE_INFLOTRT, q_sfcflx_x,q_sfcflx_y
  44
  45     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT,9):: dist
  46     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT,8)  :: SO8RT
  47     !INTEGER SO8RT_D(IXRT,JXRT,3)
  48
  49     integer  :: i,j
  50
  51
  52     !real            :: smctot2,smctot1,dsmctot
  53     !real            :: suminfxsrt,suminfxs1
  54     ! local variable
  55     real            :: chan_in1,chan_in2
  56     real            :: lake_in1,lake_in2
  57     real            :: qbdry1,qbdry2
  58     integer :: sfcrt_flag
  59
  60
  61
  62     !DJG Third, Call Overland Flow Routing Routine...
  63 #ifdef HYDRO_D
  64     print *, "Beginning OV_routing..."
  65     print *, "Routing method is ",rt_option, " direction."
  66 #endif
  67
  68     !DJG debug...OV Routing...
  69     ovrt_data%mass_balance%pre_infiltration_excess = 0.
  70     chan_in1=0.
  71     lake_in1=0.
  72     qbdry1=0.
  73     do i=1,IXRT
  74         do j=1,JXRT
  75             ovrt_data%mass_balance%pre_infiltration_excess = &
  76                 ovrt_data%mass_balance%pre_infiltration_excess + ( ovrt_data%control%infiltration_excess(I,J)/float(IXRT*JXRT) )
  77             chan_in1 = chan_in1 + ( ovrt_data%streams_and_lakes%surface_water_to_channel(I,J)/float(IXRT*JXRT) )
  78             lake_in1 = lake_in1 + ( ovrt_data%streams_and_lakes%surface_water_to_lake(I,J)/float(IXRT*JXRT) )
  79             qbdry1 = qbdry1 + ( ovrt_data%control%boundary_flux(I,J)/float(IXRT*JXRT) )
  80         end do
  81     end do
  82
  83 #ifdef MPP_LAND
  84     ! not tested
  85     CALL sum_real1(ovrt_data%mass_balance%pre_infiltration_excess)
  86     CALL sum_real1(chan_in1)
  87     CALL sum_real1(lake_in1)
  88     CALL sum_real1(qbdry1)
  89     ovrt_data%mass_balance%pre_infiltration_excess = ovrt_data%mass_balance%pre_infiltration_excess/float(numprocs)
  90     chan_in1 = chan_in1/float(numprocs)
  91     lake_in1 = lake_in1/float(numprocs)
  92     qbdry1 = qbdry1/float(numprocs)
  93 #endif
  94
  95
  96     !DJG.7.20.2007 - Global check for infxs>retdep & skip if not...(set sfcrt_flag)
  97     !DJG.7.20.2007 - this check will skip ov rtng when no flow is present...
  98
  99     sfcrt_flag = 0
 100
 101     do j=1,jxrt
 102         do i=1,ixrt
 103             if( ovrt_data%control%infiltration_excess(i,j) .gt. ovrt_data%properties%retention_depth(i,j) ) then
 104                 sfcrt_flag = 1
 105                 exit
 106             end if
 107         end do
 108         if(sfcrt_flag .eq. 1) exit
 109     end do
 110
 111 #ifdef MPP_LAND
 112     call mpp_land_max_int1(sfcrt_flag)
 113 #endif
 114     !DJG.7.20.2007 - Global check for infxs>retdep & skip if not...(IF)
 115
 116     if (sfcrt_flag.eq.1) then  !If/then for sfc_rt check...
 117 #ifdef HYDRO_D
 118         write(6,*) "calling OV_RTNG "
 119 #endif
 120         CALL OV_RTNG(           &
 121             ovrt_data, &
 122             DT,                   &
 123             DTRT_TER,             &
 124             IXRT,                 &
 125             JXRT,                 &
 126             rt_option,            &
 127             q_sfcflx_x,           &
 128             q_sfcflx_y)
 129     else
 130         ovrt_data%control%surface_water_head_routing = ovrt_data%control%infiltration_excess
 131 #ifdef HYDRO_D
 132         print *, "No water to route overland..."
 133 #endif
 134     end if  !Endif for sfc_rt check...
 135
 136     !DJG.7.20.2007 - Global check for infxs>retdep & skip if not...(ENDIF)
 137
 138 #ifdef HYDRO_D
 139     print *, "OV routing called and returned..."
 140 #endif
 141
 142     !DJG Debug...OV Routing...
 143     ovrt_data%mass_balance%post_infiltration_excess = 0.
 144     chan_in2=0.
 145     lake_in2=0.
 146     qbdry2=0.
 147     do i=1,IXRT
 148         do j=1,JXRT
 149             ovrt_data%mass_balance%post_infiltration_excess = &
 150                 ovrt_data%mass_balance%post_infiltration_excess + (ovrt_data%control%surface_water_head_routing(I,J)/float(IXRT*JXRT))
 151             chan_in2=chan_in2 + (ovrt_data%streams_and_lakes%surface_water_to_channel(I,J)/float(IXRT*JXRT))
 152             lake_in2=lake_in2 + (ovrt_data%streams_and_lakes%surface_water_to_lake(I,J)/float(IXRT*JXRT))
 153             qbdry2=qbdry2 + (ovrt_data%control%boundary_flux(I,J)/float(IXRT*JXRT))
 154         end do
 155     end do
 156 #ifdef MPP_LAND
 157     ! not tested
 158     CALL sum_real1(ovrt_data%mass_balance%post_infiltration_excess)
 159     CALL sum_real1(chan_in2)
 160     CALL sum_real1(lake_in2)
 161     CALL sum_real1(qbdry2)
 162     ovrt_data%mass_balance%post_infiltration_excess = ovrt_data%mass_balance%post_infiltration_excess/float(numprocs)
 163     chan_in2 = chan_in2/float(numprocs)
 164     lake_in2 = lake_in2/float(numprocs)
 165     qbdry2 = qbdry2/float(numprocs)
 166 #endif
 167
 168 #ifdef HYDRO_D
 169 #ifdef MPP_LAND
 170     if(my_id .eq. IO_id) then
 171 #endif
 172         print *, "OV Routing Mass Bal: "
 173         print *, "WB_OV!InfxsDiff", ovrt_data%mass_balance%post_infiltration_excess - ovrt_data%mass_balance%pre_infiltration_excess
 174         print *, "WB_OV!Infxs1", ovrt_data%mass_balance%pre_infiltration_excess
 175         print *, "WB_OV!Infxs2", ovrt_data%mass_balance%post_infiltration_excess
 176         print *, "WB_OV!ChaninDiff", chan_in2-chan_in1
 177         print *, "WB_OV!Chanin1", chan_in1
 178         print *, "WB_OV!Chanin2", chan_in2
 179         print *, "WB_OV!LakeinDiff", lake_in2-lake_in1
 180         print *, "WB_OV!Lakein1", lake_in1
 181         print *, "WB_OV!Lakein2", lake_in2
 182         print *, "WB_OV!QbdryDiff", qbdry2-qbdry1
 183         print *, "WB_OV!Qbdry1", qbdry1
 184         print *, "WB_OV!Qbdry2", qbdry2
 185         print *, "WB_OV!Residual", (ovrt_data%mass_balance%post_infiltration_excess - ovrt_data%mass_balance%pre_infiltration_excess)&
 186             -(chan_in2-chan_in1) -(lake_in2 - lake_in1) - (qbdry2 - qbdry1)
 187 #ifdef MPP_LAND
 188     endif
 189 #endif
 190 #endif
 191
 192
 193 end subroutine OverlandRouting
 194
 195 !DJG ------------------------------------------------
 196 !DJG   SUBROUTINE OV_RTNG
 197 !DJG ------------------------------------------------
 198
 199 !SUBROUTINE OV_RTNG(DT,DTRT_TER,IXRT,JXRT,INFXSUBRT,      &
 200     !  SFCHEADSUBRT,DHRT,CH_NETRT,RETDEPRT,OVROUGHRT,      &
 201     !  QSTRMVOLRT,QBDRYRT,QSTRMVOLTRT,QBDRYTRT,SOXRT,     &
 202     !  SOYRT,dist,LAKE_MSKRT,LAKE_INFLORT,LAKE_INFLOTRT,  &
 203     !  SO8RT,SO8RT_D,rt_option,q_sfcflx_x,q_sfcflx_y)
 204
 205 subroutine ov_rtng( &
 206         ovrt_data, &
 207         DT, &
 208         DTRT_TER, &
 209         IXRT, &
 210         JXRT, &
 211         rt_option, &
 212         q_sfcflx_x, &
 213         q_sfcflx_y &
 214         )
 215
 216     !yyww
 217 #ifdef MPP_LAND
 218     use module_mpp_land, only: left_id,down_id,right_id, &
 219         up_id,mpp_land_com_real, my_id, &
 220         mpp_land_sync
 221 #endif
 222     use overland_data
 223     IMPLICIT NONE
 224
 225     !DJG --------DECLARATIONS----------------------------
 226
 227     type (overland_struct), intent(inout) :: ovrt_data
 228     INTEGER, INTENT(IN)                 :: IXRT,JXRT
 229     REAL, INTENT(IN)                    :: DT,DTRT_TER
 230     integer, intent(in) :: rt_option
 231
 232     !INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: CH_NETRT
 233     !INTEGER, INTENT(IN), DIMENSION(IXRT,JXRT) :: LAKE_MSKRT
 234     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)     :: INFXSUBRT
 235     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)     :: SOXRT
 236     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)     :: SOYRT
 237     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT,9):: dist
 238     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT)  :: RETDEPRT
 239     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT)     :: OVROUGHRT
 240
 241     !REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)    :: SFCHEADSUBRT
 242     !REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)    :: DHRT
 243
 244     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: QSTRMVOLRT
 245     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: LAKE_INFLORT
 246     !REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: QBDRYRT
 247     REAL, INTENT(INOUT), DIMENSION(IXRT,JXRT) :: q_sfcflx_x,q_sfcflx_y
 248     !REAL, INTENT(INOUT)     :: QSTRMVOLTRT,QBDRYTRT,LAKE_INFLOTRT
 249     !REAL, INTENT(IN), DIMENSION(IXRT,JXRT,8)  :: SO8RT
 250
 251     !DJG Local Variables
 252
 253     INTEGER :: KRT,I,J,ct
 254
 255     REAL, DIMENSION(IXRT,JXRT)  :: INFXS_FRAC
 256     REAL        :: DT_FRAC,SUM_INFXS,sum_head
 257     !INTEGER SO8RT_D(IXRT,JXRT,3), rt_option
 258
 259     !DJG ----------------------------------------------------------------------
 260     ! DJG BEGIN 1-D or 2-D OVERLAND FLOW ROUTING LOOP
 261     !DJG ---------------------------------------------------------------------
 262     !DJG  Loop over 'routing time step'
 263     !DJG  Compute the number of time steps based on NOAH DT and routing DTRT_TER
 264
 265     DT_FRAC=INT(DT/DTRT_TER)
 266
 267 #ifdef HYDRO_D
 268     write(6,*) "OV_RTNG  DT_FRAC, DT, DTRT_TER",DT_FRAC, DT, DTRT_TER
 269     write(6,*) "IXRT, JXRT = ",ixrt,jxrt
 270 #endif
 271
 272     !DJG NOTE: Applying all infiltration excess water at once then routing
 273     !DJG       Pre-existing SFHEAD gets combined with Precip. in the
 274     !DJG       calculation of INFXS1 during subroutine SRT.f.
 275     !DJG debug
 276
 277
 278     !DJG Assign all infiltration excess to surface head...
 279     ovrt_data%control%surface_water_head_routing=ovrt_data%control%infiltration_excess
 280
 281     !DJG Divide infiltration excess over all routing time-steps
 282     !        INFXS_FRAC=INFXSUBRT/(DT/DTRT_TER)
 283
 284     !DJG Set flux accumulation fields to 0. before each loop...
 285     q_sfcflx_x = 0.
 286     q_sfcflx_y = 0.
 287     ct =0
 288
 289
 290     !DJG Execute routing time-step loop...
 291
 292
 293     DO KRT=1,DT_FRAC
 294
 295         DO J=1,JXRT
 296             DO I=1,IXRT
 297
 298                 !DJG Removed 4_29_05, sfhead now updated in route_overland subroutine...
 299                 !           SFCHEADSUBRT(I,J)=SFCHEADSUBRT(I,J)+DHRT(I,J)
 300                 !!           SFCHEADSUBRT(I,J)=SFCHEADSUBRT(I,J)+DHRT(I,J)+INFXS_FRAC(I,J)
 301                 !           DHRT(I,J)=0.
 302
 303                 !DJG ERROR Check...
 304
 305                 IF (ovrt_data%control%surface_water_head_routing(I,J).lt.0.) THEN
 306 #ifdef HYDRO_D
 307                     print *, "ywcheck 2 ERROR!!!: Neg. Surface Head Value at (i,j):",    &
 308                         i,j,ovrt_data%control%surface_water_head_routing(I,J)
 309                     print *, "RETDEPRT(I,J) = ",ovrt_data%properties%retention_depth(I,J), "KRT=",KRT
 310                     print *, "INFXSUBRT(i,j)=",ovrt_data%control%infiltration_excess(i,j)
 311                     print *, "jxrt=",jxrt," ixrt=",ixrt
 312 #endif
 313                 END IF
 314
 315                 !DJG Remove surface water from channel cells
 316                 !DJG Channel inflo cells specified as nonzeros from CH_NET
 317                 !DJG 9/16/04  Channel Extractions Removed until stream model implemented...
 318
 319
 320
 321                 !yw            IF (CH_NETRT(I,J).ne.-9999) THEN
 322                 IF (ovrt_data%streams_and_lakes%CH_NETRT(I,J).ge.0) THEN
 323                     ct = ct +1
 324
 325                     !DJG Temporary test to up the retention depth of channel grid cells to 'soak'
 326                     !more water into valleys....set retdep = retdep*100 (=5 mm)
 327
 328                     !        RETDEPRT(I,J) = RETDEPRT(I,J) * 100.0    !DJG TEMP HARDWIRE!!!!
 329                     !        RETDEPRT(I,J) = 10.0    !DJG TEMP HARDWIRE!!!!
 330
 331                     ! AD hardwire to force channel retention depth to be 5mm.
 332                     ovrt_data%properties%retention_depth(I,J) = 5.0
 333
 334                     IF (ovrt_data%control%surface_water_head_routing(I,J) .GT. ovrt_data%properties%retention_depth(I,J)) THEN
 335                         !!               QINFLO(CH_NET(I,J)=QINFLO(CH_NET(I,J)+SFCHEAD(I,J) - RETDEPRT(I,J)
 336                         ovrt_data%streams_and_lakes%accumulated_surface_water_to_channel = ovrt_data%streams_and_lakes%accumulated_surface_water_to_channel + &
 337                             (ovrt_data%control%surface_water_head_routing(I,J) - ovrt_data%properties%retention_depth(I,J))
 338                         ovrt_data%streams_and_lakes%surface_water_to_channel(I,J) = ovrt_data%streams_and_lakes%surface_water_to_channel(I,J) + &
 339                             (ovrt_data%control%surface_water_head_routing(I,J) - ovrt_data%properties%retention_depth(I,J))
 340
 341                         ! if(QSTRMVOLRT(I,J) .gt. 0) then
 342                         !     print *, "QSTRVOL GT 0", QSTRMVOLRT(I,J),I,J
 343                         !  endif
 344
 345                         ovrt_data%control%surface_water_head_routing(I,J) = ovrt_data%properties%retention_depth(I,J)
 346                     END IF
 347                 END IF
 348
 349                 !DJG Lake inflow withdrawl from surface head...(4/29/05)
 350
 351
 352                 IF (ovrt_data%streams_and_lakes%lake_mask(I,J) .gt. 0) THEN
 353                     IF (ovrt_data%control%surface_water_head_routing(I,J) .GT. ovrt_data%properties%retention_depth(I,J)) THEN
 354                         ovrt_data%streams_and_lakes%accumulated_surface_water_to_lake = ovrt_data%streams_and_lakes%accumulated_surface_water_to_lake + &
 355                             (ovrt_data%control%surface_water_head_routing(I,J) - ovrt_data%properties%retention_depth(I,J))
 356                         ovrt_data%streams_and_lakes%surface_water_to_lake(I,J) = ovrt_data%streams_and_lakes%surface_water_to_lake(I,J) + &
 357                             (ovrt_data%control%surface_water_head_routing(I,J)- ovrt_data%properties%retention_depth(I,J))
 358                         ovrt_data%control%surface_water_head_routing(I,J) = ovrt_data%properties%retention_depth(I,J)
 359                     END IF
 360                 END IF
 361
 362
 363
 364             END DO
 365         END DO
 366
 367         !yw check         call MPP_LAND_COM_REAL(QSTRMVOLRT,IXRT,JXRT,99)
 368         !DJG----------------------------------------------------------------------
 369         !DJG CALL OVERLAND FLOW ROUTING SUBROUTINE
 370         !DJG----------------------------------------------------------------------
 371
 372         !DJG Debug...
 373
 374
 375         if(rt_option .eq. 1) then
 376             CALL ROUTE_OVERLAND1(DTRT_TER, &
 377                 ovrt_data%properties%distance_to_neighbor, &
 378                 ovrt_data%control%surface_water_head_routing, &
 379                 ovrt_data%control%DHRT, &
 380                 ovrt_data%properties%surface_slope_x,   &
 381                 ovrt_data%properties%surface_slope_y, &
 382                 ovrt_data%properties%retention_depth, &
 383                 ovrt_data%properties%roughness, &
 384                 IXRT, JXRT, &
 385                 ovrt_data%control%boundary_flux, &
 386                 ovrt_data%control%boundary_flux_total,    &
 387                 ovrt_data%properties%surface_slope, &
 388                 ovrt_data%properties%max_surface_slope_index, &
 389                 q_sfcflx_x, q_sfcflx_y)
 390         else
 391             CALL ROUTE_OVERLAND2(DTRT_TER, &
 392                 ovrt_data%properties%distance_to_neighbor, &
 393                 ovrt_data%control%surface_water_head_routing, &
 394                 ovrt_data%control%DHRT, &
 395                 ovrt_data%properties%surface_slope_x, &
 396                 ovrt_data%properties%surface_slope_y, &
 397                 ovrt_data%properties%retention_depth, &
 398                 ovrt_data%properties%roughness, &
 399                 IXRT, JXRT, &
 400                 ovrt_data%control%boundary_flux, &
 401                 ovrt_data%control%boundary_flux_total,  &
 402                 q_sfcflx_x,q_sfcflx_y)
 403         end if
 404
 405     END DO          ! END routing time steps
 406
 407 #ifdef HYDRO_D
 408     print *, "End of OV_routing call..."
 409 #endif
 410
 411     !----------------------------------------------------------------------
 412     ! END OVERLAND FLOW ROUTING LOOP
 413     !     CHANNEL ROUTING TO FOLLOW
 414     !----------------------------------------------------------------------
 415
 416     !DJG ----------------------------------------------------------------
 417 END SUBROUTINE OV_RTNG
 418
 419 !DJG ----------------------------------------------------------------
 420
 421 !DJG     SUBROUTINE ROUTE_OVERLAND1
 422 !DJG ----------------------------------------------------------------
 423
 424 SUBROUTINE ROUTE_OVERLAND1(dt,                                &
 425         &          gsize,h,qsfc,sox,soy,                                   &
 426         &     retent_dep,dist_rough,XX,YY,QBDRY,QBDRYT,SO8RT,SO8RT_D,      &
 427         &     q_sfcflx_x,q_sfcflx_y)
 428     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 429     !
 430     !  Subroutine to route excess rainfall over the watershed
 431     !     using a 2d diffusion routing scheme.
 432     !
 433     !  Called from: main.f
 434     !
 435     !      Will try to formulate this to be called from NOAH
 436     !
 437     !  Returns: qsfc=DQOV   which in turn becomes DH in head calc.
 438     !
 439     !  Created:  Adaptded from CASC2D source code
 440     !  NOTE: dh from original code has been replaced by qsfc
 441     !        dhh replaced by qqsfc
 442     !
 443     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 444
 445 #ifdef MPP_LAND
 446     use module_mpp_land, only: left_id,down_id,right_id, &
 447         up_id,mpp_land_com_real, my_id, mpp_land_com_real8,&
 448         mpp_land_sync
 449 #endif
 450
 451     IMPLICIT NONE
 452
 453
 454     !! Declare Passed variables
 455
 456     INTEGER, INTENT(IN) :: XX,YY
 457     REAL, INTENT(IN) :: dt, gsize(xx,yy,9)
 458
 459     !! Declare passed arrays
 460
 461     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: h
 462     REAL, INTENT(IN), DIMENSION(XX,YY) :: qsfc
 463     REAL, INTENT(IN), DIMENSION(XX,YY) :: sox
 464     REAL, INTENT(IN), DIMENSION(XX,YY) :: soy
 465     REAL, INTENT(IN), DIMENSION(XX,YY) :: retent_dep
 466     REAL, INTENT(IN), DIMENSION(XX,YY) :: dist_rough
 467     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: QBDRY
 468     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: q_sfcflx_x, q_sfcflx_y
 469     REAL, INTENT(INOUT) :: QBDRYT
 470     REAL, INTENT(IN), DIMENSION(XX,YY,8) :: SO8RT
 471     REAL*8, DIMENSION(XX,YY) :: QBDRY_tmp, DH
 472     REAL*8, DIMENSION(XX,YY) :: DH_tmp
 473     REAL, DIMENSION(XX,YY) :: edge_adjust ! mm
 474
 475     !!! Declare Local Variables
 476
 477     REAL :: dhdx,dhdy,alfax,alfay
 478     REAL :: hh53,qqsfc,hh,dt_new,hmax
 479     REAL :: sfx,sfy
 480     REAL :: tmp_adjust
 481
 482     INTEGER :: i,j
 483     REAL IXX8,IYY8
 484     INTEGER  IXX0,JYY0,index, SO8RT_D(XX,YY,3)
 485     REAL :: tmp_gsize, hsum, tmp_gsize_arr(9)
 486
 487     !!! Initialize variables
 488
 489
 490
 491     !!! Begin Routing of Excess Rainfall over the Watershed
 492
 493     DH=0.
 494     DH_tmp=0.
 495     QBDRY_tmp =0.
 496
 497     !!! Loop to route water
 498     do j=2,YY-1
 499         do i=2,XX-1
 500             if (h(I,J).GT.retent_dep(I,J)) then
 501                 IXX0 = SO8RT_D(i,j,1)
 502                 JYY0 = SO8RT_D(i,j,2)
 503                 index = SO8RT_D(i,j,3)
 504                 tmp_gsize = 1.0/gsize(i,j,index)
 505                 sfx = so8RT(i,j,index)-(h(IXX0,JYY0)-h(i,j))*0.001*tmp_gsize
 506                 hmax = h(i,j)*0.001  !Specify max head for mass flux limit...
 507                 if(sfx .lt. 1E-20) then
 508                     tmp_gsize_arr = gsize(i,j,:)
 509                     call GETMAX8DIR(IXX0,JYY0,I,J,H,RETENT_DEP,so8rt,tmp_gsize_arr,sfx,XX,YY)
 510                 end if
 511                 if(IXX0 > 0) then  ! do the rest if the lowest grid can be found.
 512                     if(sfx .lt. 1E-20) then
 513 #ifdef HYDRO_D
 514                         print*, "Message: sfx reset to 1E-20. sfx =",sfx
 515                         print*, "i,j,index,IXX0,JYY0",i,j,index,IXX0,JYY0
 516                         print*, "so8RT(i,j,index), h(IXX0,JYY0), h(i,j), gsize(i,j,index) ", &
 517                             so8RT(i,j,index), h(IXX0,JYY0), h(i,j), gsize(i,j,index)
 518 #endif
 519                         sfx = 1E-20
 520                     end if
 521                     alfax = sqrt(sfx) / dist_rough(i,j)
 522                     hh=(h(i,j)-retent_dep(i,j)) * 0.001
 523                     hh53=hh**(5./3.)
 524
 525                     ! Calculate q-flux...
 526                     qqsfc = alfax*hh53*dt * tmp_gsize
 527
 528                     !Courant check (simple mass limit on overland flow)...
 529                     if (qqsfc.ge.(hmax*dt*tmp_gsize)) qqsfc = hmax*dt*tmp_gsize
 530
 531                     ! Accumulate directional fluxes on routing subgrid...
 532                     if (IXX0.gt.i) then
 533                         q_sfcflx_x(I,J) = q_sfcflx_x(I,J) + qqsfc * &
 534                             (1.0 - 0.5 * (ABS(j-JYY0)))
 535                     else if (IXX0.lt.i) then
 536                         q_sfcflx_x(I,J) = q_sfcflx_x(I,J) - 1.0 * &
 537                             qqsfc * (1.0 - 0.5 * (ABS(j-JYY0)))
 538                     else
 539                         q_sfcflx_x(I,J) = q_sfcflx_x(I,J) + 0.
 540                     end if
 541                     if (JYY0.gt.j) then
 542                         q_sfcflx_y(I,J) = q_sfcflx_y(I,J) + qqsfc * &
 543                             (1.0 - 0.5 * (ABS(i-IXX0)))
 544                     elseif (JYY0.lt.j) then
 545                         q_sfcflx_y(I,J) = q_sfcflx_y(I,J) - 1.0 * &
 546                             qqsfc * (1.0 - 0.5 * (ABS(i-IXX0)))
 547                     else
 548                         q_sfcflx_y(I,J) = q_sfcflx_y(I,J) + 0.
 549                     end if
 550
 551                     !DJG put adjustment in for (h) due to qqsfc
 552
 553                     !yw changed as following:
 554                     tmp_adjust=qqsfc*1000
 555
 556                     if((h(i,j) - tmp_adjust) <0 )  then
 557 #ifdef HYDRO_D
 558                         print*, "Error Warning: surface head is negative:  ",i,j,ixx0,jyy0, &
 559                             h(i,j) - tmp_adjust
 560 #endif
 561                         tmp_adjust = h(i,j)
 562                     end if
 563                     DH(i,j) = DH(i,j)-tmp_adjust
 564                     DH_tmp(ixx0,jyy0) = DH_tmp(ixx0,jyy0) + tmp_adjust
 565                     !yw end change
 566
 567                     !DG Boundary adjustments here
 568                     !DG Constant Flux Condition
 569 #ifdef MPP_LAND
 570                     if( ((ixx0.eq.XX).and.(right_id .lt. 0)) .or. &
 571                             ((ixx0.eq.1) .and.(left_id  .lt. 0)) .or. &
 572                             ((jyy0.eq.1) .and.(down_id  .lt. 0)) .or. &
 573                             ((JYY0.eq.YY).and.(up_id    .lt. 0)) ) then
 574                         !QBDRY_tmp(IXX0,JYY0)=QBDRY_tmp(IXX0,JYY0) - qqsfc*1000.
 575 #else
 576                         if ((ixx0.eq.XX).or.(ixx0.eq.1).or.(jyy0.eq.1)   &
 577                                 .or.(JYY0.eq.YY )) then
 578                             !QBDRY(IXX0,JYY0)=QBDRY(IXX0,JYY0) - qqsfc*1000.
 579 #endif
 580                             QBDRY_tmp(IXX0,JYY0)=QBDRY_tmp(IXX0,JYY0) - qqsfc*1000.
 581                             QBDRYT=QBDRYT - qqsfc
 582                             DH_tmp(IXX0,JYY0)= DH_tmp(IXX0,JYY0)-tmp_adjust
 583
 584                         end if
 585                     end if
 586
 587                     !! End loop to route sfc water
 588                 end if
 589
 590         end do
 591     end do
 592
 593 #ifdef MPP_LAND
 594     ! use double precision to solve the underflow problem.
 595     call MPP_LAND_COM_REAL8(DH_tmp,XX,YY,1)
 596     call MPP_LAND_COM_REAL8(QBDRY_tmp,XX,YY,1)
 597 #endif
 598     QBDRY = QBDRY + QBDRY_tmp
 599     DH = DH+DH_tmp
 600
 601 #ifdef MPP_LAND
 602     call MPP_LAND_COM_REAL8(DH,XX,YY,99)
 603     call MPP_LAND_COM_REAL(QBDRY,XX,YY,99)
 604 #endif
 605
 606     H = H + DH
 607
 608 !!! Scrape the outermost edges
 609 edge_adjust = 0.0
 610 do j=1,YY,YY-1
 611    do i=1,XX
 612 #ifdef MPP_LAND
 613       if( ((i.eq.XX).and.(right_id .lt. 0)) .or. &
 614           ((i.eq.1) .and.(left_id  .lt. 0)) .or. &
 615           ((j.eq.1) .and.(down_id  .lt. 0)) .or. &
 616           ((j.eq.YY).and.(up_id    .lt. 0)) ) then
 617 #else
 618           if ((i.eq.XX).or.(i.eq.1).or.(j.eq.1)   &
 619                .or.(j.eq.YY )) then
 620 #endif
 621               if (h(i,j) .GT. retent_dep(i,j)) then
 622                  edge_adjust(i,j) = h(i,j) - retent_dep(i,j) ! positive mm
 623               end if
 624
 625           end if
 626    end do
 627 end do
 628
 629 do i=1,XX,XX-1
 630    do j=1,YY
 631 #ifdef MPP_LAND
 632       if( ((i.eq.XX).and.(right_id .lt. 0)) .or. &
 633           ((i.eq.1) .and.(left_id  .lt. 0)) .or. &
 634           ((j.eq.1) .and.(down_id  .lt. 0)) .or. &
 635           ((j.eq.YY).and.(up_id    .lt. 0)) ) then
 636 #else
 637           if ((i.eq.XX).or.(i.eq.1).or.(j.eq.1)   &
 638                .or.(j.eq.YY )) then
 639 #endif
 640               if (h(i,j) .GT. retent_dep(i,j)) then
 641                  edge_adjust(i,j) = h(i,j) - retent_dep(i,j) ! positive mm
 642               end if
 643
 644           end if
 645    end do
 646 end do
 647
 648
 649 #ifdef MPP_LAND
 650     call MPP_LAND_COM_REAL(edge_adjust,XX,YY,99)
 651 #endif
 652     QBDRY = QBDRY - edge_adjust ! making this negative term more negative
 653     H = H - edge_adjust ! making this positive term less positive
 654 !!! End outermost edge scrape
 655
 656     return
 657
 658     !DJG ----------------------------------------------------------------------
 659 end subroutine ROUTE_OVERLAND1
 660
 661 !DJG ----------------------------------------------------------------------
 662
 663 !DJG     SUBROUTINE ROUTE_OVERLAND2
 664 !DJG ----------------------------------------------------------------
 665
 666 SUBROUTINE ROUTE_OVERLAND2 (dt,                               &
 667         &          dist,h,qsfc,sox,soy,                                   &
 668         &          retent_dep,dist_rough,XX,YY,QBDRY,QBDRYT,               &
 669         &          q_sfcflx_x,q_sfcflx_y)
 670
 671     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 672     !
 673     !  Subroutine to route excess rainfall over the watershed
 674     !     using a 2d diffusion routing scheme.
 675     !
 676     !  Called from: main.f
 677     !
 678     !      Will try to formulate this to be called from NOAH
 679     !
 680     !  Returns: qsfc=DQOV   which in turn becomes DH in head calc.
 681     !
 682     !  Created:  Adaptded from CASC2D source code
 683     !  NOTE: dh from original code has been replaced by qsfc
 684     !        dhh replaced by qqsfc
 685     !
 686     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 687 #ifdef MPP_LAND
 688     use module_mpp_land, only: left_id,down_id,right_id,&
 689         up_id,mpp_land_com_real,MPP_LAND_UB_COM, &
 690         MPP_LAND_LR_COM,mpp_land_com_integer
 691 #endif
 692
 693     IMPLICIT NONE
 694
 695
 696     !! Declare Passed variables
 697
 698     real :: gsize
 699     INTEGER, INTENT(IN) :: XX,YY
 700     REAL, INTENT(IN) :: dt , dist(XX,YY,9)
 701
 702     !! Declare passed arrays
 703
 704     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: h
 705     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: qsfc
 706     REAL, INTENT(IN), DIMENSION(XX,YY) :: sox
 707     REAL, INTENT(IN), DIMENSION(XX,YY) :: soy
 708     REAL, INTENT(IN), DIMENSION(XX,YY) :: retent_dep
 709     REAL, INTENT(IN), DIMENSION(XX,YY) :: dist_rough
 710     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: QBDRY
 711     REAL, INTENT(INOUT), DIMENSION(XX,YY) :: q_sfcflx_x,q_sfcflx_y
 712     REAL, INTENT(INOUT) :: QBDRYT
 713     REAL  :: DH(XX,YY)
 714
 715     !!! Declare Local Variables
 716
 717     REAL :: dhdx,dhdy,alfax,alfay
 718     REAL :: hh53,qqsfc,hh,dt_new
 719     REAL :: sfx,sfy
 720     REAL :: tmp_adjust
 721
 722     INTEGER :: i,j
 723
 724     !!! Initialize variables
 725
 726
 727
 728
 729     !!! Begin Routing of Excess Rainfall over the Watershed
 730
 731
 732     DH = 0
 733     !!! Loop to route water in x-direction
 734     do j=1,YY
 735         do i=1,XX
 736             ! check for boundary gridpoint?
 737             if (i.eq.XX) GOTO 998
 738             gsize = dist(i,j,3)
 739
 740             ! check for detention storage?
 741             if (h(i,j).lt.retent_dep(i,j).AND.     &
 742                 h(i+1,j).lt.retent_dep(i+1,j)) GOTO 998
 743
 744             dhdx = (h(i+1,j)/1000. - h(i,j)/1000.) / gsize  ! gisze-(m),h-(mm)
 745
 746             sfx = (sox(i,j)-dhdx+1E-30)
 747             if (abs(sfx).lt.1E-20) sfx=1E-20
 748             alfax = ((abs(sfx))**0.5)/dist_rough(i,j)
 749             if (sfx.lt.0.) then
 750                 hh=(h(i+1,j)-retent_dep(i+1,j))/1000.
 751             else
 752                 hh=(h(i,j)-retent_dep(i,j))/1000.
 753             end if
 754
 755             if ((retent_dep(i,j).gt.0.).AND.(hh.le.0.)) GOTO 998
 756             if (hh.lt.0.) then
 757                 GOTO 998
 758             end if
 759
 760             hh53=hh**(5./3.)
 761
 762             ! Calculate q-flux... (units (m))
 763             qqsfc = (sfx/abs(sfx))*alfax*hh53*dt/gsize
 764             q_sfcflx_x(I,J) = q_sfcflx_x(I,J) + qqsfc
 765
 766             !DJG put adjustment in for (h) due to qqsfc
 767
 768             !yw changed as following:
 769             tmp_adjust=qqsfc*1000
 770             if(tmp_adjust .le. 0 ) GOTO 998
 771             if((h(i,j) - tmp_adjust) <0 )  then
 772 #ifdef HYDRO_D
 773                 print*, "WARNING: surface head is negative:  ",i,j
 774 #endif
 775                 tmp_adjust = h(i,j)
 776             end if
 777             if((h(i+1,j) + tmp_adjust) <0) then
 778 #ifdef HYDRO_D
 779                 print*, "WARNING: surface head is negative: ",i+1,j
 780 #endif
 781                 tmp_adjust = -1*h(i+1,j)
 782             end if
 783             Dh(i,j) = Dh(i,j)-tmp_adjust
 784             Dh(i+1,j) = Dh(i+1,j) + tmp_adjust
 785             !yw end change
 786
 787
 788
 789             !DG Boundary adjustments here
 790             !DG Constant Flux Condition
 791 #ifdef MPP_LAND
 792             if ((i.eq.1).AND.(sfx.lt.0).and. &
 793                     (left_id .lt. 0) ) then
 794 #else
 795                 if ((i.eq.1).AND.(sfx.lt.0)) then
 796 #endif
 797                     Dh(i,j) = Dh(i,j) + qqsfc*1000.
 798                     QBDRY(I,J)=QBDRY(I,J) + qqsfc*1000.
 799                     QBDRYT=QBDRYT + qqsfc*1000.
 800 #ifdef MPP_LAND
 801                 else if ( (i.eq.(XX-1)).AND.(sfx.gt.0) &
 802                         .and. (right_id .lt. 0) ) then
 803 #else
 804                 else if ((i.eq.(XX-1)).AND.(sfx.gt.0)) then
 805 #endif
 806                     tmp_adjust = qqsfc*1000.
 807                     if(h(i+1,j).lt.tmp_adjust) tmp_adjust = h(i+1,j)
 808                     Dh(i+1,j) = Dh(i+1,j) - tmp_adjust
 809                     !DJG Re-assign h(i+1) = 0.0 when <0.0 (from rounding/truncation error)
 810                     QBDRY(I+1,J)=QBDRY(I+1,J) - tmp_adjust
 811                     QBDRYT=QBDRYT - tmp_adjust
 812                 end if
 813
 814
 815                 998     continue
 816
 817                 !! End loop to route sfc water in x-direction
 818         end do
 819     end do
 820
 821     H = H + DH
 822 #ifdef MPP_LAND
 823     call MPP_LAND_LR_COM(H,XX,YY,99)
 824     call MPP_LAND_LR_COM(QBDRY,XX,YY,99)
 825 #endif
 826
 827
 828     DH = 0
 829     !!!! Loop to route water in y-direction
 830     do j=1,YY
 831         do i=1,XX
 832
 833             !! check for boundary grid point?
 834             if (j.eq.YY) GOTO 999
 835             gsize = dist(i,j,1)
 836
 837
 838             !! check for detention storage?
 839             if (h(i,j).lt.retent_dep(i,j).AND.     &
 840                 h(i,j+1).lt.retent_dep(i,j+1)) GOTO 999
 841
 842             dhdy = (h(i,j+1)/1000. - h(i,j)/1000.) / gsize
 843
 844             sfy = (soy(i,j)-dhdy+1E-30)
 845             if (abs(sfy).lt.1E-20) sfy=1E-20
 846             alfay = ((abs(sfy))**0.5)/dist_rough(i,j)
 847             if (sfy.lt.0.) then
 848                 hh=(h(i,j+1)-retent_dep(i,j+1))/1000.
 849             else
 850                 hh=(h(i,j)-retent_dep(i,j))/1000.
 851             end if
 852
 853             if ((retent_dep(i,j).gt.0.).AND.(hh.le.0.)) GOTO 999
 854             if (hh.lt.0.) then
 855                 GOTO 999
 856             end if
 857
 858             hh53=hh**(5./3.)
 859
 860             ! Calculate q-flux...
 861             qqsfc = (sfy/abs(sfy))*alfay*hh53*dt / gsize
 862             q_sfcflx_y(I,J) = q_sfcflx_y(I,J) + qqsfc
 863
 864
 865             !DJG put adjustment in for (h) due to qqsfc
 866             !yw   h(i,j) = h(i,j)-qqsfc*1000.
 867             !yw          h(i,j+1) = h(i,j+1) + qqsfc*1000.
 868             !yw changed as following:
 869             tmp_adjust=qqsfc*1000
 870             if(tmp_adjust .le. 0 ) GOTO 999
 871
 872             if((h(i,j) - tmp_adjust) <0 )  then
 873 #ifdef HYDRO_D
 874                 print *, "WARNING: surface head is negative:  ",i,j
 875 #endif
 876                 tmp_adjust = h(i,j)
 877             end if
 878             if((h(i,j+1) + tmp_adjust) <0) then
 879 #ifdef HYDRO_D
 880                 print *, "WARNING: surface head is negative: ",i,j+1
 881 #endif
 882                 tmp_adjust = -1*h(i,j+1)
 883             end if
 884             Dh(i,j) = Dh(i,j)-tmp_adjust
 885             Dh(i,j+1) = Dh(i,j+1) + tmp_adjust
 886             !yw end change
 887
 888             !qsfc(i,j) = qsfc(i,j)-qqsfc
 889             !qsfc(i,j+1) = qsfc(i,j+1) + qqsfc
 890             !!DG Boundary adjustments here
 891             !!DG Constant Flux Condition
 892 #ifdef MPP_LAND
 893             if ((j.eq.1).AND.(sfy.lt.0)   &
 894                     .and. (down_id .lt. 0) ) then
 895 #else
 896                 if ((j.eq.1).AND.(sfy.lt.0)) then
 897 #endif
 898                     Dh(i,j) = Dh(i,j) + qqsfc*1000.
 899                     QBDRY(I,J)=QBDRY(I,J) + qqsfc*1000.
 900                     QBDRYT=QBDRYT + qqsfc*1000.
 901 #ifdef MPP_LAND
 902                 else if ((j.eq.(YY-1)).AND.(sfy.gt.0) &
 903                         .and. (up_id .lt. 0) ) then
 904 #else
 905                 else if ((j.eq.(YY-1)).AND.(sfy.gt.0)) then
 906 #endif
 907                     tmp_adjust = qqsfc*1000.
 908                     if(h(i,j+1).lt.tmp_adjust) tmp_adjust = h(i,j+1)
 909                     Dh(i,j+1) = Dh(i,j+1) - tmp_adjust
 910                     !DJG Re-assign h(j+1) = 0.0 when <0.0 (from rounding/truncation error)
 911                     QBDRY(I,J+1)=QBDRY(I,J+1) - tmp_adjust
 912                     QBDRYT=QBDRYT - tmp_adjust
 913                 end if
 914
 915                 999     continue
 916
 917                 !!!! End loop to route sfc water in y-direction
 918         end do
 919     end do
 920
 921     H = H +DH
 922 #ifdef MPP_LAND
 923     call MPP_LAND_UB_COM(H,XX,YY,99)
 924     call MPP_LAND_UB_COM(QBDRY,XX,YY,99)
 925 #endif
 926     return
 927
 928     !DJG ----------------------------------------------------------------------
 929 end subroutine ROUTE_OVERLAND2