hydro/HYDRO_drv/module_HYDRO_drv.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_HYDRO_drv
  22 #ifdef MPP_LAND
  23    use module_HYDRO_io, only:  output_rt, mpp_output_chrt, mpp_output_lakes, mpp_output_chrtgrd, &
  24                                restart_out_bi, restart_in_bi, mpp_output_chrt2, mpp_output_lakes2, &
  25                                hdtbl_in_nc, hdtbl_out
  26    USE module_mpp_land
  27 #else
  28    use module_HYDRO_io, only:  output_rt, output_chrt, output_chrt2, output_lakes
  29 #endif
  30    use module_NWM_io, only: output_chrt_NWM, output_rt_NWM, output_lakes_NWM,&
  31                             output_chrtout_grd_NWM, output_lsmOut_NWM, &
  32                             output_frxstPts, output_chanObs_NWM, output_gw_NWM
  33    use module_HYDRO_io, only: sub_output_gw, restart_out_nc, restart_in_nc,  &
  34         get_file_dimension , get_file_globalatts, get2d_lsm_real, get2d_lsm_vegtyp, get2d_lsm_soltyp, &
  35         output_lsm,  output_GW_Diag
  36    use module_HYDRO_io, only : output_lakes2
  37    use module_rt_data, only: rt_domain
  38    use module_GW_baseflow
  39    use module_gw_gw2d
  40    use module_gw_gw2d_data, only: gw2d
  41    use module_channel_routing, only: drive_channel, drive_channel_rsl
  42    use orchestrator_base
  43    use config_base, only: nlst, noah_lsm
  44    use module_routing, only: getChanDim, landrt_ini
  45    use module_HYDRO_utils
  46    use module_lsm_forcing, only: geth_newdate
  47 #ifdef WRF_HYDRO_NUDGING
  48    use module_stream_nudging,  only: init_stream_nudging
  49 #endif
  50    use module_hydro_stop, only: HYDRO_stop
  51    use module_UDMAP, only: get_basn_area_nhd
  52    use netcdf
  53
  54    implicit none
  55
  56
  57 #ifdef HYDRO_D
  58   real :: timeOr = 0
  59   real :: timeSr = 0
  60   real :: timeCr = 0
  61   real :: timeGW = 0
  62   integer :: clock_count_1 = 0
  63   integer :: clock_count_2 = 0
  64   integer :: clock_rate    = 0
  65 #endif
  66   integer :: rtout_factor  = 0
  67
  68   integer, parameter :: r4 = selected_real_kind(4)
  69   real,    parameter :: zeroFlt=0.0000000000000000000_r4
  70   integer, parameter :: r8 = selected_real_kind(8)
  71   real*8,  parameter :: zeroDbl=0.0000000000000000000_r8
  72
  73    contains
  74    subroutine HYDRO_rst_out(did)
  75 #ifdef WRF_HYDRO_NUDGING
  76       use module_stream_nudging, only: output_nudging_last_obs
  77 #endif
  78       implicit none
  79       integer:: rst_out
  80       integer did, outflag
  81       character(len=19) out_date
  82 #ifdef MPP_LAND
  83       character(len=19) str_tmp
  84 #endif
  85       rst_out = -99
  86 #ifdef MPP_LAND
  87    if(IO_id .eq. my_id) then
  88 #endif
  89      if(nlst(did)%dt .gt. nlst(did)%rst_dt*60) then
  90         call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%dt*rt_domain(did)%rst_counts))
  91      else
  92         call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%rst_dt*60*rt_domain(did)%rst_counts))
  93      endif
  94      if ( (nlst(did)%rst_dt .gt. 0) .and. (out_date(1:19) == nlst(did)%olddate(1:19)) ) then
  95            rst_out = 99
  96            rt_domain(did)%rst_counts = rt_domain(did)%rst_counts + 1
  97      endif
  98 ! restart every month automatically.
  99      if ( (nlst(did)%olddate(9:10) == "01") .and. (nlst(did)%olddate(12:13) == "00") .and. &
 100           (nlst(did)%olddate(15:16) == "00").and. (nlst(did)%olddate(18:19) == "00") .and. &
 101           (nlst(did)%rst_dt .le. 0)  ) then
 102            if(nlst(did)%startdate(1:16) .ne. nlst(did)%olddate(1:16) ) then
 103                rst_out = 99
 104            endif
 105      endif
 106
 107 #ifdef MPP_LAND
 108    endif
 109      call mpp_land_bcast_int1(rst_out)
 110 #endif
 111     if(rst_out .gt. 0) then
 112       write(6,*) "yw check output restart at ",nlst(did)%olddate(1:16)
 113 #ifdef MPP_LAND
 114       if(nlst(did)%rst_bi_out .eq. 1) then
 115              if(my_id .lt. 10) then
 116                   write(str_tmp,'(I1)') my_id
 117              else if(my_id .lt. 100) then
 118                   write(str_tmp,'(I2)') my_id
 119              else if(my_id .lt. 1000) then
 120                   write(str_tmp,'(I3)') my_id
 121              else if(my_id .lt. 10000) then
 122                   write(str_tmp,'(I4)') my_id
 123              else if(my_id .lt. 100000) then
 124                   write(str_tmp,'(I5)') my_id
 125              else
 126                 continue
 127              endif
 128              call mpp_land_bcast_char(16,nlst(did)%olddate(1:16))
 129              call   RESTART_OUT_bi(trim("HYDRO_RST."//nlst(did)%olddate(1:16)   &
 130                  //"_DOMAIN"//trim(nlst(did)%hgrid)//"."//trim(str_tmp)),  did)
 131       else
 132 #endif
 133              call   RESTART_OUT_nc(trim("HYDRO_RST."//nlst(did)%olddate(1:16)   &
 134                  //"_DOMAIN"//trim(nlst(did)%hgrid)),  did)
 135 #ifdef MPP_LAND
 136       endif
 137 #endif
 138
 139 #ifdef WRF_HYDRO_NUDGING
 140       call output_nudging_last_obs
 141 #endif
 142    endif
 143
 144
 145    end subroutine HYDRO_rst_out
 146
 147     subroutine HYDRO_out(did, rstflag)
 148
 149         implicit none
 150         integer did, outflag, rtflag, iniflag
 151         integer rstflag
 152         character(len=19) out_date
 153         integer :: Kt, ounit, i
 154         real, dimension(RT_DOMAIN(did)%NLINKS,2) :: str_out
 155         real, dimension(RT_DOMAIN(did)%NLINKS) :: vel_out
 156
 157         !    real, dimension(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx):: soilmx_tmp, &
 158             !           runoff1x_tmp, runoff2x_tmp, runoff3x_tmp,etax_tmp, &
 159             !           EDIRX_tmp,ECX_tmp,ETTX_tmp,RCX_tmp,HX_tmp,acrain_tmp, &
 160             !           ACSNOM_tmp, esnow2d_tmp, drip2d_tmp,dewfall_tmp, fpar_tmp, &
 161             !           qfx_tmp, prcp_out_tmp, etpndx_tmp
 162
 163         outflag = -99
 164
 165 #ifdef MPP_LAND
 166         if(IO_id .eq. my_id) then
 167 #endif
 168             if(nlst(did)%olddate(1:19) .eq. nlst(did)%startdate(1:19) .and. rt_domain(did)%his_out_counts .eq. 0) then
 169 #ifdef HYDRO_D
 170 #ifndef NCEP_WCOSS
 171                 write(6,*) "output hydrology at time : ",nlst(did)%olddate(1:19), rt_domain(did)%his_out_counts
 172 #else
 173                 write(78,*) "output hydrology at time : ",nlst(did)%olddate(1:19), rt_domain(did)%his_out_counts
 174 #endif
 175 #endif
 176                 outflag = 99
 177             else
 178                 if(nlst(did)%dt .gt. nlst(did)%out_dt*60) then
 179                     call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%dt*rt_domain(did)%out_counts))
 180                 else
 181                     call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%out_dt*60*rt_domain(did)%out_counts))
 182                 endif
 183                 if ( out_date(1:19) == nlst(did)%olddate(1:19) ) then
 184 #ifdef HYDRO_D
 185 #ifndef NCEP_WCOSS
 186                     write(6,*) "output hydrology at time : ",nlst(did)%olddate(1:19)
 187 #else
 188                     write(78,*) "output hydrology at time : ",nlst(did)%olddate(1:19)
 189 #endif
 190 #endif
 191                     outflag = 99
 192                 endif
 193             endif
 194 #ifdef MPP_LAND
 195         endif
 196         call mpp_land_bcast_int1(outflag)
 197 #endif
 198
 199         if (rstflag .eq. 1) call HYDRO_rst_out(did)
 200
 201         if (outflag .lt. 0) return
 202
 203         rt_domain(did)%out_counts = rt_domain(did)%out_counts + 1
 204         rt_domain(did)%his_out_counts = rt_domain(did)%his_out_counts + 1
 205
 206         if(nlst(did)%out_dt*60 .gt. nlst(did)%DT) then
 207             kt = rt_domain(did)%his_out_counts*nlst(did)%out_dt*60/nlst(did)%DT
 208         else
 209             kt = rt_domain(did)%his_out_counts
 210         endif
 211
 212         ! jump the ouput for the initial time when it has restart file from routing.
 213         rtflag = -99
 214         iniflag = -99
 215 #ifdef MPP_LAND
 216         if(IO_id .eq. my_id) then
 217 #endif
 218             !       if ( (trim(nlst_rt(did)%restart_file) /= "") .and. ( nlst_rt(did)%startdate(1:19) == nlst_rt(did)%olddate(1:19) ) ) then
 219             !#ifndef NCEP_WCOSS
 220             !             print*, "yyyywww restart_file = ", trim(nlst_rt(did)%restart_file)
 221             !#else
 222             !             write(78,*) "yyyywww restart_file = ", trim(nlst_rt(did)%restart_file)
 223             !#endif
 224             if ( nlst(did)%startdate(1:19) == nlst(did)%olddate(1:19) ) iniflag = 1
 225             if ( (trim(nlst(did)%restart_file) /= "") .and. ( nlst(did)%startdate(1:19) == nlst(did)%olddate(1:19) ) ) rtflag = 1
 226             !       endif
 227 #ifdef MPP_LAND
 228         endif
 229         call mpp_land_bcast_int1(rtflag)
 230         call mpp_land_bcast_int1(iniflag)
 231 #endif
 232
 233
 234         !yw keep the initial time otuput for debug
 235         if(rtflag == 1) then
 236             rt_domain(did)%restQSTRM = .false.   !!! do not reset QSTRM.. at initial time.
 237             if(nlst(did)%t0OutputFlag .eq. 0) return
 238         endif
 239
 240         if (iniflag == 1) then
 241             if(nlst(did)%t0OutputFlag .eq. 0) return
 242         endif
 243
 244       if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then
 245         if(nlst(did)%LSMOUT_DOMAIN .eq. 1)  then
 246             if(nlst(did)%io_form_outputs .eq. 0) then
 247                 call output_lsm(trim(nlst(did)%olddate(1:4)//nlst(did)%olddate(6:7)//nlst(did)%olddate(9:10)  &
 248                     //nlst(did)%olddate(12:13)//nlst(did)%olddate(15:16)//  &
 249                     ".LSMOUT_DOMAIN"//trim(nlst(did)%hgrid)),     &
 250                     did)
 251             else
 252                 call output_lsmOut_NWM(did)
 253             endif
 254         endif
 255       end if
 256
 257         if(nlst(did)%SUBRTSWCRT         .gt. 0 .or. &
 258                 nlst(did)%OVRTSWCRT          .gt. 0 .or. &
 259                 nlst(did)%GWBASESWCRT        .gt. 0 .or. &
 260    nlst(did)%CHANRTSWCRT        .gt. 0 .or. &
 261                 nlst(did)%channel_only       .gt. 0 .or. &
 262                 nlst(did)%channelBucket_only .gt. 0      ) then
 263
 264
 265             if(nlst(did)%RTOUT_DOMAIN       .eq. 1 .and. &
 266                     nlst(did)%channel_only       .eq. 0 .and. &
 267                     nlst(did)%channelBucket_only .eq. 0       ) then
 268                 if(mod(rtout_factor,3) .eq. 2 .or. &
 269                    nlst(did)%io_config_outputs .ne. 5 .and. &
 270                    nlst(did)%io_config_outputs .ne. 3) then
 271                     ! Output gridded routing variables on National Water Model
 272                     ! high-res routing grid
 273                     if(nlst(did)%io_form_outputs .ne. 0) then
 274                         call output_rt_NWM(did,nlst(did)%igrid)
 275                     else
 276                         call output_rt(    &
 277                             nlst(did)%igrid, nlst(did)%split_output_count, &
 278                             RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt, &
 279                             nlst(did)%nsoil, &
 280                             !                  nlst_rt(did)%startdate, nlst_rt(did)%olddate,
 281                         !                  rt_domain(did)%subsurface%state%qsubrt,&
 282                             nlst(did)%sincedate, nlst(did)%olddate, rt_domain(did)%subsurface%state%qsubrt,&
 283                             rt_domain(did)%subsurface%properties%zwattablrt,RT_DOMAIN(did)%subsurface%grid_transform%smcrt,&
 284                             RT_DOMAIN(did)%SUB_RESID,       &
 285                             RT_DOMAIN(did)%q_sfcflx_x,RT_DOMAIN(did)%q_sfcflx_y,&
 286                             rt_domain(did)%overland%properties%surface_slope_x,rt_domain(did)%overland%properties%surface_slope_y,&
 287                             RT_DOMAIN(did)%QSTRMVOLRT_ACC,rt_domain(did)%overland%control%surface_water_head_routing, &
 288                             nlst(did)%geo_finegrid_flnm,nlst(did)%DT,&
 289                             rt_domain(did)%subsurface%properties%sldpth,RT_DOMAIN(did)%LATVAL,&
 290                             RT_DOMAIN(did)%LONVAL,rt_domain(did)%overland%properties%distance_to_neighbor,nlst(did)%RTOUT_DOMAIN,&
 291                             rt_domain(did)%overland%control%boundary_flux, &
 292                             nlst(did)%io_config_outputs &
 293                             )
 294                     endif ! End check for io_form_outputs value
 295                 endif ! End check for rtout_factor
 296                 rtout_factor = rtout_factor + 1
 297             endif
 298             !! JLM disable GW output for NWM. Bring this line back when runtime output options avail.
 299             !! JLM This seems like a more logical place?
 300             if(nlst(did)%io_form_outputs .ne. 0) then
 301                 if(nlst(did)%GWBASESWCRT .ne. 0) then
 302                     if(nlst(did)%channel_only .eq. 0) then
 303                         if(nlst(did)%output_gw .ne. 0) then
 304                             call output_gw_NWM(did,nlst(did)%igrid)
 305                         endif
 306                     endif
 307                 endif
 308             else
 309                 if((nlst(did)%GWBASESWCRT .eq. 1) .or. (nlst(did)%GWBASESWCRT .ge. 4)) then
 310                     if(nlst(did)%channel_only .eq. 0) then
 311                         if(nlst(did)%output_gw  .eq. 1) call output_GW_Diag(did)
 312                     endif
 313                 end if
 314             endif
 315
 316
 317             if(nlst(did)%GWBASESWCRT .eq. 3) then
 318
 319                 if(nlst(did)%output_gw  .eq. 1)  &
 320                     call sub_output_gw(    &
 321                     nlst(did)%igrid, nlst(did)%split_output_count, &
 322                     RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt,          &
 323                     nlst(did)%nsoil,                               &
 324                     !               nlst(did)%startdate, nlst(did)%olddate,    &
 325                     nlst(did)%sincedate, nlst(did)%olddate,    &
 326                     gw2d(did)%h, rt_domain(did)%subsurface%grid_transform%smcrt,                   &
 327                     gw2d(did)%convgw, gw2d(did)%excess,                  &
 328                     gw2d(did)%qsgwrt, gw2d(did)%qgw_chanrt,              &
 329                     nlst(did)%geo_finegrid_flnm,nlst(did)%DT, &
 330                     rt_domain(did)%subsurface%properties%sldpth,RT_DOMAIN(did)%LATVAL,       &
 331                     RT_DOMAIN(did)%LONVAL,rt_domain(did)%overland%properties%distance_to_neighbor,           &
 332                     nlst(did)%output_gw)
 333
 334             endif
 335             ! BF end gw2d output section
 336
 337 #ifdef HYDRO_D
 338 #ifndef NCEP_WCOSS
 339             write(6,*) "before call output_chrt"
 340             call flush(6)
 341 #else
 342             write(78,*) "before call output_chrt"
 343 #endif
 344 #endif
 345
 346             if (nlst(did)%CHANRTSWCRT.eq.1.or.nlst(did)%CHANRTSWCRT.eq.2) then
 347
 348                 !ADCHANGE: Change values for within lake reaches to NA
 349                 str_out = RT_DOMAIN(did)%QLINK
 350                 vel_out = RT_DOMAIN(did)%velocity
 351
 352                 if (RT_DOMAIN(did)%NLAKES .gt. 0)  then
 353                     do i=1,RT_DOMAIN(did)%NLINKS
 354                         if (RT_DOMAIN(did)%TYPEL(i) .eq. 2) then
 355                             str_out(i,1) = -9.E15
 356                             vel_out(i) = -9.E15
 357                         endif
 358                     end do
 359                 endif
 360                 !ADCHANGE: End
 361
 362                 if(nlst(did)%io_form_outputs .ne. 0) then
 363                     ! Call National Water Model output routine for output on NHD forecast points.
 364                     if(nlst(did)%CHRTOUT_DOMAIN .ne. 0) then
 365                         call output_chrt_NWM(did)
 366                     endif
 367                     ! Call the subroutine to output frxstPts.
 368                     if(nlst(did)%frxst_pts_out .ne. 0) then
 369                         call output_frxstPts(did)
 370                     endif
 371                     ! Call the subroutine to output CHANOBS
 372                     if(nlst(did)%CHANOBS_DOMAIN .ne. 0) then
 373                         call output_chanObs_NWM(did)
 374                     endif
 375                 else
 376                     ! Call traditional output routines
 377                     !ADCHANGE: We suspect this routine is broken so default is now output_chrtout2
 378                     !             if(nlst_rt(did)%CHRTOUT_DOMAIN  .eq. 1)  then
 379                     !#ifdef MPP_LAND
 380                     !                 call mpp_output_chrt( &
 381                         !                      rt_domain(did)%gnlinks,rt_domain(did)%gnlinksl,rt_domain(did)%map_l2g, &
 382                         !#else
 383                     !                 call output_chrt(  &
 384                         !#endif
 385                     !                   nlst_rt(did)%igrid, nlst_rt(did)%split_output_count, &
 386                         !                   RT_DOMAIN(did)%NLINKS,RT_DOMAIN(did)%ORDER, &
 387                         !                   nlst_rt(did)%sincedate,nlst_rt(did)%olddate,RT_DOMAIN(did)%CHLON,&
 388                         !                   RT_DOMAIN(did)%CHLAT,                                      &
 389                         !                   RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ZELEV,                &
 390                         !                   !RT_DOMAIN(did)%QLINK,nlst_rt(did)%DT,Kt,                   &
 391                         !                   str_out, nlst_rt(did)%DT,Kt,                   &
 392                         !                   RT_DOMAIN(did)%STRMFRXSTPTS,nlst_rt(did)%order_to_write,   &
 393                         !                   RT_DOMAIN(did)%NLINKSL,nlst_rt(did)%channel_option,        &
 394                         !                   rt_domain(did)%gages, rt_domain(did)%gageMiss,             &
 395                         !                   nlst_rt(did)%dt                                            &
 396                         !#ifdef WRF_HYDRO_NUDGING
 397                     !                   , RT_DOMAIN(did)%nudge                                     &
 398                         !#endif
 399                     !                   , RT_DOMAIN(did)%accSfcLatRunoff, RT_DOMAIN(did)%accBucket &
 400                         !                   , RT_DOMAIN(did)%qSfcLatRunoff,   RT_DOMAIN(did)%qBucket   &
 401                         !                   , RT_DOMAIN(did)%qin_gwsubbas                              &
 402                         !                   , nlst_rt(did)%UDMP_OPT                                    &
 403                         !                   )
 404                     !              else
 405                     if(nlst(did)%CHRTOUT_DOMAIN  .gt. 0)  then
 406 #ifdef MPP_LAND
 407                         call mpp_output_chrt2(&
 408                             rt_domain(did)%gnlinks,rt_domain(did)%gnlinksl,rt_domain(did)%map_l2g, &
 409 #else
 410                             call output_chrt2(  &
 411 #endif
 412                             nlst(did)%igrid, nlst(did)%split_output_count, &
 413                             RT_DOMAIN(did)%NLINKS,RT_DOMAIN(did)%ORDER,          &
 414                             nlst(did)%sincedate,nlst(did)%olddate,         &
 415                             RT_DOMAIN(did)%CHLON, RT_DOMAIN(did)%CHLAT,          &
 416                             RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ZELEV,          &
 417                             !RT_DOMAIN(did)%QLINK,nlst_rt(did)%DT,Kt,             &
 418                             str_out, nlst(did)%DT,Kt,                         &
 419                             RT_DOMAIN(did)%NLINKSL,nlst(did)%channel_option,  &
 420                             rt_domain(did)%linkid                                &
 421 #ifdef WRF_HYDRO_NUDGING
 422                             , RT_DOMAIN(did)%nudge &
 423 #endif
 424                             !, RT_DOMAIN(did)%QLateral, nlst_rt(did)%io_config_outputs,
 425                         !RT_DOMAIN(did)%velocity &
 426                             , RT_DOMAIN(did)%QLateral, nlst(did)%io_config_outputs, vel_out &
 427                             , RT_DOMAIN(did)%accSfcLatRunoff, RT_DOMAIN(did)%accBucket &
 428                             , RT_DOMAIN(did)%qSfcLatRunoff,   RT_DOMAIN(did)%qBucket &
 429                             , RT_DOMAIN(did)%qin_gwsubbas,    nlst(did)%UDMP_OPT &
 430                             )
 431                     endif
 432
 433                 endif
 434             endif ! End of checking for io_form_outputs flag value
 435
 436 #ifdef MPP_LAND
 437             if(nlst(did)%CHRTOUT_GRID  .eq. 1)  then
 438                 if(nlst(did)%io_form_outputs .eq. 0) then
 439                     call mpp_output_chrtgrd(nlst(did)%igrid, nlst(did)%split_output_count, &
 440                         RT_DOMAIN(did)%ixrt,RT_DOMAIN(did)%jxrt, RT_DOMAIN(did)%NLINKS,   &
 441                         RT_DOMAIN(did)%GCH_NETLNK, &
 442                         nlst(did)%startdate, nlst(did)%olddate, &
 443                         !RT_DOMAIN(did)%qlink, nlst_rt(did)%dt, nlst_rt(did)%geo_finegrid_flnm,   &
 444                         str_out, nlst(did)%dt, nlst(did)%geo_finegrid_flnm,   &
 445                         RT_DOMAIN(did)%gnlinks,RT_DOMAIN(did)%map_l2g,                   &
 446                         RT_DOMAIN(did)%g_ixrt,RT_DOMAIN(did)%g_jxrt )
 447 #endif
 448                 else
 449                     call output_chrtout_grd_NWM(did,nlst(did)%igrid)
 450                 endif
 451             endif
 452             if (RT_DOMAIN(did)%NLAKES.gt.0)  then
 453                 if(nlst(did)%io_form_outputs .ne. 0) then
 454                     ! Output lakes in NWM format
 455                     if(nlst(did)%outlake .ne. 0) then
 456                         call output_lakes_NWM(did,nlst(did)%igrid)
 457                     endif
 458                 else
 459                     if(nlst(did)%outlake .eq. 1) then
 460 #ifdef MPP_LAND
 461                         call mpp_output_lakes( RT_DOMAIN(did)%lake_index, &
 462 #else
 463                             call output_lakes(  &
 464 #endif
 465                             nlst(did)%igrid, nlst(did)%split_output_count, &
 466                             RT_DOMAIN(did)%NLAKES, &
 467                             trim(nlst(did)%sincedate), trim(nlst(did)%olddate), &
 468                             RT_DOMAIN(did)%LATLAKE,RT_DOMAIN(did)%LONLAKE, &
 469                             RT_DOMAIN(did)%ELEVLAKE,RT_DOMAIN(did)%QLAKEI, &
 470                             RT_DOMAIN(did)%QLAKEO, &
 471                             RT_DOMAIN(did)%RESHT,nlst(did)%DT,Kt)
 472                     endif
 473                     if(nlst(did)%outlake .eq. 2) then
 474 #ifdef MPP_LAND
 475                         call mpp_output_lakes2( RT_DOMAIN(did)%lake_index, &
 476 #else
 477                             call output_lakes2(  &
 478 #endif
 479                             nlst(did)%igrid, nlst(did)%split_output_count, &
 480                             RT_DOMAIN(did)%NLAKES, &
 481                             trim(nlst(did)%sincedate), trim(nlst(did)%olddate), &
 482                             RT_DOMAIN(did)%LATLAKE,RT_DOMAIN(did)%LONLAKE, &
 483                             RT_DOMAIN(did)%ELEVLAKE,RT_DOMAIN(did)%QLAKEI, &
 484                             RT_DOMAIN(did)%QLAKEO, &
 485                             RT_DOMAIN(did)%RESHT,nlst(did)%DT,Kt,RT_DOMAIN(did)%LAKEIDM)
 486                     endif
 487
 488                 endif ! end of check for io_form_outputs value
 489             endif   ! end if block of rNLAKES .gt. 0
 490 #ifdef HYDRO_D
 491 #ifndef NCEP_WCOSS
 492             write(6,*) "end calling output functions"
 493 #else
 494             write(78,*) "end calling output functions"
 495 #endif
 496 #endif
 497
 498         endif  ! end of routing switch
 499
 500
 501     end subroutine HYDRO_out
 502
 503
 504       subroutine HYDRO_rst_in(did)
 505         integer :: did
 506         integer:: flag
 507
 508
 509
 510    flag = -1
 511 #ifdef MPP_LAND
 512    if(my_id.eq.IO_id) then
 513 #endif
 514       if (trim(nlst(did)%restart_file) /= "") then
 515           flag = 99
 516           rt_domain(did)%timestep_flag = 99   ! continue run
 517       endif
 518 #ifdef MPP_LAND
 519    endif
 520    call mpp_land_bcast_int1(flag)
 521 #endif
 522
 523    nlst(did)%sincedate = nlst(did)%startdate
 524
 525    if (flag.eq.99) then
 526
 527 #ifdef MPP_LAND
 528      if(my_id.eq.IO_id) then
 529 #endif
 530 #ifdef HYDRO_D
 531 #ifndef NCEP_WCOSS
 532         write(6,*) "*** read restart data: ",trim(nlst(did)%restart_file)
 533 #else
 534         write(78,*) "*** read restart data: ",trim(nlst(did)%restart_file)
 535 #endif
 536 #endif
 537 #ifdef MPP_LAND
 538      endif
 539 #endif
 540
 541 #ifdef MPP_LAND
 542       if(nlst(did)%rst_bi_in .eq. 1) then
 543          call RESTART_IN_bi(trim(nlst(did)%restart_file), did)
 544       else
 545 #endif
 546          call RESTART_IN_nc(trim(nlst(did)%restart_file), did)
 547 #ifdef MPP_LAND
 548       endif
 549 #endif
 550
 551 !yw  if (trim(nlst_rt(did)%restart_file) /= "") then
 552 !yw          nlst_rt(did)%restart_file = ""
 553 !yw  endif
 554
 555   endif
 556  end subroutine HYDRO_rst_in
 557
 558      subroutine HYDRO_time_adv(did)
 559         implicit none
 560         character(len = 19) :: newdate
 561         integer did
 562
 563 #ifdef MPP_LAND
 564    if(IO_id.eq.my_id) then
 565 #endif
 566          call geth_newdate(newdate, nlst(did)%olddate, nint( nlst(did)%dt))
 567          nlst(did)%olddate = newdate
 568 #ifdef HYDRO_D
 569 #ifndef NCEP_WCOSS
 570          write(6,*) "current time is ",newdate
 571 #else
 572          write(78,*) "current time is ",newdate
 573 #endif
 574 #endif
 575 #ifdef MPP_LAND
 576    endif
 577 #endif
 578      end subroutine HYDRO_time_adv
 579
 580      subroutine HYDRO_exe(did)
 581
 582
 583         implicit none
 584         integer:: did
 585         integer:: rst_out
 586
 587 !       call HYDRO_out(did)
 588
 589
 590 ! running land surface model
 591 ! cpl: 0--offline run;
 592 !      1-- coupling with WRF but running offline lsm;
 593 !      2-- coupling with WRF but do not run offline lsm
 594 !      3-- coupling with LIS and do not run offline lsm
 595 !      4:  coupling with CLM
 596 !          if(nlst_rt(did)%SYS_CPL .eq. 0 .or. nlst_rt(did)%SYS_CPL .eq. 1 )then
 597 !                  call drive_noahLSF(did,kt)
 598 !          else
 599 !              ! does not run the NOAH LASF model, only read the parameter
 600 !              call read_land_par(did,lsm(did)%ix,lsm(did)%jx)
 601 !          endif
 602
 603
 604 if (nlst(did)%GWBASESWCRT        .ne. 0 .or. &
 605     nlst(did)%SUBRTSWCRT         .ne. 0 .or. &
 606     nlst(did)%OVRTSWCRT          .ne. 0 .or. &
 607     nlst(did)%channel_only       .ne. 0 .or. &
 608     nlst(did)%channelBucket_only .ne. 0      ) then
 609
 610    ! step 1) disaggregate specific fields from LSM to Hydro grid
 611    if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then
 612
 613       RT_DOMAIN(did)%overland%streams_and_lakes%surface_water_to_channel = zeroFlt
 614       RT_DOMAIN(did)%LAKE_INFLORT_DUM = rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake
 615
 616       if(nlst(did)%SUBRTSWCRT .ne. 0 .or. nlst(did)%OVRTSWCRT .ne. 0) then
 617          call disaggregateDomain_drv(did)
 618       endif
 619       if(nlst(did)%OVRTSWCRT .eq. 0) then
 620          if(nlst(did)%UDMP_OPT .eq. 1) then
 621             call RunOffDisag(RT_DOMAIN(did)%INFXSRT, RT_DOMAIN(did)%landRunOff,        &
 622                  rt_domain(did)%dist_lsm(:,:,9),rt_domain(did)%overland%properties%distance_to_neighbor(:,:,9), &
 623                  RT_DOMAIN(did)%INFXSWGT, nlst(did)%AGGFACTRT, RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx)
 624          endif
 625       endif
 626
 627 #ifdef HYDRO_D
 628       call system_clock(count=clock_count_1, count_rate=clock_rate)
 629 #endif
 630    endif !! channel_only & channelBucket_only == 0
 631
 632
 633    ! step 2)
 634    if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then
 635       if(nlst(did)%SUBRTSWCRT .ne.0) then
 636          call SubsurfaceRouting_drv(did)
 637       endif
 638 #ifdef HYDRO_D
 639       call system_clock(count=clock_count_2, count_rate=clock_rate)
 640       timeSr = timeSr     + float(clock_count_2-clock_count_1)/float(clock_rate)
 641 #ifndef NCEP_WCOSS
 642       write(6,*) "Timing: Subsurface Routing  accumulated time--", timeSr
 643 #else
 644       write(78,*) "Timing: Subsurface Routing  accumulated time--", timeSr
 645 #endif
 646 #endif
 647    end if !! channel_only & channelBucket_only == 0
 648
 649
 650    ! step 3) todo split
 651    if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then
 652 #ifdef HYDRO_D
 653       call system_clock(count=clock_count_1, count_rate=clock_rate)
 654 #endif
 655       if(nlst(did)%OVRTSWCRT .ne. 0) then
 656          call OverlandRouting_drv(did)
 657       else
 658          !ADCHANGE: Updating landRunoff instead of surface_water_head_routing. This now allows
 659          !          landRunoff to include both infxsrt from the LSM and exfiltration (if subsfc
 660          !          is active) and prevents surface_water_head_routing from inadvertently being
 661          !          passed back to the LSM.
 662          if (nlst(did)%UDMP_OPT .eq. 1) then
 663            ! If subsurface is on, we update landRunOff to include the updated term w/ exfiltration.
 664            ! If subsurface is off, landRunOff does not change from original value so we leave as-is.
 665            if (nlst(did)%SUBRTSWCRT .ne. 0) then
 666              rt_domain(did)%landRunOff = rt_domain(did)%overland%control%infiltration_excess
 667            endif
 668          else
 669            ! If overland is off and subsurface is on, we need to update INFXSRT (LSM grid)
 670            ! since that is what gets fed through the buckets into the channels. So we aggregate
 671            ! the high-res infiltration_excess back to coarse grid.
 672            if (nlst(did)%SUBRTSWCRT .ne. 0) then
 673              call RunoffAggregate(rt_domain(did)%overland%control%infiltration_excess, &
 674                                   rt_domain(did)%INFXSRT, nlst(did)%AGGFACTRT, &
 675                                   rt_domain(did)%ix, rt_domain(did)%jx)
 676            endif
 677          endif
 678          ! In either case, if overland is off we need to zero-out surface_water_head since this
 679          ! water is being scraped into channel and should NOT be passed back to the LSM.
 680          rt_domain(did)%overland%control%infiltration_excess = 0.
 681          rt_domain(did)%overland%control%surface_water_head_routing = 0.
 682       endif
 683 #ifdef HYDRO_D
 684       call system_clock(count=clock_count_2, count_rate=clock_rate)
 685       timeOr = timeOr     + float(clock_count_2-clock_count_1)/float(clock_rate)
 686 #ifndef NCEP_WCOSS
 687       write(6,*) "Timing: Overland Routing  accumulated time--", timeOr
 688 #else
 689       write(78,*) "Timing: Overland Routing  accumulated time--", timeOr
 690 #endif
 691 #endif
 692
 693       RT_DOMAIN(did)%QSTRMVOLRT_TS = rt_domain(did)%overland%streams_and_lakes%surface_water_to_channel
 694       RT_DOMAIN(did)%QSTRMVOLRT_ACC = RT_DOMAIN(did)%QSTRMVOLRT_ACC + RT_DOMAIN(did)%QSTRMVOLRT_TS
 695
 696       RT_DOMAIN(did)%LAKE_INFLORT_TS = rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake-RT_DOMAIN(did)%LAKE_INFLORT_DUM
 697
 698 #ifdef HYDRO_D
 699       call system_clock(count=clock_count_1, count_rate=clock_rate)
 700 #endif
 701    end if !! channel_only & channelBucket_only == 0
 702
 703
 704    ! step 4) baseflow or groundwater physics
 705    !! channelBucket_only can be anything: the only time you dont run this is if channel_only=1
 706    if(nlst(did)%channel_only .eq. 0) then
 707       if (nlst(did)%GWBASESWCRT .gt. 0) then
 708          call driveGwBaseflow(did)
 709       endif
 710 #ifdef HYDRO_D
 711       call system_clock(count=clock_count_2, count_rate=clock_rate)
 712       timeGw = timeGw     + float(clock_count_2-clock_count_1)/float(clock_rate)
 713 #ifndef NCEP_WCOSS
 714       write(6,*) "Timing: GwBaseflow  accumulated time--", timeGw
 715 #else
 716       write(78,*) "Timing: GwBaseflow  accumulated time--", timeGw
 717 #endif
 718 #endif
 719 #ifdef HYDRO_D
 720       call system_clock(count=clock_count_1, count_rate=clock_rate)
 721 #endif
 722    end if !! channel_only == 0
 723
 724    ! step 5) river channel physics
 725    call driveChannelRouting(did)
 726 #ifdef HYDRO_D
 727    call system_clock(count=clock_count_2, count_rate=clock_rate)
 728    timeCr = timeCr     + float(clock_count_2-clock_count_1)/float(clock_rate)
 729 #ifndef NCEP_WCOSS
 730    write(6,*) "Timing: Channel Routing  accumulated time--", timeCr
 731 #else
 732    write(78,*) "Timing: Channel Routing  accumulated time--", timeCr
 733 #endif
 734 #endif
 735
 736    !! if not channel_only
 737    if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then
 738
 739       ! step 6) aggregate specific fields from Hydro to LSM grid
 740       if (nlst(did)%SUBRTSWCRT .ne.0  .or. nlst(did)%OVRTSWCRT .ne. 0 ) then
 741          call aggregateDomain(did)
 742       endif
 743
 744    end if !! channel_only & channelBucket_only == 0
 745
 746 end if
 747
 748
 749 !yw  if (nlst_rt(did)%sys_cpl .eq. 2) then
 750       ! advance to next time step
 751 !          call HYDRO_time_adv(did)
 752       ! output for history
 753 !          call HYDRO_out(did)
 754 !yw  endif
 755            call HYDRO_time_adv(did)
 756            call HYDRO_out(did, 1)
 757
 758
 759 !           write(90 + my_id,*) "finish calling hydro_exe"
 760 !           call flush(90+my_id)
 761 !          call mpp_land_sync()
 762
 763
 764
 765            !! Under channel-only, these variables are not allocated
 766            if(allocated(RT_DOMAIN(did)%SOLDRAIN)) RT_DOMAIN(did)%SOLDRAIN = 0
 767            if(allocated(rt_domain(did)%subsurface%state%qsubrt))   RT_DOMAIN(did)%subsurface%state%qsubrt   = 0
 768
 769
 770
 771       end subroutine HYDRO_exe
 772
 773
 774
 775 !----------------------------------------------------
 776 subroutine driveGwBaseflow(did)
 777
 778     implicit none
 779     integer, intent(in) :: did
 780
 781     integer :: i, jj, ii
 782
 783     !------------------------------------------------------------------
 784     !DJG Begin GW/Baseflow Routines
 785     !-------------------------------------------------------------------
 786
 787     if (nlst(did)%GWBASESWCRT.ge.1) then     ! Switch to activate/specify GW/Baseflow
 788
 789         !  IF (nlst(did)%GWBASESWCRT.GE.1000) THEN     ! Switch to activate/specify GW/Baseflow
 790
 791         if (nlst(did)%GWBASESWCRT.eq.1 .or. nlst(did)%GWBASESWCRT.eq.2 .or. nlst(did)%GWBASESWCRT.ge.4) then   ! Call simple bucket baseflow scheme
 792
 793 #ifdef HYDRO_D
 794 #ifndef NCEP_WCOSS
 795             write(6,*) "*****yw******start simp_gw_buck "
 796 #else
 797             write(78,*) "*****yw******start simp_gw_buck "
 798 #endif
 799 #endif
 800
 801             if(nlst(did)%UDMP_OPT .eq. 1) then
 802                 call simp_gw_buck_nhd(                                             &
 803                     RT_DOMAIN(did)%ix,            RT_DOMAIN(did)%jx,              &
 804                     RT_DOMAIN(did)%ixrt,          RT_DOMAIN(did)%jxrt,            &
 805                     RT_DOMAIN(did)%numbasns,      nlst(did)%AGGFACTRT,         &
 806                     nlst(did)%DT,              RT_DOMAIN(did)%INFXSWGT,        &
 807                     RT_DOMAIN(did)%INFXSRT,       RT_DOMAIN(did)%SOLDRAIN,        &
 808                     rt_domain(did)%overland%properties%distance_to_neighbor(:,:,9),   rt_domain(did)%dist_lsm(:,:,9), &
 809                     RT_DOMAIN(did)%gw_buck_coeff, RT_DOMAIN(did)%gw_buck_exp,     &
 810               RT_DOMAIN(did)%gw_buck_loss,                                  &
 811                     RT_DOMAIN(did)%z_max,         RT_DOMAIN(did)%z_gwsubbas,      &
 812                     RT_DOMAIN(did)%qout_gwsubbas, RT_DOMAIN(did)%qin_gwsubbas,    &
 813               RT_DOMAIN(did)%qloss_gwsubbas,                                &
 814                     nlst(did)%GWBASESWCRT,     nlst(did)%OVRTSWCRT,         &
 815 #ifdef MPP_LAND
 816                 RT_DOMAIN(did)%LNLINKSL,                                      &
 817 #else
 818                     RT_DOMAIN(did)%numbasns,                                      &
 819 #endif
 820                     rt_domain(did)%basns_area,                                    &
 821               rt_domain(did)%nhdBuckMask,   nlst(did)%bucket_loss,       &
 822               nlst(did)%channelBucket_only )
 823
 824             else
 825                 call simp_gw_buck(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,RT_DOMAIN(did)%ixrt,&
 826                     RT_DOMAIN(did)%jxrt,RT_DOMAIN(did)%numbasns,RT_DOMAIN(did)%gnumbasns,&
 827                 RT_DOMAIN(did)%basns_area,&
 828                     RT_DOMAIN(did)%basnsInd, RT_DOMAIN(did)%gw_strm_msk_lind,             &
 829                     RT_DOMAIN(did)%gwsubbasmsk, RT_DOMAIN(did)%INFXSRT, &
 830                     RT_DOMAIN(did)%SOLDRAIN, &
 831                     RT_DOMAIN(did)%z_gwsubbas,&
 832                     RT_DOMAIN(did)%qin_gwsubbas,RT_DOMAIN(did)%qout_gwsubbas,&
 833                     RT_DOMAIN(did)%qinflowbase,&
 834                     RT_DOMAIN(did)%gw_strm_msk,RT_DOMAIN(did)%gwbas_pix_ct, &
 835                     rt_domain(did)%overland%properties%distance_to_neighbor,nlst(did)%DT,&
 836                     RT_DOMAIN(did)%gw_buck_coeff,RT_DOMAIN(did)%gw_buck_exp, &
 837                     RT_DOMAIN(did)%z_max,&
 838                     nlst(did)%GWBASESWCRT,nlst(did)%OVRTSWCRT)
 839             endif
 840
 841             !! JLM: There's *perhaps* a better location for this output above.
 842             !! If above is better, remove this when runtime output options are avail.
 843             !#ifndef HYDRO_REALTIME
 844             !        call output_GW_Diag(did)
 845             !#endif
 846
 847 #ifdef HYDRO_D
 848 #ifndef NCEP_WCOSS
 849             write(6,*) "*****yw******end simp_gw_buck "
 850 #else
 851             write(78,*) "*****yw******end simp_gw_buck "
 852 #endif
 853 #endif
 854
 855             !!!For parameter setup runs output the percolation for each basin,
 856             !!!otherwise comment out this output...
 857         else if (nlst(did)%gwBaseSwCRT .eq. 3) then
 858
 859 #ifdef HYDRO_D
 860 #ifndef NCEP_WCOSS
 861             write(6,*) "*****bf******start 2d_gw_model "
 862 #else
 863             write(78,*) "*****bf******start 2d_gw_model "
 864 #endif
 865 #endif
 866
 867             !      compute qsgwrt  between lsm and gw with namelist selected coupling method
 868             !      qsgwrt is defined on the routing grid  and needs to be aggregated for SFLX
 869             if (nlst(did)%gwsoilcpl .GT. 0) THEN
 870
 871                 call gwSoilFlux(did)
 872
 873             end if
 874
 875             gw2d(did)%excess = 0.
 876
 877             call gwstep(gw2d(did)%ix, gw2d(did)%jx, gw2d(did)%dx, &
 878                 gw2d(did)%ltype, gw2d(did)%elev, gw2d(did)%bot, &
 879                 gw2d(did)%hycond, gw2d(did)%poros, gw2d(did)%compres, &
 880                 gw2d(did)%ho, gw2d(did)%h, gw2d(did)%convgw, &
 881                 gw2d(did)%excess, &
 882                 gw2d(did)%ebot, gw2d(did)%eocn, gw2d(did)%dt, &
 883                 gw2d(did)%istep)
 884
 885             gw2d(did)%ho = gw2d(did)%h
 886
 887             ! put surface exceeding groundwater to surface routing inflow
 888             rt_domain(did)%overland%control%surface_water_head_routing = rt_domain(did)%overland%control%surface_water_head_routing &
 889                 + gw2d(did)%excess*1000. ! convert to mm
 890
 891             ! aggregate  qsgw from routing to lsm grid
 892             call aggregateQsgw(did)
 893
 894             gw2d(did)%istep =  gw2d(did)%istep + 1
 895
 896 #ifdef HYDRO_D
 897 #ifndef NCEP_WCOSS
 898             write(6,*) "*****bf******end 2d_gw_model "
 899 #else
 900             write(78,*) "*****bf******end 2d_gw_model "
 901 #endif
 902 #endif
 903
 904         end if
 905
 906     end if    !DJG (End if for RTE SWC activation)
 907     !------------------------------------------------------------------
 908     !DJG End GW/Baseflow Routines
 909     !-------------------------------------------------------------------
 910 end subroutine driveGwBaseflow
 911
 912
 913
 914
 915 !-------------------------------------------
 916       subroutine driveChannelRouting(did)
 917
 918        implicit none
 919        integer, intent(in) :: did
 920
 921 !-------------------------------------------------------------------
 922 !-------------------------------------------------------------------
 923 !DJG,DNY  Begin Channel and Lake Routing Routines
 924 !-------------------------------------------------------------------
 925
 926 if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT.eq.2) then
 927
 928  if(nlst(did)%UDMP_OPT .eq. 1) then
 929      !!! for user defined Reach based Routing method.
 930
 931     call drive_CHANNEL_RSL(did, nlst(did)%UDMP_OPT,RT_DOMAIN(did)%timestep_flag, RT_DOMAIN(did)%IXRT,RT_DOMAIN(did)%JXRT,  &
 932        RT_DOMAIN(did)%LAKE_INFLORT_TS, RT_DOMAIN(did)%QSTRMVOLRT_TS, RT_DOMAIN(did)%TO_NODE, RT_DOMAIN(did)%FROM_NODE, &
 933        RT_DOMAIN(did)%TYPEL, RT_DOMAIN(did)%ORDER, RT_DOMAIN(did)%MAXORDER,   RT_DOMAIN(did)%CH_LNKRT, &
 934        rt_domain(did)%overland%streams_and_lakes%lake_mask, nlst(did)%DT, nlst(did)%DTCT, nlst(did)%DTRT_CH, &
 935        RT_DOMAIN(did)%MUSK, RT_DOMAIN(did)%MUSX, RT_DOMAIN(did)%QLINK, &
 936        RT_DOMAIN(did)%CHANLEN, RT_DOMAIN(did)%MannN, RT_DOMAIN(did)%So, RT_DOMAIN(did)%ChSSlp,RT_DOMAIN(did)%Bw, &
 937        RT_DOMAIN(did)%Tw,  RT_DOMAIN(did)%Tw_CC, RT_DOMAIN(did)%n_CC, &
 938        RT_DOMAIN(did)%ChannK,&
 939        RT_DOMAIN(did)%RESHT, &
 940        RT_DOMAIN(did)%CVOL, RT_DOMAIN(did)%QLAKEI, &
 941        RT_DOMAIN(did)%QLAKEO, RT_DOMAIN(did)%LAKENODE, &
 942        RT_DOMAIN(did)%QINFLOWBASE, RT_DOMAIN(did)%CHANXI, RT_DOMAIN(did)%CHANYJ, nlst(did)%channel_option,  &
 943        RT_DOMAIN(did)%nlinks, RT_DOMAIN(did)%NLINKSL, RT_DOMAIN(did)%LINKID, RT_DOMAIN(did)%node_area,         &
 944        RT_DOMAIN(did)%qout_gwsubbas, &
 945        RT_DOMAIN(did)%LAKEIDA, RT_DOMAIN(did)%LAKEIDM, RT_DOMAIN(did)%NLAKES, RT_DOMAIN(did)%LAKEIDX   &
 946 #ifdef MPP_LAND
 947        , RT_DOMAIN(did)%nlinks_index, RT_DOMAIN(did)%mpp_nlinks, RT_DOMAIN(did)%yw_mpp_nlinks  &
 948        , RT_DOMAIN(did)%LNLINKSL   &
 949        , RT_DOMAIN(did)%gtoNode, RT_DOMAIN(did)%toNodeInd, RT_DOMAIN(did)%nToInd      &
 950 #endif
 951        , RT_DOMAIN(did)%CH_LNKRT_SL, RT_DOMAIN(did)%landRunOff   &
 952 #ifdef WRF_HYDRO_NUDGING
 953        , RT_DOMAIN(did)%nudge &
 954 #endif
 955        , rt_domain(did)%accSfcLatRunoff,    rt_domain(did)%accBucket &
 956        , rt_domain(did)%qSfcLatRunoff,        rt_domain(did)%qBucket &
 957        , rt_domain(did)%QLateral,            rt_domain(did)%velocity &
 958        ,                                     rt_domain(did)%qloss    &
 959        ,                                     RT_DOMAIN(did)%HLINK    &
 960        , rt_domain(did)%nlinksize,            nlst(did)%OVRTSWCRT &
 961        ,                                     nlst(did)%SUBRTSWCRT &
 962        , nlst(did)%channel_only , nlst(did)%channelBucket_only &
 963        , nlst(did)%channel_bypass )
 964
 965 else
 966
 967     call drive_CHANNEL(did, RT_DOMAIN(did)%latval,RT_DOMAIN(did)%lonval, &
 968        RT_DOMAIN(did)%timestep_flag,RT_DOMAIN(did)%IXRT,RT_DOMAIN(did)%JXRT, &
 969        nlst(did)%SUBRTSWCRT, rt_domain(did)%subsurface%state%qsubrt, &
 970        RT_DOMAIN(did)%LAKE_INFLORT_TS, RT_DOMAIN(did)%QSTRMVOLRT_TS,&
 971        RT_DOMAIN(did)%TO_NODE, RT_DOMAIN(did)%FROM_NODE, RT_DOMAIN(did)%TYPEL,&
 972        RT_DOMAIN(did)%ORDER, RT_DOMAIN(did)%MAXORDER, RT_DOMAIN(did)%NLINKS,&
 973        RT_DOMAIN(did)%CH_NETLNK, rt_domain(did)%overland%streams_and_lakes%ch_netrt,RT_DOMAIN(did)%CH_LNKRT,&
 974        rt_domain(did)%overland%streams_and_lakes%lake_mask, nlst(did)%DT, nlst(did)%DTCT, nlst(did)%DTRT_CH,&
 975        RT_DOMAIN(did)%MUSK, RT_DOMAIN(did)%MUSX,  RT_DOMAIN(did)%QLINK, &
 976        RT_DOMAIN(did)%QLateral, &
 977        RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ELRT,RT_DOMAIN(did)%CHANLEN,&
 978        RT_DOMAIN(did)%MannN,RT_DOMAIN(did)%So, RT_DOMAIN(did)%ChSSlp, &
 979        RT_DOMAIN(did)%Bw,RT_DOMAIN(did)%Tw,RT_DOMAIN(did)%Tw_CC, RT_DOMAIN(did)%n_CC, &
 980        RT_DOMAIN(did)%ChannK,&
 981        RT_DOMAIN(did)%RESHT, &
 982        RT_DOMAIN(did)%ZELEV, RT_DOMAIN(did)%CVOL, &
 983        RT_DOMAIN(did)%NLAKES, RT_DOMAIN(did)%QLAKEI, RT_DOMAIN(did)%QLAKEO,&
 984        RT_DOMAIN(did)%LAKENODE, rt_domain(did)%overland%properties%distance_to_neighbor, &
 985        RT_DOMAIN(did)%QINFLOWBASE, RT_DOMAIN(did)%CHANXI, &
 986        RT_DOMAIN(did)%CHANYJ, nlst(did)%channel_option, &
 987        RT_DOMAIN(did)%RETDEP_CHAN, RT_DOMAIN(did)%NLINKSL, RT_DOMAIN(did)%LINKID, &
 988        RT_DOMAIN(did)%node_area  &
 989 #ifdef MPP_LAND
 990        ,RT_DOMAIN(did)%lake_index,RT_DOMAIN(did)%link_location,&
 991        RT_DOMAIN(did)%mpp_nlinks,RT_DOMAIN(did)%nlinks_index, &
 992        RT_DOMAIN(did)%yw_mpp_nlinks  &
 993        , RT_DOMAIN(did)%LNLINKSL,RT_DOMAIN(did)%LLINKID &
 994        , rt_domain(did)%gtoNode,rt_domain(did)%toNodeInd,rt_domain(did)%nToInd  &
 995 #endif
 996        , rt_domain(did)%CH_LNKRT_SL   &
 997        ,nlst(did)%GwBaseSwCRT, gw2d(did)%ho, gw2d(did)%qgw_chanrt, &
 998        nlst(did)%gwChanCondSw, nlst(did)%gwChanCondConstIn, &
 999        nlst(did)%gwChanCondConstOut, rt_domain(did)%velocity, rt_domain(did)%qloss &
1000        )
1001 endif
1002
1003     if((nlst(did)%gwBaseSwCRT == 3) .and. (nlst(did)%gwChanCondSw .eq. 1)) then
1004
1005            ! add/rm channel-aquifer exchange contribution
1006
1007            gw2d(did)%ho =  gw2d(did)%ho  &
1008                         +(((gw2d(did)%qgw_chanrt*(-1)) * gw2d(did)%dt / gw2d(did)%dx**2) &
1009                         /  gw2d(did)%poros)
1010
1011     endif
1012   endif
1013
1014 #ifdef HYDRO_D
1015 #ifndef NCEP_WCOSS
1016            write(6,*) "*****yw******end drive_CHANNEL "
1017 #else
1018            write(78,*) "*****yw******end drive_CHANNEL "
1019 #endif
1020 #endif
1021
1022       end subroutine driveChannelRouting
1023
1024
1025
1026 !------------------------------------------------
1027       subroutine aggregateDomain(did)
1028
1029 #ifdef MPP_LAND
1030         use module_mpp_land, only:  sum_real1, my_id, io_id, numprocs
1031 #endif
1032
1033        implicit none
1034        integer, intent(in) :: did
1035
1036        integer :: i, j, krt, ixxrt, jyyrt, &
1037                   AGGFACYRT, AGGFACXRT
1038
1039 #ifdef HYDRO_D
1040 ! ADCHANGE: Water balance variables
1041        integer :: kk
1042        real    :: smcrttot1,smctot2,sicetot2
1043        real    :: suminfxsrt1,suminfxs2
1044 #endif
1045
1046 #ifdef HYDRO_D
1047 #ifndef NCEP_WCOSS
1048         print *, "Beginning Aggregation..."
1049 #else
1050        write(78,*) "Beginning Aggregation..."
1051 #endif
1052 #endif
1053
1054 #ifdef HYDRO_D
1055 ! ADCHANGE: START Initial water balance variables
1056 ! ALL VARS in MM
1057         suminfxsrt1 = 0.
1058         smcrttot1 = 0.
1059         do i=1,RT_DOMAIN(did)%IXRT
1060          do j=1,RT_DOMAIN(did)%JXRT
1061             suminfxsrt1 = suminfxsrt1 + rt_domain(did)%overland%control%surface_water_head_routing(I,J) &
1062                                / float(RT_DOMAIN(did)%IXRT * RT_DOMAIN(did)%JXRT)
1063             do kk=1,nlst(did)%NSOIL
1064                 smcrttot1 = smcrttot1 + rt_domain(did)%subsurface%grid_transform%smcrt(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. &
1065                                / float(RT_DOMAIN(did)%IXRT * RT_DOMAIN(did)%JXRT)
1066             end do
1067          end do
1068         end do
1069 #ifdef MPP_LAND
1070 ! not tested
1071         CALL sum_real1(suminfxsrt1)
1072         CALL sum_real1(smcrttot1)
1073         suminfxsrt1 = suminfxsrt1/float(numprocs)
1074         smcrttot1 = smcrttot1/float(numprocs)
1075 #endif
1076 ! END Initial water balance variables
1077 #endif
1078
1079         do J=1,RT_DOMAIN(did)%JX
1080           do I=1,RT_DOMAIN(did)%IX
1081
1082              RT_DOMAIN(did)%SFCHEADAGGRT = 0.
1083 !DJG Subgrid weighting edit...
1084              RT_DOMAIN(did)%LSMVOL=0.
1085              do KRT=1,nlst(did)%NSOIL
1086 !                SMCAGGRT(KRT) = 0.
1087                RT_DOMAIN(did)%SH2OAGGRT(KRT) = 0.
1088              end do
1089
1090
1091              do AGGFACYRT=nlst(did)%AGGFACTRT-1,0,-1
1092               do AGGFACXRT=nlst(did)%AGGFACTRT-1,0,-1
1093
1094
1095                 IXXRT=I*nlst(did)%AGGFACTRT-AGGFACXRT
1096                 JYYRT=J*nlst(did)%AGGFACTRT-AGGFACYRT
1097 #ifdef MPP_LAND
1098        if(left_id.ge.0) IXXRT=IXXRT+1
1099        if(down_id.ge.0) JYYRT=JYYRT+1
1100 #else
1101 !yw ????
1102 !       IXXRT=IXXRT+1
1103 !       JYYRT=JYYRT+1
1104 #endif
1105
1106 !State Variables
1107                 RT_DOMAIN(did)%SFCHEADAGGRT = RT_DOMAIN(did)%SFCHEADAGGRT &
1108                                             + rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT)
1109 !DJG Subgrid weighting edit...
1110                 RT_DOMAIN(did)%LSMVOL = RT_DOMAIN(did)%LSMVOL &
1111                                       + rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) &
1112                                       * rt_domain(did)%overland%properties%distance_to_neighbor(IXXRT,JYYRT,9)
1113
1114                 do KRT=1,nlst(did)%NSOIL
1115 !DJG               SMCAGGRT(KRT)=SMCAGGRT(KRT)+SMCRT(IXXRT,JYYRT,KRT)
1116                    RT_DOMAIN(did)%SH2OAGGRT(KRT) = RT_DOMAIN(did)%SH2OAGGRT(KRT) &
1117                                                  + rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT)
1118                 end do
1119
1120               end do
1121              end do
1122
1123
1124
1125             rt_domain(did)%overland%control%surface_water_head_lsm(I,J) = RT_DOMAIN(did)%SFCHEADAGGRT &
1126                                           / (nlst(did)%AGGFACTRT**2)
1127
1128             do KRT=1,nlst(did)%NSOIL
1129 !DJG              SMC(I,J,KRT)=SMCAGGRT(KRT)/(AGGFACTRT**2)
1130                RT_DOMAIN(did)%SH2OX(I,J,KRT) = RT_DOMAIN(did)%SH2OAGGRT(KRT) &
1131                                              / (nlst(did)%AGGFACTRT**2)
1132             end do
1133
1134
1135
1136 !DJG Calculate subgrid weighting array...
1137
1138               do AGGFACYRT=nlst(did)%AGGFACTRT-1,0,-1
1139                 do AGGFACXRT=nlst(did)%AGGFACTRT-1,0,-1
1140                   IXXRT=I*nlst(did)%AGGFACTRT-AGGFACXRT
1141                   JYYRT=J*nlst(did)%AGGFACTRT-AGGFACYRT
1142 #ifdef MPP_LAND
1143        if(left_id.ge.0) IXXRT=IXXRT+1
1144        if(down_id.ge.0) JYYRT=JYYRT+1
1145 #else
1146 !yw ???
1147 !       IXXRT=IXXRT+1
1148 !       JYYRT=JYYRT+1
1149 #endif
1150                   if (RT_DOMAIN(did)%LSMVOL.gt.0.) then
1151                     RT_DOMAIN(did)%INFXSWGT(IXXRT,JYYRT) &
1152                                           = rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) &
1153                                           * rt_domain(did)%overland%properties%distance_to_neighbor(IXXRT,JYYRT,9) &
1154                                           / RT_DOMAIN(did)%LSMVOL
1155                   else
1156                     RT_DOMAIN(did)%INFXSWGT(IXXRT,JYYRT) &
1157                                           = 1./FLOAT(nlst(did)%AGGFACTRT**2)
1158                   end if
1159
1160                   do KRT=1,nlst(did)%NSOIL
1161
1162 !!!yw added for debug
1163                    if(rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) .lt. 0) then
1164 #ifndef NCEP_WCOSS
1165                       print*, "Error negative SMCRT", rt_domain(did)%SH2OWGT(IXXRT,JYYRT,KRT), RT_DOMAIN(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),RT_DOMAIN(did)%SH2OX(I,J,KRT)
1166 #else
1167                       write(78,*) "WARNING: negative SMCRT", rt_domain(did)%SH2OWGT(IXXRT,JYYRT,KRT), RT_DOMAIN(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),RT_DOMAIN(did)%SH2OX(I,J,KRT)
1168 #endif
1169                    endif
1170                    if(RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) .lt. 0) then
1171 #ifndef NCEP_WCOSS
1172                       print *, "Error negative SH2OWGT", rt_domain(did)%SH2OWGT(IXXRT,JYYRT,KRT), RT_DOMAIN(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),RT_DOMAIN(did)%SH2OX(I,J,KRT)
1173 #else
1174                       write(78,*) "WARNING: negative SH2OWGT", rt_domain(did)%SH2OWGT(IXXRT,JYYRT,KRT), RT_DOMAIN(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),RT_DOMAIN(did)%SH2OX(I,J,KRT)
1175 #endif
1176                    endif
1177
1178                     IF ( (rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) - &
1179                           rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT)) .GT. 0.000001 ) THEN
1180 #ifdef HYDRO_D
1181 #ifndef NCEP_WCOSS
1182                       print *, "SMCMAX exceeded upon aggregation...", &
1183                            rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),  &
1184                            rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT)
1185 #else
1186                      write(78,*) "FATAL ERROR: SMCMAX exceeded upon aggregation...", &
1187                            rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),  &
1188                            rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT)
1189 #endif
1190 #endif
1191                       call hydro_stop("In module_HYDRO_drv.F aggregateDomain() - "// &
1192                                       "SMCMAX exceeded upon aggregation.")
1193                     END IF
1194                     IF(RT_DOMAIN(did)%SH2OX(I,J,KRT).LT.0.) THEN
1195 #ifdef HYDRO_D
1196 #ifndef NCEP_WCOSS
1197                       print *, "Erroneous value of SH2O...", &
1198                                 RT_DOMAIN(did)%SH2OX(I,J,KRT),I,J,KRT
1199                       print *, "Error negative SH2OX", rt_domain(did)%SH2OWGT(IXXRT,JYYRT,KRT), RT_DOMAIN(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),RT_DOMAIN(did)%SH2OX(I,J,KRT)
1200 #else
1201                      write(78,*) "Erroneous value of SH2O...", &
1202                                 RT_DOMAIN(did)%SH2OX(I,J,KRT),I,J,KRT
1203                       write(78,*) "FATAL ERROR: negative SH2OX", rt_domain(did)%SH2OWGT(IXXRT,JYYRT,KRT), RT_DOMAIN(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),RT_DOMAIN(did)%SH2OX(I,J,KRT)
1204 #endif
1205 #endif
1206                       call hydro_stop("In module_HYDRO_drv.F aggregateDomain() "// &
1207                                       "- Error negative SH2OX")
1208                     END IF
1209
1210                     IF ( RT_DOMAIN(did)%SH2OX(I,J,KRT) .gt. 0 ) THEN
1211                         RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) &
1212                                  = rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) &
1213                                  / RT_DOMAIN(did)%SH2OX(I,J,KRT)
1214                     ELSE
1215 #ifdef HYDRO_D
1216                          print *, "Error zero SH2OX", rt_domain(did)%SH2OWGT(IXXRT,JYYRT,KRT), RT_DOMAIN(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT),RT_DOMAIN(did)%SH2OX(I,J,KRT)
1217 #endif
1218                          RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) = 0.0
1219                     ENDIF
1220 !?yw
1221                     RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) = max(1.0E-05, RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT))
1222                   end do
1223
1224                 end do
1225               end do
1226
1227          end do
1228         end do
1229
1230
1231 #ifdef MPP_LAND
1232         call MPP_LAND_COM_REAL(RT_DOMAIN(did)%INFXSWGT, &
1233                                RT_DOMAIN(did)%IXRT,    &
1234                                RT_DOMAIN(did)%JXRT, 99)
1235
1236         do i = 1, nlst(did)%NSOIL
1237            call MPP_LAND_COM_REAL(RT_DOMAIN(did)%SH2OWGT(:,:,i), &
1238                                   RT_DOMAIN(did)%IXRT, &
1239                                   RT_DOMAIN(did)%JXRT, 99)
1240         end do
1241 #endif
1242
1243 !DJG Update SMC with SICE (unchanged) and new value of SH2O from routing...
1244         RT_DOMAIN(did)%SMC = RT_DOMAIN(did)%SH2OX + RT_DOMAIN(did)%SICE
1245
1246 #ifdef HYDRO_D
1247 ! ADCHANGE: START Final water balance variables
1248 ! ALL VARS in MM
1249         suminfxs2 = 0.
1250         smctot2 = 0.
1251         sicetot2 = 0.
1252         do i=1,RT_DOMAIN(did)%IX
1253          do j=1,RT_DOMAIN(did)%JX
1254             suminfxs2 = suminfxs2 + rt_domain(did)%overland%control%surface_water_head_lsm(I,J) &
1255                                / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX)
1256             do kk=1,nlst(did)%NSOIL
1257                 smctot2 = smctot2 + rt_domain(did)%SMC(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. &
1258                                / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX)
1259                sicetot2 = sicetot2 + rt_domain(did)%SICE(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. &
1260                                 / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX)
1261             end do
1262          end do
1263         end do
1264
1265 #ifdef MPP_LAND
1266 ! not tested
1267         CALL sum_real1(suminfxs2)
1268         CALL sum_real1(smctot2)
1269         CALL sum_real1(sicetot2)
1270         suminfxs2 = suminfxs2/float(numprocs)
1271         smctot2 = smctot2/float(numprocs)
1272         sicetot2 = sicetot2/float(numprocs)
1273 #endif
1274
1275 #ifdef MPP_LAND
1276        if (my_id .eq. IO_id) then
1277 #endif
1278          print *, "Agg Mass Bal: "
1279          print *, "WB_AGG!InfxsDiff", suminfxs2-suminfxsrt1
1280          print *, "WB_AGG!Infxs1", suminfxsrt1
1281          print *, "WB_AGG!Infxs2", suminfxs2
1282          print *, "WB_AGG!SMCDiff", smctot2-smcrttot1-sicetot2
1283          print *, "WB_AGG!SMC1", smcrttot1
1284          print *, "WB_AGG!SMC2", smctot2
1285          print *, "WB_AGG!SICE2", sicetot2
1286          print *, "WB_AGG!Residual", (suminfxs2-suminfxsrt1) + &
1287                          (smctot2-smcrttot1-sicetot2)
1288 #ifdef MPP_LAND
1289         endif
1290 #endif
1291 ! END Final water balance variables
1292 #endif
1293
1294 #ifdef HYDRO_D
1295 #ifndef NCEP_WCOSS
1296         print *, "Finished Aggregation..."
1297 #else
1298        write(78,*) "Finished Aggregation..."
1299 #endif
1300 #endif
1301
1302
1303       end subroutine aggregateDomain
1304
1305       subroutine RunOffDisag(runoff1x_in, runoff1x, area_lsm,cellArea, infxswgt, AGGFACTRT, ix,jx)
1306         implicit none
1307         real, dimension(:,:) :: runoff1x_in, runoff1x, area_lsm, cellArea, infxswgt
1308         integer :: i,j,ix,jx,AGGFACYRT, AGGFACXRT, AGGFACTRT, IXXRT, JYYRT
1309
1310         do J=1,JX
1311         do I=1,IX
1312              do AGGFACYRT=AGGFACTRT-1,0,-1
1313              do AGGFACXRT=AGGFACTRT-1,0,-1
1314                IXXRT=I*AGGFACTRT-AGGFACXRT
1315                JYYRT=J*AGGFACTRT-AGGFACYRT
1316 #ifdef MPP_LAND
1317        if(left_id.ge.0) IXXRT=IXXRT+1
1318        if(down_id.ge.0) JYYRT=JYYRT+1
1319 #endif
1320 !DJG Implement subgrid weighting routine...
1321                if( (runoff1x_in(i,j) .lt. 0) .or. (runoff1x_in(i,j) .gt. 1000) ) then
1322                     runoff1x(IXXRT,JYYRT) = 0
1323                else
1324                     runoff1x(IXXRT,JYYRT)=runoff1x_in(i,j)*area_lsm(I,J)     &
1325                         *INFXSWGT(IXXRT,JYYRT)/cellArea(IXXRT,JYYRT)
1326                endif
1327
1328              enddo
1329              enddo
1330         enddo
1331         enddo
1332
1333       end subroutine RunOffDisag
1334
1335
1336 ! This routine was extracted from the aggregateDomain routine above to do simple depth aggregation.
1337 ! There might be a simpler way.
1338 subroutine RunoffAggregate(runoff_in, runoff_out, aggfactrt, ix, jx)
1339   implicit none
1340   ! Input variables
1341   integer, intent(in) :: aggfactrt, ix, jx
1342   real, intent(in), dimension(:,:) :: runoff_in
1343   real, intent(inout), dimension(:,:) :: runoff_out
1344   ! Local variables
1345   integer :: i, j, aggfacyrt, aggfacxrt, ixxrt, jyyrt
1346   real :: runoffagg
1347   do j=1,jx
1348     do i=1,ix
1349     runoffagg = 0.
1350     do aggfacyrt=aggfactrt-1,0,-1
1351       do aggfacxrt=aggfactrt-1,0,-1
1352         ixxrt = i * aggfactrt - aggfacxrt
1353         jyyrt = j * aggfactrt - aggfacyrt
1354 #ifdef MPP_LAND
1355         if(left_id.ge.0) ixxrt = ixxrt+1
1356         if(down_id.ge.0) jyyrt = jyyrt+1
1357 #endif
1358         runoffagg = runoffagg + runoff_in(ixxrt,jyyrt)
1359       end do
1360     end do
1361    runoff_out(i,j) = runoffagg / (aggfactrt**2)
1362   end do
1363 end do
1364 end subroutine RunoffAggregate
1365
1366 subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp)
1367 implicit none
1368 integer ntime, did
1369 integer rst_out, ix,jx
1370 !        integer, OPTIONAL:: ix0,jx0
1371 integer:: ix0,jx0
1372 integer, dimension(ix0,jx0),optional :: vegtyp, soltyp
1373 integer            :: iret, ncid, ascIndId
1374
1375
1376
1377 ! read the namelist
1378 ! the lsm namelist will be read by rtland sequentially again.
1379 !call read_rt_nlst(nlst(did) )
1380
1381 ! Some field of this structure are already initialized by the CPL component (e.g. DT)
1382 call orchestrator%config%init_nlst(did)
1383
1384 if(nlst(did)%rtFlag .eq. 0) return
1385
1386 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1387
1388 ! get the dimension
1389 call get_file_dimension(trim(nlst(did)%geo_static_flnm), ix,jx)
1390
1391
1392 #ifdef MPP_LAND
1393
1394 if (nlst(did)%sys_cpl .eq. 1 .or. nlst(did)%sys_cpl .eq. 4) then
1395    !sys_cpl: 1-- coupling with HRLDAS but running offline lsm;
1396    !         2-- coupling with WRF but do not run offline lsm
1397    !         3-- coupling with LIS and do not run offline lsm
1398    !         4:  coupling with CLM
1399
1400    ! create 2 dimensiaon logical mapping of the CPUs for coupling with CLM or HRLDAS.
1401    call log_map2d()
1402
1403    global_nx = ix  ! get from land model
1404    global_ny = jx  ! get from land model
1405
1406    call mpp_land_bcast_int1(global_nx)
1407    call mpp_land_bcast_int1(global_ny)
1408
1409 !!! temp set global_nx to ix
1410    rt_domain(did)%ix = global_nx
1411    rt_domain(did)%jx = global_ny
1412
1413    ! over write the ix and jx
1414    call MPP_LAND_PAR_INI(1,rt_domain(did)%ix,rt_domain(did)%jx,&
1415         nlst(did)%AGGFACTRT)
1416 else
1417    !  coupled with WRF, LIS
1418    numprocs = node_info(1,1)
1419
1420    call wrf_LAND_set_INIT(node_info,numprocs,nlst(did)%AGGFACTRT)
1421
1422    rt_domain(did)%ix = local_nx
1423    rt_domain(did)%jx = local_ny
1424 endif
1425
1426
1427 rt_domain(did)%g_IXRT=global_rt_nx
1428 rt_domain(did)%g_JXRT=global_rt_ny
1429 rt_domain(did)%ixrt = local_rt_nx
1430 rt_domain(did)%jxrt = local_rt_ny
1431
1432 #ifdef HYDRO_D
1433 #ifndef NCEP_WCOSS
1434 write(6,*) "rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt"
1435 write(6,*)  rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt
1436 write(6,*) "rt_domain(did)%ix, rt_domain(did)%jx "
1437 write(6,*) rt_domain(did)%ix, rt_domain(did)%jx
1438 write(6,*) "global_nx, global_ny, local_nx, local_ny"
1439 write(6,*) global_nx, global_ny, local_nx, local_ny
1440 #else
1441 write(78,*) "rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt"
1442 write(78,*)  rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt
1443 write(78,*) "rt_domain(did)%ix, rt_domain(did)%jx "
1444 write(78,*) rt_domain(did)%ix, rt_domain(did)%jx
1445 write(78,*) "global_nx, global_ny, local_nx, local_ny"
1446 write(78,*) global_nx, global_ny, local_nx, local_ny
1447 #endif
1448 #endif
1449 #else
1450 ! sequential
1451 rt_domain(did)%ix = ix
1452 rt_domain(did)%jx = jx
1453 rt_domain(did)%ixrt = ix*nlst(did)%AGGFACTRT
1454 rt_domain(did)%jxrt = jx*nlst(did)%AGGFACTRT
1455 #endif
1456
1457
1458 !      allocate rt arrays
1459
1460
1461 call getChanDim(did)
1462
1463 #ifdef HYDRO_D
1464 #ifndef NCEP_WCOSS
1465 write(6,*) "finish getChanDim "
1466 #else
1467 write(78,*) "finish getChanDim "
1468 #endif
1469 #endif
1470
1471 ! ADCHANGE: get global attributes
1472 ! need to set these after getChanDim since it allocates rt_domain vals to 0
1473      call get_file_globalatts(trim(nlst(did)%geo_static_flnm), &
1474           rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater)
1475
1476 #ifdef HYDRO_D
1477 #ifndef NCEP_WCOSS
1478       write(6,*) "hydro_ini: rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater"
1479       write(6,*) rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater
1480 #endif
1481 #endif
1482
1483 if(nlst(did)%GWBASESWCRT .eq. 3 ) then
1484    call gw2d_allocate(did,&
1485         rt_domain(did)%ixrt,&
1486         rt_domain(did)%jxrt,&
1487         nlst(did)%nsoil)
1488 #ifdef HYDRO_D
1489 #ifndef NCEP_WCOSS
1490    write(6,*) "finish gw2d_allocate"
1491 #else
1492    write(78,*) "finish gw2d_allocate"
1493 #endif
1494 #endif
1495 endif
1496
1497 ! calculate the distance between grids for routing.
1498 ! decompose the land parameter/data
1499
1500
1501 !      ix0= rt_domain(did)%ix
1502 !      jx0= rt_domain(did)%jx
1503 if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then
1504    if(present(vegtyp)) then
1505       call lsm_input(did,ix0=ix0,jx0=jx0,vegtyp0=vegtyp,soltyp0=soltyp)
1506    else
1507       call lsm_input(did,ix0=ix0,jx0=jx0)
1508    endif
1509 endif
1510
1511
1512 #ifdef HYDRO_D
1513 #ifndef NCEP_WCOSS
1514 write(6,*) "finish decomposion"
1515 #else
1516 write(78,*) "finish decomposion"
1517 #endif
1518 #endif
1519
1520 if((nlst(did)%channel_only .eq. 1 .or. nlst(did)%channelBucket_only .eq. 1) .and. &
1521      nlst(1)%io_form_outputs .ne. 0) then
1522    !! This is the "decoder ring" for reading channel-only forcing from io_form_outputs=1,2 CHRTOUT files.
1523    !! Only needed on io_id
1524    if(my_id .eq. io_id) then
1525       allocate(rt_domain(did)%ascendIndex(rt_domain(did)%gnlinksl))
1526       iret = nf90_open(trim(nlst(1)%route_link_f),NF90_NOWRITE,ncid=ncid)
1527       !if(iret .ne. 0) call hdyro_stop
1528       if(iret .ne. 0) call hydro_stop('ERROR: Unable to open RouteLink file for index extraction')
1529       iret = nf90_inq_varid(ncid,'ascendingIndex',ascIndId)
1530       if(iret .ne. 0) call hydro_stop('ERROR: Unable to find ascendingIndex from RouteLink file.')
1531       iret = nf90_get_var(ncid,ascIndId,rt_domain(did)%ascendIndex)
1532       if(iret .ne. 0) call hydro_stop('ERROR: Unable to extract ascendingIndex from RouteLink file.')
1533       iret = nf90_close(ncid)
1534       if(iret .ne. 0) call hydro_stop('ERROR: Unable to close RouteLink file.')
1535    else
1536       allocate(rt_domain(did)%ascendIndex(1))
1537       rt_domain(did)%ascendIndex(1)=-9
1538    endif
1539 endif
1540
1541
1542 call get_dist_lsm(did) !! always needed (channel_only and channelBucket_only)
1543 if(nlst(did)%channel_only .ne. 1)  call get_dist_lrt(did) !! needed forchannelBucket_only
1544
1545 ! rt model initilization
1546 call LandRT_ini(did)
1547
1548 if(nlst(did)%GWBASESWCRT .eq. 3 ) then
1549
1550    call gw2d_ini(did,&
1551         nlst(did)%dt,&
1552         nlst(did)%dxrt0)
1553 #ifdef HYDRO_D
1554 #ifndef NCEP_WCOSS
1555    write(6,*) "finish gw2d_ini"
1556 #else
1557    write(78,*) "finish gw2d_ini"
1558 #endif
1559 #endif
1560 endif
1561 #ifdef HYDRO_D
1562 #ifndef NCEP_WCOSS
1563 write(6,*) "finish LandRT_ini"
1564 #else
1565 write(78,*) "finish LandRT_ini"
1566 #endif
1567 #endif
1568
1569
1570 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1571 if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then
1572
1573    if (nlst(did)%TERADJ_SOLAR.eq.1 .and. nlst(did)%CHANRTSWCRT.ne.2) then   ! Perform ter rain adjustment of incoming solar
1574 #ifdef MPP_LAND
1575       call MPP_seq_land_SO8(rt_domain(did)%SO8LD_D,rt_domain(did)%SO8LD_Vmax,&
1576            rt_domain(did)%TERRAIN, rt_domain(did)%dist_lsm,&
1577            rt_domain(did)%ix,rt_domain(did)%jx,global_nx,global_ny)
1578 #else
1579       call seq_land_SO8(rt_domain(did)%SO8LD_D,rt_domain(did)%SO8LD_Vmax,&
1580            rt_domain(did)%TERRAIN,rt_domain(did)%dist_lsm,&
1581            rt_domain(did)%ix,rt_domain(did)%jx)
1582 #endif
1583    endif
1584 endif
1585
1586 if (nlst(did)%GWBASESWCRT .gt. 0) then
1587    if(nlst(did)%UDMP_OPT .eq. 1) then
1588       call get_basn_area_nhd(rt_domain(did)%basns_area)
1589    else
1590       call get_basn_area(did)
1591    endif
1592 endif
1593
1594 if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT .eq. 2 ) then
1595    call get_node_area(did)
1596 endif
1597
1598
1599 #ifdef WRF_HYDRO_NUDGING
1600 if(nlst(did)%CHANRTSWCRT .ne. 0) call init_stream_nudging
1601 #endif
1602
1603
1604 !    if (trim(nlst_rt(did)%restart_file) == "") then
1605 ! output at the initial time
1606 !        call HYDRO_out(did)
1607 !        return
1608 !    endif
1609
1610 ! restart the file
1611
1612         ! jummp the initial time output
1613 !        rt_domain(did)%out_counts = rt_domain(did)%out_counts + 1
1614 !        rt_domain(did)%his_out_counts = rt_domain(did)%his_out_counts + 1
1615
1616
1617 call HYDRO_rst_in(did)
1618 !#ifdef HYDRO_REALTIME
1619 if (trim(nlst(did)%restart_file) == "") then
1620   call HYDRO_out(did, 0)
1621 else
1622   call HYDRO_out(did, 1)
1623 endif
1624 !! JLM: This is only currently part 1/2 of a better accumulation tracking strategy.
1625 !! The parts:
1626 !! 1) (this) zero accumulations on restart/init after any t=0 outputs are written.
1627 !! 2) introduce a variable in the output files which indicates if accumulations are
1628 !!    reset after this output.
1629 !! Here we move zeroing to after AFTER outputing the accumulations at the initial time.
1630 !! This was previously done in HYDRO_rst_in and so output accumulations at time
1631 !! zero were getting zeroed and then writtent to file, which looses information.
1632 !! Note that nlst_rt(did)%rstrt_swc is not changed at any point in between here and the rst_in.
1633 if(nlst(did)%rstrt_swc.eq.1) then  !Switch for rest of restart accum vars...
1634    print *, "Resetting RESTART Accumulation Variables to 0...",nlst(did)%rstrt_swc
1635    ! Under channel-only , these first three variables are not allocated.
1636    if(allocated(rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake))    rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake = zeroFlt
1637    if(allocated(rt_domain(did)%QSTRMVOLRT_ACC))  rt_domain(did)%QSTRMVOLRT_ACC   = zeroFlt
1638    ! These variables are optionally allocated, if their output is requested.
1639    if(allocated(rt_domain(did)%accSfcLatRunoff)) rt_domain(did)%accSfcLatRunoff = zeroDbl
1640    if(allocated(rt_domain(did)%accBucket))       rt_domain(did)%accBucket    = zeroDbl
1641 end if
1642
1643 end subroutine HYDRO_ini
1644
1645       subroutine lsm_input(did,ix0,jx0,vegtyp0,soltyp0)
1646          implicit none
1647          integer did, leng, ncid, ierr_flg
1648          parameter(leng=100)
1649          integer :: i,j, nn
1650          integer, allocatable, dimension(:,:) :: soltyp
1651          real, dimension(leng) :: xdum1, MAXSMC,refsmc,wltsmc
1652
1653         integer :: ix0,jx0
1654         integer, dimension(ix0,jx0),OPTIONAL :: vegtyp0, soltyp0
1655         integer :: iret, istmp
1656
1657 #ifdef HYDRO_D
1658 #ifndef NCEP_WCOSS
1659          write(6,*) RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx
1660 #else
1661          write(78,*) RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx
1662 #endif
1663 #endif
1664
1665          allocate(soltyp(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx) )
1666
1667          soltyp = 0
1668          call get2d_lsm_soltyp(soltyp,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm))
1669
1670
1671          call get2d_lsm_real("HGT",RT_DOMAIN(did)%TERRAIN,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm))
1672
1673          call get2d_lsm_real("XLAT",RT_DOMAIN(did)%lat_lsm,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm))
1674          call get2d_lsm_real("XLONG",RT_DOMAIN(did)%lon_lsm,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm))
1675          call get2d_lsm_vegtyp(RT_DOMAIN(did)%VEGTYP,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm))
1676
1677
1678             if(nlst(did)%sys_cpl .eq. 2 ) then
1679               ! coupling with WRF
1680                 if(present(soltyp0) ) then
1681                   where(VEGTYP0 == rt_domain(did)%iswater .or. VEGTYP0 == rt_domain(did)%islake) soltyp0 = rt_domain(did)%isoilwater
1682                   where(soltyp0 == rt_domain(did)%isoilwater) VEGTYP0 = rt_domain(did)%iswater
1683                    soltyp = soltyp0
1684                    RT_DOMAIN(did)%VEGTYP = VEGTYP0
1685                 endif
1686             endif
1687
1688          where(RT_DOMAIN(did)%VEGTYP == rt_domain(did)%iswater .or. RT_DOMAIN(did)%VEGTYP == rt_domain(did)%islake) soltyp = rt_domain(did)%isoilwater
1689          where(soltyp == rt_domain(did)%isoilwater) RT_DOMAIN(did)%VEGTYP = rt_domain(did)%iswater
1690
1691 ! LKSAT,
1692 ! temporary set
1693        RT_DOMAIN(did)%SMCRTCHK = 0
1694        RT_DOMAIN(did)%SMCAGGRT = 0
1695        RT_DOMAIN(did)%STCAGGRT = 0
1696        RT_DOMAIN(did)%SH2OAGGRT = 0
1697
1698
1699        rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil) = nlst(did)%zsoil8(1:nlst(did)%nsoil)
1700
1701        rt_domain(did)%subsurface%properties%sldpth(1) = abs( RT_DOMAIN(did)%subsurface%properties%zsoil(1) )
1702        do i = 2, nlst(did)%nsoil
1703           rt_domain(did)%subsurface%properties%sldpth(i) = RT_DOMAIN(did)%subsurface%properties%zsoil(i-1)-RT_DOMAIN(did)%subsurface%properties%zsoil(i)
1704        enddo
1705        rt_domain(did)%subsurface%properties%soldeprt = -1.0*RT_DOMAIN(did)%subsurface%properties%zsoil(nlst(did)%NSOIL)
1706
1707        ierr_flg = 99
1708        if(trim(nlst(did)%hydrotbl_f) == "") then
1709            call hydro_stop("FATAL ERROR: hydrotbl_f is empty. Please input a netcdf file. ")
1710        endif
1711
1712 #ifdef MPP_LAND
1713        if(my_id .eq. IO_id) then
1714 #endif
1715           ierr_flg = nf90_open(trim(nlst(did)%hydrotbl_f), nf90_NOWRITE, ncid)
1716 #ifdef MPP_LAND
1717        endif
1718        call mpp_land_bcast_int1(ierr_flg)
1719 #endif
1720        if( ierr_flg .ne. 0) then
1721           ! input from HYDRO.tbl FILE
1722 !      input OV_ROUGH from OVROUGH.TBL
1723 #ifdef MPP_LAND
1724        if(my_id .eq. IO_id) then
1725 #endif
1726
1727 #ifndef NCEP_WCOSS
1728        open(71,file="HYDRO.TBL", form="formatted")
1729 !read OV_ROUGH first
1730           read(71,*) nn
1731           read(71,*)
1732           do i = 1, nn
1733              read(71,*) RT_DOMAIN(did)%OV_ROUGH(i)
1734           end do
1735 !read parameter for LKSAT
1736           read(71,*) nn
1737           read(71,*)
1738           do i = 1, nn
1739              read(71,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i)
1740           end do
1741        close(71)
1742 #else
1743        open(13, form="formatted")
1744 !read OV_ROUGH first
1745           read(13,*) nn
1746           read(13,*)
1747           do i = 1, nn
1748              read(13,*) RT_DOMAIN(did)%OV_ROUGH(i)
1749           end do
1750 !read parameter for LKSAT
1751           read(13,*) nn
1752           read(13,*)
1753           do i = 1, nn
1754              read(13,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i)
1755           end do
1756        close(13)
1757 #endif
1758
1759 #ifdef MPP_LAND
1760        endif
1761        call mpp_land_bcast_real(leng,RT_DOMAIN(did)%OV_ROUGH)
1762        call mpp_land_bcast_real(leng,xdum1)
1763        call mpp_land_bcast_real(leng,MAXSMC)
1764        call mpp_land_bcast_real(leng,refsmc)
1765        call mpp_land_bcast_real(leng,wltsmc)
1766 #endif
1767
1768        rt_domain(did)%lksat = 0.0
1769        do j = 1, RT_DOMAIN(did)%jx
1770              do i = 1, RT_DOMAIN(did)%ix
1771                 !yw rt_domain(did)%lksat(i,j) = xdum1(soltyp(i,j) ) * 1000.0
1772                 rt_domain(did)%lksat(i,j) = xdum1(soltyp(i,j) )
1773                 rt_domain(did)%SMCMAX1(i,j) = MAXSMC(soltyp(I,J))
1774                 rt_domain(did)%SMCREF1(i,j) = refsmc(soltyp(I,J))
1775                 rt_domain(did)%SMCWLT1(i,j) = wltsmc(soltyp(I,J))
1776                 !ADCHANGE: Add some sanity checks in case calibration knocks the order of these out of sequence.
1777                 !The min diffs were pulled from the existing HYDRO.TBL defaults.
1778                 !Currently water is 0, so enforcing 0 as the absolute min.
1779                 rt_domain(did)%SMCMAX1(i,j) = min(rt_domain(did)%SMCMAX1(i,j), 1.0)
1780                 rt_domain(did)%SMCREF1(i,j) = max(min(rt_domain(did)%SMCREF1(i,j), rt_domain(did)%SMCMAX1(i,j) - 0.01), 0.0)
1781                 rt_domain(did)%SMCWLT1(i,j) = max(min(rt_domain(did)%SMCWLT1(i,j), rt_domain(did)%SMCREF1(i,j) - 0.01), 0.0)
1782                 IF(rt_domain(did)%VEGTYP(i,j) > 0 ) THEN   ! created 2d ov_rough
1783                     rt_domain(did)%OV_ROUGH2d(i,j) = RT_DOMAIN(did)%OV_ROUGH(rt_domain(did)%VEGTYP(I,J))
1784                 endif
1785              end do
1786        end do
1787
1788        call hdtbl_out(did)
1789     else
1790        ! input from HYDRO.TBL.nc file
1791        print*, "reading from hydrotbl_f(HYDRO.TBL.nc)  file ...."
1792        call hdtbl_in_nc(did)
1793        if (noah_lsm%imperv_option .eq. 9) then
1794          !ADCHANGE: For consistency, mirror urban and param value checks used in table read
1795          where (rt_domain(did)%VEGTYP == rt_domain(did)%isurban)
1796            rt_domain(did)%SMCMAX1 = 0.45
1797            rt_domain(did)%SMCREF1 = 0.42
1798            rt_domain(did)%SMCWLT1 = 0.40
1799          endwhere
1800        endif
1801        where (rt_domain(did)%SMCMAX1 .gt. 1.0) rt_domain(did)%SMCMAX1 = 1.0
1802        rt_domain(did)%SMCREF1 = max(min(rt_domain(did)%SMCREF1, rt_domain(did)%SMCMAX1 - 0.01), 0.0)
1803        rt_domain(did)%SMCWLT1 = max(min(rt_domain(did)%SMCWLT1, rt_domain(did)%SMCREF1 - 0.01), 0.0)
1804     endif
1805
1806        rt_domain(did)%soiltyp = soltyp
1807
1808        if(allocated(soltyp)) deallocate(soltyp)
1809
1810
1811       end subroutine lsm_input
1812
1813
1814 end module module_HYDRO_drv
1815
1816 ! stop the job due to the fatal error.
1817 subroutine HYDRO_finish()
1818 #ifdef MPP_LAND
1819     USE module_mpp_land
1820 #endif
1821 #ifdef WRF_HYDRO_NUDGING
1822     use module_stream_nudging,  only: finish_stream_nudging
1823 #endif
1824
1825     integer :: ierr
1826
1827 #ifdef WRF_HYDRO_NUDGING
1828     call finish_stream_nudging()
1829 #endif
1830 #ifndef NCEP_WCOSS
1831     print*, "The model finished successfully......."
1832 #else
1833     write(78,*) "The model finished successfully......."
1834 #endif
1835 #ifdef MPP_LAND
1836 !         call mpp_land_abort()
1837 #ifndef NCEP_WCOSS
1838     call flush(6)
1839 #else
1840     call flush(78)
1841     close(78)
1842 #endif
1843     call mpp_land_sync()
1844     call MPI_finalize(ierr)
1845     stop
1846 #else
1847
1848 #ifndef WRF_HYDRO_NUDGING
1849     stop  !!JLM want to time at the top NoahMP level.
1850 #endif
1851
1852 #endif
1853
1854     return
1855 end  subroutine HYDRO_finish