From b95205c57949d5d4b428a06524a53df3603209f6 Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Mon, 20 May 2024 15:17:12 -0600 Subject: [PATCH] Don't reset accumulation if flag is set in namelist --- src/CPL/NoahMP_cpl/hrldas_drv_HYDRO.F | 8 + src/HYDRO_drv/module_HYDRO_drv.F90 | 1910 ++++++++--------- .../NoahMP/IO_code/main_hrldas_driver.F | 2 + .../IO_code/module_NoahMP_hrldas_driver.F | 8 +- src/utils/module_hydro_stop.F90 | 2 +- 5 files changed, 968 insertions(+), 962 deletions(-) diff --git a/src/CPL/NoahMP_cpl/hrldas_drv_HYDRO.F b/src/CPL/NoahMP_cpl/hrldas_drv_HYDRO.F index 2e92fd3b7..53e77b1d9 100644 --- a/src/CPL/NoahMP_cpl/hrldas_drv_HYDRO.F +++ b/src/CPL/NoahMP_cpl/hrldas_drv_HYDRO.F @@ -395,6 +395,14 @@ subroutine get_iocflag(did, iocflag) iocflag = nlst(did)%io_config_outputs end subroutine get_iocflag + subroutine get_rstrt_swc(did, restart_flag) + use config_base, only: nlst + implicit none + integer, intent(in) :: did + integer, intent(out) :: restart_flag + restart_flag = nlst(did)%RSTRT_SWC + end subroutine + subroutine get_t0OutputFlag(did, t0OutputFlag) use config_base, only: nlst implicit none diff --git a/src/HYDRO_drv/module_HYDRO_drv.F90 b/src/HYDRO_drv/module_HYDRO_drv.F90 index 15ad1e02c..695188cd3 100644 --- a/src/HYDRO_drv/module_HYDRO_drv.F90 +++ b/src/HYDRO_drv/module_HYDRO_drv.F90 @@ -1,128 +1,128 @@ module module_HYDRO_drv #ifdef MPP_LAND - use module_HYDRO_io, only: output_rt, mpp_output_chrt, mpp_output_lakes, mpp_output_chrtgrd, & - restart_out_bi, restart_in_bi, mpp_output_chrt2, mpp_output_lakes2, & - hdtbl_in_nc, hdtbl_out + use module_HYDRO_io, only: output_rt, mpp_output_chrt, mpp_output_lakes, mpp_output_chrtgrd, & + restart_out_bi, restart_in_bi, mpp_output_chrt2, mpp_output_lakes2, & + hdtbl_in_nc, hdtbl_out USE module_mpp_land #else - use module_HYDRO_io, only: output_rt, output_chrt, output_chrt2, output_lakes + use module_HYDRO_io, only: output_rt, output_chrt, output_chrt2, output_lakes #endif - use module_NWM_io, only: output_chrt_NWM, output_rt_NWM, output_lakes_NWM,& - output_chrtout_grd_NWM, output_lsmOut_NWM, & - output_frxstPts, output_chanObs_NWM, output_gw_NWM - use module_HYDRO_io, only: sub_output_gw, restart_out_nc, restart_in_nc, & + use module_NWM_io, only: output_chrt_NWM, output_rt_NWM, output_lakes_NWM,& + output_chrtout_grd_NWM, output_lsmOut_NWM, & + output_frxstPts, output_chanObs_NWM, output_gw_NWM + use module_HYDRO_io, only: sub_output_gw, restart_out_nc, restart_in_nc, & get_file_dimension , get_file_globalatts, get2d_lsm_real, get2d_lsm_vegtyp, get2d_lsm_soltyp, & output_lsm, output_GW_Diag - use module_HYDRO_io, only : output_lakes2 - use module_rt_data, only: rt_domain - use module_GW_baseflow - use module_gw_gw2d - use module_gw_gw2d_data, only: gw2d - use module_channel_routing, only: drive_channel, drive_channel_rsl - use orchestrator_base - use config_base, only: nlst, noah_lsm - use module_routing, only: getChanDim, landrt_ini - use module_HYDRO_utils - use module_lsm_forcing, only: geth_newdate + use module_HYDRO_io, only : output_lakes2 + use module_rt_data, only: rt_domain + use module_GW_baseflow + use module_gw_gw2d + use module_gw_gw2d_data, only: gw2d + use module_channel_routing, only: drive_channel, drive_channel_rsl + use orchestrator_base + use config_base, only: nlst, noah_lsm + use module_routing, only: getChanDim, landrt_ini + use module_HYDRO_utils + use module_lsm_forcing, only: geth_newdate #ifdef WRF_HYDRO_NUDGING - use module_stream_nudging, only: init_stream_nudging + use module_stream_nudging, only: init_stream_nudging #endif - use module_hydro_stop, only: HYDRO_stop - use module_UDMAP, only: get_basn_area_nhd - use netcdf + use module_hydro_stop, only: HYDRO_stop + use module_UDMAP, only: get_basn_area_nhd + use netcdf - implicit none + implicit none #ifdef HYDRO_D - real :: timeOr = 0 - real :: timeSr = 0 - real :: timeCr = 0 - real :: timeGW = 0 - integer :: clock_count_1 = 0 - integer :: clock_count_2 = 0 - integer :: clock_rate = 0 -#endif - integer :: rtout_factor = 0 - - integer, parameter :: r4 = selected_real_kind(4) - real, parameter :: zeroFlt=0.0000000000000000000_r4 - integer, parameter :: r8 = selected_real_kind(8) - real*8, parameter :: zeroDbl=0.0000000000000000000_r8 - - contains - subroutine HYDRO_rst_out(did) + real :: timeOr = 0 + real :: timeSr = 0 + real :: timeCr = 0 + real :: timeGW = 0 + integer :: clock_count_1 = 0 + integer :: clock_count_2 = 0 + integer :: clock_rate = 0 +#endif + integer :: rtout_factor = 0 + + integer, parameter :: r4 = selected_real_kind(4) + real, parameter :: zeroFlt=0.0000000000000000000_r4 + integer, parameter :: r8 = selected_real_kind(8) + real*8, parameter :: zeroDbl=0.0000000000000000000_r8 + +contains + subroutine HYDRO_rst_out(did) #ifdef WRF_HYDRO_NUDGING - use module_stream_nudging, only: output_nudging_last_obs + use module_stream_nudging, only: output_nudging_last_obs #endif - implicit none - integer:: rst_out - integer did, outflag - character(len=19) out_date + implicit none + integer:: rst_out + integer did, outflag + character(len=19) out_date #ifdef MPP_LAND - character(len=19) str_tmp + character(len=19) str_tmp #endif - rst_out = -99 + rst_out = -99 #ifdef MPP_LAND - if(IO_id .eq. my_id) then -#endif - if(nlst(did)%dt .gt. nlst(did)%rst_dt*60) then - call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%dt*rt_domain(did)%rst_counts)) - else - call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%rst_dt*60*rt_domain(did)%rst_counts)) - endif - if ( (nlst(did)%rst_dt .gt. 0) .and. (out_date(1:19) == nlst(did)%olddate(1:19)) ) then - rst_out = 99 - rt_domain(did)%rst_counts = rt_domain(did)%rst_counts + 1 - endif + if(IO_id .eq. my_id) then +#endif + if(nlst(did)%dt .gt. nlst(did)%rst_dt*60) then + call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%dt*rt_domain(did)%rst_counts)) + else + call geth_newdate(out_date, nlst(did)%startdate, nint(nlst(did)%rst_dt*60*rt_domain(did)%rst_counts)) + endif + if ( (nlst(did)%rst_dt .gt. 0) .and. (out_date(1:19) == nlst(did)%olddate(1:19)) ) then + rst_out = 99 + rt_domain(did)%rst_counts = rt_domain(did)%rst_counts + 1 + endif ! restart every month automatically. - if ( (nlst(did)%olddate(9:10) == "01") .and. (nlst(did)%olddate(12:13) == "00") .and. & - (nlst(did)%olddate(15:16) == "00").and. (nlst(did)%olddate(18:19) == "00") .and. & - (nlst(did)%rst_dt .le. 0) ) then - if(nlst(did)%startdate(1:16) .ne. nlst(did)%olddate(1:16) ) then - rst_out = 99 - endif - endif + if ( (nlst(did)%olddate(9:10) == "01") .and. (nlst(did)%olddate(12:13) == "00") .and. & + (nlst(did)%olddate(15:16) == "00").and. (nlst(did)%olddate(18:19) == "00") .and. & + (nlst(did)%rst_dt .le. 0) ) then + if(nlst(did)%startdate(1:16) .ne. nlst(did)%olddate(1:16) ) then + rst_out = 99 + endif + endif #ifdef MPP_LAND - endif - call mpp_land_bcast_int1(rst_out) + endif + call mpp_land_bcast_int1(rst_out) #endif - if(rst_out .gt. 0) then - write(6,*) "yw check output restart at ",nlst(did)%olddate(1:16) + if(rst_out .gt. 0) then + write(6,*) "yw check output restart at ",nlst(did)%olddate(1:16) #ifdef MPP_LAND - if(nlst(did)%rst_bi_out .eq. 1) then - if(my_id .lt. 10) then - write(str_tmp,'(I1)') my_id - else if(my_id .lt. 100) then - write(str_tmp,'(I2)') my_id - else if(my_id .lt. 1000) then - write(str_tmp,'(I3)') my_id - else if(my_id .lt. 10000) then - write(str_tmp,'(I4)') my_id - else if(my_id .lt. 100000) then - write(str_tmp,'(I5)') my_id - else - continue - endif - call mpp_land_bcast_char(16,nlst(did)%olddate(1:16)) - call RESTART_OUT_bi(trim("HYDRO_RST."//nlst(did)%olddate(1:16) & - //"_DOMAIN"//trim(nlst(did)%hgrid)//"."//trim(str_tmp)), did) - else -#endif - call RESTART_OUT_nc(trim("HYDRO_RST."//nlst(did)%olddate(1:16) & - //"_DOMAIN"//trim(nlst(did)%hgrid)), did) + if(nlst(did)%rst_bi_out .eq. 1) then + if(my_id .lt. 10) then + write(str_tmp,'(I1)') my_id + else if(my_id .lt. 100) then + write(str_tmp,'(I2)') my_id + else if(my_id .lt. 1000) then + write(str_tmp,'(I3)') my_id + else if(my_id .lt. 10000) then + write(str_tmp,'(I4)') my_id + else if(my_id .lt. 100000) then + write(str_tmp,'(I5)') my_id + else + continue + endif + call mpp_land_bcast_char(16,nlst(did)%olddate(1:16)) + call RESTART_OUT_bi(trim("HYDRO_RST."//nlst(did)%olddate(1:16) & + //"_DOMAIN"//trim(nlst(did)%hgrid)//"."//trim(str_tmp)), did) + else +#endif + call RESTART_OUT_nc(trim("HYDRO_RST."//nlst(did)%olddate(1:16) & + //"_DOMAIN"//trim(nlst(did)%hgrid)), did) #ifdef MPP_LAND - endif + endif #endif #ifdef WRF_HYDRO_NUDGING - call output_nudging_last_obs + call output_nudging_last_obs #endif - endif + endif - end subroutine HYDRO_rst_out + end subroutine HYDRO_rst_out subroutine HYDRO_out(did, rstflag) @@ -134,11 +134,11 @@ subroutine HYDRO_out(did, rstflag) real, dimension(RT_DOMAIN(did)%NLINKS,2) :: str_out real, dimension(RT_DOMAIN(did)%NLINKS) :: vel_out - ! real, dimension(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx):: soilmx_tmp, & - ! runoff1x_tmp, runoff2x_tmp, runoff3x_tmp,etax_tmp, & - ! EDIRX_tmp,ECX_tmp,ETTX_tmp,RCX_tmp,HX_tmp,acrain_tmp, & - ! ACSNOM_tmp, esnow2d_tmp, drip2d_tmp,dewfall_tmp, fpar_tmp, & - ! qfx_tmp, prcp_out_tmp, etpndx_tmp +! real, dimension(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx):: soilmx_tmp, & +! runoff1x_tmp, runoff2x_tmp, runoff3x_tmp,etax_tmp, & +! EDIRX_tmp,ECX_tmp,ETTX_tmp,RCX_tmp,HX_tmp,acrain_tmp, & +! ACSNOM_tmp, esnow2d_tmp, drip2d_tmp,dewfall_tmp, fpar_tmp, & +! qfx_tmp, prcp_out_tmp, etpndx_tmp outflag = -99 @@ -189,13 +189,13 @@ subroutine HYDRO_out(did, rstflag) kt = rt_domain(did)%his_out_counts endif - ! jump the ouput for the initial time when it has restart file from routing. +! jump the ouput for the initial time when it has restart file from routing. rtflag = -99 iniflag = -99 #ifdef MPP_LAND if(IO_id .eq. my_id) then #endif - ! if ( (trim(nlst_rt(did)%restart_file) /= "") .and. ( nlst_rt(did)%startdate(1:19) == nlst_rt(did)%olddate(1:19) ) ) then +! if ( (trim(nlst_rt(did)%restart_file) /= "") .and. ( nlst_rt(did)%startdate(1:19) == nlst_rt(did)%olddate(1:19) ) ) then !#ifndef NCEP_WCOSS ! print*, "yyyywww restart_file = ", trim(nlst_rt(did)%restart_file) !#else @@ -203,7 +203,7 @@ subroutine HYDRO_out(did, rstflag) !#endif if ( nlst(did)%startdate(1:19) == nlst(did)%olddate(1:19) ) iniflag = 1 if ( (trim(nlst(did)%restart_file) /= "") .and. ( nlst(did)%startdate(1:19) == nlst(did)%olddate(1:19) ) ) rtflag = 1 - ! endif +! endif #ifdef MPP_LAND endif call mpp_land_bcast_int1(rtflag) @@ -211,7 +211,7 @@ subroutine HYDRO_out(did, rstflag) #endif - !yw keep the initial time otuput for debug +!yw keep the initial time otuput for debug if(rtflag == 1) then rt_domain(did)%restQSTRM = .false. !!! do not reset QSTRM.. at initial time. if(nlst(did)%t0OutputFlag .eq. 0) return @@ -221,35 +221,35 @@ subroutine HYDRO_out(did, rstflag) if(nlst(did)%t0OutputFlag .eq. 0) return endif - if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then - if(nlst(did)%LSMOUT_DOMAIN .eq. 1) then - if(nlst(did)%io_form_outputs .eq. 0) then - call output_lsm(trim(nlst(did)%olddate(1:4)//nlst(did)%olddate(6:7)//nlst(did)%olddate(9:10) & - //nlst(did)%olddate(12:13)//nlst(did)%olddate(15:16)// & - ".LSMOUT_DOMAIN"//trim(nlst(did)%hgrid)), & - did) - else - call output_lsmOut_NWM(did) + if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then + if(nlst(did)%LSMOUT_DOMAIN .eq. 1) then + if(nlst(did)%io_form_outputs .eq. 0) then + call output_lsm(trim(nlst(did)%olddate(1:4)//nlst(did)%olddate(6:7)//nlst(did)%olddate(9:10) & + //nlst(did)%olddate(12:13)//nlst(did)%olddate(15:16)// & + ".LSMOUT_DOMAIN"//trim(nlst(did)%hgrid)), & + did) + else + call output_lsmOut_NWM(did) + endif endif - endif - end if + end if if(nlst(did)%SUBRTSWCRT .gt. 0 .or. & - nlst(did)%OVRTSWCRT .gt. 0 .or. & - nlst(did)%GWBASESWCRT .gt. 0 .or. & - nlst(did)%CHANRTSWCRT .gt. 0 .or. & - nlst(did)%channel_only .gt. 0 .or. & - nlst(did)%channelBucket_only .gt. 0 ) then + nlst(did)%OVRTSWCRT .gt. 0 .or. & + nlst(did)%GWBASESWCRT .gt. 0 .or. & + nlst(did)%CHANRTSWCRT .gt. 0 .or. & + nlst(did)%channel_only .gt. 0 .or. & + nlst(did)%channelBucket_only .gt. 0 ) then if(nlst(did)%RTOUT_DOMAIN .eq. 1 .and. & - nlst(did)%channel_only .eq. 0 .and. & - nlst(did)%channelBucket_only .eq. 0 ) then + nlst(did)%channel_only .eq. 0 .and. & + nlst(did)%channelBucket_only .eq. 0 ) then if(mod(rtout_factor,3) .eq. 2 .or. & - nlst(did)%io_config_outputs .ne. 5 .and. & - nlst(did)%io_config_outputs .ne. 3) then - ! Output gridded routing variables on National Water Model - ! high-res routing grid + nlst(did)%io_config_outputs .ne. 5 .and. & + nlst(did)%io_config_outputs .ne. 3) then +! Output gridded routing variables on National Water Model +! high-res routing grid if(nlst(did)%io_form_outputs .ne. 0) then call output_rt_NWM(did,nlst(did)%igrid) else @@ -257,8 +257,8 @@ subroutine HYDRO_out(did, rstflag) nlst(did)%igrid, nlst(did)%split_output_count, & RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt, & nlst(did)%nsoil, & - ! nlst_rt(did)%startdate, nlst_rt(did)%olddate, - ! rt_domain(did)%subsurface%state%qsubrt,& +! nlst_rt(did)%startdate, nlst_rt(did)%olddate, +! rt_domain(did)%subsurface%state%qsubrt,& nlst(did)%sincedate, nlst(did)%olddate, rt_domain(did)%subsurface%state%qsubrt,& rt_domain(did)%subsurface%properties%zwattablrt,RT_DOMAIN(did)%subsurface%grid_transform%smcrt,& RT_DOMAIN(did)%SUB_RESID, & @@ -275,8 +275,8 @@ subroutine HYDRO_out(did, rstflag) endif ! End check for rtout_factor rtout_factor = rtout_factor + 1 endif - !! JLM disable GW output for NWM. Bring this line back when runtime output options avail. - !! JLM This seems like a more logical place? +!! JLM disable GW output for NWM. Bring this line back when runtime output options avail. +!! JLM This seems like a more logical place? if(nlst(did)%io_form_outputs .ne. 0) then if(nlst(did)%GWBASESWCRT .ne. 0) then if(nlst(did)%channel_only .eq. 0) then @@ -301,7 +301,7 @@ subroutine HYDRO_out(did, rstflag) nlst(did)%igrid, nlst(did)%split_output_count, & RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt, & nlst(did)%nsoil, & - ! nlst(did)%startdate, nlst(did)%olddate, & +! nlst(did)%startdate, nlst(did)%olddate, & nlst(did)%sincedate, nlst(did)%olddate, & gw2d(did)%h, rt_domain(did)%subsurface%grid_transform%smcrt, & gw2d(did)%convgw, gw2d(did)%excess, & @@ -312,7 +312,7 @@ subroutine HYDRO_out(did, rstflag) nlst(did)%output_gw) endif - ! BF end gw2d output section +! BF end gw2d output section #ifdef HYDRO_D #ifndef NCEP_WCOSS @@ -325,7 +325,7 @@ subroutine HYDRO_out(did, rstflag) if (nlst(did)%CHANRTSWCRT.eq.1.or.nlst(did)%CHANRTSWCRT.eq.2) then - !ADCHANGE: Change values for within lake reaches to NA +!ADCHANGE: Change values for within lake reaches to NA str_out = RT_DOMAIN(did)%QLINK vel_out = RT_DOMAIN(did)%velocity @@ -337,51 +337,51 @@ subroutine HYDRO_out(did, rstflag) endif end do endif - !ADCHANGE: End +!ADCHANGE: End if(nlst(did)%io_form_outputs .ne. 0) then - ! Call National Water Model output routine for output on NHD forecast points. +! Call National Water Model output routine for output on NHD forecast points. if(nlst(did)%CHRTOUT_DOMAIN .ne. 0) then call output_chrt_NWM(did) endif - ! Call the subroutine to output frxstPts. +! Call the subroutine to output frxstPts. if(nlst(did)%frxst_pts_out .ne. 0) then call output_frxstPts(did) endif - ! Call the subroutine to output CHANOBS +! Call the subroutine to output CHANOBS if(nlst(did)%CHANOBS_DOMAIN .ne. 0) then call output_chanObs_NWM(did) endif else - ! Call traditional output routines - !ADCHANGE: We suspect this routine is broken so default is now output_chrtout2 - ! if(nlst_rt(did)%CHRTOUT_DOMAIN .eq. 1) then - !#ifdef MPP_LAND - ! call mpp_output_chrt( & - ! rt_domain(did)%gnlinks,rt_domain(did)%gnlinksl,rt_domain(did)%map_l2g, & - !#else - ! call output_chrt( & - !#endif - ! nlst_rt(did)%igrid, nlst_rt(did)%split_output_count, & - ! RT_DOMAIN(did)%NLINKS,RT_DOMAIN(did)%ORDER, & - ! nlst_rt(did)%sincedate,nlst_rt(did)%olddate,RT_DOMAIN(did)%CHLON,& - ! RT_DOMAIN(did)%CHLAT, & - ! RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ZELEV, & - ! !RT_DOMAIN(did)%QLINK,nlst_rt(did)%DT,Kt, & - ! str_out, nlst_rt(did)%DT,Kt, & - ! RT_DOMAIN(did)%STRMFRXSTPTS,nlst_rt(did)%order_to_write, & - ! RT_DOMAIN(did)%NLINKSL,nlst_rt(did)%channel_option, & - ! rt_domain(did)%gages, rt_domain(did)%gageMiss, & - ! nlst_rt(did)%dt & - !#ifdef WRF_HYDRO_NUDGING - ! , RT_DOMAIN(did)%nudge & - !#endif - ! , RT_DOMAIN(did)%accSfcLatRunoff, RT_DOMAIN(did)%accBucket & - ! , RT_DOMAIN(did)%qSfcLatRunoff, RT_DOMAIN(did)%qBucket & - ! , RT_DOMAIN(did)%qin_gwsubbas & - ! , nlst_rt(did)%UDMP_OPT & - ! ) - ! else +! Call traditional output routines +!ADCHANGE: We suspect this routine is broken so default is now output_chrtout2 +! if(nlst_rt(did)%CHRTOUT_DOMAIN .eq. 1) then +!#ifdef MPP_LAND +! call mpp_output_chrt( & +! rt_domain(did)%gnlinks,rt_domain(did)%gnlinksl,rt_domain(did)%map_l2g, & +!#else +! call output_chrt( & +!#endif +! nlst_rt(did)%igrid, nlst_rt(did)%split_output_count, & +! RT_DOMAIN(did)%NLINKS,RT_DOMAIN(did)%ORDER, & +! nlst_rt(did)%sincedate,nlst_rt(did)%olddate,RT_DOMAIN(did)%CHLON,& +! RT_DOMAIN(did)%CHLAT, & +! RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ZELEV, & +! !RT_DOMAIN(did)%QLINK,nlst_rt(did)%DT,Kt, & +! str_out, nlst_rt(did)%DT,Kt, & +! RT_DOMAIN(did)%STRMFRXSTPTS,nlst_rt(did)%order_to_write, & +! RT_DOMAIN(did)%NLINKSL,nlst_rt(did)%channel_option, & +! rt_domain(did)%gages, rt_domain(did)%gageMiss, & +! nlst_rt(did)%dt & +!#ifdef WRF_HYDRO_NUDGING +! , RT_DOMAIN(did)%nudge & +!#endif +! , RT_DOMAIN(did)%accSfcLatRunoff, RT_DOMAIN(did)%accBucket & +! , RT_DOMAIN(did)%qSfcLatRunoff, RT_DOMAIN(did)%qBucket & +! , RT_DOMAIN(did)%qin_gwsubbas & +! , nlst_rt(did)%UDMP_OPT & +! ) +! else if(nlst(did)%CHRTOUT_DOMAIN .gt. 0) then #ifdef MPP_LAND call mpp_output_chrt2(& @@ -394,15 +394,15 @@ subroutine HYDRO_out(did, rstflag) nlst(did)%sincedate,nlst(did)%olddate, & RT_DOMAIN(did)%CHLON, RT_DOMAIN(did)%CHLAT, & RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ZELEV, & - !RT_DOMAIN(did)%QLINK,nlst_rt(did)%DT,Kt, & +!RT_DOMAIN(did)%QLINK,nlst_rt(did)%DT,Kt, & str_out, nlst(did)%DT,Kt, & RT_DOMAIN(did)%NLINKSL,nlst(did)%channel_option, & rt_domain(did)%linkid & #ifdef WRF_HYDRO_NUDGING , RT_DOMAIN(did)%nudge & #endif - !, RT_DOMAIN(did)%QLateral, nlst_rt(did)%io_config_outputs, - !RT_DOMAIN(did)%velocity & +!, RT_DOMAIN(did)%QLateral, nlst_rt(did)%io_config_outputs, +!RT_DOMAIN(did)%velocity & , RT_DOMAIN(did)%QLateral, nlst(did)%io_config_outputs, vel_out & , RT_DOMAIN(did)%accSfcLatRunoff, RT_DOMAIN(did)%accBucket & , RT_DOMAIN(did)%qSfcLatRunoff, RT_DOMAIN(did)%qBucket & @@ -420,7 +420,7 @@ subroutine HYDRO_out(did, rstflag) RT_DOMAIN(did)%ixrt,RT_DOMAIN(did)%jxrt, RT_DOMAIN(did)%NLINKS, & RT_DOMAIN(did)%GCH_NETLNK, & nlst(did)%startdate, nlst(did)%olddate, & - !RT_DOMAIN(did)%qlink, nlst_rt(did)%dt, nlst_rt(did)%geo_finegrid_flnm, & +!RT_DOMAIN(did)%qlink, nlst_rt(did)%dt, nlst_rt(did)%geo_finegrid_flnm, & str_out, nlst(did)%dt, nlst(did)%geo_finegrid_flnm, & RT_DOMAIN(did)%gnlinks,RT_DOMAIN(did)%map_l2g, & RT_DOMAIN(did)%g_ixrt,RT_DOMAIN(did)%g_jxrt ) @@ -431,7 +431,7 @@ subroutine HYDRO_out(did, rstflag) endif if (RT_DOMAIN(did)%NLAKES.gt.0) then if(nlst(did)%io_form_outputs .ne. 0) then - ! Output lakes in NWM format +! Output lakes in NWM format if(nlst(did)%outlake .ne. 0) then call output_lakes_NWM(did,nlst(did)%igrid) endif @@ -481,83 +481,83 @@ subroutine HYDRO_out(did, rstflag) end subroutine HYDRO_out - subroutine HYDRO_rst_in(did) + subroutine HYDRO_rst_in(did) integer :: did integer:: flag - flag = -1 + flag = -1 #ifdef MPP_LAND - if(my_id.eq.IO_id) then + if(my_id.eq.IO_id) then #endif - if (trim(nlst(did)%restart_file) /= "") then - flag = 99 - rt_domain(did)%timestep_flag = 99 ! continue run - endif + if (trim(nlst(did)%restart_file) /= "") then + flag = 99 + rt_domain(did)%timestep_flag = 99 ! continue run + endif #ifdef MPP_LAND - endif - call mpp_land_bcast_int1(flag) + endif + call mpp_land_bcast_int1(flag) #endif - nlst(did)%sincedate = nlst(did)%startdate + nlst(did)%sincedate = nlst(did)%startdate - if (flag.eq.99) then + if (flag.eq.99) then #ifdef MPP_LAND - if(my_id.eq.IO_id) then + if(my_id.eq.IO_id) then #endif #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "*** read restart data: ",trim(nlst(did)%restart_file) + write(6,*) "*** read restart data: ",trim(nlst(did)%restart_file) #else - write(78,*) "*** read restart data: ",trim(nlst(did)%restart_file) + write(78,*) "*** read restart data: ",trim(nlst(did)%restart_file) #endif #endif #ifdef MPP_LAND - endif + endif #endif #ifdef MPP_LAND - if(nlst(did)%rst_bi_in .eq. 1) then - call RESTART_IN_bi(trim(nlst(did)%restart_file), did) - else + if(nlst(did)%rst_bi_in .eq. 1) then + call RESTART_IN_bi(trim(nlst(did)%restart_file), did) + else #endif - call RESTART_IN_nc(trim(nlst(did)%restart_file), did) + call RESTART_IN_nc(trim(nlst(did)%restart_file), did) #ifdef MPP_LAND - endif + endif #endif !yw if (trim(nlst_rt(did)%restart_file) /= "") then !yw nlst_rt(did)%restart_file = "" !yw endif - endif - end subroutine HYDRO_rst_in + endif + end subroutine HYDRO_rst_in - subroutine HYDRO_time_adv(did) + subroutine HYDRO_time_adv(did) implicit none character(len = 19) :: newdate integer did #ifdef MPP_LAND - if(IO_id.eq.my_id) then + if(IO_id.eq.my_id) then #endif - call geth_newdate(newdate, nlst(did)%olddate, nint( nlst(did)%dt)) - nlst(did)%olddate = newdate + call geth_newdate(newdate, nlst(did)%olddate, nint( nlst(did)%dt)) + nlst(did)%olddate = newdate #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "current time is ",newdate + write(6,*) "current time is ",newdate #else - write(78,*) "current time is ",newdate + write(78,*) "current time is ",newdate #endif #endif #ifdef MPP_LAND - endif + endif #endif - end subroutine HYDRO_time_adv + end subroutine HYDRO_time_adv - subroutine HYDRO_exe(did) + subroutine HYDRO_exe(did) implicit none @@ -581,159 +581,159 @@ subroutine HYDRO_exe(did) ! endif -if (nlst(did)%GWBASESWCRT .ne. 0 .or. & - nlst(did)%SUBRTSWCRT .ne. 0 .or. & - nlst(did)%OVRTSWCRT .ne. 0 .or. & - nlst(did)%channel_only .ne. 0 .or. & - nlst(did)%channelBucket_only .ne. 0 ) then + if (nlst(did)%GWBASESWCRT .ne. 0 .or. & + nlst(did)%SUBRTSWCRT .ne. 0 .or. & + nlst(did)%OVRTSWCRT .ne. 0 .or. & + nlst(did)%channel_only .ne. 0 .or. & + nlst(did)%channelBucket_only .ne. 0 ) then - ! step 1) disaggregate specific fields from LSM to Hydro grid - if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then +! step 1) disaggregate specific fields from LSM to Hydro grid + if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then - RT_DOMAIN(did)%overland%streams_and_lakes%surface_water_to_channel = zeroFlt - RT_DOMAIN(did)%LAKE_INFLORT_DUM = rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake + RT_DOMAIN(did)%overland%streams_and_lakes%surface_water_to_channel = zeroFlt + RT_DOMAIN(did)%LAKE_INFLORT_DUM = rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake - if(nlst(did)%SUBRTSWCRT .ne. 0 .or. nlst(did)%OVRTSWCRT .ne. 0) then - call disaggregateDomain_drv(did) - endif - if(nlst(did)%OVRTSWCRT .eq. 0) then - if(nlst(did)%UDMP_OPT .eq. 1) then - call RunOffDisag(RT_DOMAIN(did)%INFXSRT, RT_DOMAIN(did)%landRunOff, & - rt_domain(did)%dist_lsm(:,:,9),rt_domain(did)%overland%properties%distance_to_neighbor(:,:,9), & - RT_DOMAIN(did)%INFXSWGT, nlst(did)%AGGFACTRT, RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx) - endif - endif + if(nlst(did)%SUBRTSWCRT .ne. 0 .or. nlst(did)%OVRTSWCRT .ne. 0) then + call disaggregateDomain_drv(did) + endif + if(nlst(did)%OVRTSWCRT .eq. 0) then + if(nlst(did)%UDMP_OPT .eq. 1) then + call RunOffDisag(RT_DOMAIN(did)%INFXSRT, RT_DOMAIN(did)%landRunOff, & + rt_domain(did)%dist_lsm(:,:,9),rt_domain(did)%overland%properties%distance_to_neighbor(:,:,9), & + RT_DOMAIN(did)%INFXSWGT, nlst(did)%AGGFACTRT, RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx) + endif + endif #ifdef HYDRO_D - call system_clock(count=clock_count_1, count_rate=clock_rate) + call system_clock(count=clock_count_1, count_rate=clock_rate) #endif - endif !! channel_only & channelBucket_only == 0 + endif !! channel_only & channelBucket_only == 0 - ! step 2) - if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then - if(nlst(did)%SUBRTSWCRT .ne.0) then - call SubsurfaceRouting_drv(did) - endif +! step 2) + if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then + if(nlst(did)%SUBRTSWCRT .ne.0) then + call SubsurfaceRouting_drv(did) + endif #ifdef HYDRO_D - call system_clock(count=clock_count_2, count_rate=clock_rate) - timeSr = timeSr + float(clock_count_2-clock_count_1)/float(clock_rate) + call system_clock(count=clock_count_2, count_rate=clock_rate) + timeSr = timeSr + float(clock_count_2-clock_count_1)/float(clock_rate) #ifndef NCEP_WCOSS - write(6,*) "Timing: Subsurface Routing accumulated time--", timeSr + write(6,*) "Timing: Subsurface Routing accumulated time--", timeSr #else - write(78,*) "Timing: Subsurface Routing accumulated time--", timeSr + write(78,*) "Timing: Subsurface Routing accumulated time--", timeSr #endif #endif - end if !! channel_only & channelBucket_only == 0 + end if !! channel_only & channelBucket_only == 0 - ! step 3) todo split - if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then +! step 3) todo split + if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then #ifdef HYDRO_D - call system_clock(count=clock_count_1, count_rate=clock_rate) -#endif - if(nlst(did)%OVRTSWCRT .ne. 0) then - call OverlandRouting_drv(did) - else - !ADCHANGE: Updating landRunoff instead of surface_water_head_routing. This now allows - ! landRunoff to include both infxsrt from the LSM and exfiltration (if subsfc - ! is active) and prevents surface_water_head_routing from inadvertently being - ! passed back to the LSM. - if (nlst(did)%UDMP_OPT .eq. 1) then - ! If subsurface is on, we update landRunOff to include the updated term w/ exfiltration. - ! If subsurface is off, landRunOff does not change from original value so we leave as-is. - if (nlst(did)%SUBRTSWCRT .ne. 0) then - rt_domain(did)%landRunOff = rt_domain(did)%overland%control%infiltration_excess - endif - else - ! If overland is off and subsurface is on, we need to update INFXSRT (LSM grid) - ! since that is what gets fed through the buckets into the channels. So we aggregate - ! the high-res infiltration_excess back to coarse grid. - if (nlst(did)%SUBRTSWCRT .ne. 0) then - call RunoffAggregate(rt_domain(did)%overland%control%infiltration_excess, & - rt_domain(did)%INFXSRT, nlst(did)%AGGFACTRT, & - rt_domain(did)%ix, rt_domain(did)%jx) - endif - endif - ! In either case, if overland is off we need to zero-out surface_water_head since this - ! water is being scraped into channel and should NOT be passed back to the LSM. - rt_domain(did)%overland%control%infiltration_excess = 0. - rt_domain(did)%overland%control%surface_water_head_routing = 0. - endif + call system_clock(count=clock_count_1, count_rate=clock_rate) +#endif + if(nlst(did)%OVRTSWCRT .ne. 0) then + call OverlandRouting_drv(did) + else +!ADCHANGE: Updating landRunoff instead of surface_water_head_routing. This now allows +! landRunoff to include both infxsrt from the LSM and exfiltration (if subsfc +! is active) and prevents surface_water_head_routing from inadvertently being +! passed back to the LSM. + if (nlst(did)%UDMP_OPT .eq. 1) then +! If subsurface is on, we update landRunOff to include the updated term w/ exfiltration. +! If subsurface is off, landRunOff does not change from original value so we leave as-is. + if (nlst(did)%SUBRTSWCRT .ne. 0) then + rt_domain(did)%landRunOff = rt_domain(did)%overland%control%infiltration_excess + endif + else +! If overland is off and subsurface is on, we need to update INFXSRT (LSM grid) +! since that is what gets fed through the buckets into the channels. So we aggregate +! the high-res infiltration_excess back to coarse grid. + if (nlst(did)%SUBRTSWCRT .ne. 0) then + call RunoffAggregate(rt_domain(did)%overland%control%infiltration_excess, & + rt_domain(did)%INFXSRT, nlst(did)%AGGFACTRT, & + rt_domain(did)%ix, rt_domain(did)%jx) + endif + endif +! In either case, if overland is off we need to zero-out surface_water_head since this +! water is being scraped into channel and should NOT be passed back to the LSM. + rt_domain(did)%overland%control%infiltration_excess = 0. + rt_domain(did)%overland%control%surface_water_head_routing = 0. + endif #ifdef HYDRO_D - call system_clock(count=clock_count_2, count_rate=clock_rate) - timeOr = timeOr + float(clock_count_2-clock_count_1)/float(clock_rate) + call system_clock(count=clock_count_2, count_rate=clock_rate) + timeOr = timeOr + float(clock_count_2-clock_count_1)/float(clock_rate) #ifndef NCEP_WCOSS - write(6,*) "Timing: Overland Routing accumulated time--", timeOr + write(6,*) "Timing: Overland Routing accumulated time--", timeOr #else - write(78,*) "Timing: Overland Routing accumulated time--", timeOr + write(78,*) "Timing: Overland Routing accumulated time--", timeOr #endif #endif - RT_DOMAIN(did)%QSTRMVOLRT_TS = rt_domain(did)%overland%streams_and_lakes%surface_water_to_channel - RT_DOMAIN(did)%QSTRMVOLRT_ACC = RT_DOMAIN(did)%QSTRMVOLRT_ACC + RT_DOMAIN(did)%QSTRMVOLRT_TS + RT_DOMAIN(did)%QSTRMVOLRT_TS = rt_domain(did)%overland%streams_and_lakes%surface_water_to_channel + RT_DOMAIN(did)%QSTRMVOLRT_ACC = RT_DOMAIN(did)%QSTRMVOLRT_ACC + RT_DOMAIN(did)%QSTRMVOLRT_TS - RT_DOMAIN(did)%LAKE_INFLORT_TS = rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake-RT_DOMAIN(did)%LAKE_INFLORT_DUM + RT_DOMAIN(did)%LAKE_INFLORT_TS = rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake-RT_DOMAIN(did)%LAKE_INFLORT_DUM #ifdef HYDRO_D - call system_clock(count=clock_count_1, count_rate=clock_rate) + call system_clock(count=clock_count_1, count_rate=clock_rate) #endif - end if !! channel_only & channelBucket_only == 0 + end if !! channel_only & channelBucket_only == 0 - ! step 4) baseflow or groundwater physics - !! channelBucket_only can be anything: the only time you dont run this is if channel_only=1 - if(nlst(did)%channel_only .eq. 0) then - if (nlst(did)%GWBASESWCRT .gt. 0) then - call driveGwBaseflow(did) - endif +! step 4) baseflow or groundwater physics +!! channelBucket_only can be anything: the only time you dont run this is if channel_only=1 + if(nlst(did)%channel_only .eq. 0) then + if (nlst(did)%GWBASESWCRT .gt. 0) then + call driveGwBaseflow(did) + endif #ifdef HYDRO_D - call system_clock(count=clock_count_2, count_rate=clock_rate) - timeGw = timeGw + float(clock_count_2-clock_count_1)/float(clock_rate) + call system_clock(count=clock_count_2, count_rate=clock_rate) + timeGw = timeGw + float(clock_count_2-clock_count_1)/float(clock_rate) #ifndef NCEP_WCOSS - write(6,*) "Timing: GwBaseflow accumulated time--", timeGw + write(6,*) "Timing: GwBaseflow accumulated time--", timeGw #else - write(78,*) "Timing: GwBaseflow accumulated time--", timeGw + write(78,*) "Timing: GwBaseflow accumulated time--", timeGw #endif #endif #ifdef HYDRO_D - call system_clock(count=clock_count_1, count_rate=clock_rate) + call system_clock(count=clock_count_1, count_rate=clock_rate) #endif - end if !! channel_only == 0 + end if !! channel_only == 0 - ! step 5) river channel physics - call driveChannelRouting(did) +! step 5) river channel physics + call driveChannelRouting(did) #ifdef HYDRO_D - call system_clock(count=clock_count_2, count_rate=clock_rate) - timeCr = timeCr + float(clock_count_2-clock_count_1)/float(clock_rate) + call system_clock(count=clock_count_2, count_rate=clock_rate) + timeCr = timeCr + float(clock_count_2-clock_count_1)/float(clock_rate) #ifndef NCEP_WCOSS - write(6,*) "Timing: Channel Routing accumulated time--", timeCr + write(6,*) "Timing: Channel Routing accumulated time--", timeCr #else - write(78,*) "Timing: Channel Routing accumulated time--", timeCr + write(78,*) "Timing: Channel Routing accumulated time--", timeCr #endif #endif - !! if not channel_only - if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then +!! if not channel_only + if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then - ! step 6) aggregate specific fields from Hydro to LSM grid - if (nlst(did)%SUBRTSWCRT .ne.0 .or. nlst(did)%OVRTSWCRT .ne. 0 ) then - call aggregateDomain(did) - endif +! step 6) aggregate specific fields from Hydro to LSM grid + if (nlst(did)%SUBRTSWCRT .ne.0 .or. nlst(did)%OVRTSWCRT .ne. 0 ) then + call aggregateDomain(did) + endif - end if !! channel_only & channelBucket_only == 0 + end if !! channel_only & channelBucket_only == 0 -end if + end if !yw if (nlst_rt(did)%sys_cpl .eq. 2) then - ! advance to next time step +! advance to next time step ! call HYDRO_time_adv(did) - ! output for history +! output for history ! call HYDRO_out(did) !yw endif - call HYDRO_time_adv(did) - call HYDRO_out(did, 1) + call HYDRO_time_adv(did) + call HYDRO_out(did, 1) ! write(90 + my_id,*) "finish calling hydro_exe" @@ -742,292 +742,292 @@ subroutine HYDRO_exe(did) - !! Under channel-only, these variables are not allocated - if(allocated(RT_DOMAIN(did)%SOLDRAIN)) RT_DOMAIN(did)%SOLDRAIN = 0 - if(allocated(rt_domain(did)%subsurface%state%qsubrt)) RT_DOMAIN(did)%subsurface%state%qsubrt = 0 +!! Under channel-only, these variables are not allocated + if(allocated(RT_DOMAIN(did)%SOLDRAIN)) RT_DOMAIN(did)%SOLDRAIN = 0 + if(allocated(rt_domain(did)%subsurface%state%qsubrt)) RT_DOMAIN(did)%subsurface%state%qsubrt = 0 - end subroutine HYDRO_exe + end subroutine HYDRO_exe !---------------------------------------------------- -subroutine driveGwBaseflow(did) + subroutine driveGwBaseflow(did) - implicit none - integer, intent(in) :: did + implicit none + integer, intent(in) :: did - integer :: i, jj, ii + integer :: i, jj, ii - !------------------------------------------------------------------ - !DJG Begin GW/Baseflow Routines - !------------------------------------------------------------------- +!------------------------------------------------------------------ +!DJG Begin GW/Baseflow Routines +!------------------------------------------------------------------- - if (nlst(did)%GWBASESWCRT.ge.1) then ! Switch to activate/specify GW/Baseflow + if (nlst(did)%GWBASESWCRT.ge.1) then ! Switch to activate/specify GW/Baseflow - ! IF (nlst(did)%GWBASESWCRT.GE.1000) THEN ! Switch to activate/specify GW/Baseflow +! IF (nlst(did)%GWBASESWCRT.GE.1000) THEN ! Switch to activate/specify GW/Baseflow - if (nlst(did)%GWBASESWCRT.eq.1 .or. nlst(did)%GWBASESWCRT.eq.2 .or. nlst(did)%GWBASESWCRT.ge.4) then ! Call simple bucket baseflow scheme + if (nlst(did)%GWBASESWCRT.eq.1 .or. nlst(did)%GWBASESWCRT.eq.2 .or. nlst(did)%GWBASESWCRT.ge.4) then ! Call simple bucket baseflow scheme #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "*****yw******start simp_gw_buck " + write(6,*) "*****yw******start simp_gw_buck " #else - write(78,*) "*****yw******start simp_gw_buck " -#endif -#endif - - if(nlst(did)%UDMP_OPT .eq. 1) then - call simp_gw_buck_nhd( & - RT_DOMAIN(did)%ix, RT_DOMAIN(did)%jx, & - RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt, & - RT_DOMAIN(did)%numbasns, nlst(did)%AGGFACTRT, & - nlst(did)%DT, RT_DOMAIN(did)%INFXSWGT, & - RT_DOMAIN(did)%INFXSRT, RT_DOMAIN(did)%SOLDRAIN, & - rt_domain(did)%overland%properties%distance_to_neighbor(:,:,9), rt_domain(did)%dist_lsm(:,:,9), & - RT_DOMAIN(did)%gw_buck_coeff, RT_DOMAIN(did)%gw_buck_exp, & - RT_DOMAIN(did)%gw_buck_loss, & - RT_DOMAIN(did)%z_max, RT_DOMAIN(did)%z_gwsubbas, & - RT_DOMAIN(did)%qout_gwsubbas, RT_DOMAIN(did)%qin_gwsubbas, & - RT_DOMAIN(did)%qloss_gwsubbas, & - nlst(did)%GWBASESWCRT, nlst(did)%OVRTSWCRT, & + write(78,*) "*****yw******start simp_gw_buck " +#endif +#endif + + if(nlst(did)%UDMP_OPT .eq. 1) then + call simp_gw_buck_nhd( & + RT_DOMAIN(did)%ix, RT_DOMAIN(did)%jx, & + RT_DOMAIN(did)%ixrt, RT_DOMAIN(did)%jxrt, & + RT_DOMAIN(did)%numbasns, nlst(did)%AGGFACTRT, & + nlst(did)%DT, RT_DOMAIN(did)%INFXSWGT, & + RT_DOMAIN(did)%INFXSRT, RT_DOMAIN(did)%SOLDRAIN, & + rt_domain(did)%overland%properties%distance_to_neighbor(:,:,9), rt_domain(did)%dist_lsm(:,:,9), & + RT_DOMAIN(did)%gw_buck_coeff, RT_DOMAIN(did)%gw_buck_exp, & + RT_DOMAIN(did)%gw_buck_loss, & + RT_DOMAIN(did)%z_max, RT_DOMAIN(did)%z_gwsubbas, & + RT_DOMAIN(did)%qout_gwsubbas, RT_DOMAIN(did)%qin_gwsubbas, & + RT_DOMAIN(did)%qloss_gwsubbas, & + nlst(did)%GWBASESWCRT, nlst(did)%OVRTSWCRT, & #ifdef MPP_LAND - RT_DOMAIN(did)%LNLINKSL, & + RT_DOMAIN(did)%LNLINKSL, & #else - RT_DOMAIN(did)%numbasns, & + RT_DOMAIN(did)%numbasns, & #endif - rt_domain(did)%basns_area, & - rt_domain(did)%nhdBuckMask, nlst(did)%bucket_loss, & - nlst(did)%channelBucket_only ) + rt_domain(did)%basns_area, & + rt_domain(did)%nhdBuckMask, nlst(did)%bucket_loss, & + nlst(did)%channelBucket_only ) - else - call simp_gw_buck(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,RT_DOMAIN(did)%ixrt,& - RT_DOMAIN(did)%jxrt,RT_DOMAIN(did)%numbasns,RT_DOMAIN(did)%gnumbasns,& - RT_DOMAIN(did)%basns_area,& - RT_DOMAIN(did)%basnsInd, RT_DOMAIN(did)%gw_strm_msk_lind, & - RT_DOMAIN(did)%gwsubbasmsk, RT_DOMAIN(did)%INFXSRT, & - RT_DOMAIN(did)%SOLDRAIN, & - RT_DOMAIN(did)%z_gwsubbas,& - RT_DOMAIN(did)%qin_gwsubbas,RT_DOMAIN(did)%qout_gwsubbas,& - RT_DOMAIN(did)%qinflowbase,& - RT_DOMAIN(did)%gw_strm_msk,RT_DOMAIN(did)%gwbas_pix_ct, & - rt_domain(did)%overland%properties%distance_to_neighbor,nlst(did)%DT,& - RT_DOMAIN(did)%gw_buck_coeff,RT_DOMAIN(did)%gw_buck_exp, & - RT_DOMAIN(did)%z_max,& - nlst(did)%GWBASESWCRT,nlst(did)%OVRTSWCRT) - endif + else + call simp_gw_buck(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,RT_DOMAIN(did)%ixrt,& + RT_DOMAIN(did)%jxrt,RT_DOMAIN(did)%numbasns,RT_DOMAIN(did)%gnumbasns,& + RT_DOMAIN(did)%basns_area,& + RT_DOMAIN(did)%basnsInd, RT_DOMAIN(did)%gw_strm_msk_lind, & + RT_DOMAIN(did)%gwsubbasmsk, RT_DOMAIN(did)%INFXSRT, & + RT_DOMAIN(did)%SOLDRAIN, & + RT_DOMAIN(did)%z_gwsubbas,& + RT_DOMAIN(did)%qin_gwsubbas,RT_DOMAIN(did)%qout_gwsubbas,& + RT_DOMAIN(did)%qinflowbase,& + RT_DOMAIN(did)%gw_strm_msk,RT_DOMAIN(did)%gwbas_pix_ct, & + rt_domain(did)%overland%properties%distance_to_neighbor,nlst(did)%DT,& + RT_DOMAIN(did)%gw_buck_coeff,RT_DOMAIN(did)%gw_buck_exp, & + RT_DOMAIN(did)%z_max,& + nlst(did)%GWBASESWCRT,nlst(did)%OVRTSWCRT) + endif - !! JLM: There's *perhaps* a better location for this output above. - !! If above is better, remove this when runtime output options are avail. - !#ifndef HYDRO_REALTIME - ! call output_GW_Diag(did) - !#endif +!! JLM: There's *perhaps* a better location for this output above. +!! If above is better, remove this when runtime output options are avail. +!#ifndef HYDRO_REALTIME +! call output_GW_Diag(did) +!#endif #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "*****yw******end simp_gw_buck " + write(6,*) "*****yw******end simp_gw_buck " #else - write(78,*) "*****yw******end simp_gw_buck " + write(78,*) "*****yw******end simp_gw_buck " #endif #endif - !!!For parameter setup runs output the percolation for each basin, - !!!otherwise comment out this output... - else if (nlst(did)%gwBaseSwCRT .eq. 3) then +!!!For parameter setup runs output the percolation for each basin, +!!!otherwise comment out this output... + else if (nlst(did)%gwBaseSwCRT .eq. 3) then #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "*****bf******start 2d_gw_model " + write(6,*) "*****bf******start 2d_gw_model " #else - write(78,*) "*****bf******start 2d_gw_model " + write(78,*) "*****bf******start 2d_gw_model " #endif #endif - ! compute qsgwrt between lsm and gw with namelist selected coupling method - ! qsgwrt is defined on the routing grid and needs to be aggregated for SFLX - if (nlst(did)%gwsoilcpl .GT. 0) THEN +! compute qsgwrt between lsm and gw with namelist selected coupling method +! qsgwrt is defined on the routing grid and needs to be aggregated for SFLX + if (nlst(did)%gwsoilcpl .GT. 0) THEN - call gwSoilFlux(did) + call gwSoilFlux(did) - end if + end if - gw2d(did)%excess = 0. + gw2d(did)%excess = 0. - call gwstep(gw2d(did)%ix, gw2d(did)%jx, gw2d(did)%dx, & - gw2d(did)%ltype, gw2d(did)%elev, gw2d(did)%bot, & - gw2d(did)%hycond, gw2d(did)%poros, gw2d(did)%compres, & - gw2d(did)%ho, gw2d(did)%h, gw2d(did)%convgw, & - gw2d(did)%excess, & - gw2d(did)%ebot, gw2d(did)%eocn, gw2d(did)%dt, & - gw2d(did)%istep) + call gwstep(gw2d(did)%ix, gw2d(did)%jx, gw2d(did)%dx, & + gw2d(did)%ltype, gw2d(did)%elev, gw2d(did)%bot, & + gw2d(did)%hycond, gw2d(did)%poros, gw2d(did)%compres, & + gw2d(did)%ho, gw2d(did)%h, gw2d(did)%convgw, & + gw2d(did)%excess, & + gw2d(did)%ebot, gw2d(did)%eocn, gw2d(did)%dt, & + gw2d(did)%istep) - gw2d(did)%ho = gw2d(did)%h + gw2d(did)%ho = gw2d(did)%h - ! put surface exceeding groundwater to surface routing inflow - rt_domain(did)%overland%control%surface_water_head_routing = rt_domain(did)%overland%control%surface_water_head_routing & - + gw2d(did)%excess*1000. ! convert to mm +! put surface exceeding groundwater to surface routing inflow + rt_domain(did)%overland%control%surface_water_head_routing = rt_domain(did)%overland%control%surface_water_head_routing & + + gw2d(did)%excess*1000. ! convert to mm - ! aggregate qsgw from routing to lsm grid - call aggregateQsgw(did) +! aggregate qsgw from routing to lsm grid + call aggregateQsgw(did) - gw2d(did)%istep = gw2d(did)%istep + 1 + gw2d(did)%istep = gw2d(did)%istep + 1 #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "*****bf******end 2d_gw_model " + write(6,*) "*****bf******end 2d_gw_model " #else - write(78,*) "*****bf******end 2d_gw_model " + write(78,*) "*****bf******end 2d_gw_model " #endif #endif - end if + end if - end if !DJG (End if for RTE SWC activation) - !------------------------------------------------------------------ - !DJG End GW/Baseflow Routines - !------------------------------------------------------------------- -end subroutine driveGwBaseflow + end if !DJG (End if for RTE SWC activation) +!------------------------------------------------------------------ +!DJG End GW/Baseflow Routines +!------------------------------------------------------------------- + end subroutine driveGwBaseflow !------------------------------------------- - subroutine driveChannelRouting(did) + subroutine driveChannelRouting(did) - implicit none - integer, intent(in) :: did + implicit none + integer, intent(in) :: did !------------------------------------------------------------------- !------------------------------------------------------------------- !DJG,DNY Begin Channel and Lake Routing Routines !------------------------------------------------------------------- -if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT.eq.2) then - - if(nlst(did)%UDMP_OPT .eq. 1) then - !!! for user defined Reach based Routing method. - - call drive_CHANNEL_RSL(did, nlst(did)%UDMP_OPT,RT_DOMAIN(did)%timestep_flag, RT_DOMAIN(did)%IXRT,RT_DOMAIN(did)%JXRT, & - RT_DOMAIN(did)%LAKE_INFLORT_TS, RT_DOMAIN(did)%QSTRMVOLRT_TS, RT_DOMAIN(did)%TO_NODE, RT_DOMAIN(did)%FROM_NODE, & - RT_DOMAIN(did)%TYPEL, RT_DOMAIN(did)%ORDER, RT_DOMAIN(did)%MAXORDER, RT_DOMAIN(did)%CH_LNKRT, & - rt_domain(did)%overland%streams_and_lakes%lake_mask, nlst(did)%DT, nlst(did)%DTCT, nlst(did)%DTRT_CH, & - RT_DOMAIN(did)%MUSK, RT_DOMAIN(did)%MUSX, RT_DOMAIN(did)%QLINK, & - RT_DOMAIN(did)%CHANLEN, RT_DOMAIN(did)%MannN, RT_DOMAIN(did)%So, RT_DOMAIN(did)%ChSSlp,RT_DOMAIN(did)%Bw, & - RT_DOMAIN(did)%Tw, RT_DOMAIN(did)%Tw_CC, RT_DOMAIN(did)%n_CC, & - RT_DOMAIN(did)%ChannK,& - RT_DOMAIN(did)%RESHT, & - RT_DOMAIN(did)%CVOL, RT_DOMAIN(did)%QLAKEI, & - RT_DOMAIN(did)%QLAKEO, RT_DOMAIN(did)%LAKENODE, & - RT_DOMAIN(did)%QINFLOWBASE, RT_DOMAIN(did)%CHANXI, RT_DOMAIN(did)%CHANYJ, nlst(did)%channel_option, & - RT_DOMAIN(did)%nlinks, RT_DOMAIN(did)%NLINKSL, RT_DOMAIN(did)%LINKID, RT_DOMAIN(did)%node_area, & - RT_DOMAIN(did)%qout_gwsubbas, & - RT_DOMAIN(did)%LAKEIDA, RT_DOMAIN(did)%LAKEIDM, RT_DOMAIN(did)%NLAKES, RT_DOMAIN(did)%LAKEIDX & + if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT.eq.2) then + + if(nlst(did)%UDMP_OPT .eq. 1) then +!!! for user defined Reach based Routing method. + + call drive_CHANNEL_RSL(did, nlst(did)%UDMP_OPT,RT_DOMAIN(did)%timestep_flag, RT_DOMAIN(did)%IXRT,RT_DOMAIN(did)%JXRT, & + RT_DOMAIN(did)%LAKE_INFLORT_TS, RT_DOMAIN(did)%QSTRMVOLRT_TS, RT_DOMAIN(did)%TO_NODE, RT_DOMAIN(did)%FROM_NODE, & + RT_DOMAIN(did)%TYPEL, RT_DOMAIN(did)%ORDER, RT_DOMAIN(did)%MAXORDER, RT_DOMAIN(did)%CH_LNKRT, & + rt_domain(did)%overland%streams_and_lakes%lake_mask, nlst(did)%DT, nlst(did)%DTCT, nlst(did)%DTRT_CH, & + RT_DOMAIN(did)%MUSK, RT_DOMAIN(did)%MUSX, RT_DOMAIN(did)%QLINK, & + RT_DOMAIN(did)%CHANLEN, RT_DOMAIN(did)%MannN, RT_DOMAIN(did)%So, RT_DOMAIN(did)%ChSSlp,RT_DOMAIN(did)%Bw, & + RT_DOMAIN(did)%Tw, RT_DOMAIN(did)%Tw_CC, RT_DOMAIN(did)%n_CC, & + RT_DOMAIN(did)%ChannK,& + RT_DOMAIN(did)%RESHT, & + RT_DOMAIN(did)%CVOL, RT_DOMAIN(did)%QLAKEI, & + RT_DOMAIN(did)%QLAKEO, RT_DOMAIN(did)%LAKENODE, & + RT_DOMAIN(did)%QINFLOWBASE, RT_DOMAIN(did)%CHANXI, RT_DOMAIN(did)%CHANYJ, nlst(did)%channel_option, & + RT_DOMAIN(did)%nlinks, RT_DOMAIN(did)%NLINKSL, RT_DOMAIN(did)%LINKID, RT_DOMAIN(did)%node_area, & + RT_DOMAIN(did)%qout_gwsubbas, & + RT_DOMAIN(did)%LAKEIDA, RT_DOMAIN(did)%LAKEIDM, RT_DOMAIN(did)%NLAKES, RT_DOMAIN(did)%LAKEIDX & #ifdef MPP_LAND - , RT_DOMAIN(did)%nlinks_index, RT_DOMAIN(did)%mpp_nlinks, RT_DOMAIN(did)%yw_mpp_nlinks & - , RT_DOMAIN(did)%LNLINKSL & - , RT_DOMAIN(did)%gtoNode, RT_DOMAIN(did)%toNodeInd, RT_DOMAIN(did)%nToInd & + , RT_DOMAIN(did)%nlinks_index, RT_DOMAIN(did)%mpp_nlinks, RT_DOMAIN(did)%yw_mpp_nlinks & + , RT_DOMAIN(did)%LNLINKSL & + , RT_DOMAIN(did)%gtoNode, RT_DOMAIN(did)%toNodeInd, RT_DOMAIN(did)%nToInd & #endif - , RT_DOMAIN(did)%CH_LNKRT_SL, RT_DOMAIN(did)%landRunOff & + , RT_DOMAIN(did)%CH_LNKRT_SL, RT_DOMAIN(did)%landRunOff & #ifdef WRF_HYDRO_NUDGING - , RT_DOMAIN(did)%nudge & -#endif - , rt_domain(did)%accSfcLatRunoff, rt_domain(did)%accBucket & - , rt_domain(did)%qSfcLatRunoff, rt_domain(did)%qBucket & - , rt_domain(did)%QLateral, rt_domain(did)%velocity & - , rt_domain(did)%qloss & - , RT_DOMAIN(did)%HLINK & - , rt_domain(did)%nlinksize, nlst(did)%OVRTSWCRT & - , nlst(did)%SUBRTSWCRT & - , nlst(did)%channel_only , nlst(did)%channelBucket_only & - , nlst(did)%channel_bypass ) - -else - - call drive_CHANNEL(did, RT_DOMAIN(did)%latval,RT_DOMAIN(did)%lonval, & - RT_DOMAIN(did)%timestep_flag,RT_DOMAIN(did)%IXRT,RT_DOMAIN(did)%JXRT, & - nlst(did)%SUBRTSWCRT, rt_domain(did)%subsurface%state%qsubrt, & - RT_DOMAIN(did)%LAKE_INFLORT_TS, RT_DOMAIN(did)%QSTRMVOLRT_TS,& - RT_DOMAIN(did)%TO_NODE, RT_DOMAIN(did)%FROM_NODE, RT_DOMAIN(did)%TYPEL,& - RT_DOMAIN(did)%ORDER, RT_DOMAIN(did)%MAXORDER, RT_DOMAIN(did)%NLINKS,& - RT_DOMAIN(did)%CH_NETLNK, rt_domain(did)%overland%streams_and_lakes%ch_netrt,RT_DOMAIN(did)%CH_LNKRT,& - rt_domain(did)%overland%streams_and_lakes%lake_mask, nlst(did)%DT, nlst(did)%DTCT, nlst(did)%DTRT_CH,& - RT_DOMAIN(did)%MUSK, RT_DOMAIN(did)%MUSX, RT_DOMAIN(did)%QLINK, & - RT_DOMAIN(did)%QLateral, & - RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ELRT,RT_DOMAIN(did)%CHANLEN,& - RT_DOMAIN(did)%MannN,RT_DOMAIN(did)%So, RT_DOMAIN(did)%ChSSlp, & - RT_DOMAIN(did)%Bw,RT_DOMAIN(did)%Tw,RT_DOMAIN(did)%Tw_CC, RT_DOMAIN(did)%n_CC, & - RT_DOMAIN(did)%ChannK,& - RT_DOMAIN(did)%RESHT, & - RT_DOMAIN(did)%ZELEV, RT_DOMAIN(did)%CVOL, & - RT_DOMAIN(did)%NLAKES, RT_DOMAIN(did)%QLAKEI, RT_DOMAIN(did)%QLAKEO,& - RT_DOMAIN(did)%LAKENODE, rt_domain(did)%overland%properties%distance_to_neighbor, & - RT_DOMAIN(did)%QINFLOWBASE, RT_DOMAIN(did)%CHANXI, & - RT_DOMAIN(did)%CHANYJ, nlst(did)%channel_option, & - RT_DOMAIN(did)%RETDEP_CHAN, RT_DOMAIN(did)%NLINKSL, RT_DOMAIN(did)%LINKID, & - RT_DOMAIN(did)%node_area & + , RT_DOMAIN(did)%nudge & +#endif + , rt_domain(did)%accSfcLatRunoff, rt_domain(did)%accBucket & + , rt_domain(did)%qSfcLatRunoff, rt_domain(did)%qBucket & + , rt_domain(did)%QLateral, rt_domain(did)%velocity & + , rt_domain(did)%qloss & + , RT_DOMAIN(did)%HLINK & + , rt_domain(did)%nlinksize, nlst(did)%OVRTSWCRT & + , nlst(did)%SUBRTSWCRT & + , nlst(did)%channel_only , nlst(did)%channelBucket_only & + , nlst(did)%channel_bypass ) + + else + + call drive_CHANNEL(did, RT_DOMAIN(did)%latval,RT_DOMAIN(did)%lonval, & + RT_DOMAIN(did)%timestep_flag,RT_DOMAIN(did)%IXRT,RT_DOMAIN(did)%JXRT, & + nlst(did)%SUBRTSWCRT, rt_domain(did)%subsurface%state%qsubrt, & + RT_DOMAIN(did)%LAKE_INFLORT_TS, RT_DOMAIN(did)%QSTRMVOLRT_TS,& + RT_DOMAIN(did)%TO_NODE, RT_DOMAIN(did)%FROM_NODE, RT_DOMAIN(did)%TYPEL,& + RT_DOMAIN(did)%ORDER, RT_DOMAIN(did)%MAXORDER, RT_DOMAIN(did)%NLINKS,& + RT_DOMAIN(did)%CH_NETLNK, rt_domain(did)%overland%streams_and_lakes%ch_netrt,RT_DOMAIN(did)%CH_LNKRT,& + rt_domain(did)%overland%streams_and_lakes%lake_mask, nlst(did)%DT, nlst(did)%DTCT, nlst(did)%DTRT_CH,& + RT_DOMAIN(did)%MUSK, RT_DOMAIN(did)%MUSX, RT_DOMAIN(did)%QLINK, & + RT_DOMAIN(did)%QLateral, & + RT_DOMAIN(did)%HLINK, RT_DOMAIN(did)%ELRT,RT_DOMAIN(did)%CHANLEN,& + RT_DOMAIN(did)%MannN,RT_DOMAIN(did)%So, RT_DOMAIN(did)%ChSSlp, & + RT_DOMAIN(did)%Bw,RT_DOMAIN(did)%Tw,RT_DOMAIN(did)%Tw_CC, RT_DOMAIN(did)%n_CC, & + RT_DOMAIN(did)%ChannK,& + RT_DOMAIN(did)%RESHT, & + RT_DOMAIN(did)%ZELEV, RT_DOMAIN(did)%CVOL, & + RT_DOMAIN(did)%NLAKES, RT_DOMAIN(did)%QLAKEI, RT_DOMAIN(did)%QLAKEO,& + RT_DOMAIN(did)%LAKENODE, rt_domain(did)%overland%properties%distance_to_neighbor, & + RT_DOMAIN(did)%QINFLOWBASE, RT_DOMAIN(did)%CHANXI, & + RT_DOMAIN(did)%CHANYJ, nlst(did)%channel_option, & + RT_DOMAIN(did)%RETDEP_CHAN, RT_DOMAIN(did)%NLINKSL, RT_DOMAIN(did)%LINKID, & + RT_DOMAIN(did)%node_area & #ifdef MPP_LAND - ,RT_DOMAIN(did)%lake_index,RT_DOMAIN(did)%link_location,& - RT_DOMAIN(did)%mpp_nlinks,RT_DOMAIN(did)%nlinks_index, & - RT_DOMAIN(did)%yw_mpp_nlinks & - , RT_DOMAIN(did)%LNLINKSL & - , rt_domain(did)%gtoNode,rt_domain(did)%toNodeInd,rt_domain(did)%nToInd & -#endif - , rt_domain(did)%CH_LNKRT_SL & - ,nlst(did)%GwBaseSwCRT, gw2d(did)%ho, gw2d(did)%qgw_chanrt, & - nlst(did)%gwChanCondSw, nlst(did)%gwChanCondConstIn, & - nlst(did)%gwChanCondConstOut, rt_domain(did)%velocity, rt_domain(did)%qloss & - ) -endif + ,RT_DOMAIN(did)%lake_index,RT_DOMAIN(did)%link_location,& + RT_DOMAIN(did)%mpp_nlinks,RT_DOMAIN(did)%nlinks_index, & + RT_DOMAIN(did)%yw_mpp_nlinks & + , RT_DOMAIN(did)%LNLINKSL & + , rt_domain(did)%gtoNode,rt_domain(did)%toNodeInd,rt_domain(did)%nToInd & +#endif + , rt_domain(did)%CH_LNKRT_SL & + ,nlst(did)%GwBaseSwCRT, gw2d(did)%ho, gw2d(did)%qgw_chanrt, & + nlst(did)%gwChanCondSw, nlst(did)%gwChanCondConstIn, & + nlst(did)%gwChanCondConstOut, rt_domain(did)%velocity, rt_domain(did)%qloss & + ) + endif - if((nlst(did)%gwBaseSwCRT == 3) .and. (nlst(did)%gwChanCondSw .eq. 1)) then + if((nlst(did)%gwBaseSwCRT == 3) .and. (nlst(did)%gwChanCondSw .eq. 1)) then - ! add/rm channel-aquifer exchange contribution +! add/rm channel-aquifer exchange contribution - gw2d(did)%ho = gw2d(did)%ho & - +(((gw2d(did)%qgw_chanrt*(-1)) * gw2d(did)%dt / gw2d(did)%dx**2) & - / gw2d(did)%poros) + gw2d(did)%ho = gw2d(did)%ho & + +(((gw2d(did)%qgw_chanrt*(-1)) * gw2d(did)%dt / gw2d(did)%dx**2) & + / gw2d(did)%poros) - endif - endif + endif + endif #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "*****yw******end drive_CHANNEL " + write(6,*) "*****yw******end drive_CHANNEL " #else - write(78,*) "*****yw******end drive_CHANNEL " + write(78,*) "*****yw******end drive_CHANNEL " #endif #endif - end subroutine driveChannelRouting + end subroutine driveChannelRouting !------------------------------------------------ - subroutine aggregateDomain(did) + subroutine aggregateDomain(did) #ifdef MPP_LAND use module_mpp_land, only: sum_real1, my_id, io_id, numprocs #endif - implicit none - integer, intent(in) :: did + implicit none + integer, intent(in) :: did - integer :: i, j, krt, ixxrt, jyyrt, & - AGGFACYRT, AGGFACXRT + integer :: i, j, krt, ixxrt, jyyrt, & + AGGFACYRT, AGGFACXRT #ifdef HYDRO_D ! ADCHANGE: Water balance variables - integer :: kk - real :: smcrttot1,smctot2,sicetot2 - real :: suminfxsrt1,suminfxs2 + integer :: kk + real :: smcrttot1,smctot2,sicetot2 + real :: suminfxsrt1,suminfxs2 #endif #ifdef HYDRO_D #ifndef NCEP_WCOSS - print *, "Beginning Aggregation..." + print *, "Beginning Aggregation..." #else - write(78,*) "Beginning Aggregation..." + write(78,*) "Beginning Aggregation..." #endif #endif @@ -1037,14 +1037,14 @@ subroutine aggregateDomain(did) suminfxsrt1 = 0. smcrttot1 = 0. do i=1,RT_DOMAIN(did)%IXRT - do j=1,RT_DOMAIN(did)%JXRT - suminfxsrt1 = suminfxsrt1 + rt_domain(did)%overland%control%surface_water_head_routing(I,J) & - / float(RT_DOMAIN(did)%IXRT * RT_DOMAIN(did)%JXRT) - do kk=1,nlst(did)%NSOIL - smcrttot1 = smcrttot1 + rt_domain(did)%subsurface%grid_transform%smcrt(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. & - / float(RT_DOMAIN(did)%IXRT * RT_DOMAIN(did)%JXRT) + do j=1,RT_DOMAIN(did)%JXRT + suminfxsrt1 = suminfxsrt1 + rt_domain(did)%overland%control%surface_water_head_routing(I,J) & + / float(RT_DOMAIN(did)%IXRT * RT_DOMAIN(did)%JXRT) + do kk=1,nlst(did)%NSOIL + smcrttot1 = smcrttot1 + rt_domain(did)%subsurface%grid_transform%smcrt(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. & + / float(RT_DOMAIN(did)%IXRT * RT_DOMAIN(did)%JXRT) + end do end do - end do end do #ifdef MPP_LAND ! not tested @@ -1057,26 +1057,26 @@ subroutine aggregateDomain(did) #endif do J=1,RT_DOMAIN(did)%JX - do I=1,RT_DOMAIN(did)%IX + do I=1,RT_DOMAIN(did)%IX - RT_DOMAIN(did)%SFCHEADAGGRT = 0. + RT_DOMAIN(did)%SFCHEADAGGRT = 0. !DJG Subgrid weighting edit... - RT_DOMAIN(did)%LSMVOL=0. - do KRT=1,nlst(did)%NSOIL + RT_DOMAIN(did)%LSMVOL=0. + do KRT=1,nlst(did)%NSOIL ! SMCAGGRT(KRT) = 0. - RT_DOMAIN(did)%SH2OAGGRT(KRT) = 0. - end do + RT_DOMAIN(did)%SH2OAGGRT(KRT) = 0. + end do - do AGGFACYRT=nlst(did)%AGGFACTRT-1,0,-1 - do AGGFACXRT=nlst(did)%AGGFACTRT-1,0,-1 + do AGGFACYRT=nlst(did)%AGGFACTRT-1,0,-1 + do AGGFACXRT=nlst(did)%AGGFACTRT-1,0,-1 - IXXRT=I*nlst(did)%AGGFACTRT-AGGFACXRT - JYYRT=J*nlst(did)%AGGFACTRT-AGGFACYRT + IXXRT=I*nlst(did)%AGGFACTRT-AGGFACXRT + JYYRT=J*nlst(did)%AGGFACTRT-AGGFACYRT #ifdef MPP_LAND - if(left_id.ge.0) IXXRT=IXXRT+1 - if(down_id.ge.0) JYYRT=JYYRT+1 + if(left_id.ge.0) IXXRT=IXXRT+1 + if(down_id.ge.0) JYYRT=JYYRT+1 #else !yw ???? ! IXXRT=IXXRT+1 @@ -1084,144 +1084,144 @@ subroutine aggregateDomain(did) #endif !State Variables - RT_DOMAIN(did)%SFCHEADAGGRT = RT_DOMAIN(did)%SFCHEADAGGRT & - + rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) + RT_DOMAIN(did)%SFCHEADAGGRT = RT_DOMAIN(did)%SFCHEADAGGRT & + + rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) !DJG Subgrid weighting edit... - RT_DOMAIN(did)%LSMVOL = RT_DOMAIN(did)%LSMVOL & - + rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) & - * rt_domain(did)%overland%properties%distance_to_neighbor(IXXRT,JYYRT,9) + RT_DOMAIN(did)%LSMVOL = RT_DOMAIN(did)%LSMVOL & + + rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) & + * rt_domain(did)%overland%properties%distance_to_neighbor(IXXRT,JYYRT,9) - do KRT=1,nlst(did)%NSOIL + do KRT=1,nlst(did)%NSOIL !DJG SMCAGGRT(KRT)=SMCAGGRT(KRT)+SMCRT(IXXRT,JYYRT,KRT) - RT_DOMAIN(did)%SH2OAGGRT(KRT) = RT_DOMAIN(did)%SH2OAGGRT(KRT) & - + rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) - end do + RT_DOMAIN(did)%SH2OAGGRT(KRT) = RT_DOMAIN(did)%SH2OAGGRT(KRT) & + + rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) + end do - end do - end do + end do + end do - rt_domain(did)%overland%control%surface_water_head_lsm(I,J) = RT_DOMAIN(did)%SFCHEADAGGRT & - / (nlst(did)%AGGFACTRT**2) + rt_domain(did)%overland%control%surface_water_head_lsm(I,J) = RT_DOMAIN(did)%SFCHEADAGGRT & + / (nlst(did)%AGGFACTRT**2) - do KRT=1,nlst(did)%NSOIL + do KRT=1,nlst(did)%NSOIL !DJG SMC(I,J,KRT)=SMCAGGRT(KRT)/(AGGFACTRT**2) - RT_DOMAIN(did)%SH2OX(I,J,KRT) = RT_DOMAIN(did)%SH2OAGGRT(KRT) & - / (nlst(did)%AGGFACTRT**2) - end do + RT_DOMAIN(did)%SH2OX(I,J,KRT) = RT_DOMAIN(did)%SH2OAGGRT(KRT) & + / (nlst(did)%AGGFACTRT**2) + end do !DJG Calculate subgrid weighting array... - do AGGFACYRT=nlst(did)%AGGFACTRT-1,0,-1 - do AGGFACXRT=nlst(did)%AGGFACTRT-1,0,-1 - IXXRT=I*nlst(did)%AGGFACTRT-AGGFACXRT - JYYRT=J*nlst(did)%AGGFACTRT-AGGFACYRT + do AGGFACYRT=nlst(did)%AGGFACTRT-1,0,-1 + do AGGFACXRT=nlst(did)%AGGFACTRT-1,0,-1 + IXXRT=I*nlst(did)%AGGFACTRT-AGGFACXRT + JYYRT=J*nlst(did)%AGGFACTRT-AGGFACYRT #ifdef MPP_LAND - if(left_id.ge.0) IXXRT=IXXRT+1 - if(down_id.ge.0) JYYRT=JYYRT+1 + if(left_id.ge.0) IXXRT=IXXRT+1 + if(down_id.ge.0) JYYRT=JYYRT+1 #else !yw ??? ! IXXRT=IXXRT+1 ! JYYRT=JYYRT+1 #endif - if (RT_DOMAIN(did)%LSMVOL.gt.0.) then - RT_DOMAIN(did)%INFXSWGT(IXXRT,JYYRT) & - = rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) & - * rt_domain(did)%overland%properties%distance_to_neighbor(IXXRT,JYYRT,9) & - / RT_DOMAIN(did)%LSMVOL - else - RT_DOMAIN(did)%INFXSWGT(IXXRT,JYYRT) & - = 1./FLOAT(nlst(did)%AGGFACTRT**2) - end if + if (RT_DOMAIN(did)%LSMVOL.gt.0.) then + RT_DOMAIN(did)%INFXSWGT(IXXRT,JYYRT) & + = rt_domain(did)%overland%control%surface_water_head_routing(IXXRT,JYYRT) & + * rt_domain(did)%overland%properties%distance_to_neighbor(IXXRT,JYYRT,9) & + / RT_DOMAIN(did)%LSMVOL + else + RT_DOMAIN(did)%INFXSWGT(IXXRT,JYYRT) & + = 1./FLOAT(nlst(did)%AGGFACTRT**2) + end if - do KRT=1,nlst(did)%NSOIL + do KRT=1,nlst(did)%NSOIL !!!yw added for debug - if(rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) .lt. 0) then + if(rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) .lt. 0) then #ifndef NCEP_WCOSS - 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) + 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) #else - 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) + 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) #endif - endif - if(RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) .lt. 0) then + endif + if(RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) .lt. 0) then #ifndef NCEP_WCOSS - 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) + 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) #else - 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) + 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) #endif - endif + endif - IF ( (rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) - & - rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT)) .GT. 0.000001 ) THEN + IF ( (rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) - & + rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT)) .GT. 0.000001 ) THEN #ifdef HYDRO_D #ifndef NCEP_WCOSS - print *, "SMCMAX exceeded upon aggregation...", & - rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT), & - rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT) + print *, "SMCMAX exceeded upon aggregation...", & + rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT), & + rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT) #else - write(78,*) "FATAL ERROR: SMCMAX exceeded upon aggregation...", & - rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT), & - rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT) + write(78,*) "FATAL ERROR: SMCMAX exceeded upon aggregation...", & + rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT), & + rt_domain(did)%subsurface%grid_transform%smcmaxrt(IXXRT,JYYRT,KRT) #endif #endif - call hydro_stop("In module_HYDRO_drv.F aggregateDomain() - "// & - "SMCMAX exceeded upon aggregation.") - END IF - IF(RT_DOMAIN(did)%SH2OX(I,J,KRT).LT.0.) THEN + call hydro_stop("In module_HYDRO_drv.F aggregateDomain() - "// & + "SMCMAX exceeded upon aggregation.") + END IF + IF(RT_DOMAIN(did)%SH2OX(I,J,KRT).LT.0.) THEN #ifdef HYDRO_D #ifndef NCEP_WCOSS - print *, "Erroneous value of SH2O...", & - RT_DOMAIN(did)%SH2OX(I,J,KRT),I,J,KRT - 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) + print *, "Erroneous value of SH2O...", & + RT_DOMAIN(did)%SH2OX(I,J,KRT),I,J,KRT + 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) #else - write(78,*) "Erroneous value of SH2O...", & - RT_DOMAIN(did)%SH2OX(I,J,KRT),I,J,KRT - 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) + write(78,*) "Erroneous value of SH2O...", & + RT_DOMAIN(did)%SH2OX(I,J,KRT),I,J,KRT + 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) #endif #endif - call hydro_stop("In module_HYDRO_drv.F aggregateDomain() "// & - "- Error negative SH2OX") - END IF + call hydro_stop("In module_HYDRO_drv.F aggregateDomain() "// & + "- Error negative SH2OX") + END IF - IF ( RT_DOMAIN(did)%SH2OX(I,J,KRT) .gt. 0 ) THEN - RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) & - = rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) & - / RT_DOMAIN(did)%SH2OX(I,J,KRT) - ELSE + IF ( RT_DOMAIN(did)%SH2OX(I,J,KRT) .gt. 0 ) THEN + RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) & + = rt_domain(did)%subsurface%grid_transform%smcrt(IXXRT,JYYRT,KRT) & + / RT_DOMAIN(did)%SH2OX(I,J,KRT) + ELSE #ifdef HYDRO_D - 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) + 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) #endif - RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) = 0.0 - ENDIF + RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) = 0.0 + ENDIF !?yw - RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) = max(1.0E-05, RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT)) - end do + RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT) = max(1.0E-05, RT_DOMAIN(did)%SH2OWGT(IXXRT,JYYRT,KRT)) + end do + end do end do - end do - end do + end do end do #ifdef MPP_LAND call MPP_LAND_COM_REAL(RT_DOMAIN(did)%INFXSWGT, & - RT_DOMAIN(did)%IXRT, & - RT_DOMAIN(did)%JXRT, 99) + RT_DOMAIN(did)%IXRT, & + RT_DOMAIN(did)%JXRT, 99) do i = 1, nlst(did)%NSOIL - call MPP_LAND_COM_REAL(RT_DOMAIN(did)%SH2OWGT(:,:,i), & - RT_DOMAIN(did)%IXRT, & - RT_DOMAIN(did)%JXRT, 99) + call MPP_LAND_COM_REAL(RT_DOMAIN(did)%SH2OWGT(:,:,i), & + RT_DOMAIN(did)%IXRT, & + RT_DOMAIN(did)%JXRT, 99) end do #endif !DJG Update SMC with SICE (unchanged) and new value of SH2O from routing... - RT_DOMAIN(did)%SMC = RT_DOMAIN(did)%SH2OX + RT_DOMAIN(did)%SICE + RT_DOMAIN(did)%SMC = RT_DOMAIN(did)%SH2OX + RT_DOMAIN(did)%SICE #ifdef HYDRO_D ! ADCHANGE: START Final water balance variables @@ -1230,16 +1230,16 @@ subroutine aggregateDomain(did) smctot2 = 0. sicetot2 = 0. do i=1,RT_DOMAIN(did)%IX - do j=1,RT_DOMAIN(did)%JX - suminfxs2 = suminfxs2 + rt_domain(did)%overland%control%surface_water_head_lsm(I,J) & - / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX) - do kk=1,nlst(did)%NSOIL - smctot2 = smctot2 + rt_domain(did)%SMC(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. & - / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX) - sicetot2 = sicetot2 + rt_domain(did)%SICE(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. & - / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX) + do j=1,RT_DOMAIN(did)%JX + suminfxs2 = suminfxs2 + rt_domain(did)%overland%control%surface_water_head_lsm(I,J) & + / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX) + do kk=1,nlst(did)%NSOIL + smctot2 = smctot2 + rt_domain(did)%SMC(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. & + / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX) + sicetot2 = sicetot2 + rt_domain(did)%SICE(I,J,KK)*RT_DOMAIN(did)%subsurface%properties%sldpth(KK)*1000. & + / float(RT_DOMAIN(did)%IX * RT_DOMAIN(did)%JX) + end do end do - end do end do #ifdef MPP_LAND @@ -1253,104 +1253,104 @@ subroutine aggregateDomain(did) #endif #ifdef MPP_LAND - if (my_id .eq. IO_id) then -#endif - print *, "Agg Mass Bal: " - print *, "WB_AGG!InfxsDiff", suminfxs2-suminfxsrt1 - print *, "WB_AGG!Infxs1", suminfxsrt1 - print *, "WB_AGG!Infxs2", suminfxs2 - print *, "WB_AGG!SMCDiff", smctot2-smcrttot1-sicetot2 - print *, "WB_AGG!SMC1", smcrttot1 - print *, "WB_AGG!SMC2", smctot2 - print *, "WB_AGG!SICE2", sicetot2 - print *, "WB_AGG!Residual", (suminfxs2-suminfxsrt1) + & - (smctot2-smcrttot1-sicetot2) + if (my_id .eq. IO_id) then +#endif + print *, "Agg Mass Bal: " + print *, "WB_AGG!InfxsDiff", suminfxs2-suminfxsrt1 + print *, "WB_AGG!Infxs1", suminfxsrt1 + print *, "WB_AGG!Infxs2", suminfxs2 + print *, "WB_AGG!SMCDiff", smctot2-smcrttot1-sicetot2 + print *, "WB_AGG!SMC1", smcrttot1 + print *, "WB_AGG!SMC2", smctot2 + print *, "WB_AGG!SICE2", sicetot2 + print *, "WB_AGG!Residual", (suminfxs2-suminfxsrt1) + & + (smctot2-smcrttot1-sicetot2) #ifdef MPP_LAND - endif + endif #endif ! END Final water balance variables #endif #ifdef HYDRO_D #ifndef NCEP_WCOSS - print *, "Finished Aggregation..." + print *, "Finished Aggregation..." #else - write(78,*) "Finished Aggregation..." + write(78,*) "Finished Aggregation..." #endif #endif - end subroutine aggregateDomain + end subroutine aggregateDomain - subroutine RunOffDisag(runoff1x_in, runoff1x, area_lsm,cellArea, infxswgt, AGGFACTRT, ix,jx) + subroutine RunOffDisag(runoff1x_in, runoff1x, area_lsm,cellArea, infxswgt, AGGFACTRT, ix,jx) implicit none real, dimension(:,:) :: runoff1x_in, runoff1x, area_lsm, cellArea, infxswgt integer :: i,j,ix,jx,AGGFACYRT, AGGFACXRT, AGGFACTRT, IXXRT, JYYRT do J=1,JX - do I=1,IX - do AGGFACYRT=AGGFACTRT-1,0,-1 - do AGGFACXRT=AGGFACTRT-1,0,-1 - IXXRT=I*AGGFACTRT-AGGFACXRT - JYYRT=J*AGGFACTRT-AGGFACYRT + do I=1,IX + do AGGFACYRT=AGGFACTRT-1,0,-1 + do AGGFACXRT=AGGFACTRT-1,0,-1 + IXXRT=I*AGGFACTRT-AGGFACXRT + JYYRT=J*AGGFACTRT-AGGFACYRT #ifdef MPP_LAND - if(left_id.ge.0) IXXRT=IXXRT+1 - if(down_id.ge.0) JYYRT=JYYRT+1 + if(left_id.ge.0) IXXRT=IXXRT+1 + if(down_id.ge.0) JYYRT=JYYRT+1 #endif !DJG Implement subgrid weighting routine... - if( (runoff1x_in(i,j) .lt. 0) .or. (runoff1x_in(i,j) .gt. 1000) ) then - runoff1x(IXXRT,JYYRT) = 0 - else - runoff1x(IXXRT,JYYRT)=runoff1x_in(i,j)*area_lsm(I,J) & - *INFXSWGT(IXXRT,JYYRT)/cellArea(IXXRT,JYYRT) - endif - - enddo - enddo - enddo + if( (runoff1x_in(i,j) .lt. 0) .or. (runoff1x_in(i,j) .gt. 1000) ) then + runoff1x(IXXRT,JYYRT) = 0 + else + runoff1x(IXXRT,JYYRT)=runoff1x_in(i,j)*area_lsm(I,J) & + *INFXSWGT(IXXRT,JYYRT)/cellArea(IXXRT,JYYRT) + endif + + enddo + enddo + enddo enddo - end subroutine RunOffDisag + end subroutine RunOffDisag ! This routine was extracted from the aggregateDomain routine above to do simple depth aggregation. ! There might be a simpler way. -subroutine RunoffAggregate(runoff_in, runoff_out, aggfactrt, ix, jx) - implicit none - ! Input variables - integer, intent(in) :: aggfactrt, ix, jx - real, intent(in), dimension(:,:) :: runoff_in - real, intent(inout), dimension(:,:) :: runoff_out - ! Local variables - integer :: i, j, aggfacyrt, aggfacxrt, ixxrt, jyyrt - real :: runoffagg - do j=1,jx - do i=1,ix - runoffagg = 0. - do aggfacyrt=aggfactrt-1,0,-1 - do aggfacxrt=aggfactrt-1,0,-1 - ixxrt = i * aggfactrt - aggfacxrt - jyyrt = j * aggfactrt - aggfacyrt + subroutine RunoffAggregate(runoff_in, runoff_out, aggfactrt, ix, jx) + implicit none +! Input variables + integer, intent(in) :: aggfactrt, ix, jx + real, intent(in), dimension(:,:) :: runoff_in + real, intent(inout), dimension(:,:) :: runoff_out +! Local variables + integer :: i, j, aggfacyrt, aggfacxrt, ixxrt, jyyrt + real :: runoffagg + do j=1,jx + do i=1,ix + runoffagg = 0. + do aggfacyrt=aggfactrt-1,0,-1 + do aggfacxrt=aggfactrt-1,0,-1 + ixxrt = i * aggfactrt - aggfacxrt + jyyrt = j * aggfactrt - aggfacyrt #ifdef MPP_LAND - if(left_id.ge.0) ixxrt = ixxrt+1 - if(down_id.ge.0) jyyrt = jyyrt+1 -#endif - runoffagg = runoffagg + runoff_in(ixxrt,jyyrt) - end do - end do - runoff_out(i,j) = runoffagg / (aggfactrt**2) - end do -end do -end subroutine RunoffAggregate - -subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp) -implicit none -integer ntime, did -integer rst_out, ix,jx + if(left_id.ge.0) ixxrt = ixxrt+1 + if(down_id.ge.0) jyyrt = jyyrt+1 +#endif + runoffagg = runoffagg + runoff_in(ixxrt,jyyrt) + end do + end do + runoff_out(i,j) = runoffagg / (aggfactrt**2) + end do + end do + end subroutine RunoffAggregate + + subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp) + implicit none + integer ntime, did + integer rst_out, ix,jx ! integer, OPTIONAL:: ix0,jx0 -integer:: ix0,jx0 -integer, dimension(ix0,jx0),optional :: vegtyp, soltyp -integer :: iret, ncid, ascIndId + integer:: ix0,jx0 + integer, dimension(ix0,jx0),optional :: vegtyp, soltyp + integer :: iret, ncid, ascIndId @@ -1359,120 +1359,120 @@ subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp) !call read_rt_nlst(nlst(did) ) ! Some field of this structure are already initialized by the CPL component (e.g. DT) -call orchestrator%config%init_nlst(did) + call orchestrator%config%init_nlst(did) -if(nlst(did)%rtFlag .eq. 0) return + if(nlst(did)%rtFlag .eq. 0) return !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! get the dimension -call get_file_dimension(trim(nlst(did)%geo_static_flnm), ix,jx) + call get_file_dimension(trim(nlst(did)%geo_static_flnm), ix,jx) #ifdef MPP_LAND -if (nlst(did)%sys_cpl .eq. 1 .or. nlst(did)%sys_cpl .eq. 4) then - !sys_cpl: 1-- coupling with HRLDAS but running offline lsm; - ! 2-- coupling with WRF but do not run offline lsm - ! 3-- coupling with LIS and do not run offline lsm - ! 4: coupling with CLM + if (nlst(did)%sys_cpl .eq. 1 .or. nlst(did)%sys_cpl .eq. 4) then +!sys_cpl: 1-- coupling with HRLDAS but running offline lsm; +! 2-- coupling with WRF but do not run offline lsm +! 3-- coupling with LIS and do not run offline lsm +! 4: coupling with CLM - ! create 2 dimensiaon logical mapping of the CPUs for coupling with CLM or HRLDAS. - call log_map2d() +! create 2 dimensiaon logical mapping of the CPUs for coupling with CLM or HRLDAS. + call log_map2d() - global_nx = ix ! get from land model - global_ny = jx ! get from land model + global_nx = ix ! get from land model + global_ny = jx ! get from land model - call mpp_land_bcast_int1(global_nx) - call mpp_land_bcast_int1(global_ny) + call mpp_land_bcast_int1(global_nx) + call mpp_land_bcast_int1(global_ny) !!! temp set global_nx to ix - rt_domain(did)%ix = global_nx - rt_domain(did)%jx = global_ny + rt_domain(did)%ix = global_nx + rt_domain(did)%jx = global_ny - ! over write the ix and jx - call MPP_LAND_PAR_INI(1,rt_domain(did)%ix,rt_domain(did)%jx,& - nlst(did)%AGGFACTRT) -else - ! coupled with WRF, LIS - numprocs = node_info(1,1) +! over write the ix and jx + call MPP_LAND_PAR_INI(1,rt_domain(did)%ix,rt_domain(did)%jx,& + nlst(did)%AGGFACTRT) + else +! coupled with WRF, LIS + numprocs = node_info(1,1) - call wrf_LAND_set_INIT(node_info,numprocs,nlst(did)%AGGFACTRT) + call wrf_LAND_set_INIT(node_info,numprocs,nlst(did)%AGGFACTRT) - rt_domain(did)%ix = local_nx - rt_domain(did)%jx = local_ny -endif + rt_domain(did)%ix = local_nx + rt_domain(did)%jx = local_ny + endif -rt_domain(did)%g_IXRT=global_rt_nx -rt_domain(did)%g_JXRT=global_rt_ny -rt_domain(did)%ixrt = local_rt_nx -rt_domain(did)%jxrt = local_rt_ny + rt_domain(did)%g_IXRT=global_rt_nx + rt_domain(did)%g_JXRT=global_rt_ny + rt_domain(did)%ixrt = local_rt_nx + rt_domain(did)%jxrt = local_rt_ny #ifdef HYDRO_D #ifndef NCEP_WCOSS -write(6,*) "rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt" -write(6,*) rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt -write(6,*) "rt_domain(did)%ix, rt_domain(did)%jx " -write(6,*) rt_domain(did)%ix, rt_domain(did)%jx -write(6,*) "global_nx, global_ny, local_nx, local_ny" -write(6,*) global_nx, global_ny, local_nx, local_ny + write(6,*) "rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt" + write(6,*) rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt + write(6,*) "rt_domain(did)%ix, rt_domain(did)%jx " + write(6,*) rt_domain(did)%ix, rt_domain(did)%jx + write(6,*) "global_nx, global_ny, local_nx, local_ny" + write(6,*) global_nx, global_ny, local_nx, local_ny #else -write(78,*) "rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt" -write(78,*) rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt -write(78,*) "rt_domain(did)%ix, rt_domain(did)%jx " -write(78,*) rt_domain(did)%ix, rt_domain(did)%jx -write(78,*) "global_nx, global_ny, local_nx, local_ny" -write(78,*) global_nx, global_ny, local_nx, local_ny + write(78,*) "rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt" + write(78,*) rt_domain(did)%g_IXRT, rt_domain(did)%g_JXRT, rt_domain(did)%ixrt, rt_domain(did)%jxrt + write(78,*) "rt_domain(did)%ix, rt_domain(did)%jx " + write(78,*) rt_domain(did)%ix, rt_domain(did)%jx + write(78,*) "global_nx, global_ny, local_nx, local_ny" + write(78,*) global_nx, global_ny, local_nx, local_ny #endif #endif #else ! sequential -rt_domain(did)%ix = ix -rt_domain(did)%jx = jx -rt_domain(did)%ixrt = ix*nlst(did)%AGGFACTRT -rt_domain(did)%jxrt = jx*nlst(did)%AGGFACTRT + rt_domain(did)%ix = ix + rt_domain(did)%jx = jx + rt_domain(did)%ixrt = ix*nlst(did)%AGGFACTRT + rt_domain(did)%jxrt = jx*nlst(did)%AGGFACTRT #endif ! allocate rt arrays -call getChanDim(did) + call getChanDim(did) #ifdef HYDRO_D #ifndef NCEP_WCOSS -write(6,*) "finish getChanDim " + write(6,*) "finish getChanDim " #else -write(78,*) "finish getChanDim " + write(78,*) "finish getChanDim " #endif #endif ! ADCHANGE: get global attributes ! need to set these after getChanDim since it allocates rt_domain vals to 0 - call get_file_globalatts(trim(nlst(did)%geo_static_flnm), & - rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater) + call get_file_globalatts(trim(nlst(did)%geo_static_flnm), & + rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater) #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "hydro_ini: rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater" - write(6,*) rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater + write(6,*) "hydro_ini: rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater" + write(6,*) rt_domain(did)%iswater, rt_domain(did)%islake, rt_domain(did)%isurban, rt_domain(did)%isoilwater #endif #endif -if(nlst(did)%GWBASESWCRT .eq. 3 ) then - call gw2d_allocate(did,& - rt_domain(did)%ixrt,& - rt_domain(did)%jxrt,& - nlst(did)%nsoil) + if(nlst(did)%GWBASESWCRT .eq. 3 ) then + call gw2d_allocate(did,& + rt_domain(did)%ixrt,& + rt_domain(did)%jxrt,& + nlst(did)%nsoil) #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "finish gw2d_allocate" + write(6,*) "finish gw2d_allocate" #else - write(78,*) "finish gw2d_allocate" + write(78,*) "finish gw2d_allocate" #endif #endif -endif + endif ! calculate the distance between grids for routing. ! decompose the land parameter/data @@ -1480,104 +1480,104 @@ subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp) ! ix0= rt_domain(did)%ix ! jx0= rt_domain(did)%jx -if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then - if(present(vegtyp)) then - call lsm_input(did,ix0=ix0,jx0=jx0,vegtyp0=vegtyp,soltyp0=soltyp) - else - call lsm_input(did,ix0=ix0,jx0=jx0) - endif -endif + if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then + if(present(vegtyp)) then + call lsm_input(did,ix0=ix0,jx0=jx0,vegtyp0=vegtyp,soltyp0=soltyp) + else + call lsm_input(did,ix0=ix0,jx0=jx0) + endif + endif #ifdef HYDRO_D #ifndef NCEP_WCOSS -write(6,*) "finish decomposion" + write(6,*) "finish decomposion" #else -write(78,*) "finish decomposion" -#endif -#endif - -if((nlst(did)%channel_only .eq. 1 .or. nlst(did)%channelBucket_only .eq. 1) .and. & - nlst(1)%io_form_outputs .ne. 0) then - !! This is the "decoder ring" for reading channel-only forcing from io_form_outputs=1,2 CHRTOUT files. - !! Only needed on io_id - if(my_id .eq. io_id) then - allocate(rt_domain(did)%ascendIndex(rt_domain(did)%gnlinksl)) - iret = nf90_open(trim(nlst(1)%route_link_f),NF90_NOWRITE,ncid=ncid) - !if(iret .ne. 0) call hdyro_stop - if(iret .ne. 0) call hydro_stop('ERROR: Unable to open RouteLink file for index extraction') - iret = nf90_inq_varid(ncid,'ascendingIndex',ascIndId) - if(iret .ne. 0) call hydro_stop('ERROR: Unable to find ascendingIndex from RouteLink file.') - iret = nf90_get_var(ncid,ascIndId,rt_domain(did)%ascendIndex) - if(iret .ne. 0) call hydro_stop('ERROR: Unable to extract ascendingIndex from RouteLink file.') - iret = nf90_close(ncid) - if(iret .ne. 0) call hydro_stop('ERROR: Unable to close RouteLink file.') - else - allocate(rt_domain(did)%ascendIndex(1)) - rt_domain(did)%ascendIndex(1)=-9 - endif -endif - - -call get_dist_lsm(did) !! always needed (channel_only and channelBucket_only) -if(nlst(did)%channel_only .ne. 1) call get_dist_lrt(did) !! needed forchannelBucket_only + write(78,*) "finish decomposion" +#endif +#endif + + if((nlst(did)%channel_only .eq. 1 .or. nlst(did)%channelBucket_only .eq. 1) .and. & + nlst(1)%io_form_outputs .ne. 0) then +!! This is the "decoder ring" for reading channel-only forcing from io_form_outputs=1,2 CHRTOUT files. +!! Only needed on io_id + if(my_id .eq. io_id) then + allocate(rt_domain(did)%ascendIndex(rt_domain(did)%gnlinksl)) + iret = nf90_open(trim(nlst(1)%route_link_f),NF90_NOWRITE,ncid=ncid) +!if(iret .ne. 0) call hdyro_stop + if(iret .ne. 0) call hydro_stop('ERROR: Unable to open RouteLink file for index extraction') + iret = nf90_inq_varid(ncid,'ascendingIndex',ascIndId) + if(iret .ne. 0) call hydro_stop('ERROR: Unable to find ascendingIndex from RouteLink file.') + iret = nf90_get_var(ncid,ascIndId,rt_domain(did)%ascendIndex) + if(iret .ne. 0) call hydro_stop('ERROR: Unable to extract ascendingIndex from RouteLink file.') + iret = nf90_close(ncid) + if(iret .ne. 0) call hydro_stop('ERROR: Unable to close RouteLink file.') + else + allocate(rt_domain(did)%ascendIndex(1)) + rt_domain(did)%ascendIndex(1)=-9 + endif + endif + + + call get_dist_lsm(did) !! always needed (channel_only and channelBucket_only) + if(nlst(did)%channel_only .ne. 1) call get_dist_lrt(did) !! needed forchannelBucket_only ! rt model initilization -call LandRT_ini(did) + call LandRT_ini(did) -if(nlst(did)%GWBASESWCRT .eq. 3 ) then + if(nlst(did)%GWBASESWCRT .eq. 3 ) then - call gw2d_ini(did,& - nlst(did)%dt,& - nlst(did)%dxrt0) + call gw2d_ini(did,& + nlst(did)%dt,& + nlst(did)%dxrt0) #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) "finish gw2d_ini" + write(6,*) "finish gw2d_ini" #else - write(78,*) "finish gw2d_ini" + write(78,*) "finish gw2d_ini" #endif #endif -endif + endif #ifdef HYDRO_D #ifndef NCEP_WCOSS -write(6,*) "finish LandRT_ini" + write(6,*) "finish LandRT_ini" #else -write(78,*) "finish LandRT_ini" + write(78,*) "finish LandRT_ini" #endif #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then + if(nlst(did)%channel_only .eq. 0 .and. nlst(did)%channelBucket_only .eq. 0) then - if (nlst(did)%TERADJ_SOLAR.eq.1 .and. nlst(did)%CHANRTSWCRT.ne.2) then ! Perform ter rain adjustment of incoming solar + if (nlst(did)%TERADJ_SOLAR.eq.1 .and. nlst(did)%CHANRTSWCRT.ne.2) then ! Perform ter rain adjustment of incoming solar #ifdef MPP_LAND - call MPP_seq_land_SO8(rt_domain(did)%SO8LD_D,rt_domain(did)%SO8LD_Vmax,& - rt_domain(did)%TERRAIN, rt_domain(did)%dist_lsm,& - rt_domain(did)%ix,rt_domain(did)%jx,global_nx,global_ny) + call MPP_seq_land_SO8(rt_domain(did)%SO8LD_D,rt_domain(did)%SO8LD_Vmax,& + rt_domain(did)%TERRAIN, rt_domain(did)%dist_lsm,& + rt_domain(did)%ix,rt_domain(did)%jx,global_nx,global_ny) #else - call seq_land_SO8(rt_domain(did)%SO8LD_D,rt_domain(did)%SO8LD_Vmax,& - rt_domain(did)%TERRAIN,rt_domain(did)%dist_lsm,& - rt_domain(did)%ix,rt_domain(did)%jx) + call seq_land_SO8(rt_domain(did)%SO8LD_D,rt_domain(did)%SO8LD_Vmax,& + rt_domain(did)%TERRAIN,rt_domain(did)%dist_lsm,& + rt_domain(did)%ix,rt_domain(did)%jx) #endif - endif -endif + endif + endif -if (nlst(did)%GWBASESWCRT .gt. 0) then - if(nlst(did)%UDMP_OPT .eq. 1) then - call get_basn_area_nhd(rt_domain(did)%basns_area) - else - call get_basn_area(did) - endif -endif + if (nlst(did)%GWBASESWCRT .gt. 0) then + if(nlst(did)%UDMP_OPT .eq. 1) then + call get_basn_area_nhd(rt_domain(did)%basns_area) + else + call get_basn_area(did) + endif + endif -if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT .eq. 2 ) then - call get_node_area(did) -endif + if (nlst(did)%CHANRTSWCRT.eq.1 .or. nlst(did)%CHANRTSWCRT .eq. 2 ) then + call get_node_area(did) + endif #ifdef WRF_HYDRO_NUDGING -if(nlst(did)%CHANRTSWCRT .ne. 0) call init_stream_nudging + if(nlst(did)%CHANRTSWCRT .ne. 0) call init_stream_nudging #endif @@ -1589,18 +1589,18 @@ subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp) ! restart the file - ! jummp the initial time output +! jummp the initial time output ! rt_domain(did)%out_counts = rt_domain(did)%out_counts + 1 ! rt_domain(did)%his_out_counts = rt_domain(did)%his_out_counts + 1 -call HYDRO_rst_in(did) + call HYDRO_rst_in(did) !#ifdef HYDRO_REALTIME -if (trim(nlst(did)%restart_file) == "") then - call HYDRO_out(did, 0) -else - call HYDRO_out(did, 1) -endif + if (trim(nlst(did)%restart_file) == "") then + call HYDRO_out(did, 0) + else + call HYDRO_out(did, 1) + endif !! JLM: This is only currently part 1/2 of a better accumulation tracking strategy. !! The parts: !! 1) (this) zero accumulations on restart/init after any t=0 outputs are written. @@ -1610,25 +1610,25 @@ subroutine HYDRO_ini(ntime, did,ix0,jx0, vegtyp,soltyp) !! This was previously done in HYDRO_rst_in and so output accumulations at time !! zero were getting zeroed and then writtent to file, which looses information. !! Note that nlst_rt(did)%rstrt_swc is not changed at any point in between here and the rst_in. -if(nlst(did)%rstrt_swc.eq.1) then !Switch for rest of restart accum vars... - print *, "Resetting RESTART Accumulation Variables to 0...",nlst(did)%rstrt_swc - ! Under channel-only , these first three variables are not allocated. - 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 - if(allocated(rt_domain(did)%QSTRMVOLRT_ACC)) rt_domain(did)%QSTRMVOLRT_ACC = zeroFlt - ! These variables are optionally allocated, if their output is requested. - if(allocated(rt_domain(did)%accSfcLatRunoff)) rt_domain(did)%accSfcLatRunoff = zeroDbl - if(allocated(rt_domain(did)%accBucket)) rt_domain(did)%accBucket = zeroDbl -end if - -end subroutine HYDRO_ini - - subroutine lsm_input(did,ix0,jx0,vegtyp0,soltyp0) - implicit none - integer did, leng, ncid, ierr_flg - parameter(leng=100) - integer :: i,j, nn - integer, allocatable, dimension(:,:) :: soltyp - real, dimension(leng) :: xdum1, MAXSMC,refsmc,wltsmc + if(nlst(did)%rstrt_swc.eq.1) then !Switch for reset of restart accum vars... + print *, "Resetting RESTART Accumulation Variables to 0...",nlst(did)%rstrt_swc +! Under channel-only , these first three variables are not allocated. + 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 + if(allocated(rt_domain(did)%QSTRMVOLRT_ACC)) rt_domain(did)%QSTRMVOLRT_ACC = zeroFlt +! These variables are optionally allocated, if their output is requested. + if(allocated(rt_domain(did)%accSfcLatRunoff)) rt_domain(did)%accSfcLatRunoff = zeroDbl + if(allocated(rt_domain(did)%accBucket)) rt_domain(did)%accBucket = zeroDbl + end if + + end subroutine HYDRO_ini + + subroutine lsm_input(did,ix0,jx0,vegtyp0,soltyp0) + implicit none + integer did, leng, ncid, ierr_flg + parameter(leng=100) + integer :: i,j, nn + integer, allocatable, dimension(:,:) :: soltyp + real, dimension(leng) :: xdum1, MAXSMC,refsmc,wltsmc integer :: ix0,jx0 integer, dimension(ix0,jx0),OPTIONAL :: vegtyp0, soltyp0 @@ -1636,200 +1636,194 @@ subroutine lsm_input(did,ix0,jx0,vegtyp0,soltyp0) #ifdef HYDRO_D #ifndef NCEP_WCOSS - write(6,*) RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx + write(6,*) RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx #else - write(78,*) RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx + write(78,*) RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx #endif #endif - allocate(soltyp(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx) ) + allocate(soltyp(RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx) ) - soltyp = 0 - call get2d_lsm_soltyp(soltyp,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) + soltyp = 0 + call get2d_lsm_soltyp(soltyp,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) - call get2d_lsm_real("HGT",RT_DOMAIN(did)%TERRAIN,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) + call get2d_lsm_real("HGT",RT_DOMAIN(did)%TERRAIN,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) - call get2d_lsm_real("XLAT",RT_DOMAIN(did)%lat_lsm,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) - call get2d_lsm_real("XLONG",RT_DOMAIN(did)%lon_lsm,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) - call get2d_lsm_vegtyp(RT_DOMAIN(did)%VEGTYP,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) + call get2d_lsm_real("XLAT",RT_DOMAIN(did)%lat_lsm,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) + call get2d_lsm_real("XLONG",RT_DOMAIN(did)%lon_lsm,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) + call get2d_lsm_vegtyp(RT_DOMAIN(did)%VEGTYP,RT_DOMAIN(did)%ix,RT_DOMAIN(did)%jx,trim(nlst(did)%geo_static_flnm)) - if(nlst(did)%sys_cpl .eq. 2 ) then - ! coupling with WRF - if(present(soltyp0) ) then - where(VEGTYP0 == rt_domain(did)%iswater .or. VEGTYP0 == rt_domain(did)%islake) soltyp0 = rt_domain(did)%isoilwater - where(soltyp0 == rt_domain(did)%isoilwater) VEGTYP0 = rt_domain(did)%iswater - soltyp = soltyp0 - RT_DOMAIN(did)%VEGTYP = VEGTYP0 - endif + if(nlst(did)%sys_cpl .eq. 2 ) then +! coupling with WRF + if(present(soltyp0) ) then + where(VEGTYP0 == rt_domain(did)%iswater .or. VEGTYP0 == rt_domain(did)%islake) soltyp0 = rt_domain(did)%isoilwater + where(soltyp0 == rt_domain(did)%isoilwater) VEGTYP0 = rt_domain(did)%iswater + soltyp = soltyp0 + RT_DOMAIN(did)%VEGTYP = VEGTYP0 endif + endif - where(RT_DOMAIN(did)%VEGTYP == rt_domain(did)%iswater .or. RT_DOMAIN(did)%VEGTYP == rt_domain(did)%islake) soltyp = rt_domain(did)%isoilwater - where(soltyp == rt_domain(did)%isoilwater) RT_DOMAIN(did)%VEGTYP = rt_domain(did)%iswater + where(RT_DOMAIN(did)%VEGTYP == rt_domain(did)%iswater .or. RT_DOMAIN(did)%VEGTYP == rt_domain(did)%islake) soltyp = rt_domain(did)%isoilwater + where(soltyp == rt_domain(did)%isoilwater) RT_DOMAIN(did)%VEGTYP = rt_domain(did)%iswater ! LKSAT, ! temporary set - RT_DOMAIN(did)%SMCRTCHK = 0 - RT_DOMAIN(did)%SMCAGGRT = 0 - RT_DOMAIN(did)%STCAGGRT = 0 - RT_DOMAIN(did)%SH2OAGGRT = 0 + RT_DOMAIN(did)%SMCRTCHK = 0 + RT_DOMAIN(did)%SMCAGGRT = 0 + RT_DOMAIN(did)%STCAGGRT = 0 + RT_DOMAIN(did)%SH2OAGGRT = 0 - rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil) = nlst(did)%zsoil8(1:nlst(did)%nsoil) + rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil) = nlst(did)%zsoil8(1:nlst(did)%nsoil) - rt_domain(did)%subsurface%properties%sldpth(1) = abs( RT_DOMAIN(did)%subsurface%properties%zsoil(1) ) - do i = 2, nlst(did)%nsoil - rt_domain(did)%subsurface%properties%sldpth(i) = RT_DOMAIN(did)%subsurface%properties%zsoil(i-1)-RT_DOMAIN(did)%subsurface%properties%zsoil(i) - enddo - rt_domain(did)%subsurface%properties%soldeprt = -1.0*RT_DOMAIN(did)%subsurface%properties%zsoil(nlst(did)%NSOIL) + rt_domain(did)%subsurface%properties%sldpth(1) = abs( RT_DOMAIN(did)%subsurface%properties%zsoil(1) ) + do i = 2, nlst(did)%nsoil + rt_domain(did)%subsurface%properties%sldpth(i) = RT_DOMAIN(did)%subsurface%properties%zsoil(i-1)-RT_DOMAIN(did)%subsurface%properties%zsoil(i) + enddo + rt_domain(did)%subsurface%properties%soldeprt = -1.0*RT_DOMAIN(did)%subsurface%properties%zsoil(nlst(did)%NSOIL) - ierr_flg = 99 - if(trim(nlst(did)%hydrotbl_f) == "") then - call hydro_stop("FATAL ERROR: hydrotbl_f is empty. Please input a netcdf file. ") - endif + ierr_flg = 99 + if(trim(nlst(did)%hydrotbl_f) == "") then + call hydro_stop("FATAL ERROR: hydrotbl_f is empty. Please input a netcdf file. ") + endif #ifdef MPP_LAND - if(my_id .eq. IO_id) then + if(my_id .eq. IO_id) then #endif - ierr_flg = nf90_open(trim(nlst(did)%hydrotbl_f), nf90_NOWRITE, ncid) + ierr_flg = nf90_open(trim(nlst(did)%hydrotbl_f), nf90_NOWRITE, ncid) #ifdef MPP_LAND - endif - call mpp_land_bcast_int1(ierr_flg) + endif + call mpp_land_bcast_int1(ierr_flg) #endif - if( ierr_flg .ne. 0) then - ! input from HYDRO.tbl FILE + if( ierr_flg .ne. 0) then +! input from HYDRO.tbl FILE ! input OV_ROUGH from OVROUGH.TBL #ifdef MPP_LAND - if(my_id .eq. IO_id) then + if(my_id .eq. IO_id) then #endif #ifndef NCEP_WCOSS - open(71,file="HYDRO.TBL", form="formatted") + open(71,file="HYDRO.TBL", form="formatted") !read OV_ROUGH first - read(71,*) nn - read(71,*) - do i = 1, nn - read(71,*) RT_DOMAIN(did)%OV_ROUGH(i) - end do + read(71,*) nn + read(71,*) + do i = 1, nn + read(71,*) RT_DOMAIN(did)%OV_ROUGH(i) + end do !read parameter for LKSAT - read(71,*) nn - read(71,*) - do i = 1, nn - read(71,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i) - end do - close(71) + read(71,*) nn + read(71,*) + do i = 1, nn + read(71,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i) + end do + close(71) #else - open(13, form="formatted") + open(13, form="formatted") !read OV_ROUGH first - read(13,*) nn - read(13,*) - do i = 1, nn - read(13,*) RT_DOMAIN(did)%OV_ROUGH(i) - end do + read(13,*) nn + read(13,*) + do i = 1, nn + read(13,*) RT_DOMAIN(did)%OV_ROUGH(i) + end do !read parameter for LKSAT - read(13,*) nn - read(13,*) - do i = 1, nn - read(13,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i) - end do - close(13) + read(13,*) nn + read(13,*) + do i = 1, nn + read(13,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i) + end do + close(13) #endif #ifdef MPP_LAND - endif - call mpp_land_bcast_real(leng,RT_DOMAIN(did)%OV_ROUGH) - call mpp_land_bcast_real(leng,xdum1) - call mpp_land_bcast_real(leng,MAXSMC) - call mpp_land_bcast_real(leng,refsmc) - call mpp_land_bcast_real(leng,wltsmc) -#endif - - rt_domain(did)%lksat = 0.0 - do j = 1, RT_DOMAIN(did)%jx - do i = 1, RT_DOMAIN(did)%ix - !yw rt_domain(did)%lksat(i,j) = xdum1(soltyp(i,j) ) * 1000.0 - rt_domain(did)%lksat(i,j) = xdum1(soltyp(i,j) ) - rt_domain(did)%SMCMAX1(i,j) = MAXSMC(soltyp(I,J)) - rt_domain(did)%SMCREF1(i,j) = refsmc(soltyp(I,J)) - rt_domain(did)%SMCWLT1(i,j) = wltsmc(soltyp(I,J)) - !ADCHANGE: Add some sanity checks in case calibration knocks the order of these out of sequence. - !The min diffs were pulled from the existing HYDRO.TBL defaults. - !Currently water is 0, so enforcing 0 as the absolute min. - rt_domain(did)%SMCMAX1(i,j) = min(rt_domain(did)%SMCMAX1(i,j), 1.0) - rt_domain(did)%SMCREF1(i,j) = max(min(rt_domain(did)%SMCREF1(i,j), rt_domain(did)%SMCMAX1(i,j) - 0.01), 0.0) - rt_domain(did)%SMCWLT1(i,j) = max(min(rt_domain(did)%SMCWLT1(i,j), rt_domain(did)%SMCREF1(i,j) - 0.01), 0.0) - IF(rt_domain(did)%VEGTYP(i,j) > 0 ) THEN ! created 2d ov_rough - rt_domain(did)%OV_ROUGH2d(i,j) = RT_DOMAIN(did)%OV_ROUGH(rt_domain(did)%VEGTYP(I,J)) - endif - end do - end do - - call hdtbl_out(did) - else - ! input from HYDRO.TBL.nc file - print*, "reading from hydrotbl_f(HYDRO.TBL.nc) file ...." - call hdtbl_in_nc(did) - if (noah_lsm%imperv_option .eq. 9) then - !ADCHANGE: For consistency, mirror urban and param value checks used in table read - where (rt_domain(did)%VEGTYP == rt_domain(did)%isurban) - rt_domain(did)%SMCMAX1 = 0.45 - rt_domain(did)%SMCREF1 = 0.42 - rt_domain(did)%SMCWLT1 = 0.40 - endwhere - endif - where (rt_domain(did)%SMCMAX1 .gt. 1.0) rt_domain(did)%SMCMAX1 = 1.0 - rt_domain(did)%SMCREF1 = max(min(rt_domain(did)%SMCREF1, rt_domain(did)%SMCMAX1 - 0.01), 0.0) - rt_domain(did)%SMCWLT1 = max(min(rt_domain(did)%SMCWLT1, rt_domain(did)%SMCREF1 - 0.01), 0.0) - endif - - rt_domain(did)%soiltyp = soltyp - - if(allocated(soltyp)) deallocate(soltyp) - + endif + call mpp_land_bcast_real(leng,RT_DOMAIN(did)%OV_ROUGH) + call mpp_land_bcast_real(leng,xdum1) + call mpp_land_bcast_real(leng,MAXSMC) + call mpp_land_bcast_real(leng,refsmc) + call mpp_land_bcast_real(leng,wltsmc) +#endif + + rt_domain(did)%lksat = 0.0 + do j = 1, RT_DOMAIN(did)%jx + do i = 1, RT_DOMAIN(did)%ix +!yw rt_domain(did)%lksat(i,j) = xdum1(soltyp(i,j) ) * 1000.0 + rt_domain(did)%lksat(i,j) = xdum1(soltyp(i,j) ) + rt_domain(did)%SMCMAX1(i,j) = MAXSMC(soltyp(I,J)) + rt_domain(did)%SMCREF1(i,j) = refsmc(soltyp(I,J)) + rt_domain(did)%SMCWLT1(i,j) = wltsmc(soltyp(I,J)) +!ADCHANGE: Add some sanity checks in case calibration knocks the order of these out of sequence. +!The min diffs were pulled from the existing HYDRO.TBL defaults. +!Currently water is 0, so enforcing 0 as the absolute min. + rt_domain(did)%SMCMAX1(i,j) = min(rt_domain(did)%SMCMAX1(i,j), 1.0) + rt_domain(did)%SMCREF1(i,j) = max(min(rt_domain(did)%SMCREF1(i,j), rt_domain(did)%SMCMAX1(i,j) - 0.01), 0.0) + rt_domain(did)%SMCWLT1(i,j) = max(min(rt_domain(did)%SMCWLT1(i,j), rt_domain(did)%SMCREF1(i,j) - 0.01), 0.0) + IF(rt_domain(did)%VEGTYP(i,j) > 0 ) THEN ! created 2d ov_rough + rt_domain(did)%OV_ROUGH2d(i,j) = RT_DOMAIN(did)%OV_ROUGH(rt_domain(did)%VEGTYP(I,J)) + endif + end do + end do - end subroutine lsm_input + call hdtbl_out(did) + else +! input from HYDRO.TBL.nc file + print*, "reading from hydrotbl_f(HYDRO.TBL.nc) file ...." + call hdtbl_in_nc(did) + if (noah_lsm%imperv_option .eq. 9) then +!ADCHANGE: For consistency, mirror urban and param value checks used in table read + where (rt_domain(did)%VEGTYP == rt_domain(did)%isurban) + rt_domain(did)%SMCMAX1 = 0.45 + rt_domain(did)%SMCREF1 = 0.42 + rt_domain(did)%SMCWLT1 = 0.40 + endwhere + endif + where (rt_domain(did)%SMCMAX1 .gt. 1.0) rt_domain(did)%SMCMAX1 = 1.0 + rt_domain(did)%SMCREF1 = max(min(rt_domain(did)%SMCREF1, rt_domain(did)%SMCMAX1 - 0.01), 0.0) + rt_domain(did)%SMCWLT1 = max(min(rt_domain(did)%SMCWLT1, rt_domain(did)%SMCREF1 - 0.01), 0.0) + endif + rt_domain(did)%soiltyp = soltyp -end module module_HYDRO_drv + if(allocated(soltyp)) deallocate(soltyp) + end subroutine lsm_input -! stop the job due to the fatal error. -subroutine HYDRO_finish() + subroutine HYDRO_finish() #ifdef MPP_LAND - USE module_mpp_land + use module_mpp_land #endif #ifdef WRF_HYDRO_NUDGING - use module_stream_nudging, only: finish_stream_nudging + use module_stream_nudging, only: finish_stream_nudging #endif - integer :: ierr + integer :: ierr #ifdef WRF_HYDRO_NUDGING - call finish_stream_nudging() + call finish_stream_nudging() #endif #ifndef NCEP_WCOSS - print*, "The model finished successfully......." + print*, "The model finished successfully......." #else - write(78,*) "The model finished successfully......." + write(78,*) "The model finished successfully......." #endif #ifdef MPP_LAND -! call mpp_land_abort() #ifndef NCEP_WCOSS - call flush(6) + call flush(6) #else - call flush(78) - close(78) + call flush(78) + close(78) #endif - call mpp_land_sync() - call MPI_Finalize(ierr) - stop + call mpp_land_sync() + call MPI_Finalize(ierr) + stop #else #ifndef WRF_HYDRO_NUDGING - stop !!JLM want to time at the top NoahMP level. + stop !!JLM want to time at the top NoahMP level. #endif #endif + return + end subroutine HYDRO_finish - return -end subroutine HYDRO_finish +end module module_HYDRO_drv diff --git a/src/Land_models/NoahMP/IO_code/main_hrldas_driver.F b/src/Land_models/NoahMP/IO_code/main_hrldas_driver.F index d479a7f0d..72a4fe8bb 100644 --- a/src/Land_models/NoahMP/IO_code/main_hrldas_driver.F +++ b/src/Land_models/NoahMP/IO_code/main_hrldas_driver.F @@ -9,6 +9,8 @@ program noah_hrldas_driver use module_noahmp_hrldas_driver, only: land_driver_ini, land_driver_exe #endif + use module_HYDRO_drv, only: HYDRO_finish + ! NEW MODULE WITH DT INCLUDING STATE, PARAMETER,FORCING AND GEOMETRY. Different modules use state_module, only: state_type use parameters diff --git a/src/Land_models/NoahMP/IO_code/module_NoahMP_hrldas_driver.F b/src/Land_models/NoahMP/IO_code/module_NoahMP_hrldas_driver.F index a584b0510..93f8fae9f 100644 --- a/src/Land_models/NoahMP/IO_code/module_NoahMP_hrldas_driver.F +++ b/src/Land_models/NoahMP/IO_code/module_NoahMP_hrldas_driver.F @@ -41,8 +41,9 @@ module module_NoahMP_hrldas_driver REAL, ALLOCATABLE, DIMENSION(:,:) :: ACCECAN ! accumulated canopy evap [mm] REAL, ALLOCATABLE, DIMENSION(:,:) :: ACCETRAN ! accumulated transpiration [mm] REAL, ALLOCATABLE, DIMENSION(:,:) :: ACCEDIR ! accumulated direct soil evap [mm] - integer :: io_config_outputs=0 - integer :: t0OutputFlag=0 + integer :: io_config_outputs = 0 + integer :: reset_flag = 0 + integer :: t0OutputFlag = 0 REAL, ALLOCATABLE, DIMENSION(:,:) :: SOILSAT_TOP ! top 2 layer soil saturation [fraction] REAL, ALLOCATABLE, DIMENSION(:,:) :: SOILSAT ! column integrated soil saturation [fraction] REAL, ALLOCATABLE, DIMENSION(:,:) :: SOILICE ! fraction of soil moisture that is ice [fraction] @@ -1533,9 +1534,10 @@ subroutine land_driver_ini(NTIME_out, state, wrfits,wrfite,wrfjts,wrfjte) olddate,zsoil(1:NSOIL)) call get_iocflag(1, io_config_outputs) +call get_rstrt_swc(1, reset_flag) #ifdef WRF_HYDRO -if (io_config_outputs .gt. 0) then +if (reset_flag == 1) then ACCPRCP = 0.0 ACCECAN = 0.0 ACCEDIR = 0.0 diff --git a/src/utils/module_hydro_stop.F90 b/src/utils/module_hydro_stop.F90 index 47fe9c182..c4b25a21a 100644 --- a/src/utils/module_hydro_stop.F90 +++ b/src/utils/module_hydro_stop.F90 @@ -14,7 +14,7 @@ subroutine HYDRO_stop(msg) ierr = 1 #ifndef NCEP_WCOSS !#ifdef HYDRO_D !! PLEASE NEVER UNCOMMENT THIS IFDEF, it's just one incredibly useful string. - write(6,*) "The job is stopped due to the fatal error. ", trim(msg) + write(6,'(a)') "The job has stopped due to a fatal error: ", trim(msg) call flush(6) !#endif #else