hydro/Routing/module_GW_baseflow.F

   1 !  Program Name:
   2 !  Author(s)/Contact(s):
   3 !  Abstract:
   4 !  History Log:
   5 !
   6 !  Usage:
   7 !  Parameters: <Specify typical arguments passed>
   8 !  Input Files:
   9 !        <list file names and briefly describe the data they include>
  10 !  Output Files:
  11 !        <list file names and briefly describe the information they include>
  12 !
  13 !  Condition codes:
  14 !        <list exit condition or error codes returned >
  15 !        If appropriate, descriptive troubleshooting instructions or
  16 !        likely causes for failures could be mentioned here with the
  17 !        appropriate error code
  18 !
  19 !  User controllable options: <if applicable>
  20
  21 module module_GW_baseflow
  22
  23 !   use overland_data
  24 #ifdef MPP_LAND
  25    use module_mpp_land
  26    use MODULE_mpp_GWBUCKET, only: gw_sum_real, gw_write_io_real
  27    use MODULE_mpp_ReachLS, only : updatelinkv
  28 #endif
  29    use iso_fortran_env, only: int64
  30    implicit none
  31
  32 !#include "rt_include.inc"
  33 !yw #include "namelist.inc"
  34 contains
  35
  36 !------------------------------------------------------------------------------
  37 !DJG   Simple GW Bucket Model
  38 !      for NHDPLUS mapping
  39 !------------------------------------------------------------------------------
  40
  41    subroutine simp_gw_buck_nhd(           &
  42         ix,            jx,                &
  43         ixrt,          jxrt,              &
  44         numbasns,      AGGFACTRT,         &
  45         DT,            INFXSWGT,          &
  46         runoff1x_in,   runoff2x_in,       &
  47         cellArea,      area_lsm,          &
  48         c,             ex,                &
  49         loss_fraction,                    &
  50         z_mx,          z_gwsubbas_tmp,    &
  51         qout_gwsubbas, qin_gwsubbas,      &
  52         qloss_gwsubbas,                   &
  53         GWBASESWCRT,   OVRTSWCRT,         &
  54         LNLINKSL,                         &
  55         basns_area,                       &
  56         nhdBuckMask,   bucket_loss,       &
  57         channelBucket_only )
  58
  59    use module_UDMAP, only: LNUMRSL, LUDRSL
  60
  61    implicit none
  62
  63 !!!Declarations...
  64    integer, intent(in)                               :: ix,jx,ixrt,jxrt
  65    integer, intent(in)                               :: numbasns, lnlinksl
  66    real, intent(in), dimension(ix,jx)                :: runoff2x_in
  67    real, dimension(ixrt,jxrt)                            :: runoff2x , runoff1x
  68    real, intent(in), dimension(ix,jx)                :: runoff1x_in, area_lsm
  69    real, intent(in)                                  :: cellArea(ixrt,jxrt),DT
  70    real, intent(in),dimension(numbasns)              :: C,ex,loss_fraction
  71    real, intent(inout),dimension(numbasns)              :: z_mx
  72    real, intent(out),dimension(numbasns)             :: qout_gwsubbas
  73    !! intent inout for channelBucket_only .eq. 1
  74    real, intent(inout),dimension(numbasns)           :: qin_gwsubbas
  75    real, intent(out),dimension(numbasns)             :: qloss_gwsubbas
  76    real*8                                            :: z_gwsubbas(numbasns)
  77    real                                              :: qout_max, qout_spill, z_gw_spill
  78    real, intent(inout),dimension(:)                  :: z_gwsubbas_tmp
  79    real, intent(in),dimension(ixrt,jxrt)             :: INFXSWGT
  80    integer, intent(in)                               :: GWBASESWCRT
  81    integer, intent(in)                               :: OVRTSWCRT
  82    real, intent(in), dimension(numbasns)             :: basns_area
  83    integer, intent(in)                               :: channelBucket_only
  84    integer, intent(in)                               :: bucket_loss
  85
  86    real, dimension(numbasns)                         :: net_perc
  87    integer, dimension(numbasns)                      :: nhdBuckMask
  88
  89    integer                                           :: i,j,bas, k, m, ii,jj
  90
  91    integer :: AGGFACYRT, AGGFACTRT, AGGFACXRT, IXXRT, JYYRT
  92    real*8,  dimension(LNLINKSL) :: LQLateral
  93
  94    real, parameter :: one         = 1.000000000000000
  95    real, parameter :: m_to_mm     = 1000.000000000000
  96    real, parameter :: mm_to_m     = one/m_to_mm
  97    real, parameter :: sqm_to_sqkm = one/1000000.0000000000
  98
  99
 100 !!!Initialize variables...
 101    net_perc = 0.
 102    qout_gwsubbas = 0.
 103    if (bucket_loss .eq. 1) then
 104       qloss_gwsubbas = 0.
 105    endif
 106    z_gwsubbas(1:numbasns) = z_gwsubbas_tmp(1:numbasns)
 107
 108    if(channelBucket_only .eq. 0) then
 109
 110       !! Initialize if not passed in
 111       qin_gwsubbas = 0.
 112
 113 !Assign local value of runoff2 (drainage) for flux caluclation to buckets...
 114
 115         do J=1,JX
 116         do I=1,IX
 117              do AGGFACYRT=AGGFACTRT-1,0,-1
 118              do AGGFACXRT=AGGFACTRT-1,0,-1
 119                IXXRT=I*AGGFACTRT-AGGFACXRT
 120                JYYRT=J*AGGFACTRT-AGGFACYRT
 121 #ifdef MPP_LAND
 122        if(left_id.ge.0) IXXRT=IXXRT+1
 123        if(down_id.ge.0) JYYRT=JYYRT+1
 124 !              if(AGGFACTRT .eq. 1) then
 125 !                  IXXRT=I
 126 !                  JYYRT=J
 127 !             endif
 128 #endif
 129 !DJG Implement subgrid weighting routine...
 130                if( (runoff1x_in(i,j) .lt. 0) .or. (runoff1x_in(i,j) .gt. 1000) ) then
 131                     runoff1x(IXXRT,JYYRT) = 0
 132                else
 133                     runoff1x(IXXRT,JYYRT)=runoff1x_in(i,j)*area_lsm(I,J)     &
 134                         *INFXSWGT(IXXRT,JYYRT)/cellArea(IXXRT,JYYRT)
 135                endif
 136
 137                if( (runoff2x_in(i,j) .lt. 0) .or. (runoff2x_in(i,j) .gt. 1000) ) then
 138                     runoff2x(IXXRT,JYYRT) = 0
 139                else
 140                   runoff2x(IXXRT,JYYRT)=runoff2x_in(i,j)*area_lsm(I,J)     &
 141                       *INFXSWGT(IXXRT,JYYRT)/cellArea(IXXRT,JYYRT)
 142                endif
 143              enddo
 144              enddo
 145         enddo
 146         enddo
 147
 148
 149        LQLateral = 0
 150        do k = 1, LNUMRSL
 151               ! get from land grid runoff
 152                do m = 1, LUDRSL(k)%ncell
 153                    ii =  LUDRSL(k)%cell_i(m)
 154                    jj =  LUDRSL(k)%cell_j(m)
 155                    if(ii .gt. 0 .and. jj .gt. 0) then
 156                            LQLateral(k) = LQLateral(k)+runoff2x(ii,jj)*LUDRSL(k)%cellWeight(m) * mm_to_m &
 157                                *cellArea(ii,jj)
 158                    endif
 159                end do
 160        end do
 161
 162
 163 #ifdef MPP_LAND
 164        call updateLinkV(LQLateral, net_perc)      ! m^3
 165
 166 #else
 167        net_perc = LQLateral        ! m^3
 168 #endif
 169
 170     endif !! if channelBucket_only .eq. 0 else
 171
 172 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 173 !!!Loop through GW basins to adjust for inflow/outflow
 174 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 175
 176
 177
 178    DO bas=1,numbasns     ! Loop for GW bucket calcs...
 179       if(nhdBuckMask(bas) .eq. 1) then     ! if the basn is masked
 180
 181          if(channelBucket_only .eq. 0) then
 182             !! If not using channelBucket_only, save qin_gwsubbas
 183             qin_gwsubbas(bas) = net_perc(bas)             !units (m^3)
 184          else
 185             !! If using channelBucket_only, get net_perc from the passed qin_gwsubbas
 186             net_perc(bas)     = qin_gwsubbas(bas)         !units (m^3)
 187          end if
 188
 189          ! !Adjust level of GW depth...(conceptual GW bucket units (mm))
 190          z_gwsubbas(bas) = z_gwsubbas(bas) + net_perc(bas) / basns_area(bas)   ! m
 191
 192          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 193          !Calculate baseflow as a function of GW bucket depth...
 194          ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 195          if( (GWBASESWCRT.eq.1) .or. (GWBASESWCRT.eq.4) ) then  !active bucket model
 196
 197             !DJG...Estimation of bucket 'overflow' (qout_spill) if/when bucket gets filled...
 198             qout_spill = 0.
 199             z_gw_spill = 0.
 200
 201             !!DJG...convert z_mx to millimeters...for v2 and later...
 202             !yw  added by Wei Yu...If block is to accomodate old parameter file...
 203             !                    if(z_mx(bas) .gt. 5) then
 204             !                         z_mx(bas) = z_mx(bas) * mm_to_m    ! change from mm to meters
 205             !                    endif
 206
 207                if (z_gwsubbas(bas).gt.z_mx(bas) * mm_to_m) then  !If/then for bucket overflow case...
 208
 209                     z_gw_spill = z_gwsubbas(bas) - z_mx(bas) * mm_to_m    ! meters
 210                     z_gwsubbas(bas) = z_mx(bas) * mm_to_m    ! meters
 211
 212                else
 213                       z_gw_spill = 0.
 214                end if   ! End if for bucket overflow case...
 215
 216                qout_spill = z_gw_spill*(basns_area(bas))/DT  !amount spilled from bucket overflow...units (m^3/s)
 217
 218 !DJG...Maximum estimation of bucket outlfow that is limited by total quantity in bucket...
 219                qout_max = z_gwsubbas(bas)*(basns_area(bas))/DT   ! (m^3/s)   ! Estimate max bucket disharge limit to total volume in bucket...(m^3/s)
 220
 221
 222 ! Assume exponential relation between z/zmax and Q...
 223 !DJG force asymptote to zero to prevent 'overdraft'...
 224                if(GWBASESWCRT.eq.1) then ! Exponential model
 225                    qout_gwsubbas(bas) = C(bas)*(exp(ex(bas)*z_gwsubbas(bas)/(z_mx(bas) * mm_to_m))-1) ! q_out (m^3/s)
 226
 227                elseif(GWBASESWCRT.eq.4) then ! Updated exponential model normalized by area
 228                    qout_gwsubbas(bas) = C(bas)*(basns_area(bas)*sqm_to_sqkm)*(exp(ex(bas)*z_gwsubbas(bas)/(z_mx(bas)*mm_to_m))-1) ! q_out (m^3/s)
 229
 230                end if
 231
 232 !DJG...Calculation of max bucket outlfow that is limited by total quantity in bucket...
 233                qout_gwsubbas(bas) = MIN(qout_max,qout_gwsubbas(bas))   ! Limit bucket discharge to max. bucket limit   (m^3/s)
 234
 235                if (bucket_loss .eq. 1) then
 236                     qloss_gwsubbas(bas) = qout_gwsubbas(bas)*loss_fraction(bas)
 237                     qout_gwsubbas(bas) = qout_gwsubbas(bas)-qloss_gwsubbas(bas)
 238                endif
 239
 240           elseif (GWBASESWCRT.eq.2) then  !Pass through/steady-state bucket
 241
 242 ! Assuming a steady-state (inflow=outflow) model...
 243 !DJG convert input and output units to cms...       qout_gwsubbas(bas) = qin_gwsubbas(bas)  !steady-state model...(m^3)
 244                qout_gwsubbas(bas) = qin_gwsubbas(bas)/DT  !steady-state model...(m^3/s)
 245
 246           end if    ! End if for bucket model discharge type....
 247
 248
 249 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 250 !Adjust level of GW depth in bucket...
 251 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 252
 253           if (bucket_loss .eq. 1) then
 254              z_gwsubbas(bas) = z_gwsubbas(bas) - (qout_gwsubbas(bas)+qloss_gwsubbas(bas))*DT/( basns_area(bas) )   ! units (meters)
 255           else
 256              z_gwsubbas(bas) = z_gwsubbas(bas) - (qout_gwsubbas(bas))*DT/( basns_area(bas) )   ! units (meters)
 257           endif
 258
 259 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 260 !Combine calculated bucket discharge and amount spilled from bucket...
 261 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 262
 263           qout_gwsubbas(bas) = qout_gwsubbas(bas) + qout_spill   ! units (m^3/s)
 264       else
 265           qout_gwsubbas(bas) = 0.0
 266       endif   ! the basns is masked
 267
 268
 269    END DO                 ! End loop for GW bucket calcs...
 270
 271    z_gwsubbas_tmp(1:numbasns) = z_gwsubbas(1:numbasns)     ! units (meters)
 272
 273    return
 274
 275 !------------------------------------------------------------------------------
 276    End subroutine simp_gw_buck_nhd
 277 !------------------------------------------------------------------------------
 278 !------------------------------------------------------------------------------
 279 !DJG   Simple GW Bucket Model
 280 !------------------------------------------------------------------------------
 281
 282    subroutine simp_gw_buck(ix,jx,ixrt,jxrt,numbasns,gnumbasns,basns_area,basnsInd,gw_strm_msk_lind,&
 283                             gwsubbasmsk, runoff1x_in, runoff2x_in, z_gwsubbas_tmp, qin_gwsubbas,&
 284                             qout_gwsubbas,qinflowbase,gw_strm_msk,gwbas_pix_ct,dist,DT,&
 285                             C,ex,z_mx,GWBASESWCRT,OVRTSWCRT)
 286    implicit none
 287
 288 !!!Declarations...
 289    integer, intent(in)                               :: ix,jx,ixrt,jxrt
 290    integer, intent(in)                               :: numbasns, gnumbasns
 291    integer, intent(in), dimension(ix,jx)             :: gwsubbasmsk
 292    real, intent(in), dimension(ix,jx)                :: runoff2x_in
 293    real, dimension(ix,jx)                            :: runoff2x
 294    real, intent(in), dimension(ix,jx)                :: runoff1x_in
 295    real, dimension(ix,jx)                            :: runoff1x
 296    real, intent(in)                                  :: basns_area(numbasns),dist(ixrt,jxrt,9),DT
 297    integer(kind=int64), intent(in)                   :: basnsInd(numbasns)
 298    real, intent(in),dimension(numbasns)              :: C,ex,z_mx
 299    real, intent(out),dimension(numbasns)             :: qout_gwsubbas
 300    real, intent(out),dimension(numbasns)             :: qin_gwsubbas
 301    real, dimension(numbasns)                         :: qin_gwsubbas_surf
 302    real*8                                            :: z_gwsubbas(numbasns)
 303    real                                              :: qout_max, qout_spill, z_gw_spill
 304    real, intent(inout),dimension(numbasns)           :: z_gwsubbas_tmp
 305    real, intent(out),dimension(ixrt,jxrt)            :: qinflowbase
 306    integer, intent(in),dimension(ixrt,jxrt)          :: gw_strm_msk, gw_strm_msk_lind
 307    integer, intent(in)                               :: GWBASESWCRT
 308    integer, intent(in)                               :: OVRTSWCRT
 309
 310
 311    real*8, dimension(numbasns)                      :: sum_perc8, ct_bas8, sum_perc8_surf
 312    real, dimension(numbasns)                        :: sum_perc, sum_perc_surf
 313    real, dimension(numbasns)                        :: net_perc, net_perc_surf
 314
 315    real, dimension(numbasns)                        :: ct_bas
 316    real, dimension(numbasns)                        :: gwbas_pix_ct
 317    integer                                          :: i,j,bas, k
 318    character(len=19)                                :: header
 319    character(len=1)                                 :: jnk
 320
 321
 322 !!!Initialize variables...
 323    ct_bas8 = 0
 324    sum_perc8 = 0.
 325    sum_perc8_surf = 0.
 326    net_perc = 0.
 327    net_perc_surf = 0.
 328    qout_gwsubbas = 0.
 329    qin_gwsubbas = 0.
 330    qin_gwsubbas_surf = 0.
 331    z_gwsubbas = z_gwsubbas_tmp
 332
 333 !Assign local value of runoff2 (drainage) for flux caluclation to buckets...
 334    runoff2x = runoff2x_in
 335    runoff1x = runoff1x_in
 336
 337
 338
 339
 340 !!!Calculate aggregated percolation from deep runoff into GW basins...
 341    do i=1,ix
 342      do j=1,jx
 343
 344 !!DJG 4/15/2015...reset runoff2x, runoff1x, values to 0 where extreme values exist...(<0 or
 345 !> 1000)
 346        if((runoff2x(i,j).lt.0.).OR.(runoff2x(i,j).gt.1000.)) then
 347          runoff2x(i,j)=0.
 348        end if
 349        if((runoff1x(i,j).lt.0.).OR.(runoff1x(i,j).gt.1000.)) then
 350          runoff1x(i,j)=0.
 351        end if
 352
 353        do bas=1,numbasns
 354          if(gwsubbasmsk(i,j).eq.basnsInd(bas) ) then
 355            !ADCHANGE: Separate surface input from subsurface
 356            if(OVRTSWCRT.eq.0) then
 357              sum_perc8_surf(bas) = sum_perc8_surf(bas)+runoff1x(i,j)
 358            end if
 359            sum_perc8(bas) = sum_perc8(bas)+runoff2x(i,j)  !Add only drainage to bucket...runoff2x in (mm)
 360            ct_bas8(bas) = ct_bas8(bas) + 1
 361          end if
 362        end do
 363      end do
 364    end do
 365
 366 #ifdef MPP_LAND
 367     call gw_sum_real(sum_perc8,numbasns,gnumbasns,basnsInd)
 368     call gw_sum_real(sum_perc8_surf,numbasns,gnumbasns,basnsInd)
 369     call gw_sum_real(ct_bas8,numbasns,gnumbasns,basnsInd)
 370 #endif
 371    sum_perc = sum_perc8
 372    sum_perc_surf = sum_perc8_surf
 373    ct_bas = ct_bas8
 374
 375
 376
 377
 378 !!!Loop through GW basins to adjust for inflow/outflow
 379
 380    DO bas=1,numbasns     ! Loop for GW bucket calcs...
 381 ! #ifdef MPP_LAND
 382 !      if(ct_bas(bas) .gt. 0) then
 383 ! #endif
 384
 385      net_perc(bas) = sum_perc(bas) / ct_bas(bas)   !units (mm)
 386 !DJG...old change to cms     qin_gwsubbas(bas) = net_perc(bas)/1000. * ct_bas(bas) * basns_area(bas) !units (m^3)
 387      qin_gwsubbas(bas) = net_perc(bas)/1000.* &
 388                          ct_bas(bas)*basns_area(bas)/DT    !units (m^3/s)
 389
 390 !ADCHANGE: Separate tracking for surface runoff term to push straight to pass-through
 391      net_perc_surf(bas) = sum_perc_surf(bas) / ct_bas(bas)   !units (mm)
 392      qin_gwsubbas_surf(bas) = net_perc_surf(bas)/1000.* &
 393                          ct_bas(bas)*basns_area(bas)/DT    !units (m^3/s)
 394
 395 !Adjust level of GW depth...(conceptual GW bucket units (mm))
 396 !DJG...old change to cms inflow...     z_gwsubbas(bas) = z_gwsubbas(bas) + net_perc(bas) / 1000.0   ! (m)
 397
 398 !DJG...debug    write (6,*) "DJG...before",C(bas),ex(bas),z_gwsubbas(bas),z_mx(bas),z_gwsubbas(bas)/z_mx(bas)
 399
 400      z_gwsubbas(bas) = z_gwsubbas(bas) + qin_gwsubbas(bas)*DT/( &
 401                        ct_bas(bas)*basns_area(bas))*1000.   !  units (mm)
 402
 403
 404
 405
 406
 407 !Calculate baseflow as a function of GW bucket depth...
 408
 409      if(GWBASESWCRT.eq.1) then  !active exponential bucket... if/then for bucket model discharge type...
 410
 411 !DJG...Estimation of bucket 'overflow' (qout_spill) if/when bucket gets filled...
 412      qout_spill = 0.
 413      z_gw_spill = 0.
 414      if (z_gwsubbas(bas).gt.z_mx(bas)) then  !If/then for bucket overflow case...
 415        z_gw_spill = z_gwsubbas(bas) - z_mx(bas)
 416        z_gwsubbas(bas) = z_mx(bas)
 417 #ifdef HYDRO_D
 418        write (6,*) "Bucket spilling...", bas, z_gwsubbas(bas), z_mx(bas), z_gw_spill
 419 #endif
 420      else
 421        z_gw_spill = 0.
 422      end if   ! End if for bucket overflow case...
 423
 424      qout_spill = z_gw_spill/1000.*(ct_bas(bas)*basns_area(bas))/DT  !amount spilled from bucket overflow...units (cms)
 425
 426
 427 !DJG...Maximum estimation of bucket outlfow that is limited by total quantity in bucket...
 428      qout_max = z_gwsubbas(bas)/1000.*(ct_bas(bas)*basns_area(bas))/DT   ! Estimate max bucket disharge limit to total volume in bucket...(m^3/s)
 429
 430
 431 ! Assume exponential relation between z/zmax and Q...
 432 !DJG...old...creates non-asymptotic flow...   qout_gwsubbas(bas) = C(bas)*EXP(ex(bas)*z_gwsubbas(bas)/z_mx(bas)) !Exp.model. q_out (m^3/s)
 433 !DJG force asymptote to zero to prevent 'overdraft'...
 434 !DJG debug hardwire test...       qout_gwsubbas(bas) = 1*(EXP(7.0*10./100.)-1) !Exp.model. q_out (m^3/s)
 435      qout_gwsubbas(bas) = C(bas)*(EXP(ex(bas)*z_gwsubbas(bas)/z_mx(bas))-1) !Exp.model. q_out (m^3/s)
 436
 437 !DJG...Calculation of max bucket outlfow that is limited by total quantity in bucket...
 438      qout_gwsubbas(bas) = MIN(qout_max,qout_gwsubbas(bas))   ! Limit bucket discharge to max. bucket limit
 439
 440 !DJG...debug...     write (6,*) "DJG-exp bucket...during",C(bas),ex(bas),z_gwsubbas(bas),qin_gwsubbas(bas),z_mx(bas),z_gwsubbas(bas)/z_mx(bas), qout_gwsubbas(bas), qout_max, qout_spill
 441
 442
 443
 444      elseif (GWBASESWCRT.eq.2) then  !Pass through/steady-state bucket
 445
 446 ! Assuming a steady-state (inflow=outflow) model...
 447 !DJG convert input and output units to cms...       qout_gwsubbas(bas) = qin_gwsubbas(bas)  !steady-state model...(m^3)
 448        qout_gwsubbas(bas) = qin_gwsubbas(bas)  !steady-state model...(m^3/s)
 449
 450 !DJG...debug       write (6,*) "DJG-pass through...during",C(bas),ex(bas),qin_gwsubbas(bas), z_gwsubbas(bas),z_mx(bas),z_gwsubbas(bas)/z_mx(bas), qout_gwsubbas(bas), qout_max
 451
 452      end if    ! End if for bucket model discharge type....
 453
 454
 455
 456
 457 !Adjust level of GW depth...
 458 !DJG bug adjust output to be mm and correct area bug...       z_gwsubbas(bas) = z_gwsubbas(bas) - qout_gwsubbas(bas)*DT &
 459 !DJG bug adjust output to be mm and correct area bug...                       / (ct_bas(bas)*basns_area(bas))   !units(m)
 460
 461      z_gwsubbas(bas) = z_gwsubbas(bas) - qout_gwsubbas(bas)*DT/( &
 462                        ct_bas(bas)*basns_area(bas))*1000.   ! units (mm)
 463
 464 !DJG...Combine calculated bucket discharge and amount spilled from bucket...
 465 !ADCHANGE: Add in surface runoff as direct pass-through
 466      qout_gwsubbas(bas) = qout_gwsubbas(bas) + qout_spill + qin_gwsubbas_surf(bas)  ! units (cms)
 467
 468
 469 !DJG...debug     write (6,*) "DJG...after",C(bas),ex(bas),z_gwsubbas(bas),z_mx(bas),z_gwsubbas(bas)/z_mx(bas), qout_gwsubbas(bas), qout_spill
 470 !DJG...debug     write (6,*) "DJG...after...calc",bas,ct_bas(bas),ct_bas(bas)*basns_area(bas),basns_area(bas),DT
 471
 472
 473
 474
 475 ! #ifdef MPP_LAND
 476 !      endif
 477 ! #endif
 478    END DO                 ! End loop for GW bucket calcs...
 479
 480    z_gwsubbas_tmp = z_gwsubbas
 481
 482
 483 !!!Distribute basin integrated baseflow to stream pixels as stream 'inflow'...
 484
 485       qinflowbase = 0.
 486
 487
 488       do i=1,ixrt
 489         do j=1,jxrt
 490 !!!    -simple uniform disaggregation (8.31.06)
 491            if (gw_strm_msk_lind(i,j).gt.0) then
 492
 493              qinflowbase(i,j) = qout_gwsubbas(gw_strm_msk_lind(i,j))*1000.*DT/ &
 494                 gwbas_pix_ct(gw_strm_msk_lind(i,j))/dist(i,j,9)     ! units (mm) that gets passed into chan routing as stream inflow
 495
 496            end if
 497         end do
 498       end do
 499
 500
 501 !!!    - weighted redistribution...(need to pass accum weights (slope) in...)
 502 !        NOT FINISHED just BASIC framework...
 503 !         do bas=1,numbasns
 504 !           do k=1,gwbas_pix_ct(bas)
 505 !             qinflowbase(i,j) = k*slope
 506 !           end do
 507 !         end do
 508
 509         z_gwsubbas = z_gwsubbas_tmp
 510
 511    return
 512
 513 !------------------------------------------------------------------------------
 514    End subroutine simp_gw_buck
 515 !------------------------------------------------------------------------------
 516
 517
 518
 519
 520 #ifdef MPP_LAND
 521    subroutine pix_ct_1(in_gw_strm_msk,ixrt,jxrt,gwbas_pix_ct,numbasns,gnumbasns,basnsInd)
 522       USE module_mpp_land
 523       implicit none
 524       integer ::    i,j,ixrt,jxrt,numbasns, bas, gnumbasns, k
 525       integer,dimension(ixrt,jxrt) :: in_gw_strm_msk
 526       integer,dimension(global_rt_nx,global_rt_ny) :: gw_strm_msk
 527       real,dimension(numbasns) :: gwbas_pix_ct
 528       real,dimension(gnumbasns) :: tmp_gwbas_pix_ct
 529       integer(kind=int64), intent(in), dimension(:) :: basnsInd
 530
 531       gw_strm_msk = 0
 532
 533
 534       call write_IO_rt_int(in_gw_strm_msk, gw_strm_msk)
 535
 536       call mpp_land_sync()
 537
 538       if(my_id .eq. IO_id) then
 539 !        tmp_gwbas_pix_ct = 0.0
 540 !         do bas = 1,gnumbasns
 541 !         do i=1,global_rt_nx
 542 !           do j=1,global_rt_ny
 543 !             if(gw_strm_msk(i,j) .eq. bas) then
 544 !                tmp_gwbas_pix_ct(bas) = tmp_gwbas_pix_ct(bas) + 1.0
 545 !             endif
 546 !           end do
 547 !         end do
 548 !         end do
 549
 550             tmp_gwbas_pix_ct = 0.0
 551             do i=1,global_rt_nx
 552               do j=1,global_rt_ny
 553                 if(gw_strm_msk(i,j) .gt. 0) then
 554                    bas = gw_strm_msk(i,j)
 555                    tmp_gwbas_pix_ct(bas) = tmp_gwbas_pix_ct(bas) + 1.0
 556                endif
 557               end do
 558             end do
 559       end if
 560
 561       call mpp_land_sync()
 562
 563       if(gnumbasns .gt. 0) then
 564          call mpp_land_bcast_real(gnumbasns,tmp_gwbas_pix_ct)
 565       endif
 566       do k = 1, numbasns
 567          bas = basnsInd(k)
 568          gwbas_pix_ct(k) = tmp_gwbas_pix_ct(bas)
 569       end do
 570
 571
 572       return
 573    end subroutine pix_ct_1
 574 #endif
 575
 576
 577
 578
 579
 580 end module module_GW_baseflow