hydro/Routing/module_RT.F90

   1 !  Program Name:
   2 !  Author(s)/Contact(s):
   3 !  Abstract:
   4 !  History Log:
   5 !
   6 !  Usage:
   7 !  Parameters: <Specify typical arguments passed>
   8 !  Input Files:
   9 !        <list file names and briefly describe the data they include>
  10 !  Output Files:
  11 !        <list file names and briefly describe the information they include>
  12 !
  13 !  Condition codes:
  14 !        <list exit condition or error codes returned >
  15 !        If appropriate, descriptive troubleshooting instructions or
  16 !        likely causes for failures could be mentioned here with the
  17 !        appropriate error code
  18 !
  19 !  User controllable options: <if applicable>
  20
  21 MODULE module_Routing
  22 #ifdef MPP_LAND
  23    use module_gw_baseflow, only: pix_ct_1
  24    use module_HYDRO_io, only: mpp_read_routedim, read_routing_seq, mpp_read_chrouting_new, &
  25                               mpp_read_simp_gw, read_routelink, get_nlinksl
  26    use MODULE_mpp_ReachLS, only: ReachLS_ini, getlocalindx,  getToInd
  27    USE module_mpp_land, only : left_id, up_id, right_id, down_id, mpp_land_com_integer, &
  28                                mpp_land_bcast_int, mpp_land_bcast_int1, &
  29                                updateLake_seq
  30    use module_mpp_GWBUCKET, only : collectSizeInd
  31 #else
  32    !yw use module_HYDRO_io, only: read_routedim, read_routing_old, read_chrouting,read_simp_gw
  33    use module_HYDRO_io, only: read_routedim, read_routing_seq, read_chrouting1,read_simp_gw, get_nlinksl
  34 #endif
  35    use module_HYDRO_io, only: readgw2d, simp_gw_ind,read_GWBUCKPARM, get_gw_strm_msk_lind, readBucket_nhd, read_NSIMLAKES
  36    use module_HYDRO_utils
  37
  38    use module_UDMAP, only: LNUMRSL, LUDRSL, UDMP_ini
  39    use module_levelpool, only: levelpool
  40    use module_persistence_levelpool_hybrid, only: persistence_levelpool_hybrid
  41    use module_rfc_forecasts, only: rfc_forecasts
  42    use module_hydro_stop, only: HYDRO_stop
  43    use hashtable
  44    use iso_fortran_env, only: int64
  45 #ifdef OUTPUT_CHAN_CONN
  46 #ifdef MPP_LAND
  47    use mpi
  48 #endif
  49 #endif
  50
  51    IMPLICIT NONE
  52
  53
  54    integer, parameter :: r8 = selected_real_kind(8)
  55    real*8,  parameter :: zeroDbl=0.0000000000000000000_r8
  56    integer, parameter :: r4 = selected_real_kind(4)
  57    real  ,  parameter :: zeroFlt=0.0000000000000000000_r4
  58
  59 CONTAINS
  60
  61    subroutine rt_allocate(did,ix,jx,ixrt,jxrt,nsoil,CHANRTSWCRT)
  62       use module_RT_data, only: rt_domain
  63       use config_base, only: nlst
  64
  65       implicit none
  66       integer ixrt,jxrt, ix,jx,nsoil,NLINKS, CHANRTSWCRT, NLAKES, NLINKSL
  67       integer istatus, did, nsizes
  68
  69       if(rt_domain(did)%allo_status .eq. 1) return
  70       rt_domain(did)%allo_status = 1
  71
  72       rt_domain(did)%ix = ix
  73       rt_domain(did)%jx = jx
  74       rt_domain(did)%ixrt = ixrt
  75       rt_domain(did)%jxrt = jxrt
  76 !     ixrt = rt_domain(did)%ixrt
  77 !     jxrt = rt_domain(did)%jxrt
  78
  79       ! allocate overland data structures
  80       call rt_domain(did)%overland%init(ix, jx, ixrt, jxrt)
  81
  82       ! allocate subsurface data structures
  83       call rt_domain(did)%subsurface%init(ixrt, jxrt, nsoil, rt_domain(did)%overland)
  84
  85       ! allocate susurface static data structure
  86       call rt_domain(did)%subsurface_static%init(ixrt, jxrt, nsoil, nlst(did)%DT, nlst(did)%rt_option)
  87
  88       ! allocate subsurface input structure
  89       call rt_domain(did)%subsurface_output%init(ixrt, jxrt, rt_domain(did)%overland%control%infiltration_excess)
  90
  91       !allocate subsurface output structure
  92       call rt_domain(did)%subsurface_input%init(ixrt, jxrt, rt_domain(did)%overland%control%infiltration_excess)
  93
  94 !     if( nlst_rt(did)%channel_option .eq. 1  .or. nlst_rt(did)%channel_option .eq. 2 ) then
  95 !         rt_domain(did)%NLINKS = rt_domain(did)%NLINKSL
  96 !     endif
  97       if(nlst(did)%UDMP_OPT .eq. 1) then
  98           if(rt_domain(did)%NLINKS .lt. rt_domain(did)%NLINKSL) then
  99               rt_domain(did)%NLINKS = rt_domain(did)%NLINKSL
 100           endif
 101       endif
 102
 103
 104       NLINKS = rt_domain(did)%NLINKS
 105       NLAKES = rt_domain(did)%NLAKES
 106       NLINKSL = rt_domain(did)%NLINKSL
 107
 108       if(NLINKSL .gt. NLINKS) then
 109          nsizes = nlinksl
 110       else
 111          nsizes = nlinks
 112 !           write(6,*) "Fatal Error: NLINKSL .gt. NLINKS .. "
 113 !           call hydro_stop("not solved, contact WRF-Hydro group. ")
 114       endif
 115       rt_domain(did)%nlinksize = nsizes
 116
 117
 118       if(rt_domain(did)%NLINKS .eq. 0) NLINKS = 1
 119       if(rt_domain(did)%NLAKES .eq. 0) NLAKES = 1
 120       if(rt_domain(did)%NLINKSL .eq. 0) NLINKSL = 1
 121
 122       rt_domain(did)%iswater = 0
 123       rt_domain(did)%isurban = 0
 124       rt_domain(did)%isoilwater = 0
 125
 126 !DJG Allocate routing and disaggregation arrays
 127
 128 #ifdef HYDRO_D
 129   write(6,*) "  rt_allocate ***** ixrt,jxrt, nsoil", ixrt,jxrt, nsoil
 130 #endif
 131
 132   if(nlst(did)%channel_only       .eq. 0 .and. &
 133      nlst(did)%channelBucket_only .eq. 0        ) then
 134
 135      allocate( rt_domain(did)%DSMC      (NSOIL) )
 136      rt_domain(did)%dsmc = 0
 137      allocate( rt_domain(did)%SMCRTCHK          (NSOIL) )
 138      rt_domain(did)%SMCRTCHK = 0
 139      allocate( rt_domain(did)%SH2OAGGRT         (NSOIL) )
 140      rt_domain(did)%SH2OAGGRT = 0
 141      allocate( rt_domain(did)%STCAGGRT          (NSOIL) )
 142      rt_domain(did)%STCAGGRT = 0
 143      allocate( rt_domain(did)%SMCAGGRT          (NSOIL) )
 144      rt_domain(did)%SMCAGGRT = 0
 145
 146      if(nlst(did)%UDMP_OPT .eq. 1) then
 147         allocate ( RT_DOMAIN(did)%landRunOff (ixrt,jxrt) )
 148      endif
 149
 150      !allocate( rt_domain(did)%subsurface%grid_transform%smcrt          (IXRT,JXRT,NSOIL) )
 151      !rt_domain(did)%subsurface%grid_transform%smcrt    = 0.0
 152      allocate( rt_domain(did)%soiltypRT         (IXRT,JXRT) )
 153      !!
 154
 155      !allocate( rt_domain(did)%overland%properties%surface_slope_x      (IXRT,JXRT) )
 156      !rt_domain(did)%overland%properties%surface_slope_x        = 0.0
 157      !allocate( rt_domain(did)%overland%properties%surface_slope_y      (IXRT,JXRT) )
 158      !rt_domain(did)%overland%properties%surface_slope_y        = 0.0
 159      !allocate( rt_domain(did)%overland%properties%surface_slope        (IXRT,JXRT,8) )
 160      !rt_domain(did)%overland%properties%surface_slope          = -999
 161      !allocate( rt_domain(did)%overland%properties%max_surface_slope_index      (IXRT,JXRT,3) )
 162      !rt_domain(did)%overland%properties%max_surface_slope_index        = 0.0
 163      !allocate( rt_domain(did)%overland%properties%roughness   (IXRT,JXRT) )
 164      !
 165
 166      !allocate( rt_domain(did)%QSUBBDRYTRT   (IXRT,JXRT) )
 167      !rt_domain(did)%QSUBBDRYTRT = 0.0
 168      allocate( rt_domain(did)%OVROUGHRTFAC   (IXRT,JXRT) )
 169      !rt_domain(did)%overland%properties%roughness   = 0.0
 170      !allocate( rt_domain(did)%overland%properties%retention_depth    (IXRT,JXRT) )
 171      !
 172
 173      allocate( rt_domain(did)%RETDEPRTFAC    (IXRT,JXRT) )
 174      !
 175
 176      !allocate( rt_domain(did)%overland%control%surface_water_head_routing(IXRT,JXRT) )
 177      !rt_domain(did)%overland%control%surface_water_head_routing= 0.0
 178      !allocate( rt_domain(did)%overland%control%infiltration_excess   (IXRT,JXRT) )
 179      !rt_domain(did)%overland%control%infiltration_excess   = 0.0
 180      allocate( rt_domain(did)%INFXSWGT    (IXRT,JXRT) )
 181      rt_domain(did)%INFXSWGT    = 0.0
 182      !allocate( rt_domain(did)%LKSATRT     (IXRT,JXRT) )
 183      !rt_domain(did)%LKSATRT     = 0.0
 184      allocate( rt_domain(did)%LKSATFAC    (IXRT,JXRT) )
 185      rt_domain(did)%LKSATFAC    = 0.0
 186
 187      allocate( rt_domain(did)%IMPERVFRAC    (IXRT,JXRT) )
 188      rt_domain(did)%IMPERVFRAC    = 0.0
 189
 190      !allocate( rt_domain(did)%subsurface%state%qsubrt      (IXRT,JXRT) )
 191      !rt_domain(did)%subsurface%state%qsubrt      = 0.0
 192      !allocate( rt_domain(did)%subsurface%properties%zwattablrt  (IXRT,JXRT) )
 193      !rt_domain(did)%subsurface%properties%zwattablrt  = 0.0
 194      !allocate( rt_domain(did)%QSUBBDRYRT  (IXRT,JXRT) )
 195      !rt_domain(did)%QSUBBDRYRT  = 0.0
 196      !allocate( rt_domain(did)%subsurface%properties%soldeprt    (IXRT,JXRT) )
 197      !rt_domain(did)%subsurface%properties%soldeprt    = 0.0
 198      allocate( rt_domain(did)%q_sfcflx_x  (IXRT,JXRT) )
 199      rt_domain(did)%q_sfcflx_x  = 0.0
 200      allocate( rt_domain(did)%q_sfcflx_y  (IXRT,JXRT) )
 201      rt_domain(did)%q_sfcflx_y  = 0.0
 202      !allocate( rt_domain(did)%subsurface%grid_transform%smcmaxrt       (IXRT,JXRT,NSOIL) )
 203      !rt_domain(did)%subsurface%grid_transform%smcmaxrt         = 0.0
 204      !allocate( rt_domain(did)%subsurface%grid_transform%smcwltrt       (IXRT,JXRT,NSOIL) )
 205      !rt_domain(did)%subsurface%grid_transform%smcwltrt         = 0.0
 206      allocate( rt_domain(did)%SH2OWGT           (IXRT,JXRT,NSOIL) )
 207      rt_domain(did)%SH2OWGT     = 0.0
 208      allocate( rt_domain(did)%INFXSAGGRT        (IXRT,JXRT) )
 209      rt_domain(did)%INFXSAGGRT  = 0.0
 210      !allocate( rt_domain(did)%overland%control%dhrt    (IXRT,JXRT) ) ! moved to overland control
 211      !rt_domain(did)%overland%control%dhrt      = 0.0                 ! moved to overland control
 212      !allocate( rt_domain(did)%overland%streams_and_lakes%surface_water_to_channel (IXRT,JXRT) )              ! moved to overland streams and lakes
 213      !rt_domain(did)%overland%streams_and_lakes%surface_water_to_channel = 0.0                                ! moved to overland streams and lakes
 214      allocate( rt_domain(did)%QSTRMVOLRT_TS  (IXRT,JXRT) )
 215      rt_domain(did)%QSTRMVOLRT_TS  = 0.0
 216      allocate( rt_domain(did)%QSTRMVOLRT_ACC  (IXRT,JXRT) )
 217      rt_domain(did)%QSTRMVOLRT_ACC  = 0.0
 218      !allocate( rt_domain(did)%overland%control%boundary_flux           (IXRT,JXRT) )
 219      !rt_domain(did)%overland%control%boundary_flux     = 0.0
 220      allocate( rt_domain(did)%SUB_RESID (ixrt,jxrt) )
 221      rt_domain(did)%SUB_RESID = 0.0
 222
 223      ! tmp array
 224      !allocate( rt_domain(did)%subsurface%grid_transform%smcrefrt       (IXRT,JXRT,NSOIL) )
 225      ! tmp
 226
 227      !! Variables (formerly?) needed for channel_only
 228      allocate( rt_domain(did)%ELRT      (IXRT,JXRT) )
 229      rt_domain(did)%ELRT        = 0.0
 230      !allocate( rt_domain(did)%overland%streams_and_lakes%lake_mask     (IXRT,JXRT) ) ! moved to overland%stream_and_lakes
 231      !rt_domain(did)%overland%streams_and_lakes%lake_mask       = -9999               ! moved to overland%streams_and_lakes
 232      !allocate( rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake(IXRT,JXRT) )                              ! moved to overland%streams_and_lakes
 233      !!rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake= 0.0                                               ! moved to overland%streams_and_lakes
 234      allocate( rt_domain(did)%LAKE_INFLORT_TS(IXRT,JXRT) )
 235      allocate( rt_domain(did)%LAKE_INFLORT_DUM(IXRT,JXRT) )
 236      rt_domain(did)%LAKE_INFLORT_DUM= 0.0
 237      allocate( rt_domain(did)%LATVAL (ixrt,jxrt) )
 238      allocate( rt_domain(did)%LONVAL (ixrt,jxrt) )
 239      rt_domain(did)%LONVAL = 0.0
 240      rt_domain(did)%LATVAL = 0.0
 241
 242
 243   !DJG Allocate routing and disaggregation arrays
 244   allocate(rt_domain(did)%qinflowbase  (IXRT,JXRT) )
 245   rt_domain(did)%qinflowbase = 0.0
 246
 247   allocate(rt_domain(did)%gw_strm_msk  (IXRT,JXRT) )
 248            rt_domain(did)%gw_strm_msk   = 0
 249   allocate(rt_domain(did)%gw_strm_msk_lind  (IXRT,JXRT) )
 250
 251   ! allocate land surface grid variables
 252   allocate( rt_domain(did)%SMC  (IX,JX,NSOIL) )
 253             rt_domain(did)%SMC   = 0.25
 254   allocate( rt_domain(did)%SICE (IX,JX,NSOIL) )
 255             rt_domain(did)%SICE  = 0.
 256   ! allocate( rt_domain(did)%dist_lsm (ixrt,jxrt,9) )
 257   ! allocate( rt_domain(did)%lat_lsm (ixrt,jxrt) )
 258   ! allocate( rt_domain(did)%lon_lsm (ixrt,jxrt) )
 259
 260   ! allocate( rt_domain(did)%SICE  (IX,JX,NSOIL) )
 261   allocate( rt_domain(did)%SMCMAX1  (IX,JX) )
 262             rt_domain(did)%SMCMAX1   = 0.0
 263            !rt_domain(did)%SMCMAX1   = 0.434
 264   allocate( rt_domain(did)%STC  (IX,JX,NSOIL) )
 265             rt_domain(did)%STC   = 282.0
 266   allocate( rt_domain(did)%SH2OX(IX,JX,NSOIL) )
 267             rt_domain(did)%SH2OX = rt_domain(did)%SMC
 268   allocate( rt_domain(did)%SMCWLT1  (IX,JX) )
 269             rt_domain(did)%SMCWLT1   = 0.0
 270   allocate( rt_domain(did)%SMCREF1  (IX,JX) )
 271             rt_domain(did)%SMCREF1   = 0.0
 272   allocate( rt_domain(did)%VEGTYP   (IX,JX) )
 273             rt_domain(did)%VEGTYP    = 0
 274
 275   allocate( rt_domain(did)%OV_ROUGH2d   (IX,JX) )
 276
 277   allocate( rt_domain(did)%SOILTYP   (IX,JX) )
 278
 279   allocate( rt_domain(did)%GWSUBBASMSK   (IX,JX) )
 280             rt_domain(did)%GWSUBBASMSK    = 0
 281   !allocate( rt_domain(did)%subsurface%properties%sldpth(NSOIL) )
 282   !          rt_domain(did)%subsurface%properties%sldpth = 0.0
 283   allocate( rt_domain(did)%SO8LD_D   (IX,JX,3) )
 284             rt_domain(did)%SO8LD_D    = 0.0
 285   allocate( rt_domain(did)%SO8LD_Vmax   (IX,JX) )
 286             rt_domain(did)%SO8LD_Vmax    = 0.0
 287   !allocate( rt_domain(did)%sfcheadrt   (IX,JX) ) !moved to overland control structure
 288   !          rt_domain(did)%sfcheadrt    = 0.0    !moved to overland control structure
 289   allocate( rt_domain(did)%INFXSRT   (IX,JX) )
 290             rt_domain(did)%INFXSRT    = 0.0
 291   allocate( rt_domain(did)%TERRAIN   (IX,JX) )
 292             rt_domain(did)%TERRAIN    = 0.0
 293   allocate( rt_domain(did)%LKSAT   (IX,JX) )
 294             rt_domain(did)%LKSAT    = 0.0
 295   allocate( rt_domain(did)%SOLDRAIN   (IX,JX) )
 296             rt_domain(did)%SOLDRAIN    = 0.0
 297
 298   allocate( rt_domain(did)%NEXP   (IX,JX) )
 299             rt_domain(did)%NEXP    = 1.0
 300
 301   end if ! neither channel_only nor channelBucket_only
 302
 303
 304   !! needed regardless
 305   allocate( rt_domain(did)%dist_lsm (ix,jx,9) )
 306             rt_domain(did)%dist_lsm = 0.0
 307   allocate( rt_domain(did)%lat_lsm (ix,jx) )
 308   allocate( rt_domain(did)%lon_lsm (ix,jx) )
 309   rt_domain(did)%timestep_flag = 1    ! default is cold start
 310   !allocate( rt_domain(did)%overland%properties%distance_to_neighbor (ixrt,jxrt,9) ) ! moved to overland%properties
 311   !rt_domain(did)%overland%properties%distance_to_neighbor = -999                    ! moved to overland%properties
 312
 313   !! This is needed for channelBucket_only
 314   !! because the bucket area (basns_area) depends on the initialization of the
 315   !! UDMP code, this is a required variable.
 316   !! JLM: could these be deallocated under channel_only
 317   if(nlst(did)%channel_only       .eq. 0) then
 318      !allocate( rt_domain(did)%overland%streams_and_lakes%ch_netrt      (IXRT,JXRT) ) !moved to overland%streams_and_lakes
 319      !rt_domain(did)%overland%streams_and_lakes%ch_netrt        = 0.0                 !moved to overland%streams_and_lakes
 320      allocate( rt_domain(did)%CH_LNKRT (IXRT,JXRT) )
 321      rt_domain(did)%CH_LNKRT = 0.0
 322   endif
 323
 324
 325   if (CHANRTSWCRT.eq.1 .or. CHANRTSWCRT .eq. 2) then  !IF/then for channel routing
 326
 327      !! JLM TODO: clean up this section for routing options, group 2D variables.
 328
 329      allocate( rt_domain(did)%CH_NETLNK (IXRT,JXRT) )
 330      rt_domain(did)%CH_NETLNK = 0.0
 331      allocate( rt_domain(did)%GCH_NETLNK (IXRT,JXRT) )
 332      rt_domain(did)%GCH_NETLNK = 0.0
 333
 334 #ifdef MPP_LAND
 335      allocate( rt_domain(did)%LAKE_INDEX(NLAKES) )
 336      rt_domain(did)%lake_index = -99
 337      allocate( rt_domain(did)%nlinks_INDEX(nsizes) )
 338      allocate( rt_domain(did)%Link_location(ixrt,jxrt) )
 339 #endif
 340
 341      allocate( rt_domain(did)%CH_LNKRT_SL (IXRT,JXRT) )
 342      rt_domain(did)%CH_LNKRT_SL = -99
 343      rt_domain(did)%MAXORDER = -9999
 344
 345 !tmp  if( nlst_rt(did)%channel_option .eq. 1  .or. nlst_rt(did)%channel_option .eq. 3 ) then
 346 !tmp       NLINKS = rt_domain(did)%NLINKSL
 347 !tmp       NLAKES = rt_domain(did)%NLINKSL
 348 !tmp  endif
 349
 350      allocate( rt_domain(did)%LINKID(nsizes) )
 351      allocate( rt_domain(did)%gages(nsizes) )
 352      allocate( rt_domain(did)%TO_NODE(nsizes) )
 353      allocate( rt_domain(did)%FROM_NODE(nsizes) )
 354      allocate( rt_domain(did)%CHLAT(nsizes) )   !-latitutde of channel grid point
 355      allocate( rt_domain(did)%CHLON(nsizes) )   !-longitude of channel grid point
 356      allocate( rt_domain(did)%ZELEV(nsizes) )
 357      allocate( rt_domain(did)%TYPEL(nsizes) )
 358      allocate( rt_domain(did)%ORDER(nsizes) )
 359      allocate( rt_domain(did)%QLINK(nsizes,2) )
 360
 361 #ifdef WRF_HYDRO_NUDGING
 362      allocate( rt_domain(did)%nudge(nsizes) )
 363 #endif
 364      allocate( rt_domain(did)%MUSK(nsizes) )
 365      allocate( rt_domain(did)%MUSX(nsizes) )
 366      allocate( rt_domain(did)%CHANLEN(nsizes) )
 367      allocate( rt_domain(did)%MannN(nsizes))
 368      allocate( rt_domain(did)%So(nsizes) )
 369      allocate( rt_domain(did)%ChSSlp(nsizes) )
 370      allocate( rt_domain(did)%Bw(nsizes) )
 371      allocate( rt_domain(did)%Tw(nsizes) )
 372      allocate( rt_domain(did)%Tw_CC(nsizes) )
 373      allocate( rt_domain(did)%n_CC(nsizes) )
 374      allocate( rt_domain(did)%ChannK(nsizes) )
 375      allocate( rt_domain(did)%LAKEIDA(nsizes) )
 376      allocate( rt_domain(did)%LAKEIDX(nsizes) )
 377
 378      if(NLAKES .gt. 0) then
 379
 380         allocate ( rt_domain(did)%reservoirs(NLAKES) )    ! allocate array of pointers to reservoirs
 381         allocate ( rt_domain(did)%reservoir_type(NLAKES) )        ! allocate array to specify type of reservoir
 382         allocate ( rt_domain(did)%final_reservoir_type(NLAKES) )  ! allocate array to specify final type of reservoir
 383         allocate ( rt_domain(did)%reservoir_assimilated_value(NLAKES) )  ! allocate array to specify assimilated value to reservoir discharge
 384         allocate ( rt_domain(did)%reservoir_assimilated_source_file(NLAKES) )  ! allocate array to specify assimilated source file
 385         allocate( rt_domain(did)%LAKEIDM(NLAKES) )
 386         allocate( rt_domain(did)%HRZAREA(NLAKES) )
 387         allocate( rt_domain(did)%LAKEMAXH(NLAKES) )
 388         allocate( rt_domain(did)%ELEVLAKE(NLAKES) )
 389         allocate( rt_domain(did)%WEIRH(NLAKES) )
 390         allocate( rt_domain(did)%WEIRC(NLAKES) )
 391         allocate( rt_domain(did)%WEIRL(NLAKES) )
 392         allocate( rt_domain(did)%DAML(NLAKES) )
 393         allocate( rt_domain(did)%ORIFICEC(NLAKES) )
 394         allocate( rt_domain(did)%ORIFICEA(NLAKES) )
 395         allocate( rt_domain(did)%ORIFICEE(NLAKES) )
 396
 397          rt_domain(did)%HRZAREA = 0.0
 398          rt_domain(did)%WEIRH = 0.0
 399          rt_domain(did)%WEIRC = 0.0
 400          rt_domain(did)%WEIRL = 0.0
 401          rt_domain(did)%DAML  = 0.0
 402          rt_domain(did)%LAKEMAXH = 0.0
 403          rt_domain(did)%ELEVLAKE= 0.0
 404          rt_domain(did)%ORIFICEC = 0.0
 405          rt_domain(did)%ORIFICEA = 0.0
 406          rt_domain(did)%ORIFICEE = 0.0
 407          rt_domain(did)%reservoir_type = 1
 408          rt_domain(did)%final_reservoir_type = 1
 409          rt_domain(did)%reservoir_assimilated_value = -9999.0
 410          rt_domain(did)%reservoir_assimilated_source_file = repeat(char(0), 256)
 411      endif
 412
 413
 414 !    allocate( rt_domain(did)%LAKEMAXH(nsizes) )
 415 !    allocate( rt_domain(did)%WEIRC(nsizes) )
 416 !    allocate( rt_domain(did)%WEIRL(nsizes) )
 417 !    allocate( rt_domain(did)%ORIFICEC(nsizes) )
 418 !    allocate( rt_domain(did)%ORIFICEA(nsizes) )
 419 !    allocate( rt_domain(did)%ORIFICEE(nsizes) )
 420
 421      if(nsizes .gt. 0) then
 422         if(nlst(did)%output_channelBucket_influx .eq. 1 .or. &
 423            nlst(did)%output_channelBucket_influx .eq. 2      ) then
 424            allocate( rt_domain(did)%accSfcLatRunoff(1) )
 425            allocate( rt_domain(did)%accBucket(      1) )
 426            allocate( rt_domain(did)%qSfcLatRunoff(  nsizes) )
 427            allocate( rt_domain(did)%qBucket(        nsizes) )
 428         endif
 429
 430         if(nlst(did)%output_channelBucket_influx .eq. 1 .or. &
 431            nlst(did)%output_channelBucket_influx .eq. 3      ) &
 432            allocate( rt_domain(did)%qBtmVertRunoff(     1) )
 433         if(nlst(did)%output_channelBucket_influx .eq. 2) then
 434              allocate( rt_domain(did)%qBtmVertRunoff(nsizes) )
 435              rt_domain(did)%qBtmVertRunoff  = zeroFlt
 436         endif
 437
 438         if(nlst(did)%output_channelBucket_influx .eq. 3) then
 439            allocate( rt_domain(did)%accSfcLatRunoff(nsizes) )
 440            allocate( rt_domain(did)%accBucket(      nsizes) )
 441            allocate( rt_domain(did)%qSfcLatRunoff(       1) )
 442            allocate( rt_domain(did)%qBucket(             1) )
 443            rt_domain(did)%accSfcLatRunoff = zeroDbl
 444            rt_domain(did)%accBucket       = zeroDbl
 445            rt_domain(did)%qSfcLatRunoff   = zeroFlt
 446            rt_domain(did)%qBucket         = zeroFlt
 447         endif
 448
 449         allocate( rt_domain(did)%QLateral(nsizes) )
 450         allocate( rt_domain(did)%velocity(nsizes) )
 451         rt_domain(did)%QLateral  = zeroFlt
 452         rt_domain(did)%velocity  = zeroFlt
 453         if( nlst(did)%channel_option .eq. 2 ) then
 454            allocate( rt_domain(did)%qloss(nsizes) )
 455            rt_domain(did)%qloss = zeroFlt
 456         endif
 457      endif
 458
 459   if( nlst(did)%channel_option .eq. 1  .or. nlst(did)%channel_option .eq. 2 ) then
 460        NLINKS = rt_domain(did)%NLINKS
 461        NLAKES = rt_domain(did)%NLAKES
 462   endif
 463
 464      allocate( rt_domain(did)%LINK(nsizes) )
 465      allocate( rt_domain(did)%STRMFRXSTPTS(nsizes) )
 466      allocate( rt_domain(did)%CHANXI(nsizes) )
 467      allocate( rt_domain(did)%CHANYJ(nsizes) )
 468      allocate( rt_domain(did)%CVOL(nsizes) )
 469      allocate( rt_domain(did)%LATLAKE(NLAKES) )
 470      allocate( rt_domain(did)%LONLAKE(NLAKES) )
 471 !    allocate( rt_domain(did)%ELEVLAKE(NLAKES) )
 472      allocate( rt_domain(did)%LAKENODE(nsizes) )
 473      allocate( rt_domain(did)%RESHT(NLAKES),STAT=istatus )
 474      allocate( rt_domain(did)%QLAKEI(NLAKES),STAT=istatus )
 475      allocate( rt_domain(did)%QLAKEO(NLAKES),STAT=istatus )
 476
 477      allocate( rt_domain(did)%HLINK(nsizes) )  !--used for diffusion only
 478
 479      allocate( rt_domain(did)%node_area(nsizes) )
 480
 481 !!!! tmp
 482       if(nsizes .gt. 0) then
 483       rt_domain(did)%LINK = 0.0
 484       rt_domain(did)%gages = rt_domain(did)%gageMiss
 485       rt_domain(did)%TO_NODE = 0.0
 486       rt_domain(did)%FROM_NODE = 0
 487       rt_domain(did)%TYPEL = -999
 488       rt_domain(did)%ORDER = 0.0
 489       rt_domain(did)%STRMFRXSTPTS = 0.0
 490       rt_domain(did)%MUSK = 0.0
 491       rt_domain(did)%MUSX = 0.0
 492       rt_domain(did)%CHANXI = 0.0
 493       rt_domain(did)%CHANYJ = 0.0
 494       rt_domain(did)%CHLAT = 0.0         !-latitutde of channel grid point
 495       rt_domain(did)%CHLON = 0.0         !-longitude of channel grid point
 496       rt_domain(did)%CHANLEN = 0.0
 497       rt_domain(did)%ChSSlp = 0.0
 498       rt_domain(did)%Bw = 0.0
 499       rt_domain(did)%Tw = 0.0
 500       rt_domain(did)%Tw_CC = 0.0
 501       rt_domain(did)%n_CC = 0.0
 502       rt_domain(did)%ChannK = 0.0
 503
 504
 505       rt_domain(did)%ZELEV = 0.0
 506       rt_domain(did)%CVOL = 0.0
 507       rt_domain(did)%LAKEIDA = 0
 508       rt_domain(did)%LAKEIDX = 0
 509
 510       rt_domain(did)%LATLAKE = 0.0
 511       rt_domain(did)%LONLAKE = 0.0
 512 !     rt_domain(did)%ELEVLAKE = 0.0
 513       rt_domain(did)%LAKENODE = 0.0
 514       rt_domain(did)%RESHT = 0.0
 515       rt_domain(did)%QLAKEI = 0.0
 516       rt_domain(did)%QLAKEO = 0.0
 517       rt_domain(did)%QLINK = 0
 518 #ifdef WRF_HYDRO_NUDGING
 519       rt_domain(did)%nudge = 0
 520 #endif
 521       rt_domain(did)%HLINK = -999.  !--default to -999 if not found in the restart.
 522       rt_domain(did)%MannN = 0.0
 523       rt_domain(did)%LINKID = 0.0
 524
 525       rt_domain(did)%So = 0.01
 526      endif
 527
 528      rt_domain(did)%restQSTRM = .true.
 529
 530   END IF   !IF/then for channel routing
 531
 532   rt_domain(did)%out_counts = 0
 533   rt_domain(did)%his_out_counts = 0
 534   rt_domain(did)%rst_counts = 1
 535
 536 #ifdef HYDRO_D
 537   write(6,*) "***** finish rt_allocate "
 538 #endif
 539
 540 end subroutine rt_allocate
 541
 542
 543 subroutine getChanDim(did)
 544   use config_base, only: nlst
 545   use module_RT_data, only: rt_domain
 546   implicit none
 547
 548   integer ixrt,jxrt, ix,jx, did, i,j
 549   integer, allocatable,dimension(:,:) :: CH_NETLNK, GCH_NETLNK
 550   !INTEGER, dimension( rt_domain(did)%ixrt,GCH_NETLNK(ixrt,jxrt)) :: GCH_NETLNK, CH_NETLNK
 551
 552   real :: Vmax
 553
 554   ix = rt_domain(did)%ix
 555   jx = rt_domain(did)%jx
 556   ixrt = rt_domain(did)%ixrt
 557   jxrt = rt_domain(did)%jxrt
 558
 559   if(nlst(did)%rtFlag .eq. 0) return
 560
 561   if(nlst(did)%channel_only       .eq. 1 .or. &
 562        nlst(did)%channelBucket_only .eq. 1        ) then
 563
 564    !! Try to avoid some of the 2-d initialization.
 565    !! if this is successful, it most likely will not work for gridded channel (opt 3)
 566
 567    if(my_id .eq. io_id) then
 568       call get_NLINKSL(rt_domain(did)%NLINKSL, nlst(did)%channel_option, nlst(did)%route_link_f)
 569    end if
 570    call mpp_land_bcast_int1(rt_domain(did)%NLINKSL)
 571
 572
 573    if(nlst(did)%channel_option .eq. 1 .or. nlst(did)%channel_option .eq. 2) then
 574       rt_domain(did)%GNLINKSL = rt_domain(did)%NLINKSL
 575
 576       call ReachLS_ini(rt_domain(did)%GNLINKSL,rt_domain(did)%nlinksl,   &
 577            rt_domain(did)%linklsS, rt_domain(did)%linklsE )
 578    else
 579       rt_domain(did)%GNLINKSL = 1
 580       rt_domain(did)%NLINKSL = 1
 581    endif
 582    if(nlst(did)%UDMP_OPT .eq. 1) &
 583         call read_NSIMLAKES(rt_domain(did)%NLAKES,nlst(did)%route_lake_f)
 584
 585    call rt_allocate(did,rt_domain(did)%ix,rt_domain(did)%jx,&
 586         rt_domain(did)%ixrt,rt_domain(did)%jxrt, nlst(did)%nsoil,nlst(did)%CHANRTSWCRT)
 587
 588    return
 589
 590 endif
 591
 592
 593 allocate(CH_NETLNK(ixrt,jxrt))
 594 allocate(GCH_NETLNK(ixrt,jxrt))
 595
 596 if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT .eq. 2) then  !IF/then for channel routing
 597 #ifdef MPP_LAND
 598    call MPP_READ_ROUTEDIM(did, rt_domain(did)%g_IXRT,rt_domain(did)%g_JXRT, &
 599                           GCH_NETLNK, rt_domain(did)%GNLINKS, &
 600 #else
 601    call READ_ROUTEDIM( &
 602 #endif
 603               IXRT, JXRT, nlst(did)%route_chan_f, nlst(did)%route_link_f, &
 604               nlst(did)%route_direction_f, &
 605               rt_domain(did)%NLINKS, &
 606               CH_NETLNK, nlst(did)%channel_option, nlst(did)%geo_finegrid_flnm, &
 607               rt_domain(did)%NLINKSL, nlst(did)%udmp_opt , rt_domain(did)%nlakes)
 608 #ifndef MPP_LAND
 609    call get_NLINKSL(rt_domain(did)%NLINKSL, nlst(did)%channel_option, nlst(did)%route_link_f)
 610 #endif
 611
 612 #ifdef HYDRO_D
 613    write(6,*) "before rt_allocate after READ_ROUTEDIM"
 614 #endif
 615
 616    if(nlst(did)%channel_option .eq. 1 .or. nlst(did)%channel_option .eq. 2) then
 617
 618       rt_domain(did)%GNLINKSL = rt_domain(did)%NLINKSL
 619
 620 #ifdef MPP_LAND
 621       call ReachLS_ini(rt_domain(did)%GNLINKSL,rt_domain(did)%nlinksl,   &
 622            rt_domain(did)%linklsS, rt_domain(did)%linklsE )
 623 #else
 624       rt_domain(did)%linklsS = 1
 625       rt_domain(did)%linklsE = rt_domain(did)%NLINKSL
 626 #endif
 627    else
 628       rt_domain(did)%GNLINKSL = 1
 629       rt_domain(did)%NLINKSL = 1
 630    endif
 631
 632 #ifndef MPP_LAND
 633    GCH_NETLNK = CH_NETLNK
 634 #endif
 635
 636 endif
 637
 638 if(nlst(did)%UDMP_OPT .eq. 1) then
 639    call read_NSIMLAKES(rt_domain(did)%NLAKES,nlst(did)%route_lake_f)
 640 endif
 641
 642 call rt_allocate(did,rt_domain(did)%ix,rt_domain(did)%jx,&
 643      rt_domain(did)%ixrt,rt_domain(did)%jxrt, nlst(did)%nsoil,nlst(did)%CHANRTSWCRT)
 644
 645
 646 if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT .eq. 2) then  !IF/then for channel routing
 647    rt_domain(did)%CH_NETLNK = CH_NETLNK
 648    rt_domain(did)%GCH_NETLNK = GCH_NETLNK
 649 endif
 650
 651 if(allocated(CH_NETLNK)) deallocate(CH_NETLNK)
 652 if(allocated(GCH_NETLNK)) deallocate(GCH_NETLNK)
 653
 654 end subroutine getChanDim
 655
 656 !===================================================================================================
 657 subroutine LandRT_ini(did)
 658
 659   use module_noah_chan_param_init_rt
 660   use config_base, only: nlst, noah_lsm
 661   use module_RT_data, only: rt_domain
 662   use module_gw_gw2d_data, only: gw2d
 663   use module_HYDRO_io, only: regrid_lowres_to_highres
 664 #ifdef HYDRO_D
 665   use module_HYDRO_io, only: output_lake_types
 666 #endif
 667
 668 #ifdef OUTPUT_CHAN_CONN
 669   use module_nudging_io,        only: output_chan_connectivity
 670 #endif
 671
 672   implicit none
 673
 674   integer :: did
 675   real    :: Vmax
 676
 677   integer :: bas
 678   character(len=19)                      :: header
 679   character(len=1)                       :: jnk
 680
 681   integer :: lake_index, NLAKES_total
 682
 683   real,  dimension(50)     :: BOTWID, TOPWID, HLINK_INIT, CHAN_SS, CHMann !Channel parms from table
 684   real,  dimension(50)     :: TOPWIDCC, NCC    !channnel params of compound
 685   real,  dimension(50)     :: CHANN_K ! Channel infiltration param
 686   integer :: i,j,k, ll, count
 687
 688   integer, allocatable, dimension(:) :: tmp_int
 689   real, allocatable, dimension(:) :: tmp_real
 690   integer, allocatable, dimension(:) :: buf
 691   real, allocatable, dimension(:) :: tmpRESHT
 692   integer :: new_start_i, new_start_j, new_end_i, new_end_j
 693
 694 #ifdef OUTPUT_CHAN_CONN
 695   real :: connCalcTimeStart, connCalcTimeEnd
 696 #endif
 697 !------------------------------------------------------------------------
 698 !DJG Routing Processing
 699 !------------------------------------------------------------------------
 700 !DJG IF/then to get routing terrain fields if either routing module is
 701 !DJG   activated
 702
 703   if(nlst(did)%rtFlag .eq. 0) return
 704
 705   if (nlst(did)%SUBRTSWCRT  .eq.1 .or. &
 706        nlst(did)%OVRTSWCRT   .eq.1 .or. &
 707        nlst(did)%GWBASESWCRT .ne. 0) then
 708
 709      if(nlst(did)%channel_only       .eq. 0 .and. &
 710           nlst(did)%channelBucket_only .eq. 0        ) then
 711
 712 #ifdef HYDRO_D
 713         print *, "Terrain routing initialization..."
 714 #endif
 715
 716         call READ_ROUTING_seq  (  &
 717              rt_domain(did)%IXRT,rt_domain(did)%JXRT,rt_domain(did)%ELRT,rt_domain(did)%overland%streams_and_lakes%ch_netrt, &
 718              rt_domain(did)%CH_LNKRT, &
 719              rt_domain(did)%LKSATFAC,trim(nlst(did)%route_topo_f),&
 720              nlst(did)%route_chan_f,nlst(did)%geo_finegrid_flnm  ,  &
 721              rt_domain(did)%OVROUGHRTFAC,rt_domain(did)%RETDEPRTFAC, rt_domain(did)%IMPERVFRAC, &
 722              nlst(did)%channel_option, nlst(did)%udmp_opt, nlst(did)%imperv_adj)
 723
 724
 725    !yw CALL READ_ROUTING_old(rt_domain(did)%IXRT,rt_domain(did)%JXRT,rt_domain(did)%ELRT,rt_domain(did)%overland%streams_and_lakes%ch_netrt, &
 726
 727         if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT .eq. 2) then  !IF/then for channel routing
 728
 729 #ifdef MPP_LAND
 730            call MPP_READ_CHROUTING_new(    &
 731 #else
 732                 call READ_CHROUTING1( &  !! NOT TESTED
 733 #endif
 734              rt_domain(did)%IXRT,         rt_domain(did)%JXRT,       &
 735              rt_domain(did)%ELRT,         rt_domain(did)%overland%streams_and_lakes%ch_netrt,        &
 736              rt_domain(did)%CH_LNKRT,     rt_domain(did)%overland%streams_and_lakes%lake_mask, &
 737              rt_domain(did)%FROM_NODE,    rt_domain(did)%TO_NODE, &
 738              rt_domain(did)%TYPEL,        rt_domain(did)%ORDER, &
 739              rt_domain(did)%MAXORDER,     rt_domain(did)%NLINKS, &
 740              rt_domain(did)%NLAKES,       rt_domain(did)%CHANLEN, &
 741              rt_domain(did)%MannN,        rt_domain(did)%So, &
 742              rt_domain(did)%ChSSlp,       rt_domain(did)%Bw, &
 743              rt_domain(did)%Tw,                              &
 744              rt_domain(did)%Tw_CC,                           &
 745              rt_domain(did)%n_CC,                            &
 746              rt_domain(did)%ChannK,                          &
 747              rt_domain(did)%HRZAREA,      rt_domain(did)%LAKEMAXH, &
 748              rt_domain(did)%WEIRH,        rt_domain(did)%WEIRC, &
 749              rt_domain(did)%WEIRL,        rt_domain(did)%DAML, &
 750              rt_domain(did)%ORIFICEC,     rt_domain(did)%ORIFICEA, &
 751              rt_domain(did)%ORIFICEE,                           &
 752              nlst(did)%reservoir_type_specified, rt_domain(did)%reservoir_type, &
 753              nlst(did)%reservoir_parameter_file,    &
 754              rt_domain(did)%LATLAKE,      rt_domain(did)%LONLAKE, &
 755              rt_domain(did)%ELEVLAKE,     rt_domain(did)%overland%properties%distance_to_neighbor, &
 756              rt_domain(did)%ZELEV,        rt_domain(did)%LAKENODE,        &
 757              rt_domain(did)%CH_NETLNK,    rt_domain(did)%CHANXI,          &
 758              rt_domain(did)%CHANYJ,       rt_domain(did)%CHLAT,           &
 759              rt_domain(did)%CHLON,        nlst(did)%channel_option,    &
 760              rt_domain(did)%latval,       rt_domain(did)%lonval,          &
 761              rt_domain(did)%STRMFRXSTPTS, nlst(did)%geo_finegrid_flnm, &
 762              nlst(did)%route_lake_f, rt_domain(did)%LAKEIDM,nlst(did)%UDMP_OPT   & !! no comma
 763 #ifdef MPP_LAND
 764                 ,rt_domain(did)%g_IXRT,      rt_domain(did)%g_JXRT      &
 765                 ,rt_domain(did)%gnlinks,     rt_domain(did)%GCH_NETLNK  &
 766                 ,rt_domain(did)%map_l2g,     rt_domain(did)%link_location, &
 767                 rt_domain(did)%yw_mpp_nlinks,rt_domain(did)%lake_index, &
 768                 rt_domain(did)%nlinks_index &
 769 #endif
 770              )
 771
 772       end if  !! CHANRTSWCRT 1 or 2
 773
 774    end if  !! neither channel_only nor channelBucket_only
 775
 776
 777    if((nlst(did)%CHANRTSWCRT    .eq. 1 .or. nlst(did)%CHANRTSWCRT    .eq. 2) .and. &
 778       (nlst(did)%channel_option .eq. 1 .or. nlst(did)%channel_option .eq. 2)        ) then
 779       call read_routelink( &
 780            rt_domain(did)%TO_NODE,       rt_domain(did)%TYPEL,      &
 781            rt_domain(did)%ORDER,         rt_domain(did)%MAXORDER,   &
 782            rt_domain(did)%NLAKES,        rt_domain(did)%MUSK,       &
 783            rt_domain(did)%MUSX,                                     &
 784            rt_domain(did)%QLINK,         rt_domain(did)%CHANLEN,    &
 785            rt_domain(did)%MannN,         rt_domain(did)%So,         &
 786            rt_domain(did)%ChSSlp,        rt_domain(did)%Bw,         &
 787            rt_domain(did)%Tw,                                       &
 788            rt_domain(did)%Tw_CC,                                    &
 789            rt_domain(did)%n_CC,                                     &
 790            rt_domain(did)%ChannK,                                   &
 791            rt_domain(did)%LAKEIDA,       rt_domain(did)%HRZAREA,    &
 792            rt_domain(did)%LAKEMAXH,      rt_domain(did)%WEIRH,      &
 793            rt_domain(did)%WEIRC,         rt_domain(did)%WEIRL,      &
 794            rt_domain(did)%DAML,                                     &
 795            rt_domain(did)%ORIFICEC,      rt_domain(did)%ORIFICEA,   &
 796            rt_domain(did)%ORIFICEE,                                 &
 797            nlst(did)%reservoir_type_specified,                      &
 798            rt_domain(did)%reservoir_type,                           &
 799            nlst(did)%reservoir_parameter_file,                      &
 800            rt_domain(did)%LATLAKE,                                  &
 801            rt_domain(did)%LONLAKE,       rt_domain(did)%ELEVLAKE,   &
 802            rt_domain(did)%LAKEIDM,       rt_domain(did)%LAKEIDX,    &
 803            nlst(did)%route_link_f,       nlst(did)%route_lake_f, &
 804            rt_domain(did)%ZELEV,         rt_domain(did)%CHLAT,      &
 805            rt_domain(did)%CHLON,         rt_domain(did)%NLINKSL,    &
 806            rt_domain(did)%LINKID,        rt_domain(did)%GNLINKSL,   &
 807            rt_domain(did)%NLINKS,        rt_domain(did)%gages,      &
 808            rt_domain(did)%gageMiss                                   )
 809    end if
 810
 811    NLAKES_total = rt_domain(did)%NLAKES
 812
 813    do lake_index = 1, NLAKES_total
 814      nullify(rt_domain(did)%reservoirs(lake_index)%ptr)
 815
 816    end do
 817
 818    ! For loop to cycle array of pointers to reservoirs
 819    do lake_index = 1, NLAKES_total
 820
 821        if (rt_domain(did)%reservoir_type(lake_index) == 1) then
 822           allocate( levelpool :: rt_domain(did)%reservoirs(lake_index)%ptr)
 823
 824        else if (rt_domain(did)%reservoir_type(lake_index) == 2 .and. nlst(did)%reservoir_persistence_usgs) then
 825           allocate(persistence_levelpool_hybrid :: rt_domain(did)%reservoirs(lake_index)%ptr)
 826
 827        else if (rt_domain(did)%reservoir_type(lake_index) == 3 .and. nlst(did)%reservoir_persistence_usace) then
 828           allocate(persistence_levelpool_hybrid :: rt_domain(did)%reservoirs(lake_index)%ptr)
 829
 830        else if (rt_domain(did)%reservoir_type(lake_index) == 4 .and. nlst(did)%reservoir_rfc_forecasts) then
 831           allocate(rfc_forecasts :: rt_domain(did)%reservoirs(lake_index)%ptr)
 832
 833        else if (rt_domain(did)%reservoir_type(lake_index) == 5 .and. nlst(did)%reservoir_rfc_forecasts) then
 834           allocate(rfc_forecasts :: rt_domain(did)%reservoirs(lake_index)%ptr)
 835
 836        else
 837           allocate( levelpool :: rt_domain(did)%reservoirs(lake_index)%ptr)
 838
 839        end if
 840
 841        ! Dynamically type reservoir to initialize.
 842        select type (reservoir => rt_domain(did)%reservoirs(lake_index)%ptr)
 843
 844           type is (levelpool)
 845                call reservoir%init(                               &
 846                      rt_domain(did)%RESHT(lake_index),            &
 847                      rt_domain(did)%HRZAREA(lake_index),          &
 848                      rt_domain(did)%WEIRH(lake_index),            &
 849                      rt_domain(did)%WEIRC(lake_index),            &
 850                      rt_domain(did)%WEIRL(lake_index),            &
 851                      rt_domain(did)%DAML(lake_index),             &
 852                      rt_domain(did)%ORIFICEE(lake_index),         &
 853                      rt_domain(did)%ORIFICEC(lake_index),         &
 854                      rt_domain(did)%ORIFICEA(lake_index),         &
 855                      rt_domain(did)%LAKEMAXH(lake_index),         &
 856                      rt_domain(did)%LAKEIDM(lake_index)            )
 857
 858           type is (persistence_levelpool_hybrid)
 859               call reservoir%init(                                     &
 860                     rt_domain(did)%RESHT(lake_index),                  &
 861                     rt_domain(did)%HRZAREA(lake_index),                &
 862                     rt_domain(did)%WEIRH(lake_index),                  &
 863                     rt_domain(did)%WEIRC(lake_index),                  &
 864                     rt_domain(did)%WEIRL(lake_index),                  &
 865                     rt_domain(did)%DAML(lake_index),                   &
 866                     rt_domain(did)%ORIFICEE(lake_index),               &
 867                     rt_domain(did)%ORIFICEC(lake_index),               &
 868                     rt_domain(did)%ORIFICEA(lake_index),               &
 869                     rt_domain(did)%LAKEMAXH(lake_index),               &
 870                     rt_domain(did)%ELEVLAKE(lake_index),               &
 871                     rt_domain(did)%LAKEIDM(lake_index),                &
 872                     rt_domain(did)%reservoir_type(lake_index),         &
 873                     nlst(did)%reservoir_parameter_file,                &
 874                     nlst(did)%startdate(1:19),                         &
 875                     nlst(did)%reservoir_usgs_timeslice_path,           &
 876                     nlst(did)%reservoir_usace_timeslice_path,          &
 877                     nlst(did)%reservoir_observation_lookback_hours,    &
 878                     nlst(did)%reservoir_observation_update_time_interval_seconds)
 879
 880           type is (rfc_forecasts)
 881               call reservoir%init(                                      &
 882                     rt_domain(did)%RESHT(lake_index),                   &
 883                     rt_domain(did)%HRZAREA(lake_index),                 &
 884                     rt_domain(did)%WEIRH(lake_index),                   &
 885                     rt_domain(did)%WEIRC(lake_index),                   &
 886                     rt_domain(did)%WEIRL(lake_index),                   &
 887                     rt_domain(did)%DAML(lake_index),                    &
 888                     rt_domain(did)%ORIFICEE(lake_index),                &
 889                     rt_domain(did)%ORIFICEC(lake_index),                &
 890                     rt_domain(did)%ORIFICEA(lake_index),                &
 891                     rt_domain(did)%LAKEMAXH(lake_index),                &
 892                     rt_domain(did)%ELEVLAKE(lake_index),                &
 893                     rt_domain(did)%LAKEIDM(lake_index),                 &
 894                     rt_domain(did)%reservoir_type(lake_index),          &
 895                     nlst(did)%reservoir_parameter_file,                 &
 896                     nlst(did)%startdate(1:19),                          &
 897                     nlst(did)%reservoir_rfc_forecasts_time_series_path, &
 898                     nlst(did)%reservoir_rfc_forecasts_lookback_hours)
 899
 900            end select
 901
 902    end do
 903
 904    !ADCHANGE: Add lake reach output
 905 #ifdef HYDRO_D
 906      if(nlst(did)%UDMP_OPT .eq. 1) then
 907         call output_lake_types( rt_domain(did)%GNLINKSL, rt_domain(did)%LINKID, rt_domain(did)%TYPEL )
 908      endif
 909 #endif
 910
 911 #ifdef OUTPUT_CHAN_CONN
 912 #ifdef MPP_LAND
 913      connCalcTimeEnd = MPI_Wtime()
 914 #else
 915      call cpu_time(connCalcTimeEnd)
 916 #endif
 917      if ((nlst(did)%CHANRTSWCRT .eq. 1) .and. (nlst(did)%channel_option .eq. 3)) then
 918         call output_chan_connectivity(       &
 919              rt_domain(did)%CHLAT,     &   !! Channel grid lat
 920              rt_domain(did)%CHLON,     &   !! Channel grid lat
 921              rt_domain(did)%CHANLEN,   &   !! The distance between channel grid centers in m.
 922              rt_domain(did)%FROM_NODE, &   !! Index of a given cell and ...
 923              rt_domain(did)%TO_NODE,   &   !!   ... the index which it flows to.
 924              rt_domain(did)%CHANXI,    &   !! Index on fine/routing
 925              rt_domain(did)%CHANYJ,    &   !!   grid of grid cells.
 926              rt_domain(did)%TYPEL,     &   !! Link type
 927              rt_domain(did)%LAKENODE   &   !! Lake indexing
 928              )
 929      end if
 930
 931      !if(my_id .eq. io_id) &
 932      print '("Time to calculate channel connectivity= ",f6.3," seconds.")', &
 933           connCalcTimeEnd-connCalcTimeStart
 934      call exit(17)  !! bail if you're just calculating output connectivity.
 935 #endif
 936      ! end OUTPUT_CHAN_CONN
 937
 938
 939      ! The UDMP_ini effectively sets the nhd gw bucket area (that field is not used from the file)
 940      !   this may be the only dependence of the nhd_routing on the UDMAPING in channelBucket_only
 941      if(nlst(did)%channel_only .eq. 0) then
 942
 943         if(nlst(did)%UDMP_OPT .eq. 1) then
 944            ! get NHDPLUS mapping function.
 945            !          call UDMP_ini(rt_domain(did)%GNLINKSL,rt_domain(did)%ixrt,rt_domain(did)%jxrt,rt_domain(did)%CH_LNKRT , &
 946            call UDMP_ini( rt_domain(did)%GNLINKSL, rt_domain(did)%ixrt,      &
 947                 rt_domain(did)%jxrt,     rt_domain(did)%overland%streams_and_lakes%ch_netrt , &
 948                 nlst(did)%OVRTSWCRT,  nlst(did)%SUBRTSWCRT,  &
 949                 rt_domain(did)%overland%properties%distance_to_neighbor(:,:,9)                         )
 950 #ifdef HYDRO_D
 951            write(6,*) "after UDMP_ini "
 952            call flush(6)
 953 #endif
 954         endif
 955
 956      end if ! end not channel_only
 957
 958      if ( (nlst(did)%CHANRTSWCRT .eq. 1) .and. &
 959           (nlst(did)%channel_option .eq. 1 .or. nlst(did)%channel_option .eq. 2) ) then
 960 #ifdef MPP_LAND
 961
 962         if(nlst(did)%UDMP_OPT .eq. 1) then
 963            ! NHDPLUS
 964            rt_domain(did)%LNLINKSL = LNUMRSL
 965            allocate(rt_domain(did)%LLINKID(rt_domain(did)%LNLINKSL))
 966            do k = 1,LNUMRSL
 967               rt_domain(did)%LLINKID(k) = LUDRSL(k)%myid
 968            end do
 969
 970         else
 971
 972            allocate (buf(rt_domain(did)%GNLINKS) )
 973            buf = -99
 974            do j = 1, rt_domain(did)%jxrt
 975               do i = 1, rt_domain(did)%ixrt
 976                  if( .not. ( (i .eq. 1 .and. left_id .ge. 0) .or. (i .eq. rt_domain(did)%ixrt .and. right_id .ge. 0) .or.  &
 977                       (j .eq. 1 .and. down_id .ge. 0) .or. (j .eq. rt_domain(did)%jxrt .and. up_id .ge. 0)    )   ) then
 978                     if(rt_domain(did)%CH_LNKRT(i,j) .gt. 0) then
 979                        k = rt_domain(did)%CH_LNKRT(i,j)
 980                        buf(k) = k
 981                     endif
 982                  endif
 983               end do
 984            end do
 985
 986            rt_domain(did)%LNLINKSL = 0
 987            do k = 1, rt_domain(did)%GNLINKS
 988               if(buf(k) .gt. 0) then
 989                  rt_domain(did)%LNLINKSL = rt_domain(did)%LNLINKSL + 1
 990               endif
 991            end do
 992
 993 #ifdef HYDRO_D
 994            write(6,*) "LNLINKSL, NLINKS, GNLINKS =",rt_domain(did)%LNLINKSL,rt_domain(did)%NLINKSL,rt_domain(did)%GNLINKSL
 995            call flush(6)
 996 #endif
 997
 998            allocate(rt_domain(did)%LLINKID(rt_domain(did)%LNLINKSL))
 999
1000            k = 0
1001            do i = 1, rt_domain(did)%GNLINKS
1002               if(buf(i) .gt. 0) then
1003                  k = k + 1
1004                  rt_domain(did)%LLINKID(k) = buf(i)
1005               endif
1006            end do
1007
1008            if(allocated(buf)) deallocate(buf)
1009
1010         endif  ! end if block for UDMP_OPT
1011
1012         !-------------------------------------------
1013         ! Z.Cui: Changed new_start_i to 0, otherwise, it crashes with
1014         !        an out-of-bounds error when the '-check all' is enabled
1015         !        for the ifort compiler.
1016         !        The reason is that CH_LNKRT and CH_LNKRT_SL is defined as
1017         !        1:IXRT and 1:JXRT. Now new_start_i is changed to start from 1.
1018         !-------------------------------------------
1019         new_start_i = 1; new_start_j = 1
1020         new_end_i = rt_domain(did)%ixrt; new_end_j = rt_domain(did)%jxrt
1021
1022         if(left_id .ge. 0) new_start_i = 1
1023         if(right_id .ge. 0) new_end_i = rt_domain(did)%ixrt - 1
1024         if(down_id .ge. 0) new_start_j = 2
1025         if(up_id .ge. 0) new_end_j = rt_domain(did)%jxrt - 1
1026
1027                 ! The following loops are replaced by a hashtable-based algorithm
1028         ! do k = 1, rt_domain(did)%LNLINKSL
1029         !    do j = new_start_j, new_end_j
1030         !       do i = new_start_i, new_end_i
1031         !          if(rt_domain(did)%CH_LNKRT(i,j) .eq. rt_domain(did)%LLINKID(k) ) then
1032         !             rt_domain(did)%CH_LNKRT_SL(i,j) = k   !! mapping
1033         !          endif
1034         !       end  do
1035         !    end do
1036         ! end do
1037
1038         block
1039           type(hash_t) :: hash_table
1040           integer(kind=int64) :: val,ii,jj
1041           logical :: found
1042
1043           call hash_table%set_all_idx(rt_domain(did)%LLINKID, rt_domain(did)%LNLINKSL)
1044           do jj = new_start_j, new_end_j
1045              do ii = new_start_i, new_end_i
1046                 call hash_table%get(rt_domain(did)%CH_LNKRT(ii,jj), val, found)
1047                 if(found .eqv. .true.) then
1048                    rt_domain(did)%CH_LNKRT_SL(ii,jj) = val
1049                 end if
1050              end do
1051           end do
1052           call hash_table%clear()
1053         end block
1054
1055         call getLocalIndx(rt_domain(did)%gnlinksl,rt_domain(did)%LINKID, rt_domain(did)%LLINKID)
1056
1057         call getToInd(rt_domain(did)%LINKID,rt_domain(did)%to_node,rt_domain(did)%toNodeInd,rt_domain(did)%nToInd,rt_domain(did)%gtoNode)
1058 #else
1059         do k = 1, rt_domain(did)%NLINKSL
1060            do j = 1, rt_domain(did)%jxrt
1061               do i = 1, rt_domain(did)%ixrt
1062                  if(rt_domain(did)%CH_LNKRT(i,j) .eq. rt_domain(did)%LINKID(k) ) then
1063                     rt_domain(did)%CH_LNKRT_SL(i,j) = k   !! mapping
1064                  endif
1065               end do
1066            end do
1067         end do
1068 #endif
1069
1070 !!$        ! use gage information in RouteLink like strmfrxstpts
1071 !!$        rt_domain(did)%STRMFRXSTPTS = -9999  !! existing info useless for link-based routing
1072 !!$        count = 1
1073 !!$        do ll=1,rt_domain(did)%NLINKSL
1074 !!$           if(trim(rt_domain(did)%gages(ll)) .ne. trim(rt_domain(did)%gageMiss)) then
1075 !!$              rt_domain(did)%STRMFRXSTPTS(count) = ll
1076 !!$              count = count + 1
1077 !!$           end if
1078 !!$        end do
1079
1080      endif ! end of if (nlst_rt(did)%channel_option .eq. 1 .or. nlst_rt(did)%channel_option .eq. 2)
1081
1082   end if ! end of if (nlst_rt(did)%SUBRTSWCRT  .eq.1 .or. &    nlst_rt(did)%OVRTSWCRT   .eq.1 .or. &
1083 !            nlst_rt(did)%GWBASESWCRT .ne. 0) then
1084
1085
1086
1087 !yw       allocate(tmp_int(rt_domain(did)%GNLINKS))
1088 !yw       allocate(tmp_real(rt_domain(did)%GNLINKS))
1089
1090
1091   if(nlst(did)%channel_only       .eq. 0 .and. &
1092        nlst(did)%channelBucket_only .eq. 0        ) then
1093
1094      !DJG Temporary hardwire of RETDEPRT,RETDEP_CHAN
1095      !DJG    will later make this a function of SOLTYP and VEGTYP
1096      !            OVROUGHRT(i,j) = 0.01
1097
1098      rt_domain(did)%overland%properties%retention_depth = 0.001   ! units (mm)
1099      rt_domain(did)%RETDEP_CHAN = 0.001
1100
1101
1102      !DJG Need to insert call for acquiring routing fields here...
1103      !DJG     include as a subroutine in module module_Noahlsm_wrfcode_input.F
1104      !DJG  Calculate terrain slopes 'SOXRT,SOYRT' from subgrid elevation 'ELRT'
1105
1106      rt_domain(did)%overland%properties%surface_slope = -999
1107      Vmax = 0.0
1108      do j=2,rt_domain(did)%JXRT-1
1109         do i=2,rt_domain(did)%IXRT-1
1110            rt_domain(did)%overland%properties%surface_slope_x(i,j)=(rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i+1,j))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,3)
1111            rt_domain(did)%overland%properties%surface_slope_y(i,j)=(rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i,j+1))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,1)
1112            !DJG Introduce reduction in retention depth as a linear function of terrain slope
1113            if (nlst(did)%RT_OPTION.eq.2) then
1114               if (rt_domain(did)%overland%properties%surface_slope_x(i,j).gt.rt_domain(did)%overland%properties%surface_slope_y(i,j)) then
1115                  Vmax=rt_domain(did)%overland%properties%surface_slope_x(i,j)
1116               else
1117                  Vmax=rt_domain(did)%overland%properties%surface_slope_y(i,j)
1118               end if
1119
1120               if (Vmax.gt.0.1) then
1121                  rt_domain(did)%overland%properties%retention_depth(i,j)=0.
1122               else
1123                  rt_domain(did)%RETDEPFRAC=Vmax/0.1
1124                  rt_domain(did)%overland%properties%retention_depth(i,j)=rt_domain(did)%overland%properties%retention_depth(i,j)*(1.-rt_domain(did)%RETDEPFRAC)
1125                  if (rt_domain(did)%overland%properties%retention_depth(i,j).lt.0.) rt_domain(did)%overland%properties%retention_depth(i,j)=0.
1126               end if
1127            end if
1128
1129            rt_domain(did)%overland%properties%surface_slope(i,j,1) = &
1130                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i,j+1))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,1)
1131            rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i
1132            rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j + 1
1133            rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 1
1134            Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,1)
1135
1136            rt_domain(did)%overland%properties%surface_slope(i,j,2) = &
1137                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i+1,j+1))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,2)
1138            if(rt_domain(did)%overland%properties%surface_slope(i,j,2) .gt. Vmax ) then
1139               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i + 1
1140               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j + 1
1141               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 2
1142               Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,2)
1143            end if
1144
1145            rt_domain(did)%overland%properties%surface_slope(i,j,3) = &
1146                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i+1,j))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,3)
1147            if(rt_domain(did)%overland%properties%surface_slope(i,j,3) .gt. Vmax ) then
1148               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i + 1
1149               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j
1150               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 3
1151               Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,3)
1152            end if
1153
1154            rt_domain(did)%overland%properties%surface_slope(i,j,4) = &
1155                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i+1,j-1))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,4)
1156            if(rt_domain(did)%overland%properties%surface_slope(i,j,4) .gt. Vmax ) then
1157               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i + 1
1158               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j - 1
1159               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 4
1160               Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,4)
1161            end if
1162
1163            rt_domain(did)%overland%properties%surface_slope(i,j,5) = &
1164                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i,j-1))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,5)
1165            if(rt_domain(did)%overland%properties%surface_slope(i,j,5) .gt. Vmax ) then
1166               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i
1167               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j - 1
1168               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 5
1169               Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,5)
1170            end if
1171
1172            rt_domain(did)%overland%properties%surface_slope(i,j,6) = &
1173                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i-1,j-1))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,6)
1174            if(rt_domain(did)%overland%properties%surface_slope(i,j,6) .gt. Vmax ) then
1175               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i - 1
1176               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j - 1
1177               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 6
1178               Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,6)
1179            end if
1180
1181            rt_domain(did)%overland%properties%surface_slope(i,j,7) = &
1182                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i-1,j))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,7)
1183            if(rt_domain(did)%overland%properties%surface_slope(i,j,7) .gt. Vmax ) then
1184               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i - 1
1185               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j
1186               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 7
1187               Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,7)
1188            end if
1189
1190            rt_domain(did)%overland%properties%surface_slope(i,j,8) = &
1191                 (rt_domain(did)%ELRT(i,j)-rt_domain(did)%ELRT(i-1,j+1))/rt_domain(did)%overland%properties%distance_to_neighbor(i,j,8)
1192            if(rt_domain(did)%overland%properties%surface_slope(i,j,8) .gt. Vmax ) then
1193               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,1) = i - 1
1194               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,2) = j + 1
1195               rt_domain(did)%overland%properties%max_surface_slope_index(i,j,3) = 8
1196               Vmax = rt_domain(did)%overland%properties%surface_slope(i,j,8)
1197            end if
1198
1199            !DJG Introduce reduction in retention depth as a linear function of terrain slope
1200            if (nlst(did)%RT_OPTION.eq.1) then
1201               if (Vmax.gt.0.75) then
1202                  rt_domain(did)%overland%properties%retention_depth(i,j)=0.
1203               else
1204                  rt_domain(did)%RETDEPFRAC=Vmax/0.75
1205                  rt_domain(did)%overland%properties%retention_depth(i,j)=rt_domain(did)%overland%properties%retention_depth(i,j)*(1.-rt_domain(did)%RETDEPFRAC)
1206                  if (rt_domain(did)%overland%properties%retention_depth(i,j).lt.0.) rt_domain(did)%overland%properties%retention_depth(i,j)=0.
1207               end if
1208            end if
1209         end do
1210      end do
1211
1212      !Apply impervious adjustment to retdeprt (AD)
1213      if (nlst(did)%imperv_adj .ne. 0) then
1214        rt_domain(did)%overland%properties%retention_depth = rt_domain(did)%overland%properties%retention_depth*(1.-rt_domain(did)%impervfrac)
1215      end if
1216
1217      !Apply calibration scaling factors to sfc roughness and retention depth here...
1218      rt_domain(did)%overland%properties%retention_depth = rt_domain(did)%overland%properties%retention_depth * rt_domain(did)%RETDEPRTFAC
1219      rt_domain(did)%overland%properties%roughness = rt_domain(did)%overland%properties%roughness * rt_domain(did)%OVROUGHRTFAC
1220
1221    !ADCHANGE: Moved this channel cell setting from OV_RTNG so it is outside
1222    !of overland routine (frequently called) and time loop.
1223    !Force channel retention depth to be 5mm.
1224    ! JLM: for DWJ I'm leaving this next line for you to verify as
1225    !      it's the one I translated to the following line,
1226    !where (rt_domain(did)%CH_NETRT .ge. 0) rt_domain(did)%RETDEPRT = 5.0
1227      where (rt_domain(did)%overland%streams_and_lakes%ch_netrt .ge. 0) &
1228           rt_domain(did)%overland%properties%retention_depth = 5.0
1229
1230    ! calculate the slope for boundary
1231 #ifdef MPP_LAND
1232      if(right_id .lt. 0) rt_domain(did)%overland%properties%surface_slope_x(rt_domain(did)%IXRT,:)= &
1233         rt_domain(did)%overland%properties%surface_slope_x(rt_domain(did)%IXRT-1,:)
1234      if(left_id  .lt. 0) rt_domain(did)%overland%properties%surface_slope_x(1,:)=rt_domain(did)%overland%properties%surface_slope_x(2,:)
1235      if(up_id    .lt. 0) rt_domain(did)%overland%properties%surface_slope_y(:,rt_domain(did)%JXRT)= &
1236           rt_domain(did)%overland%properties%surface_slope_y(:,rt_domain(did)%JXRT-1)
1237      if(down_id  .lt. 0) rt_domain(did)%overland%properties%surface_slope_y(:,1)=rt_domain(did)%overland%properties%surface_slope_y(:,2)
1238 #else
1239      rt_domain(did)%overland%properties%surface_slope_x(rt_domain(did)%IXRT,:)=rt_domain(did)%overland%properties%surface_slope_x(rt_domain(did)%IXRT-1,:)
1240      rt_domain(did)%overland%properties%surface_slope_x(1,:)=rt_domain(did)%overland%properties%surface_slope_x(2,:)
1241      rt_domain(did)%overland%properties%surface_slope_y(:,rt_domain(did)%JXRT)=rt_domain(did)%overland%properties%surface_slope_y(:,rt_domain(did)%JXRT-1)
1242      rt_domain(did)%overland%properties%surface_slope_y(:,1)=rt_domain(did)%overland%properties%surface_slope_y(:,2)
1243 #endif
1244
1245 #ifdef MPP_LAND
1246    ! communicate the value to
1247      call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%retention_depth,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
1248      call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%roughness,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
1249      call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%surface_slope_x,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
1250      call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%surface_slope_y,rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
1251      do i = 1, 8
1252         call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%surface_slope(:,:,i),rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
1253      end do
1254      do i = 1, 3
1255         call MPP_LAND_COM_INTEGER(rt_domain(did)%overland%properties%max_surface_slope_index(:,:,i),rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
1256      end do
1257 #endif
1258
1259   end if ! end of neither channel_only nor channelBucket_only
1260
1261   if(nlst(did)%UDMP_OPT .eq. 1) then
1262
1263      allocate (rt_domain(did)%qout_gwsubbas (rt_domain(did)%nlinksL))
1264      rt_domain(did)%qout_gwsubbas = 0
1265      ! use different baseflow for NHDPlus
1266      if (nlst(did)%GWBASESWCRT.ge.1) then
1267         rt_domain(did)%numbasns = rt_domain(did)%NLINKSL
1268         RT_DOMAIN(did)%gnumbasns = rt_domain(did)%gNLINKSL
1269
1270         allocate (rt_domain(did)%z_gwsubbas (rt_domain(did)%numbasns  ))
1271         allocate (rt_domain(did)%nhdBuckMask(rt_domain(did)%numbasns  ))  ! default is -999
1272
1273         allocate (rt_domain(did)%qin_gwsubbas (rt_domain(did)%numbasns))
1274         allocate (rt_domain(did)%gwbas_pix_ct (rt_domain(did)%numbasns))
1275         allocate (rt_domain(did)%ct2_bas (rt_domain(did)%numbasns))
1276         allocate (rt_domain(did)%bas_pcp (rt_domain(did)%numbasns))
1277         allocate (rt_domain(did)%gw_buck_coeff (rt_domain(did)%numbasns))
1278         allocate (rt_domain(did)%bas_id (rt_domain(did)%numbasns))
1279         allocate (rt_domain(did)%gw_buck_exp(rt_domain(did)%numbasns))
1280         allocate (rt_domain(did)%z_max (rt_domain(did)%numbasns))
1281         allocate (rt_domain(did)%basns_area (rt_domain(did)%numbasns))
1282
1283         if(nlst(did)%bucket_loss .eq. 1) then
1284            allocate(rt_domain(did)%qloss_gwsubbas(rt_domain(did)%numbasns))
1285            allocate(rt_domain(did)%gw_buck_loss(rt_domain(did)%numbasns))
1286
1287            rt_domain(did)%qloss_gwsubbas = 0
1288            rt_domain(did)%gw_buck_loss = 0
1289         endif
1290
1291         rt_domain(did)%qin_gwsubbas = 0
1292         rt_domain(did)%z_gwsubbas = 0
1293         rt_domain(did)%gwbas_pix_ct = 0
1294         rt_domain(did)%bas_pcp = 0
1295
1296         rt_domain(did)%gw_buck_coeff = 0.04
1297         rt_domain(did)%gw_buck_exp  = 0.2
1298         rt_domain(did)%z_max = 0.1
1299
1300         !Temporary hardwire...
1301         rt_domain(did)%z_gwsubbas = 0.05   ! This gets updated with spun-up GW level in GWBUCKPARM.TBL
1302
1303         call readBucket_nhd(trim(nlst(did)%GWBUCKPARM_file), rt_domain(did)%numbasns, &
1304              rt_domain(did)%gw_buck_coeff, rt_domain(did)%gw_buck_exp, rt_domain(did)%gw_buck_loss, &
1305              rt_domain(did)%z_max, rt_domain(did)%z_gwsubbas, rt_domain(did)%LINKID(1:rt_domain(did)%numbasns),  &
1306              rt_domain(did)%nhdBuckMask )
1307
1308       !ADCHANGE: Added in read for z_init from GWBUCKPARM. Units coming in are mm but z_gwsubbas is in m
1309       !for UDMP=1 so we convert units.
1310         rt_domain(did)%z_gwsubbas = rt_domain(did)%z_gwsubbas/1000.
1311
1312 #ifdef HYDRO_D
1313         write(6,*) "finish readBucket_nhd "
1314         call flush(6)
1315 #endif
1316      endif
1317
1318   else
1319
1320    !---------------------------------------------------------------------
1321    !DJG  If GW/Baseflow activated...Read in req'd fields...
1322    !----------------------------------------------------------------------
1323      if (nlst(did)%GWBASESWCRT.ge.1) then
1324         if (nlst(did)%GWBASESWCRT.eq.1.or.nlst(did)%GWBASESWCRT.eq.2 .or. nlst(did)%GWBASESWCRT.ge.4) then
1325 #ifdef HYDRO_D
1326            print *, "new Simple GW-Bucket Scheme selected, retrieving files..."
1327 #endif
1328 #ifdef MPP_LAND
1329            call MPP_READ_SIMP_GW(              &
1330 #else
1331                 call READ_SIMP_GW(                  &
1332 #endif
1333                 rt_domain(did)%IX,rt_domain(did)%JX,rt_domain(did)%IXRT,&
1334                 rt_domain(did)%JXRT,rt_domain(did)%GWSUBBASMSK,nlst(did)%gwbasmskfil,&
1335                 rt_domain(did)%gw_strm_msk,rt_domain(did)%numbasns,rt_domain(did)%overland%streams_and_lakes%ch_netrt,nlst(did)%AGGFACTRT)
1336
1337
1338            call SIMP_GW_IND(rt_domain(did)%ix,rt_domain(did)%jx,rt_domain(did)%GWSUBBASMSK,  &
1339                 rt_domain(did)%numbasns,rt_domain(did)%gnumbasns,rt_domain(did)%basnsInd)
1340
1341 #ifdef HYDRO_D
1342            write(6,*) "rt_domain(did)%gnumbasns, rt_domain(did)%numbasns, ", rt_domain(did)%gnumbasns , rt_domain(did)%numbasns
1343
1344 #endif
1345 #ifdef MPP_LAND
1346            call collectSizeInd(rt_domain(did)%numbasns)
1347 #endif
1348
1349            call get_gw_strm_msk_lind (rt_domain(did)%IXRT, rt_domain(did)%JXRT, rt_domain(did)%gw_strm_msk,&
1350                 rt_domain(did)%numbasns,rt_domain(did)%basnsInd,rt_domain(did)%gw_strm_msk_lind)
1351
1352            allocate (rt_domain(did)%qout_gwsubbas (rt_domain(did)%numbasns))
1353            allocate (rt_domain(did)%qin_gwsubbas (rt_domain(did)%numbasns))
1354            allocate (rt_domain(did)%z_gwsubbas (rt_domain(did)%numbasns))
1355            allocate (rt_domain(did)%gwbas_pix_ct (rt_domain(did)%numbasns))
1356            allocate (rt_domain(did)%ct2_bas (rt_domain(did)%numbasns))
1357            allocate (rt_domain(did)%bas_pcp (rt_domain(did)%numbasns))
1358            allocate (rt_domain(did)%gw_buck_coeff (rt_domain(did)%numbasns))
1359            allocate (rt_domain(did)%bas_id (rt_domain(did)%numbasns))
1360            allocate (rt_domain(did)%gw_buck_exp(rt_domain(did)%numbasns))
1361            allocate (rt_domain(did)%z_max (rt_domain(did)%numbasns))
1362            allocate (rt_domain(did)%basns_area (rt_domain(did)%numbasns))
1363
1364 #ifdef HYDRO_D
1365            write(6,*)  "end Simple GW-Bucket ..."
1366            print *, "Simple GW-Bucket Scheme selected, retrieving files..."
1367 #endif
1368
1369            !Temporary hardwire...
1370            rt_domain(did)%z_gwsubbas = 1.     ! This gets updated with spun-up GW level in GWBUCKPARM.TBL
1371
1372            call read_GWBUCKPARM(trim(nlst(did)%GWBUCKPARM_file),rt_domain(did)%numbasns,   &
1373                 rt_domain(did)%gnumbasns, rt_domain(did)%basnsInd, &
1374                 rt_domain(did)%gw_buck_coeff, rt_domain(did)%gw_buck_exp, rt_domain(did)%gw_buck_loss, &
1375                 rt_domain(did)%z_max, rt_domain(did)%z_gwsubbas, rt_domain(did)%bas_id, &
1376                 rt_domain(did)%basns_area)
1377
1378 !!! Determine number of stream pixels per GW basin for distribution...
1379 #ifdef MPP_LAND
1380            call pix_ct_1(rt_domain(did)%gw_strm_msk,rt_domain(did)%ixrt,rt_domain(did)%jxrt,rt_domain(did)%gwbas_pix_ct,rt_domain(did)%numbasns, &
1381                 rt_domain(did)%gnumbasns,rt_domain(did)%basnsInd)
1382 #else
1383            rt_domain(did)%gwbas_pix_ct = 0.
1384          !         do k = 1, rt_domain(did)%numbasns
1385          !            bas = rt_domain(did)%basnsInd(k)
1386            do i=1,rt_domain(did)%ixrt
1387               do j=1,rt_domain(did)%jxrt
1388                  if (rt_domain(did)%gw_strm_msk(i,j).gt.0) then
1389                     bas = rt_domain(did)%gw_strm_msk(i,j)
1390                     rt_domain(did)%gwbas_pix_ct(bas) = &
1391                          rt_domain(did)%gwbas_pix_ct(bas)  + 1.0
1392                  endif
1393               end do
1394            end do
1395            !         end do
1396 #endif
1397
1398 #ifdef HYDRO_D
1399            print *, "Starting GW basin levels...",rt_domain(did)%z_gwsubbas
1400 #endif
1401
1402          ! BF gw2d model
1403         elseif (nlst(did)%GWBASESWCRT.eq.3) then
1404
1405            call readGW2d(gw2d(did)%ix, gw2d(did)%jx,     &
1406                 gw2d(did)%hycond, gw2d(did)%ho, &
1407                 gw2d(did)%bot, gw2d(did)%poros, &
1408                 gw2d(did)%ltype, nlst(did)%gwIhShift)
1409
1410            gw2d(did)%elev = rt_domain(did)%elrt
1411
1412         end if ! end if (nlst_rt(did)%GWBASESWCRT.eq.1.or.nlst_rt(did)%GWBASESWCRT.eq.2)
1413
1414      end if ! end if (nlst_rt(did)%GWBASESWCRT.ge.1)
1415
1416 !---------------------------------------------------------------------
1417 !DJG  End if GW/Baseflow activated...
1418 !----------------------------------------------------------------------
1419   endif   !!! end if block for UDMP_OPT .eq. 1
1420
1421
1422
1423 !---------------------------------------------------------------------
1424 !DJG,DNY  If channel routing activated...
1425 !----------------------------------------------------------------------
1426
1427   if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT .eq. 2) then
1428
1429    !---------------------------------------------------------------------
1430    !DJG,DNY  Initalize lake and channel heights, this may be overwritten by RESTART
1431    !--------------------------------------------------------------------
1432      if (nlst(did)%channel_option .eq. 3) then
1433 #ifdef MPP_LAND
1434       ! JLM: Currently compound channel does not work for diffusive wave/gridded channel.
1435       ! This conflict of options is caught in Data_Rec/module_namelist.F
1436       ! Some of the code for reading/using top-width-necessary params from CHANPARM.TBL are available
1437       ! but commented out in the code, since they were a bridge to nowhere.
1438         call mpp_CHAN_PARM_INIT (BOTWID,CHANN_K,HLINK_INIT,CHAN_SS,CHMann)  !Read chan parms from table...
1439       !call mpp_CHAN_PARM_INIT (BOTWID,TOPWID,HLINK_INIT,CHAN_SS,CHMann,TOPWIDCC,NCC)  !Read chan parms from table...
1440 #else
1441         call CHAN_PARM_INIT (BOTWID,CHANN_K,HLINK_INIT,CHAN_SS,CHMann) !,TOPWIDCC,NCC)  !Read chan parms from table...
1442 #endif
1443      end if
1444
1445      if (nlst(did)%channel_option .ne. 3) then
1446
1447 #ifdef MPP_LAND
1448         if(my_id .eq. io_id) then
1449 #endif
1450            allocate(tmpRESHT(rt_domain(did)%nlakes))
1451            tmpRESHT = rt_domain(did)%RESHT
1452 #ifdef MPP_LAND
1453         endif
1454 #endif
1455
1456 #ifdef MPP_LAND
1457         call updateLake_seq(rt_domain(did)%RESHT, rt_domain(did)%NLAKES,tmpRESHT)
1458         if(my_id .eq. io_id) then
1459            if(allocated(tmpRESHT)) deallocate(tmpRESHT)
1460         endif
1461 #endif
1462
1463      else       !-- parameterize according to order of diffusion scheme, or if read from hi res file, use its value
1464                 !--  put condition within the if/then structure, which will assign a value if something is missing in hi res
1465
1466         do j=1,rt_domain(did)%NLINKS
1467
1468            if (rt_domain(did)%ORDER(j) .eq. 1) then    !-- smallest stream reach
1469               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1470                  rt_domain(did)%Bw(j) = BOTWID(rt_domain(did)%ORDER(j))
1471               endif
1472               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1473                  rt_domain(did)%Tw(j) = TOPWID(rt_domain(did)%ORDER(j))
1474               endif
1475               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1476                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(rt_domain(did)%ORDER(j))
1477               endif
1478               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1479                  rt_domain(did)%n_CC(j) = NCC(rt_domain(did)%ORDER(j))
1480               endif
1481               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1482                  rt_domain(did)%ChSSlp(j) = CHAN_SS(rt_domain(did)%ORDER(j))
1483               endif
1484               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1485                  rt_domain(did)%MannN(j) = CHMann(rt_domain(did)%ORDER(j))
1486               endif
1487               rt_domain(did)%HLINK(j) = HLINK_INIT(rt_domain(did)%ORDER(j))
1488            elseif (rt_domain(did)%ORDER(j) .eq. 2) then
1489               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1490                  rt_domain(did)%Bw(j) = BOTWID(rt_domain(did)%ORDER(j))
1491               endif
1492               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1493                  rt_domain(did)%Tw(j) = TOPWID(rt_domain(did)%ORDER(j))
1494               endif
1495               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1496                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(rt_domain(did)%ORDER(j))
1497               endif
1498               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1499                  rt_domain(did)%n_CC(j) = NCC(rt_domain(did)%ORDER(j))
1500               endif
1501               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1502                  rt_domain(did)%ChSSlp(j) = CHAN_SS(rt_domain(did)%ORDER(j))
1503               endif
1504               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1505                  rt_domain(did)%MannN(j) = CHMann(rt_domain(did)%ORDER(j))
1506               endif
1507               rt_domain(did)%HLINK(j) = HLINK_INIT(rt_domain(did)%ORDER(j))
1508            elseif (rt_domain(did)%ORDER(j) .eq. 3) then
1509               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1510                  rt_domain(did)%Bw(j) = BOTWID(rt_domain(did)%ORDER(j))
1511               endif
1512               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1513                  rt_domain(did)%Tw(j) = TOPWID(rt_domain(did)%ORDER(j))
1514               endif
1515               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1516                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(rt_domain(did)%ORDER(j))
1517               endif
1518               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1519                  rt_domain(did)%n_CC(j) = NCC(rt_domain(did)%ORDER(j))
1520               endif
1521               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1522                  rt_domain(did)%ChSSlp(j) = CHAN_SS(rt_domain(did)%ORDER(j))
1523               endif
1524               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1525                  rt_domain(did)%MannN(j) = CHMann(rt_domain(did)%ORDER(j))
1526               endif
1527               rt_domain(did)%HLINK(j) = HLINK_INIT(rt_domain(did)%ORDER(j))
1528            elseif (rt_domain(did)%ORDER(j) .eq. 4) then
1529               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1530                  rt_domain(did)%Bw(j) = BOTWID(rt_domain(did)%ORDER(j))
1531               endif
1532               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1533                  rt_domain(did)%Tw(j) = TOPWID(rt_domain(did)%ORDER(j))
1534               endif
1535               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1536                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(rt_domain(did)%ORDER(j))
1537               endif
1538               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1539                  rt_domain(did)%n_CC(j) = NCC(rt_domain(did)%ORDER(j))
1540               endif
1541               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1542                  rt_domain(did)%ChSSlp(j) = CHAN_SS(rt_domain(did)%ORDER(j))
1543               endif
1544               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1545                  rt_domain(did)%MannN(j) = CHMann(rt_domain(did)%ORDER(j))
1546               endif
1547               rt_domain(did)%HLINK(j) = HLINK_INIT(rt_domain(did)%ORDER(j))
1548            elseif (rt_domain(did)%ORDER(j) .eq. 5) then
1549               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1550                  rt_domain(did)%Bw(j) = BOTWID(rt_domain(did)%ORDER(j))
1551               endif
1552               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1553                  rt_domain(did)%Tw(j) = TOPWID(rt_domain(did)%ORDER(j))
1554               endif
1555               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1556                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(rt_domain(did)%ORDER(j))
1557               endif
1558               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1559                  rt_domain(did)%n_CC(j) = NCC(rt_domain(did)%ORDER(j))
1560               endif
1561               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1562                  rt_domain(did)%ChSSlp(j) = CHAN_SS(rt_domain(did)%ORDER(j))
1563               endif
1564               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1565                  rt_domain(did)%MannN(j) = CHMann(rt_domain(did)%ORDER(j))
1566               endif
1567               rt_domain(did)%HLINK(j) = HLINK_INIT(rt_domain(did)%ORDER(j))
1568            elseif (rt_domain(did)%ORDER(j) .eq. 6) then
1569               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1570                  rt_domain(did)%Bw(j) = BOTWID(rt_domain(did)%ORDER(j))
1571               endif
1572               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1573                  rt_domain(did)%Tw(j) = TOPWID(rt_domain(did)%ORDER(j))
1574               endif
1575               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1576                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(rt_domain(did)%ORDER(j))
1577               endif
1578               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1579                  rt_domain(did)%n_CC(j) = NCC(rt_domain(did)%ORDER(j))
1580               endif
1581               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1582                  rt_domain(did)%ChSSlp(j) = CHAN_SS(rt_domain(did)%ORDER(j))
1583               endif
1584               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1585                  rt_domain(did)%MannN(j) = CHMann(rt_domain(did)%ORDER(j))
1586               endif
1587               rt_domain(did)%HLINK(j) = HLINK_INIT(rt_domain(did)%ORDER(j))
1588            elseif (rt_domain(did)%ORDER(j) .ge. 7) then
1589               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1590                  rt_domain(did)%Bw(j) = BOTWID(rt_domain(did)%ORDER(j))
1591               endif
1592               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1593                  rt_domain(did)%Tw(j) = TOPWID(rt_domain(did)%ORDER(j))
1594               endif
1595               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1596                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(rt_domain(did)%ORDER(j))
1597               endif
1598               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1599                  rt_domain(did)%n_CC(j) = NCC(rt_domain(did)%ORDER(j))
1600               endif
1601               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1602                  rt_domain(did)%ChSSlp(j) = CHAN_SS(rt_domain(did)%ORDER(j))
1603               endif
1604               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1605                  rt_domain(did)%MannN(j) = CHMann(rt_domain(did)%ORDER(j))
1606               endif
1607               rt_domain(did)%HLINK(j) = HLINK_INIT(rt_domain(did)%ORDER(j))
1608            else   !-- the outlets won't have orders since there's no nodes, so
1609               !-- assign the order 5 values
1610
1611               if(rt_domain(did)%Bw(j) .eq. 0.0) then
1612                  rt_domain(did)%Bw(j) = BOTWID(5)
1613               endif
1614
1615               if(rt_domain(did)%Tw(j) .eq. 0.0) then
1616                  rt_domain(did)%Tw(j) = TOPWID(5)
1617               endif
1618               if(rt_domain(did)%Tw_CC(j) .eq. 0.0) then
1619                  rt_domain(did)%Tw_CC(j) = TOPWIDCC(5)
1620               endif
1621               if(rt_domain(did)%n_CC(j) .eq. 0.0) then
1622                  rt_domain(did)%n_CC(j) = NCC(5)
1623               endif
1624               if(rt_domain(did)%ChSSlp(j) .eq. 0.0) then  !if id didn't get set from the hi res file, use the  CHANPARAM
1625                  rt_domain(did)%ChSSlp(j) = CHAN_SS(5)
1626               endif
1627               if(rt_domain(did)%MannN(j) .eq. 0.0) then
1628                  rt_domain(did)%MannN(j) = CHMann(5)
1629               endif
1630               rt_domain(did)%HLINK(j) = HLINK_INIT(5)
1631            endif
1632
1633            rt_domain(did)%CVOL(j) = (rt_domain(did)%Bw(j)+ 1/rt_domain(did)%ChSSLP(j)*rt_domain(did)%HLINK(j))*rt_domain(did)%HLINK(j)*rt_domain(did)%CHANLEN(j) !-- initalize channel volume
1634         end do
1635
1636      endif  !End if channel option eq 3; else;
1637
1638
1639      ! Initialize Lake Elevations for Gridded and NWM routing.
1640      do j=1,rt_domain(did)%NLAKES
1641         rt_domain(did)%RESHT(j) = rt_domain(did)%ORIFICEE(j) + &
1642              ((rt_domain(did)%LAKEMAXH(j) - rt_domain(did)%ORIFICEE(j) )* rt_domain(did)%ELEVLAKE(j))
1643      end do
1644
1645   end if     ! Endif for channel routing setup
1646 !-----------------------------------------------------------------------
1647
1648   if(nlst(did)%channel_only       .eq. 0 .and. &
1649        nlst(did)%channelBucket_only .eq. 0        ) then
1650
1651      rt_domain(did)%INFXSWGT = 1./(nlst(did)%AGGFACTRT*nlst(did)%AGGFACTRT)
1652      rt_domain(did)%SH2OWGT = 1.
1653      rt_domain(did)%subsurface%properties%soldeprt = -1.0 * nlst(did)%ZSOIL8(nlst(did)%NSOIL)
1654      rt_domain(did)%subsurface%state%qsubrt = 0.0
1655      rt_domain(did)%subsurface%properties%zwattablrt = 0.0
1656      rt_domain(did)%subsurface%state%qsubbdryrt = 0.0
1657      rt_domain(did)%overland%streams_and_lakes%surface_water_to_channel = 0.0
1658      rt_domain(did)%overland%control%boundary_flux = 0.0
1659      rt_domain(did)%overland%control%surface_water_head_routing = 0.0
1660      rt_domain(did)%overland%control%infiltration_excess = 0.0
1661      rt_domain(did)%overland%control%dhrt = 0.0
1662      rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake = 0.0
1663      rt_domain(did)%LAKE_CT = 0
1664      rt_domain(did)%STRM_CT = 0
1665      rt_domain(did)%SOLDRAIN = 0.0
1666      rt_domain(did)%qinflowbase = 0.0
1667
1668    !  rt_domain(did)%BASIN_MSK = 1
1669    ! !DJG Initialize mass balance check variables...
1670      rt_domain(did)%SMC_INIT=0.
1671      rt_domain(did)%DSMC=0.
1672      rt_domain(did)%DACRAIN=0.
1673      rt_domain(did)%DSFCEVP=0.
1674      rt_domain(did)%DCANEVP=0.
1675      rt_domain(did)%DEDIR=0.
1676      rt_domain(did)%DETT=0.
1677      rt_domain(did)%DEPND=0.
1678      rt_domain(did)%DESNO=0.
1679      rt_domain(did)%DSFCRNFF=0.
1680      rt_domain(did)%DQBDRY=0.
1681      rt_domain(did)%overland%mass_balance%pre_infiltration_excess=0.
1682
1683   end if ! end of neither channel_only nor channelBucket_only
1684
1685    ! Deallocate existing arrays that contain reservoir parameters/properties
1686    if(allocated(rt_domain(did)%HRZAREA)) deallocate(rt_domain(did)%HRZAREA)
1687    if(allocated(rt_domain(did)%WEIRH)) deallocate(rt_domain(did)%WEIRH)
1688    if(allocated(rt_domain(did)%WEIRC)) deallocate(rt_domain(did)%WEIRC)
1689    if(allocated(rt_domain(did)%WEIRL)) deallocate(rt_domain(did)%WEIRL)
1690    if(allocated(rt_domain(did)%DAML)) deallocate(rt_domain(did)%DAML)
1691    if(allocated(rt_domain(did)%ORIFICEE)) deallocate(rt_domain(did)%ORIFICEE)
1692    if(allocated(rt_domain(did)%ORIFICEC)) deallocate(rt_domain(did)%ORIFICEC)
1693    if(allocated(rt_domain(did)%ORIFICEA)) deallocate(rt_domain(did)%ORIFICEA)
1694    if(allocated(rt_domain(did)%LAKEMAXH)) deallocate(rt_domain(did)%LAKEMAXH)
1695
1696 end subroutine LandRT_ini
1697
1698 subroutine deriveFromNode(did)
1699   implicit none
1700   integer :: did
1701   integer :: i,j, kk, maxv
1702   integer :: tmp(rt_domain(did)%nlinks)
1703   tmp = 0
1704   maxv = 1
1705   do i = 1, rt_domain(did)%nlinks
1706      if(rt_domain(did)%to_node(i) .gt. 0) then
1707         kk = rt_domain(did)%to_node(i)
1708         tmp(kk) = tmp(kk) + 1
1709         if(maxv .lt. tmp(kk)) maxv = tmp(kk)
1710      end if
1711   end do
1712   allocate(rt_domain(did)%pnode(rt_domain(did)%nlinks,maxv+1) )
1713   rt_domain(did)%maxv_p = maxv+1
1714   rt_domain(did)%pnode = -99
1715   rt_domain(did)%pnode(:,1) = 1
1716   do i = 1, rt_domain(did)%nlinks
1717      if(rt_domain(did)%to_node(i) .gt. 0) then
1718         j = rt_domain(did)%to_node(i)
1719         rt_domain(did)%pnode(j,1) = rt_domain(did)%pnode(j,1) + 1
1720         kk = rt_domain(did)%pnode(j,1)
1721         rt_domain(did)%pnode(j,kk) = i
1722      end if
1723   end do
1724
1725 end subroutine deriveFromNode
1726
1727 END MODULE module_Routing