From 1b1340a206f371a5b7d7ba59b4dc99b28a519c9f Mon Sep 17 00:00:00 2001 From: David Yates Date: Tue, 8 Sep 2020 13:23:30 -0600 Subject: [PATCH 1/5] made lake variables double precision --- src/Data_Rec/module_RT_data.F | 2 +- src/Data_Rec/rt_include.inc | 6 +- src/MPP/mpp_land.F | 78 +++++- .../Reservoirs/Level_Pool/module_levelpool.F | 252 ++++++++++-------- .../Level_Pool/module_levelpool_state.F | 4 +- .../module_persistence_levelpool_hybrid.F | 14 +- ...odule_persistence_levelpool_hybrid_state.F | 4 +- .../RFC_Forecasts/module_rfc_forecasts.F | 14 +- .../module_rfc_forecasts_state.F | 6 +- src/Routing/Reservoirs/module_reservoir.F | 24 +- .../Reservoirs/module_reservoir_utilities.F | 20 +- src/Routing/module_HYDRO_io.F | 58 ++-- src/Routing/module_NWM_io.F | 59 +++- src/Routing/module_RT.F | 2 +- src/Routing/module_channel_routing.F | 136 +++++----- 15 files changed, 413 insertions(+), 266 deletions(-) diff --git a/src/Data_Rec/module_RT_data.F b/src/Data_Rec/module_RT_data.F index a1fa61b4d..58e16edee 100644 --- a/src/Data_Rec/module_RT_data.F +++ b/src/Data_Rec/module_RT_data.F @@ -25,7 +25,7 @@ Module module_RT_data use module_subsurface_input use module_subsurface_output use module_reservoir, only: reservoir_container - use iso_fortran_env, only: int64 + use iso_fortran_env, only: int64, real64 IMPLICIT NONE INTEGER, PARAMETER :: max_domain=5 diff --git a/src/Data_Rec/rt_include.inc b/src/Data_Rec/rt_include.inc index 1ccc7c793..51f960301 100644 --- a/src/Data_Rec/rt_include.inc +++ b/src/Data_Rec/rt_include.inc @@ -129,8 +129,8 @@ INTEGER, allocatable, DIMENSION(:) :: TYPEL !type of link Muskingum: 0 strm 1 lake !-- Diffusion: 0 edge or pour; 1 interior; 2 lake INTEGER, allocatable, DIMENSION(:) :: TYPEN !type of link 0 strm 1 lake - REAL, allocatable, DIMENSION(:) :: QLAKEI !lake inflow in difussion scheme - REAL, allocatable, DIMENSION(:) :: QLAKEO !lake outflow in difussion scheme + REAL(kind=real64), allocatable, DIMENSION(:) :: QLAKEI !lake inflow in difussion scheme + REAL(kind=real64), allocatable, DIMENSION(:) :: QLAKEO !lake outflow in difussion scheme INTEGER(kind=int64), allocatable, DIMENSION(:) :: LAKENODE !which nodes flow into which lakes integer(kind=int64), allocatable, dimension(:) :: LINKID ! id of links on linked routing REAL, allocatable, DIMENSION(:) :: CVOL ! channel volume @@ -159,7 +159,7 @@ REAL, DIMENSION(50) :: TOPWIDCC, NCC !topwidth and mannings n of compund ! VARIABLES FOR RESERVOIRS - REAL, allocatable, DIMENSION(:) :: RESHT !reservoir height + REAL*8, allocatable, DIMENSION(:) :: RESHT !reservoir height !-- lake params integer(kind=int64), allocatable, dimension(:) :: LAKEIDA !id of lakes in routlink file integer(kind=int64), allocatable, dimension(:) :: LAKEIDM !id of LAKES Modeled in LAKEPARM.nc or tbl diff --git a/src/MPP/mpp_land.F b/src/MPP/mpp_land.F index 9de75f0bf..cbcd3f61e 100644 --- a/src/MPP/mpp_land.F +++ b/src/MPP/mpp_land.F @@ -2092,9 +2092,52 @@ subroutine write_chanel_int8(v,map_l2g,gnlinks,nlinks,g_v) end subroutine write_chanel_int8 - subroutine write_lake_real(v,nodelist_in,nlakes) + subroutine write_lake_real8(v,nodelist_in,nlakes) implicit none - real recv(nlakes), v(nlakes) + real*8 recv(nlakes) + !real recv(nlakes) + real*8 v(nlakes) + integer nodelist(nlakes), nlakes, nodelist_in(nlakes) + integer i, ierr, tag, k + integer length, node + + nodelist = nodelist_in + if(my_id .eq. IO_id) then + do i = 0, numprocs - 1 + if(i .ne. my_id) then + !block receive from other node. + tag = 129 + call mpi_recv(nodelist,nlakes,MPI_INTEGER,i, & + tag,HYDRO_COMM_WORLD,mpp_status,ierr) + tag = 139 + !call mpi_recv(recv(:),nlakes,MPI_REAL,i, & + call mpi_recv(recv(:),nlakes,MPI_DOUBLE_PRECISION,i, & + tag,HYDRO_COMM_WORLD,mpp_status,ierr) + + do k = 1,nlakes + if(nodelist(k) .gt. -99) then + node = nodelist(k) + v(node) = recv(node) + endif + enddo + end if + end do + else + tag = 129 + call mpi_send(nodelist,nlakes,MPI_INTEGER, IO_id, & + tag,HYDRO_COMM_WORLD,ierr) + tag = 139 + !call mpi_send(v,nlakes,MPI_REAL,IO_id, & + call mpi_send(v,nlakes,MPI_DOUBLE_PRECISION,IO_id, & + tag,HYDRO_COMM_WORLD,ierr) + end if + end subroutine write_lake_real8 + + + subroutine write_lake_real(v,nodelist_in,nlakes) + implicit none + real recv(nlakes) + real v(nlakes) integer nodelist(nlakes), nlakes, nodelist_in(nlakes) integer i, ierr, tag, k integer length, node @@ -2604,29 +2647,38 @@ end subroutine updateLake_seqInt8 subroutine updateLake_seq(in,nsize,in0) implicit none integer :: nsize - real, dimension(nsize) :: in - real, dimension(nsize) :: tmp - real, dimension(:) :: in0 + real*8, dimension(nsize) :: in + real*8, dimension(nsize) :: tmp + real*8, dimension(:) :: in0 integer tag, i, status, ierr, k if(nsize .le. 0) return tag = 29 if(my_id .ne. IO_id) then - call mpi_send(in,nsize,MPI_REAL, IO_id, & + + ! call mpi_send(in,nsize,MPI_REAL, IO_id, & + ! tag,HYDRO_COMM_WORLD,ierr) + call mpi_send(in,nsize,MPI_double, IO_id, & tag,HYDRO_COMM_WORLD,ierr) else do i = 0, numprocs - 1 if(i .ne. IO_id) then + + ! call mpi_recv(tmp,nsize,& + ! MPI_REAL,i,tag,HYDRO_COMM_WORLD,mpp_status,ierr) call mpi_recv(tmp,nsize,& - MPI_REAL,i,tag,HYDRO_COMM_WORLD,mpp_status,ierr) - do k = 1, nsize + MPI_double,i,tag,HYDRO_COMM_WORLD,mpp_status,ierr) + + do k = 1, nsize if(in0(k) .ne. tmp(k)) in(k) = tmp(k) end do end if end do end if - call mpp_land_bcast_real_1d(in) + !call mpp_land_bcast_real_1d(in) + call mpp_land_bcast_real8_1d(in) + end subroutine updateLake_seq @@ -2664,7 +2716,8 @@ end subroutine updateLake_seq_char subroutine updateLake_grid(in,nsize,lake_index) implicit none integer :: nsize - real, dimension(nsize) :: in + !real, dimension(nsize) :: in + real*8, dimension(nsize) :: in integer, dimension(nsize) :: lake_index real, dimension(nsize) :: tmp integer tag, i, status, ierr, k @@ -2692,8 +2745,11 @@ subroutine updateLake_grid(in,nsize,lake_index) end if end do end if - call mpp_land_bcast_real_1d(in) + !call mpp_land_bcast_real_1d(in) + call mpp_land_bcast_real8_1d(in) + + end subroutine updateLake_grid diff --git a/src/Routing/Reservoirs/Level_Pool/module_levelpool.F b/src/Routing/Reservoirs/Level_Pool/module_levelpool.F index 6e8236751..9fadf8c93 100644 --- a/src/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/src/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -59,7 +59,7 @@ subroutine levelpool_init(this, water_elevation, & orifice_area, max_depth, lake_number) implicit none class(levelpool), intent(inout) :: this ! object being initialized - real, intent(inout) :: water_elevation ! meters AMSL + real*8, intent(inout) :: water_elevation ! meters AMSL real, intent(in) :: lake_area ! area of lake (km^2) real, intent(in) :: weir_elevation ! bottom of weir elevation (meters AMSL) real, intent(in) :: weir_coeffecient ! weir coefficient @@ -150,14 +150,14 @@ subroutine run_levelpool_reservoir(this, previous_timestep_inflow, inflow, & assimilated_value, assimilated_source_file) implicit none class(levelpool), intent(inout) :: this - real, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) - real, intent(in) :: inflow ! cubic meters per second (cms) - real, intent(in) :: lateral_inflow ! cubic meters per second (cms) - real, intent(inout) :: water_elevation ! meters - real, intent(out) :: outflow ! cubic meters per second (cms) - real, intent(in) :: routing_period ! seconds - integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files - real, intent(out) :: assimilated_value ! value assimilated from observation or forecast + real*8, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) + real*8, intent(in) :: inflow ! cubic meters per second (cms) + real, intent(in) :: lateral_inflow ! cubic meters per second (cms) + real*8, intent(inout) :: water_elevation ! meters + real*8, intent(out) :: outflow ! cubic meters per second (cms) + real, intent(in) :: routing_period ! seconds + integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files + real, intent(out) :: assimilated_value ! value assimilated from observation or forecast character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value @@ -165,6 +165,9 @@ subroutine run_levelpool_reservoir(this, previous_timestep_inflow, inflow, & this%input%inflow = inflow this%input%lateral_inflow = lateral_inflow + ! write(6,29) previous_timestep_inflow,inflow,outflow + ! 29 format("BLP ",f10.6,2x,f10.6,2x,f10.6) + ! Update state variables this%state%water_elevation = water_elevation @@ -211,6 +214,9 @@ subroutine run_levelpool_reservoir(this, previous_timestep_inflow, inflow, & call log_levelpool_diagnostic_data(this%properties%lake_number, inflow, water_elevation, outflow) #endif + ! write(6,30) inflow,previous_timestep_inflow,outflow + ! 30 format("ALP ",f8.6,2x,f8.6,2x,f8.6) + end subroutine run_levelpool_reservoir ! ------------------------------------------------ @@ -222,11 +228,11 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa !! ---------------------------- argument variables !! All elevations should be relative to a common base (often belev(k)) - real, intent(INOUT) :: H ! water elevation height (m) + real*8, intent(INOUT) :: H ! water elevation height (m) real, intent(IN) :: dt ! routing period [s] - real, intent(IN) :: qi0 ! inflow at previous timestep (cms) - real, intent(IN) :: qi1 ! inflow at current timestep (cms) - real, intent(OUT) :: qo1 ! outflow at current timestep + real*8, intent(IN) :: qi0 ! inflow at previous timestep (cms) + real*8, intent(IN) :: qi1 ! inflow at current timestep (cms) + real*8, intent(OUT) :: qo1 ! outflow at current timestep real, intent(IN) :: ql ! lateral inflow real, intent(IN) :: ar ! area of reservoir (km^2) real, intent(IN) :: we ! bottom of weir elevation @@ -241,35 +247,38 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa !!DJG Add lake option switch here...move up to namelist in future versions... integer :: LAKE_OPT ! Lake model option (move to namelist later) - real :: Htmp ! Temporary assign of incoming lake el. (m) !! ---------------------------- local variables - real :: sap ! local surface area values - real :: discharge ! storage discharge m^3/s - real :: tmp1, tmp2 - real :: dh, dh1, dh2, dh3 ! Depth in weir, and height function for 3 order RK - real :: It, Itdt_3, Itdt_2_3 ! inflow hydrographs - real :: maxWeirDepth !maximum capacity of weir - !real :: hdiff_vol, qdiff_vol ! water balance check variables + real*8 :: sap ! local surface area values + real*8 :: discharge ! storage discharge m^3/s + real*8 :: tmp1, tmp2,tmp3 + real*8 :: dh, dh1, dh2, dh3 ! Depth in weir, and height function for 3 order RK + real*8 :: It, Itdt_3, Itdt_2_3 ! inflow hydrographs + real*8 :: maxWeirDepth !maximum capacity of weir + real*8 :: hdiff_vol, qdiff_vol ! water balance check variables + real*8 :: Htmp ! temporary store of WSE + !! ---------------------------- subroutine body: from chow, mad mays. pg. 252 !! -- determine from inflow hydrograph !!DJG Set hardwire for LAKE_OPT...move specification of this to namelist in !future versions... - LAKE_OPT = 2 - Htmp = H !temporary set of incoming lake water elevation... - !hdiff_vol = 0.0 - !qdiff_vol = 0.0 - +! LAKE_OPT = 2 + hdiff_vol = 0.0 + qdiff_vol = 0.0 + Htmp = H + !!DJG IF-block for lake model option 1 - outflow=inflow, 2 - Chow et al level !pool, ..... - if (LAKE_OPT == 1) then ! If-block for simple pass through scheme.... +! if (LAKE_OPT == 1) then ! If-block for simple pass through scheme.... + +! qo1 = qi1 ! Set outflow equal to inflow at current time - qo1 = qi1 ! Set outflow equal to inflow at current time - H = Htmp ! Set new lake water elevation to incoming lake el. +! else if (LAKE_OPT == 2) then ! If-block for Chow et al level pool scheme - else if (LAKE_OPT == 2) then ! If-block for Chow et al level pool scheme + open (unit=67, & + file='lake_massbalance_out.txt', status='unknown',position='append') It = qi0 Itdt_3 = qi0 + ((qi1 + ql - qi0) * 0.33) @@ -279,134 +288,151 @@ subroutine LEVELPOOL_PHYSICS(ln,qi0,qi1,qo1,ql,dt,H,ar,we,maxh,wc,wl,dl,oe,oc,oa !assume vertically walled reservoir !remove this when moving to a variable head area volume sap = ar * 1.0E6 + !sap = ar * 1.0E4 !-- determine Q(dh) from elevation-discharge relationship !-- and dh1 - dh = H - we + dh = Htmp - we if (dh > maxWeirDepth) then dh = maxWeirDepth endif - tmp1 = oc * oa * sqrt(2. * 9.81 * ( H - oe )) !orifice at capacity - tmp2 = wc * wl * (dh ** (3./2.)) !weir flows at capacity + if(Htmp - oe .gt.0) then + !tmp1 = oc * (H-oe)/(we-oe)*oa * sqrt(2. * 9.81 * ( Htmp - oe )) !orifice flow is a function of storage + tmp1 = oc * oa * sqrt(2. * 9.81 * ( Htmp - oe )) !orifice flow is a function of storage + else + tmp1 =0.0 + endif + + if(dh .gt. 0.0) then + tmp2 = wc * wl * (dh ** (3./2.)) !spillway discharge + else + tmp2 =0.0 + endif !determine the discharge based on current height if(H > maxh) then - discharge = tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop - else if (dh > 0.0 ) then !! orifice and weir discharge + discharge = tmp1 + tmp2 + (wc* dl * (Htmp - maxh)**(3./2.)) !overtop + else discharge = tmp1 + tmp2 - else if ( H > oe ) then !! only orifice flow - discharge = oc * oa * sqrt(2. * 9.81 * ( H - oe ) ) - else - discharge = 0.0 !in the dead pool - endif - - if (sap > 0) then - dh1 = ((It - discharge)/sap)*dt - else - dh1 = 0.0 endif + dh1 = dble(((It - discharge)/sap)*dt) !-- determine Q(H + dh1/3) from elevation-discharge relationship !-- dh2 - dh = (H+dh1/3) - we + dh = (H+(0.33*dh1)) - we if (dh > maxWeirDepth) then dh = maxWeirDepth endif - tmp1 = oc * oa * sqrt(2. * 9.81 * ( (H+dh1/3.) - oe ) ) - tmp2 = wc * wl * (dh ** (3./2.)) - - !determine the discharge based on current height - if(H > maxh) then - discharge = tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop - else if (dh > 0.0 ) then !! orifice and weir discharge - discharge = tmp1 + tmp2 - else if ( (H+dh1/3) > oe ) then !! only orifice flow,not full - discharge = oc * oa * sqrt(2. * 9.81 * ( (H+dh1/3.) - oe ) ) + if((H + 0.33*dh1 ) - oe .gt. 0.0) then + !tmp1 = oc * ((H+0.33*dh1)-oe)/(we-oe)*oa * sqrt(2. * 9.81 * ( (H+0.33*dh1) - oe ) ) + tmp1 = oc * oa * sqrt(2. * 9.81 * ((Htmp + 0.33*dh1) - oe) ) else - discharge = 0.0 - endif - + tmp1 =0.0 + endif - if (sap > 0.0) then - dh2 = ((Itdt_3 - discharge)/sap)*dt + if(dh .gt. 0.0) then + tmp2 = wc * wl * (dh ** (3./2.)) else - dh2 = 0.0 + tmp2 =0.0 + endif + + !determine the discharge based on current height + if((H + 0.33*dh1 ) > maxh) then + discharge = tmp1 + tmp2 + (wc* dl * ((Htmp + 0.33*dh1)-maxh)**(3./2.)) !overtop + else + discharge = tmp1 + tmp2 endif + dh2 = dble(((Itdt_3 - discharge)/sap)*dt) + !-- determine Q(H + 2/3 dh2) from elevation-discharge relationship !-- dh3 - dh = (H + (0.667*dh2)) - we + dh = (Htmp + (0.667*dh2)) - we if (dh > maxWeirDepth) then dh = maxWeirDepth endif - tmp1 = oc * oa * sqrt(2. * 9.81 * ( (H+dh2*0.667) - oe ) ) - tmp2 = wc * wl * (dh ** (3./2.)) - - !determine the discharge based on current height - if(H > maxh) then ! overtop condition, not good! - discharge = tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop - else if (dh > 0.0 ) then !! orifice and weir discharge - discharge = tmp1 + tmp2 - else if ( (H+dh2*0.667) > oe ) then !! only orifice flow,not full - discharge = oc * oa * sqrt(2. * 9.81 * ( (H+dh2*0.667) - oe ) ) + if( (Htmp + 0.667*dh2 ) - oe .gt. 0.0) then + !tmp1 = oc * ((Htmp + 0.667*dh2)-oe)/(we-oe)*oa * sqrt(2. * 9.81 * ( (Htmp + 0.667*dh2) - oe ) ) + tmp1 = oc * oa * sqrt(2. * 9.81 * ( (Htmp + 0.667*dh2) - oe) ) else - discharge = 0.0 + tmp1 =0.0 endif - if (sap > 0.0) then - dh3 = ((Itdt_2_3 - discharge)/sap)*dt + if(dh .gt. 0.0) then + tmp2 = wc * wl * (dh ** (3./2.)) else - dh3 = 0.0 + tmp2 =0.0 + endif + + ! find dh3 + if((H+0.667*dh2) > maxh) then ! overtop condition, not good! + discharge = tmp1 + tmp2 + (wc* dl * ((Htmp + 0.667*dh2)-maxh)**(3./2.)) !overtop + else + discharge = tmp1 + tmp2 endif - !-- determine dh and H - dh = (dh1/4.) + (0.75*dh3) - H = H + dh + dh3 = dble(((Itdt_2_3 - discharge)/sap)*dt) - !-- compute final discharge - dh = H - we + H = Htmp + ((0.25*dh1) + (0.75*dh3)) + + if (H-Htmp .eq. 0.0 ) then + write(6,23) H,Htmp,((0.25*dh1) + (0.75*dh3)) + 23 format("DHZERO ",f16.8,2x,f16.8,2x,f12.8) + endif + + !-- compute final discharge + dh = dble(H - we) if (dh > maxWeirDepth) then - dh = maxWeirDepth + dh = maxWeirDepth + endif + + if ((dble(H)-oe) .gt. 0.0) then + !tmp1 = oc * (H-oe)/(we-oe)*oa * sqrt(2. * 9.81 * (H - oe) ) + tmp1 = oc * oa * sqrt(2. * 9.81 * ( dble(H) - oe ) ) + else + tmp1 = 0.0 endif - tmp1 = oc * oa * sqrt(2. * 9.81 * ( H - oe ) ) - tmp2 = wc * wl * (dh ** (3./2.)) - - !determine the discharge based on current height + if(dh .gt. 0.0) then + tmp2 = wc * wl * (dh ** (3./2.)) + else + tmp2 =0.0 + endif + + !determine the discharge based on current height if(H > maxh) then ! overtop condition, not good! - discharge = tmp1 + tmp2 + (wc* (wl*dl) * (H-maxh)**(3./2.)) !overtop - else if (dh > 0.0 ) then !! orifice and overtop discharge - discharge = tmp1 + tmp2 - else if ( H > oe ) then !! only orifice flow,not full - discharge = oc * oa * sqrt(2. * 9.81 * ( H - oe ) ) + tmp3 = (wc* dl * (dble(H)-maxh)**(3./2.)) !overtop + discharge = tmp1 + tmp2 + tmp3 else - discharge = 0.0 + discharge = tmp1 + tmp2 endif - + qo1 = discharge ! return the flow rate from reservoir - !#ifdef HYDRO_D - !#ifndef NCEP_WCOSS - ! ! Water balance check - ! qdiff_vol = (qi1+ql-qo1)*dt !m3 - ! hdiff_vol = (H-Htmp)*sap !m3 - !22 format(f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f8.4,2x,f6.0,2x,f20.1,2x,f20.1) - ! open (unit=67, & - ! file='lake_massbalance_out.txt', status='unknown',position='append') - ! write(67,22) Htmp, H, qi1, ql, qo1, dt, qdiff_vol, hdiff_vol - ! close(67) - !#endif - !#endif - - !23 format('botof H dh orf wr Q',f8.4,2x,f8.4,2x,f8.3,2x,f8.3,2x,f8.2) - !24 format('ofonl H dh sap Q ',f8.4,2x,f8.4,2x,f8.0,2x,f8.2) - - - else ! ELSE for LAKE_OPT.... - endif ! ENDIF for LAKE_OPT.... + ! #ifdef HYDRO_D + ! #ifndef NCEP_WCOSS + + ! Water balance check + qdiff_vol = (qo1-qi1)*dt !m3 over time interval + hdiff_vol = (dble(H)-Htmp)*sap !m3 + + ! open (unit=67, & + ! file='lake_massbalance_out.txt', status='unknown',position='append') + ! write(67,22) H,Htmp,qi1,qo1,ql,qdiff_vol,hdiff_vol,((dh1/4.) + (0.75*dh3)) + ! close(67) + + write(67,22) H,Htmp,qi1,qo1 + 22 format(f12.6,2x,f12.6,2x,f9.6,2x,f9.6) + close(67) + + ! #endif + ! #endif + +! else ! ELSE for LAKE_OPT.... +! endif ! ENDIF for LAKE_OPT.... return diff --git a/src/Routing/Reservoirs/Level_Pool/module_levelpool_state.F b/src/Routing/Reservoirs/Level_Pool/module_levelpool_state.F index 98724de8a..6d5cbdd6c 100644 --- a/src/Routing/Reservoirs/Level_Pool/module_levelpool_state.F +++ b/src/Routing/Reservoirs/Level_Pool/module_levelpool_state.F @@ -12,7 +12,7 @@ module module_levelpool_state ! Extend/derive level pool state from the abstract base ! type for reservoir state. type, extends(reservoir_state) :: levelpool_state_interface - real :: water_elevation ! meters AMSL + real*8 :: water_elevation ! meters AMSL contains @@ -27,7 +27,7 @@ module module_levelpool_state subroutine levelpool_state_init(this, water_elevation) implicit none class(levelpool_state_interface), intent(inout) :: this ! the type object being initialized - real, intent(inout) :: water_elevation ! meters AMSL + real*8, intent(inout) :: water_elevation ! meters AMSL ! Assign the water elevation value passed in to a particular level pool reservoir ! state object's variable for water elevation diff --git a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F index 6d4b599b7..22612447a 100644 --- a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F +++ b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F @@ -66,7 +66,7 @@ subroutine hybrid_init(this, water_elevation, & observation_update_time_interval_seconds) implicit none class(persistence_levelpool_hybrid), intent(inout) :: this ! object being initialized - real, intent(inout) :: water_elevation ! meters AMSL + real*8, intent(inout) :: water_elevation ! meters AMSL real, intent(in) :: lake_area ! area of lake (km^2) real, intent(in) :: weir_elevation ! bottom of weir elevation (meters AMSL) real, intent(in) :: weir_coeffecient ! weir coefficient @@ -220,18 +220,18 @@ subroutine run_hybrid_reservoir(this, previous_timestep_inflow, inflow, & assimilated_value, assimilated_source_file) implicit none class(persistence_levelpool_hybrid), intent(inout) :: this - real, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) - real, intent(in) :: inflow ! cubic meters per second (cms) + real*8, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) + real*8, intent(in) :: inflow ! cubic meters per second (cms) real, intent(in) :: lateral_inflow ! cubic meters per second (cms) - real, intent(inout) :: water_elevation ! meters - real, intent(out) :: outflow ! cubic meters per second (cms) + real*8, intent(inout) :: water_elevation ! meters + real*8, intent(out) :: outflow ! cubic meters per second (cms) real, intent(in) :: routing_period ! seconds integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files real, intent(out) :: assimilated_value ! value assimilated from observation character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value real :: delta_storage ! timestep change in storage (cubic meters) - real :: local_water_elevation ! water elevation passed to levelpool (meters AMSL) - real :: levelpool_outflow ! cubic meters per second (cms) + real*8 :: local_water_elevation ! water elevation passed to levelpool (meters AMSL) + real*8 :: levelpool_outflow ! cubic meters per second (cms) logical :: max_storage_reached ! flag for when max storage is reached integer :: gage_lookback_seconds ! seconds between current model ! time and the time a gage diff --git a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F index 4fcf1327d..424776b94 100644 --- a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F +++ b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F @@ -14,7 +14,7 @@ module module_persistence_levelpool_hybrid_state ! Extend/derive hybrid state from the abstract base ! type for reservoir state. type, extends(reservoir_state) :: hybrid_state_interface - real :: water_elevation ! meters AMSL + real*8 :: water_elevation ! meters AMSL real*4 :: current_storage ! cubic meters real :: gage_discharge ! cubic meters per second (cms) real :: persisted_outflow ! cubic meters per second (cms) @@ -49,7 +49,7 @@ subroutine hybrid_state_init(this, water_elevation, lake_area, lake_max_water_el initial_fractional_depth, reservoir_type) implicit none class(hybrid_state_interface), intent(inout) :: this ! the type object being initialized - real, intent(in) :: water_elevation ! meters AMSL + real*8, intent(in) :: water_elevation ! meters AMSL real, intent(in) :: lake_area ! area of lake (km^2) real, intent(in) :: lake_max_water_elevation ! max water elevation (meters) real, intent(in) :: orifice_elevation ! orifice elevation (meters AMSL) diff --git a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F index d78d056f8..7aa568260 100644 --- a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F +++ b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F @@ -61,7 +61,7 @@ subroutine rfc_forecasts_init(this, water_elevation, & implicit none class(rfc_forecasts), intent(inout) :: this ! object being initialized - real, intent(inout) :: water_elevation ! meters AMSL + real*8, intent(inout) :: water_elevation ! meters AMSL real, intent(in) :: lake_area ! area of lake (km^2) real, intent(in) :: weir_elevation ! bottom of weir elevation (meters AMSL) real, intent(in) :: weir_coeffecient ! weir coefficient @@ -206,16 +206,16 @@ subroutine run_rfc_forecasts_reservoir(this, previous_timestep_inflow, inflow, & assimilated_value, assimilated_source_file) implicit none class(rfc_forecasts), intent(inout) :: this - real, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) - real, intent(in) :: inflow ! cubic meters per second (cms) + real*8, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) + real*8, intent(in) :: inflow ! cubic meters per second (cms) real, intent(in) :: lateral_inflow ! cubic meters per second (cms) - real, intent(inout) :: water_elevation ! meters - real, intent(out) :: outflow ! cubic meters per second (cms) + real*8, intent(inout) :: water_elevation ! meters + real*8, intent(out) :: outflow ! cubic meters per second (cms) real, intent(in) :: routing_period ! seconds integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files - real, intent(out) :: assimilated_value ! value assimilated from observation or forecast + real, intent(out) :: assimilated_value ! value assimilated from observation or forecast character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value - real :: levelpool_outflow ! cubic meters per second (cms) + real*8 :: levelpool_outflow ! cubic meters per second (cms) integer :: missing_outflow_index character(len=15) :: lake_number_string diff --git a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F index e71d6c1f0..779dc1a04 100644 --- a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F +++ b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F @@ -14,8 +14,8 @@ module module_rfc_forecasts_state ! Extend/derive rfc forecasts state from the abstract base ! type for reservoir state. type, extends(reservoir_state) :: rfc_forecasts_state_interface - real :: water_elevation ! meters AMSL - real :: levelpool_water_elevation ! meters AMSL + real*8 :: water_elevation ! meters AMSL + real*8 :: levelpool_water_elevation ! meters AMSL real, allocatable, dimension(:) :: discharges logical, allocatable, dimension(:) :: synthetic_values integer :: time_series_update_time ! seconds @@ -45,7 +45,7 @@ module module_rfc_forecasts_state subroutine rfc_forecasts_state_init(this, water_elevation, lake_area, lake_max_water_elevation, orifice_elevation, initial_fractional_depth) implicit none class(rfc_forecasts_state_interface), intent(inout) :: this ! the type object being initialized - real, intent(in) :: water_elevation ! meters AMSL + real*8, intent(in) :: water_elevation ! meters AMSL real, intent(in) :: lake_area ! area of lake (km^2) real, intent(in) :: lake_max_water_elevation ! max water elevation (meters) real, intent(in) :: orifice_elevation ! orifice elevation (meters AMSL) diff --git a/src/Routing/Reservoirs/module_reservoir.F b/src/Routing/Reservoirs/module_reservoir.F index d48d3347c..d0517e192 100644 --- a/src/Routing/Reservoirs/module_reservoir.F +++ b/src/Routing/Reservoirs/module_reservoir.F @@ -33,9 +33,9 @@ module module_reservoir ! object's methods/subroutines at runtime. type :: reservoir_input - real :: inflow ! cubic meters per second (cms) - real :: lateral_inflow ! cubic meters per second (cms) - real :: previous_timestep_inflow ! cubic meters per second (cms) + real*8 :: inflow ! cubic meters per second (cms) + real :: lateral_inflow ! cubic meters per second (cms) + real*8 :: previous_timestep_inflow ! cubic meters per second (cms) contains @@ -51,7 +51,7 @@ module module_reservoir ! object's methods/subroutines at runtime. type :: reservoir_output - real :: outflow ! cubic meters per second (cms) + real*8 :: outflow ! cubic meters per second (cms) contains @@ -96,14 +96,14 @@ subroutine run_reservoir_interface(this, previous_timestep_inflow, inflow, & assimilated_value, assimilated_source_file) import reservoir class(reservoir), intent(inout) :: this - real, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) - real, intent(in) :: inflow ! cubic meters per second (cms) - real, intent(in) :: lateral_inflow ! cubic meters per second (cms) - real, intent(inout) :: water_elevation ! meters AMSL - real, intent(out) :: outflow ! cubic meters per second (cms) - real, intent(in) :: routing_period ! seconds - integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files - real, intent(out) :: assimilated_value ! value assimilated from observation or forecast + real*8, intent(in) :: previous_timestep_inflow ! cubic meters per second (cms) + real*8, intent(in) :: inflow ! cubic meters per second (cms) + real, intent(in) :: lateral_inflow ! cubic meters per second (cms) + real*8, intent(inout) :: water_elevation ! meters AMSL + real*8, intent(out) :: outflow ! cubic meters per second (cms) + real, intent(in) :: routing_period ! seconds + integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files + real, intent(out) :: assimilated_value ! value assimilated from observation or forecast character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value end subroutine run_reservoir_interface diff --git a/src/Routing/Reservoirs/module_reservoir_utilities.F b/src/Routing/Reservoirs/module_reservoir_utilities.F index b6f59499d..f8c298616 100644 --- a/src/Routing/Reservoirs/module_reservoir_utilities.F +++ b/src/Routing/Reservoirs/module_reservoir_utilities.F @@ -17,7 +17,7 @@ module module_reservoir_utilities ! Checks and warns for boundary conditions of negative inflow subroutine warn_negative_inflow(inflow, lake_number, current_time) - real, intent(in) :: inflow + real*8, intent(in) :: inflow integer, intent(in) :: lake_number, current_time if (inflow < 0.0) then @@ -31,7 +31,8 @@ end subroutine warn_negative_inflow ! Checks and warns for boundary conditions of exceeding max or min storage subroutine warn_storage_boundaries(inflow, end_time_storage, min_storage, max_storage, lake_number, current_time) - real, intent(in) :: inflow, end_time_storage, min_storage, max_storage + real, intent(in) :: end_time_storage, min_storage, max_storage + real*8, intent(in) :: inflow integer, intent(in) :: lake_number, current_time ! Also have option to set the end_time_storage to the max_storage and change the outflow accordingly. @@ -58,9 +59,12 @@ end subroutine warn_storage_boundaries ! Checks and modifies for boundary conditions of negative storage subroutine modify_negative_storage(inflow, previous_time_storage, routing_period, lake_number, current_time, & update_time, end_time_storage, outflow) - real, intent(in) :: inflow, previous_time_storage, routing_period + real, intent(in) :: previous_time_storage, routing_period + real*8, intent(in) :: inflow + integer, intent(in) :: lake_number, current_time, update_time - real, intent(inout) :: end_time_storage, outflow + real, intent(inout) :: end_time_storage + real*8, intent(inout) :: outflow ! Note to user: ! In order to conserve mass, when a reservoir's storage would become 0.0 or less, the release @@ -107,7 +111,7 @@ end subroutine modify_negative_storage ! Checks and modifies for boundary conditions of negative inflow subroutine modify_negative_inflow(inflow) - real, intent(inout) :: inflow + real*8, intent(inout) :: inflow if (inflow < 0.0) inflow = 0.0 @@ -117,9 +121,11 @@ end subroutine modify_negative_inflow ! Checks and modifies for boundary conditions of negative outflow or exceeding min or max storage subroutine modify_for_projected_storage(inflow, previous_time_storage, min_storage, max_storage, & lake_number, current_time, time_interval, outflow, max_storage_reached) - real, intent(in) :: inflow, previous_time_storage, min_storage, max_storage + real, intent(in) :: previous_time_storage, min_storage, max_storage + real*8, intent(in) :: inflow + integer, intent(in) :: lake_number, current_time, time_interval - real, intent(inout) :: outflow + real*8, intent(inout) :: outflow logical, intent(inout) :: max_storage_reached real :: excess, projected_next_time_interval_storage diff --git a/src/Routing/module_HYDRO_io.F b/src/Routing/module_HYDRO_io.F index b4312a88e..bb8bb3b9a 100644 --- a/src/Routing/module_HYDRO_io.F +++ b/src/Routing/module_HYDRO_io.F @@ -4751,21 +4751,24 @@ subroutine mpp_output_lakes(lake_index,igrid, split_output_count, NLAKES, & integer, intent(in) :: igrid, K integer, intent(in) :: split_output_count integer, intent(in) :: NLAKES - real, dimension(NLAKES), intent(in) :: latlake,lonlake,elevlake,resht - real, dimension(NLAKES), intent(in) :: qlakei,qlakeo !-- inflow and outflow of lake + real, dimension(NLAKES), intent(in) :: latlake,lonlake,elevlake !resht + real*8, dimension(NLAKES), intent(in) :: RESHT + real*8, dimension(NLAKES), intent(in) :: qlakei,qlakeo !-- inflow and outflow of lake real, intent(in) :: dtrt_ch character(len=*), intent(in) :: startdate character(len=*), intent(in) :: date integer lake_index(nlakes) + ! print *, "QLAKEIOa", qlakei,qlakeo call write_lake_real(latlake,lake_index,nlakes) call write_lake_real(lonlake,lake_index,nlakes) call write_lake_real(elevlake,lake_index,nlakes) - call write_lake_real(resht,lake_index,nlakes) - call write_lake_real(qlakei,lake_index,nlakes) - call write_lake_real(qlakeo,lake_index,nlakes) + call write_lake_real8(resht,lake_index,nlakes) + call write_lake_real8(qlakei,lake_index,nlakes) + call write_lake_real8(qlakeo,lake_index,nlakes) + if(my_id.eq. IO_id) then call output_lakes(igrid, split_output_count, NLAKES, & startdate, date, latlake, lonlake, elevlake, & @@ -4785,8 +4788,10 @@ subroutine mpp_output_lakes2(lake_index,igrid, split_output_count, NLAKES, & integer, intent(in) :: igrid, K integer, intent(in) :: split_output_count integer, intent(in) :: NLAKES - real, dimension(NLAKES), intent(inout) :: latlake,lonlake,elevlake,resht - real, dimension(NLAKES), intent(inout) :: qlakei,qlakeo !-- inflow and outflow of lake + real, dimension(NLAKES), intent(inout) :: latlake,lonlake,elevlake!resht + real*8, dimension(NLAKES), intent(inout) :: RESHT + + real*8, dimension(NLAKES), intent(inout) :: qlakei,qlakeo !-- inflow and outflow of lake real, intent(in) :: dtrt_ch integer(kind=int64), dimension(NLAKES), intent(in) :: LAKEIDM ! lake id @@ -4797,9 +4802,9 @@ subroutine mpp_output_lakes2(lake_index,igrid, split_output_count, NLAKES, & call write_lake_real(latlake,lake_index,nlakes) call write_lake_real(lonlake,lake_index,nlakes) call write_lake_real(elevlake,lake_index,nlakes) - call write_lake_real(resht,lake_index,nlakes) - call write_lake_real(qlakei,lake_index,nlakes) - call write_lake_real(qlakeo,lake_index,nlakes) + call write_lake_real8(resht,lake_index,nlakes) + call write_lake_real8(qlakei,lake_index,nlakes) + call write_lake_real8(qlakeo,lake_index,nlakes) if(my_id.eq. IO_id) then call output_lakes2(igrid, split_output_count, NLAKES, & @@ -4821,8 +4826,9 @@ subroutine output_lakes(igrid, split_output_count, NLAKES, & integer, intent(in) :: igrid, K integer, intent(in) :: split_output_count integer, intent(in) :: NLAKES - real, dimension(NLAKES), intent(in) :: latlake,lonlake,elevlake,resht - real, dimension(NLAKES), intent(in) :: qlakei,qlakeo !-- inflow and outflow of lake + real, dimension(NLAKES), intent(in) :: latlake,lonlake,elevlake!resht + real*8, dimension(NLAKES), intent(in) :: RESHT + real*8, dimension(NLAKES), intent(in) :: qlakei,qlakeo !-- inflow and outflow of lake real, intent(in) :: dtrt_ch character(len=*), intent(in) :: startdate @@ -5079,9 +5085,11 @@ subroutine output_lakes2(igrid, split_output_count, NLAKES, & integer, intent(in) :: igrid, K integer, intent(in) :: split_output_count integer, intent(in) :: NLAKES - real, dimension(NLAKES), intent(in) :: latlake,lonlake,elevlake,resht - real, dimension(NLAKES), intent(in) :: qlakei,qlakeo !-- inflow and outflow of lake - integer(kind=int64), dimension(NLAKES), intent(in) :: LAKEIDM !-- LAKE ID + real, dimension(NLAKES), intent(in) :: latlake,lonlake,elevlake + real*8, dimension(NLAKES), intent(in) :: RESHT + + real*8, dimension(NLAKES), intent(in) :: qlakei,qlakeo !-- inflow and outflow of lake + integer(kind=int64), dimension(NLAKES), intent(in) :: LAKEIDM !-- LAKE ID real, intent(in) :: dtrt_ch character(len=*), intent(in) :: startdate @@ -6153,13 +6161,13 @@ subroutine RESTART_OUT_nc(outFile,did) #endif ) - call w_rst_crt_nc1_lake(ncid,rt_domain(did)%nlakes,rt_domain(did)%qlakeo,"qlakeo" & + call w_rst_crt_nc1_lake(ncid,rt_domain(did)%nlakes,dble(rt_domain(did)%qlakeo),"qlakeo" & #ifdef MPP_LAND ,rt_domain(did)%lake_index & #endif ) - call w_rst_crt_nc1_lake(ncid,rt_domain(did)%nlakes,rt_domain(did)%qlakei,"qlakei" & + call w_rst_crt_nc1_lake(ncid,rt_domain(did)%nlakes,dble(rt_domain(did)%qlakei),"qlakei" & #ifdef MPP_LAND ,rt_domain(did)%lake_index & #endif @@ -6527,12 +6535,13 @@ subroutine w_rst_crt_nc1_lake(ncid,n,inVar,varName & implicit none integer:: ncid,n,varid , iret character(len=*) varName - real inVar(n) + real*8 inVar(n) #ifdef MPP_LAND integer:: nodelist(n) if(n .eq. 0) return - call write_lake_real(inVar,nodelist,n) + call write_lake_real8(inVar,nodelist,n) + if(my_id .eq. IO_id) then #endif iret = nf90_inq_varid(ncid,varName, varid) @@ -6691,7 +6700,7 @@ subroutine read_rst_gwbucket_real(ncid,outV,numbasns,& character(len=*) :: vName integer i, j,k real, dimension(gnumbasns) :: buf - call read_rst_crt_nc(ncid,buf,gnumbasns,vName) + call read_rst_crt_nc(ncid,dble(buf),gnumbasns,vName) do k = 1, numbasns outV(k) = buf(basnsInd(k)) end do @@ -6866,7 +6875,7 @@ subroutine RESTART_IN_NC(inFile,did) endif if(rt_domain(did)%NLAKES .gt. 0) then - call read_rst_crt_nc(ncid,rt_domain(did)%RESHT,rt_domain(did)%NLAKES,"resht") + ! call read_rst_crt_nc(ncid,rt_domain(did)%RESHT,rt_domain(did)%NLAKES,"resht") call read_rst_crt_nc(ncid,rt_domain(did)%QLAKEO,rt_domain(did)%NLAKES,"qlakeo") call read_rst_crt_nc(ncid,rt_domain(did)%QLAKEI,rt_domain(did)%NLAKES,"qlakei") endif @@ -7142,7 +7151,8 @@ end subroutine read_rt_nc2 subroutine read_rst_crt_nc(ncid,var,n,varStr) implicit none integer :: ireg, ncid, varid, n, iret - real,dimension(n) :: var + !real,dimension(n) :: var + real*8,dimension(n) :: var character(len=*) :: varStr if( n .le. 0) return @@ -7169,7 +7179,9 @@ subroutine read_rst_crt_nc(ncid,var,n,varStr) #ifdef MPP_LAND endif if(n .gt. 0) then - call mpp_land_bcast_real(n,var) + call mpp_land_bcast_real(n,real(var)) + !call mpp_land_bcast_real(n,var) + endif #endif return diff --git a/src/Routing/module_NWM_io.F b/src/Routing/module_NWM_io.F index ffce2e60c..7b5605cff 100644 --- a/src/Routing/module_NWM_io.F +++ b/src/Routing/module_NWM_io.F @@ -2480,18 +2480,23 @@ subroutine output_lakes_NWM(domainId,iGrid) integer :: ierr, myId ! MPI return status, process ID integer :: ascFlag ! Flag for resorting timeseries output by feature_id. ! Allocatable arrays to hold output variables. - real, allocatable, dimension(:) :: g_lakeLat,g_lakeLon,g_lakeElev - real, allocatable, dimension(:) :: g_lakeInflow,g_lakeOutflow + real, allocatable, dimension(:) :: g_lakeLat,g_lakeLon!g_lakeElev + real*8, allocatable, dimension(:) :: g_lakeElev + + real*8, allocatable, dimension(:) :: g_lakeInflow,g_lakeOutflow real, allocatable, dimension(:) :: g_lakeType real, allocatable, dimension(:) :: g_lake_assimilated_value character(len=256), allocatable, dimension(:) :: g_lake_assimilated_source_file integer(kind=int64), allocatable, dimension(:) :: g_lakeid - real, allocatable, dimension(:) :: g_lakeLatOut,g_lakeLonOut,g_lakeElevOut - real, allocatable, dimension(:) :: g_lakeInflowOut,g_lakeOutflowOut + real, allocatable, dimension(:) :: g_lakeLatOut,g_lakeLonOut + real*8, allocatable, dimension(:) :: g_lakeElevOut + real*8, allocatable, dimension(:) :: g_lakeInflowOut,g_lakeOutflowOut integer(kind=int64), allocatable, dimension(:) :: g_lakeTypeOut, g_lakeidOut real, allocatable, dimension(:) :: g_lake_assimilated_valueOut character(len=256), allocatable, dimension(:) :: g_lake_assimilated_source_fileOut real, allocatable, dimension(:,:) :: varOutReal ! Array holding output variables in real format + real*8, allocatable, dimension(:,:) :: varOutReal8 ! Array holding output variables in real format + integer, allocatable, dimension(:) :: varOutInt ! Array holding output variables after ! scale_factor/add_offset ! have been applied. @@ -2623,16 +2628,18 @@ subroutine output_lakes_NWM(domainId,iGrid) ! Collect arrays from various processors through MPI, and ! assemble into global arrays previously allocated. #ifdef MPP_LAND + call write_lake_real(g_lakeLat,RT_DOMAIN(domainId)%lake_index,gsize) call write_lake_real(g_lakeLon,RT_DOMAIN(domainId)%lake_index,gsize) - call write_lake_real(g_lakeElev,RT_DOMAIN(domainId)%lake_index,gsize) - call write_lake_real(g_lakeInflow,RT_DOMAIN(domainId)%lake_index,gsize) - call write_lake_real(g_lakeOutflow,RT_DOMAIN(domainId)%lake_index,gsize) + call write_lake_real8(g_lakeElev,RT_DOMAIN(domainId)%lake_index,gsize) + call write_lake_real8(g_lakeInflow,RT_DOMAIN(domainId)%lake_index,gsize) + call write_lake_real8(g_lakeOutflow,RT_DOMAIN(domainId)%lake_index,gsize) call write_lake_real(g_lakeType,RT_DOMAIN(domainId)%lake_index,gsize) call write_lake_real(g_lake_assimilated_value,RT_DOMAIN(domainId)%lake_index,gsize) !call write_lake_char(g_lake_assimilated_source_file,RT_DOMAIN(domainId)%lake_index,gsize) #endif + else gSize = rt_domain(domainId)%NLAKES ! No MPI - single processor @@ -2663,6 +2670,8 @@ subroutine output_lakes_NWM(domainId,iGrid) g_lake_assimilated_source_file = RT_DOMAIN(domainID)%reservoir_assimilated_source_file endif + print *, "LAKEINOUTAGAIN",g_lakeInflow,g_lakeOutflow,gsize + ! Sync all processes up. if(mppFlag .eq. 1) then #ifdef MPP_LAND @@ -2670,6 +2679,9 @@ subroutine output_lakes_NWM(domainId,iGrid) #endif endif + print *, "AFTERSYNCin1",g_lakeInflow,gsize + print *, "AFTERSYNCout1",g_lakeOutflow,gsize + ! Calculate datetime information. ! First compose strings of EPOCH and simulation start date. epochDate = trim("1970-01-01 00:00") @@ -2735,10 +2747,18 @@ subroutine output_lakes_NWM(domainId,iGrid) iret = nf90_close(ftnRt) call nwmCheck(diagFlag,iret,'ERROR: Unable to close LAKEPARM file.') + print *, "AFTERSYNCin2",g_lakeInflow,gsize + print *, "AFTERSYNCout2",g_lakeOutflow,gsize + ! Place all output arrays into one real array that will be looped over ! during conversion to compressed integer format. allocate(varOutReal(fileMeta%numVars,gSize)) + allocate(varOutReal8(fileMeta%numVars,gSize)) allocate(varOutInt(gSize)) + + print *, "AFTERSYNCin3",g_lakeInflow,gsize + print *, "AFTERSYNCout3",g_lakeOutflow,gsize + if(ascFlag .eq. 1) then ! Sort feature_id values by ascending values using the index array ! extracted from the RouteLink file. @@ -2754,6 +2774,11 @@ subroutine output_lakes_NWM(domainId,iGrid) g_lake_assimilated_valueOut(iTmp) = g_lake_assimilated_value(indTmp) g_lake_assimilated_source_fileOut(iTmp) = g_lake_assimilated_source_file(indTmp) g_lakeidOut(iTmp) = g_lakeid(indTmp) + + print *, "AGAINIn",g_lakeInflow + print *, "AGAINInOut",g_lakeInflowOut + print *, "AGAINlke",g_lakeElevOut + end do else g_lakeInflowOut = g_lakeInflow @@ -2766,11 +2791,15 @@ subroutine output_lakes_NWM(domainId,iGrid) g_lake_assimilated_source_fileOut = g_lake_assimilated_source_file g_lakeidOut = g_lakeid endif - varOutReal(1,:) = g_lakeInflowOut(:) - varOutReal(2,:) = g_lakeOutflowOut(:) + + varOutReal8(1,:) = g_lakeInflowOut(:) + varOutReal8(2,:) = g_lakeOutflowOut(:) ! Mask out missing values - where ( varOutReal == fileMeta%modelNdv ) varOutReal = -9999.0 + where ( varOutReal == fileMeta%modelNdv ) varOutReal = -9999.0 !! might not need this here, replaced below + where ( varOutReal8 == fileMeta%modelNdv ) varOutReal8 = -9999.0 + + print *, "LAKEINOUTVAR",g_lakeInflowOut,g_lakeOutflowOut,varOutReal8(1,:) iret = nf90_create(trim(output_flnm),cmode=NF90_NETCDF4,ncid = ftn) call nwmCheck(diagFlag,iret,'ERROR: Unable to create LAKEOUT NetCDF file.') @@ -2924,7 +2953,8 @@ subroutine output_lakes_NWM(domainId,iGrid) call nwmCheck(diagFlag,iret,'ERROR: Unable to place units attribute into longitude variable') ! Create channel elevation variable - iret = nf90_def_var(ftn,"water_sfc_elev",nf90_float,dimId(1),elevVarId) + !iret = nf90_def_var(ftn,"water_sfc_elev",nf90_float,dimId(1),elevVarId) + iret = nf90_def_var(ftn,"water_sfc_elev",nf90_double,dimId(1),elevVarId) call nwmCheck(diagFlag,iret,'ERROR: Unable to create water_sfc_elev variable.') iret = nf90_put_att(ftn,elevVarId,'long_name',trim(fileMeta%elevLName)) call nwmCheck(diagFlag,iret,'ERROR: Unable to place long_name attribute into water_sfc_elev variable') @@ -3028,6 +3058,7 @@ subroutine output_lakes_NWM(domainId,iGrid) if(minSinceSim .eq. 0 .and. fileMeta%timeZeroFlag(iTmp) .eq. 0) then varOutInt(:) = fileMeta%fillComp(iTmp) varOutReal(iTmp,:) = fileMeta%fillReal(iTmp) + varOutReal8(iTmp,:) = fileMeta%fillReal(iTmp) else varOutInt(:) = NINT((varOutReal(iTmp,:)-fileMeta%addOffset(iTmp))/fileMeta%scaleFactor(iTmp)) endif @@ -3038,7 +3069,12 @@ subroutine output_lakes_NWM(domainId,iGrid) if((nlst(1)%io_form_outputs .eq. 1) .or. (nlst(1)%io_form_outputs .eq. 2)) then iret = nf90_put_var(ftn,varId,varOutInt) else + if(iTmp .gt. 2) then iret = nf90_put_var(ftn,varId,varOutReal(iTmp,:)) + else + print *, "THIS GOT HERE", iTmp,varOutReal8(iTmp,:) + iret = nf90_put_var(ftn,varId,varOutReal8(iTmp,:)) + endif endif call nwmCheck(diagFlag,iret,'ERROR: Unable to place data into output variable: '//trim(fileMeta%varNames(iTmp))) endif @@ -3112,6 +3148,7 @@ subroutine output_lakes_NWM(domainId,iGrid) ! Deallocate all memory if(myId .eq. 0) then deallocate(varOutReal) + deallocate(varOutReal8) deallocate(varOutInt) endif deallocate(g_lakeLon) diff --git a/src/Routing/module_RT.F b/src/Routing/module_RT.F index fa3ef00c0..6bd43fa91 100644 --- a/src/Routing/module_RT.F +++ b/src/Routing/module_RT.F @@ -688,7 +688,7 @@ subroutine LandRT_ini(did) integer, allocatable, dimension(:) :: tmp_int real, allocatable, dimension(:) :: tmp_real integer, allocatable, dimension(:) :: buf - real, allocatable, dimension(:) :: tmpRESHT + real*8, allocatable, dimension(:) :: tmpRESHT integer :: new_start_i, new_start_j, new_end_i, new_end_j #ifdef OUTPUT_CHAN_CONN diff --git a/src/Routing/module_channel_routing.F b/src/Routing/module_channel_routing.F index ff26557e4..6674936d2 100644 --- a/src/Routing/module_channel_routing.F +++ b/src/Routing/module_channel_routing.F @@ -163,7 +163,7 @@ subroutine SUBMUSKINGCUNGE( & REAL, intent(IN) :: qup ! flow upstream previous timestep REAL, intent(IN) :: quc ! flow upstream current timestep REAL, intent(IN) :: qdp ! flow downstream previous timestep - REAL, intent(INOUT) :: qdc ! flow downstream current timestep + REAL*8, intent(INOUT) :: qdc ! flow downstream current timestep REAL, intent(IN) :: ql ! lateral inflow through reach (m^3/sec) REAL, intent(IN) :: Bw ! bottom width (meters) REAL, intent(IN) :: Tw ! top width before bankfull (meters) @@ -613,7 +613,8 @@ Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & REAL, INTENT(IN), DIMENSION(NLINKS) :: ChannK !--Channel Infiltration REAL :: Km, X REAL , INTENT(INOUT), DIMENSION(:,:) :: QLINK - REAL , DIMENSION(NLINKS,2) :: tmpQLINK + REAL*8, DIMENSION(NLINKS,2) :: tmpQLINK + REAL , INTENT(INOUT), DIMENSION(NLINKS) :: HLINK REAL, dimension(NLINKS), intent(inout) :: QLateral !--lateral flow REAL, INTENT(IN) :: DT !-- model timestep @@ -639,11 +640,11 @@ Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & !-- lake params - REAL, INTENT(INOUT), DIMENSION(NLAKES) :: RESHT !-- reservoir height (m) + REAL*8, INTENT(INOUT), DIMENSION(NLAKES) :: RESHT !-- reservoir height (m) REAL*8, DIMENSION(NLAKES) :: QLAKEI8 !-- lake inflow (cms) - REAL, INTENT(INOUT), DIMENSION(NLAKES) :: QLAKEI !-- lake inflow (cms) - REAL, DIMENSION(NLAKES) :: QLAKEIP !-- lake inflow previous timestep (cms) - REAL, INTENT(INOUT), DIMENSION(NLAKES) :: QLAKEO !-- outflow from lake used in diffusion scheme + REAL*8, INTENT(INOUT), DIMENSION(NLAKES) :: QLAKEI !-- lake inflow (cms) + REAL*8, DIMENSION(NLAKES) :: QLAKEIP !-- lake inflow previous timestep (cms) + REAL*8, INTENT(INOUT), DIMENSION(NLAKES) :: QLAKEO !-- outflow from lake used in diffusion scheme INTEGER(kind=int64), INTENT(IN), DIMENSION(NLINKS) :: LAKENODE !-- outflow from lake used in diffusion scheme INTEGER(kind=int64), INTENT(IN), DIMENSION(NLINKS) :: LINKID !-- id of channel elements for linked scheme !REAL, DIMENSION(NLINKS) :: QLateral !--lateral flow @@ -685,15 +686,15 @@ Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & integer :: n, kk2, nt, nsteps ! tmp - QLAKEIP = 0 - QLAKEI8 = 0 - HLINKTMP = 0 - CVOLTMP = 0 - CD = 0 + QLAKEIP = 0. + QLAKEI8 = 0. + HLINKTMP = 0. + CVOLTMP = 0. + CD = 0. node = 1 - QLateral = 0 - QSUM = 0 - QLLAKE = 0 + QLateral = 0. + QSUM = 0. + QLLAKE = 0. !yw print *, "DRIVE_channel,option,nlinkl,nlinks!!", channel_option,NLINKSL,NLINKS @@ -802,22 +803,21 @@ Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & do n = 1, gtoNODE(k,1) m = gtoNODE(k,n+1) !yw if (LINKID(k) .eq. m) then - Quc = Quc + gQLINK(m,2) !--accum of upstream inflow of current timestep (2) - Qup = Qup + gQLINK(m,1) !--accum of upstream inflow of previous timestep (1) - - ! if(LINKID(k) .eq. 3259 .or. LINKID(k) .eq. 3316 .or. LINKID(k) .eq. 3219) then - ! write(6,*) "id,Uc,Up",LINKID(k),Quc,Qup - ! call flush(6) - ! endif + !Quc = Quc + gQLINK(m,2) !--accum of upstream inflow of current timestep (2) + !Qup = Qup + gQLINK(m,1) !--accum of upstream inflow of previous timestep (1) + Quc = Quc + gQLINK(m,1) !--accum of upstream inflow of current timestep (2) + Qup = Qup + gQLINK(m,2) !--accum of upstream inflow of previous timestep (1) !yw endif end do ! do i #else do m = 1, NLINKSL if (LINKID(k) .eq. TO_NODE(m)) then - Quc = Quc + QLINK(m,2) !--accum of upstream inflow of current timestep (2) - Qup = Qup + QLINK(m,1) !--accum of upstream inflow of previous timestep (1) + ! Quc = Quc + QLINK(m,2) !--accum of upstream inflow of current timestep (2) + ! Qup = Qup + QLINK(m,1) !--accum of upstream inflow of previous timestep (1) + Quc = Quc + QLINK(m,1) !--accum of upstream inflow of current timestep (2) + Qup = Qup + QLINK(m,2) !--accum of upstream inflow of previous timestep (1) endif end do ! do m #endif @@ -831,7 +831,7 @@ Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & elseif (channel_option .eq. 1) then !muskingum routing Km = MUSK(k) X = MUSX(k) - tmpQLINK(k,2) = MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X) !upstream plust lateral inflow + tmpQLINK(k,2) = dble(MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X)) !upstream plust lateral inflow elseif (channel_option .eq. 2) then ! muskingum cunge call SUBMUSKINGCUNGE(tmpQLINK(k,2), velocity(k), qloss(k), LINKID(k), & @@ -1222,6 +1222,7 @@ Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & ! in whether it uses level pool, machine learning, or persistence. ! Inflow, lateral inflow, water elevation, and the routing period ! are passed in. Updated outflow and water elevation are returned. + call rt_domain(did)%reservoirs(i)%ptr%run(QLAKEIP(i), QLAKEI(i), & QLLAKE(i), RESHT(i), QLAKEO(i), DTCT, rt_domain(did)%final_reservoir_type(i), & rt_domain(did)%reservoir_assimilated_value(i), rt_domain(did)%reservoir_assimilated_source_file(i)) @@ -1247,11 +1248,11 @@ Subroutine drive_CHANNEL(did, latval,lonval,KT, IXRT,JXRT, SUBRTSWCRT, & ! TODO: Change updateLake_grid calls below to updated distributed reservoir ! objects and not arrays as currently implemented. - call updateLake_grid(QLLAKE, nlakes,lake_index) - call updateLake_grid(RESHT, nlakes,lake_index) - call updateLake_grid(QLAKEO, nlakes,lake_index) - call updateLake_grid(QLAKEI, nlakes,lake_index) - call updateLake_grid(QLAKEIP,nlakes,lake_index) + call updateLake_grid(dble(QLLAKE), nlakes,lake_index) + call updateLake_grid(dble(RESHT), nlakes,lake_index) + call updateLake_grid(dble(QLAKEO), nlakes,lake_index) + call updateLake_grid(dble(QLAKEI), nlakes,lake_index) + call updateLake_grid(dble(QLAKEIP),nlakes,lake_index) endif #endif @@ -1649,7 +1650,7 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & real, intent(IN), dimension(:) :: ChSSlp,Bw,Tw !--properties of nodes or links real, intent(IN), dimension(:) :: Tw_CC, n_CC ! properties of compound channel real, intent(IN), dimension(:) :: ChannK !--channel infiltration - real :: Km, X + real :: Km, X real , intent(INOUT), dimension(:,:) :: QLINK real , intent(INOUT), dimension(:) :: HLINK @@ -1662,29 +1663,29 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & real, dimension(:), intent(inout) :: QLateral !--lateral flow real, dimension(:), intent(out) :: velocity, qloss - real*8, dimension(:), intent(inout) :: accSfcLatRunoff, accBucket + real*8, dimension(:), intent(inout) :: accSfcLatRunoff, accBucket real , dimension(:), intent(out) :: qSfcLatRunoff , qBucket + real*8, dimension(NLINKSL,2) :: tmpQLINK - real , dimension(NLINKSL,2) :: tmpQLINK - real, intent(IN) :: DT !-- model timestep - real, intent(IN) :: DTRT_CH !-- routing timestep - real, intent(INOUT) :: DTCT - real :: minDTCT !BF minimum routing timestep - integer, intent(IN) :: MAXORDER - real , intent(IN), dimension(:) :: node_area + real, intent(IN) :: DT !-- model timestep + real, intent(IN) :: DTRT_CH !-- routing timestep + real, intent(INOUT) :: DTCT + real :: minDTCT !BF minimum routing timestep + integer, intent(IN) :: MAXORDER + real , intent(IN), dimension(:) :: node_area !DJG GW-chan coupling variables... - real, dimension(NLINKS) :: dzGwChanHead - real, dimension(NLINKS) :: Q_GW_CHAN_FLUX !DJG !!! Change 'INTENT' to 'OUT' when ready to update groundwater state... - real, dimension(IXRT,JXRT) :: ZWATTBLRT !DJG !!! Match with subsfce/gw routing & Change 'INTENT' to 'INOUT' when ready to update groundwater state... + real, dimension(NLINKS) :: dzGwChanHead + real, dimension(NLINKS) :: Q_GW_CHAN_FLUX !DJG !!! Change 'INTENT' to 'OUT' when ready to update groundwater state... + real, dimension(IXRT,JXRT) :: ZWATTBLRT !DJG !!! Match with subsfce/gw routing & Change 'INTENT' to 'INOUT' when ready to update groundwater state... !-- lake params integer(kind=int64), intent(IN), dimension(:) :: LAKEIDM !-- NHDPLUS lakeid for lakes to be modeled - real, intent(INOUT), dimension(:) :: RESHT !-- reservoir height (m) - real, intent(INOUT), dimension(:) :: QLAKEI !-- lake inflow (cms) - real, dimension(NLAKES) :: QLAKEIP !-- lake inflow previous timestep (cms) - real, intent(INOUT), dimension(NLAKES) :: QLAKEO !-- outflow from lake used in diffusion scheme + real*8, intent(INOUT), dimension(:) :: RESHT !-- reservoir height (m) + real*8, intent(INOUT), dimension(:) :: QLAKEI !-- lake inflow (cms) + real*8, dimension(NLAKES) :: QLAKEIP !-- lake inflow previous timestep (cms) + real*8, intent(INOUT), dimension(NLAKES) :: QLAKEO !-- outflow from lake used in diffusion scheme integer(kind=int64), intent(IN), dimension(:) :: LAKENODE !-- outflow from lake used in diffusion scheme integer(kind=int64), intent(IN), dimension(:) :: LINKID !-- id of channel elements for linked scheme @@ -1729,7 +1730,8 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & integer :: n, kk2, nt, nsteps ! tmp real, intent(in), dimension(:) :: qout_gwsubbas - real, allocatable,dimension(:) :: tmpQLAKEO, tmpQLAKEI, tmpRESHT + real*8, allocatable,dimension(:) :: tmpQLAKEO, tmpQLAKEI + real*8, allocatable,dimension(:) :: tmpRESHT integer, allocatable, dimension(:) :: tmpFinalResType real, allocatable,dimension(:) :: tmpAssimilatedValue character(len=256), allocatable,dimension(:) :: tmpAssimilatedSourceFile @@ -1745,11 +1747,11 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & endif #endif - QLAKEIP = 0 - CD = 0 + QLAKEIP = 0. + CD = 0. node = 1 - QSUM = 0 - QLLAKE = 0 + QSUM = 0. + QLLAKE = 0. dzGwChanHead = 0. nsteps = (DT+0.5)/DTRT_CH @@ -1789,11 +1791,11 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & else - QLateral = 0 !! the variable solved in this section. Channel only knows this already. - LQLateral = 0 !-- initial lateral flow to 0 for this reach + QLateral = 0. !! the variable solved in this section. Channel only knows this already. + LQLateral = 0. !-- initial lateral flow to 0 for this reach - tmpQLateral = 0 !! WHY DOES THIS tmp variable EXIST?? Only for accumulations?? - tmpLQLateral = 0 + tmpQLateral = 0. !! WHY DOES THIS tmp variable EXIST?? Only for accumulations?? + tmpLQLateral = 0. ! NHDPLUS maping if(OVRTSWCRT .eq. 0) then @@ -1893,7 +1895,8 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & do nt = 1, nsteps #ifdef MPP_LAND - gQLINK = 0 + +! gQLINK = 0 call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,2), NLINKSL, gQLINK(:,2)) call gbcastReal2(toNodeInd,nToNodeInd,QLINK(1:NLINKSL,1), NLINKSL, gQLINK(:,1)) !---------- route other reaches, with upstream inflow @@ -1916,8 +1919,8 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & do k = 1,NLINKSL - Quc = 0 - Qup = 0 + Quc = 0. + Qup = 0. !process as standard link or a lake inflow link, or lake outflow link ! link flowing out of lake, accumulate all the inflows with the revised TO_NODEs @@ -1933,8 +1936,10 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & !! current would require global communication at the end of each loop through k !! (=kth reach). Additionally, how do you synchronize to make sure the upstream are all !! done before doing the downstream? + if(gQLINK(m,2) .gt. 0) Quc = Quc + gQLINK(m,2) !--accum of upstream inflow of current timestep (2) if(gQLINK(m,1) .gt. 0) Qup = Qup + gQLINK(m,1) !--accum of upstream inflow of previous timestep (1) + end do ! do i #else do m = 1, NLINKSL @@ -1954,12 +1959,14 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & lakeid = LAKEIDX(k) if(lakeid .ge. 0) then + + Qup = QLAKEI(lakeid) ! Calls run for a single reservoir depending on its type as ! in whether it uses level pool, machine learning, or persistence. ! Inflow, lateral inflow, water elevation, and the routing period ! are passed in. Updated outflow and water elevation are returned. - call rt_domain(did)%reservoirs(lakeid)%ptr%run(Qup, Quc, 0.0, & + call rt_domain(did)%reservoirs(lakeid)%ptr%run(dble(Qup), dble(Quc), 0.0, & RESHT(lakeid), tmpQLINK(k,2), DTRT_CH, rt_domain(did)%final_reservoir_type(lakeid), & rt_domain(did)%reservoir_assimilated_value(lakeid), rt_domain(did)%reservoir_assimilated_source_file(lakeid)) @@ -1968,7 +1975,10 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & ! MPI communication model. QLAKEO(lakeid) = tmpQLINK(k,2) !save outflow to lake - QLAKEI(lakeid) = Quc !save inflow to lake + QLAKEI(lakeid) = dble(Quc) !save inflow to lake + + !print *, "QLAKEOI", QLAKEO(lakeid),Quc + endif 105 continue @@ -1976,7 +1986,7 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & elseif (channel_option .eq. 1) then !muskingum routing Km = MUSK(k) X = MUSX(k) - tmpQLINK(k,2) = MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X) !upstream plust lateral inflow + tmpQLINK(k,2) = dble(MUSKING(k,Qup,(Quc+QLateral(k)),QLINK(k,1),DTRT_CH,Km,X)) !upstream plust lateral inflow elseif (channel_option .eq. 2) then ! muskingum cunge, don't process internal lake nodes TYP=2 ! tmpQLINK(k,2) = MUSKINGCUNGE(k,Qup, Quc, QLINK(k,1), & @@ -2003,8 +2013,8 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & #ifdef MPP_LAND - call updateLake_seq(QLAKEO,nlakes,tmpQLAKEO) - call updateLake_seq(QLAKEI,nlakes,tmpQLAKEI) + call updateLake_seq(dble(QLAKEO),nlakes,dble(tmpQLAKEO)) + call updateLake_seq(dble(QLAKEI),nlakes,dble(tmpQLAKEI)) call updateLake_seq(RESHT,nlakes,tmpRESHT) call updateLake_seqInt(rt_domain(did)%final_reservoir_type, nlakes, tmpFinalResType) call updateLake_seq(rt_domain(did)%reservoir_assimilated_value, nlakes, tmpAssimilatedValue) @@ -2013,7 +2023,7 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & do k = 1, NLINKSL !tmpQLINK? if(TYPEL(k) .ne. 2) then !only the internal lake nodes don't have info.. but need to save QLINK of lake out too - QLINK(k,2) = tmpQLINK(k,2) + QLINK(k,2) = real(tmpQLINK(k,2)) endif QLINK(k,1) = QLINK(k,2) !assigng link flow of current to be previous for next time step end do From 2efcfb1bddb48d3ccdea2ef50bd79421e4b6bbe6 Mon Sep 17 00:00:00 2001 From: David Yates Date: Wed, 16 Sep 2020 08:32:27 -0600 Subject: [PATCH 2/5] double precision on lake variables --- src/Routing/module_channel_routing.F | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Routing/module_channel_routing.F b/src/Routing/module_channel_routing.F index 6674936d2..ff2d1b299 100644 --- a/src/Routing/module_channel_routing.F +++ b/src/Routing/module_channel_routing.F @@ -2011,10 +2011,9 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & end do !--k links - #ifdef MPP_LAND - call updateLake_seq(dble(QLAKEO),nlakes,dble(tmpQLAKEO)) - call updateLake_seq(dble(QLAKEI),nlakes,dble(tmpQLAKEI)) + call updateLake_seq(QLAKEO,nlakes,tmpQLAKEO) + call updateLake_seq(QLAKEI,nlakes,tmpQLAKEI) call updateLake_seq(RESHT,nlakes,tmpRESHT) call updateLake_seqInt(rt_domain(did)%final_reservoir_type, nlakes, tmpFinalResType) call updateLake_seq(rt_domain(did)%reservoir_assimilated_value, nlakes, tmpAssimilatedValue) From 2864a14c85b830de15ba971de89dbde29f4a96e9 Mon Sep 17 00:00:00 2001 From: Bahram Khazaei Date: Mon, 21 Sep 2020 08:51:23 -0600 Subject: [PATCH 3/5] few changes made to the code to pass the perfect restarts test --- src/Routing/module_HYDRO_io.F | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Routing/module_HYDRO_io.F b/src/Routing/module_HYDRO_io.F index bb8bb3b9a..a03b1f8f4 100644 --- a/src/Routing/module_HYDRO_io.F +++ b/src/Routing/module_HYDRO_io.F @@ -5961,9 +5961,9 @@ subroutine RESTART_OUT_nc(outFile,did) if(nlst(did)%channel_option .eq. 3) & iret = nf90_def_var(ncid, "cvol", NF90_FLOAT, (/dimid_links/), varid) if(rt_domain(did)%nlakes .gt. 0) then - iret = nf90_def_var(ncid, "resht", NF90_FLOAT, (/dimid_lakes/), varid) - iret = nf90_def_var(ncid, "qlakeo", NF90_FLOAT, (/dimid_lakes/), varid) - iret = nf90_def_var(ncid, "qlakei", NF90_FLOAT, (/dimid_lakes/), varid) + iret = nf90_def_var(ncid, "resht", NF90_DOUBLE, (/dimid_lakes/), varid) + iret = nf90_def_var(ncid, "qlakeo", NF90_DOUBLE, (/dimid_lakes/), varid) + iret = nf90_def_var(ncid, "qlakei", NF90_DOUBLE, (/dimid_lakes/), varid) endif if( nlst(did)%channel_only .eq. 0 .and. & @@ -6875,7 +6875,7 @@ subroutine RESTART_IN_NC(inFile,did) endif if(rt_domain(did)%NLAKES .gt. 0) then - ! call read_rst_crt_nc(ncid,rt_domain(did)%RESHT,rt_domain(did)%NLAKES,"resht") + call read_rst_crt_nc(ncid,rt_domain(did)%RESHT,rt_domain(did)%NLAKES,"resht") call read_rst_crt_nc(ncid,rt_domain(did)%QLAKEO,rt_domain(did)%NLAKES,"qlakeo") call read_rst_crt_nc(ncid,rt_domain(did)%QLAKEI,rt_domain(did)%NLAKES,"qlakei") endif @@ -7179,7 +7179,7 @@ subroutine read_rst_crt_nc(ncid,var,n,varStr) #ifdef MPP_LAND endif if(n .gt. 0) then - call mpp_land_bcast_real(n,real(var)) + call mpp_land_bcast(var) !call mpp_land_bcast_real(n,var) endif From 2f91cb7881a8d605c8a9c6ca2cb6df2a9078966f Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Tue, 7 Dec 2021 15:21:06 -0700 Subject: [PATCH 4/5] Update reservoir assimilated values to doubles --- src/Data_Rec/rt_include.inc | 2 +- src/Routing/Reservoirs/Level_Pool/module_levelpool.F | 2 +- .../module_persistence_levelpool_hybrid.F | 2 +- .../module_persistence_levelpool_hybrid_state.F | 2 +- src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F | 2 +- .../Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F | 2 +- src/Routing/Reservoirs/module_reservoir.F | 2 +- src/Routing/module_channel_routing.F | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Data_Rec/rt_include.inc b/src/Data_Rec/rt_include.inc index 51f960301..484519e4e 100644 --- a/src/Data_Rec/rt_include.inc +++ b/src/Data_Rec/rt_include.inc @@ -8,7 +8,7 @@ class (reservoir_container), allocatable, dimension(:) :: reservoirs integer, allocatable, dimension(:) :: reservoir_type ! specifying type of reservoir integer, allocatable, dimension(:) :: final_reservoir_type ! resolved reservoir type (since it can change) - real, allocatable, dimension(:) :: reservoir_assimilated_value ! observation or forecast assimilated to reservoir discharge + real*8, allocatable, dimension(:) :: reservoir_assimilated_value ! observation or forecast assimilated to reservoir discharge character(len=256), allocatable, dimension(:) :: reservoir_assimilated_source_file ! source file of assimilated value INTEGER :: IX, JX diff --git a/src/Routing/Reservoirs/Level_Pool/module_levelpool.F b/src/Routing/Reservoirs/Level_Pool/module_levelpool.F index 9fadf8c93..d6c2d6ce6 100644 --- a/src/Routing/Reservoirs/Level_Pool/module_levelpool.F +++ b/src/Routing/Reservoirs/Level_Pool/module_levelpool.F @@ -157,7 +157,7 @@ subroutine run_levelpool_reservoir(this, previous_timestep_inflow, inflow, & real*8, intent(out) :: outflow ! cubic meters per second (cms) real, intent(in) :: routing_period ! seconds integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files - real, intent(out) :: assimilated_value ! value assimilated from observation or forecast + real*8, intent(out) :: assimilated_value ! value assimilated from observation or forecast character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value diff --git a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F index 22612447a..82171a731 100644 --- a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F +++ b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid.F @@ -227,7 +227,7 @@ subroutine run_hybrid_reservoir(this, previous_timestep_inflow, inflow, & real*8, intent(out) :: outflow ! cubic meters per second (cms) real, intent(in) :: routing_period ! seconds integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files - real, intent(out) :: assimilated_value ! value assimilated from observation + real*8, intent(out) :: assimilated_value ! value assimilated from observation character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value real :: delta_storage ! timestep change in storage (cubic meters) real*8 :: local_water_elevation ! water elevation passed to levelpool (meters AMSL) diff --git a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F index 424776b94..4d7faa9f9 100644 --- a/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F +++ b/src/Routing/Reservoirs/Persistence_Level_Pool_Hybrid/module_persistence_levelpool_hybrid_state.F @@ -28,7 +28,7 @@ module module_persistence_levelpool_hybrid_state real :: assimilated_value character(len=256) :: assimilated_source_file integer :: levelpool_reservoir_type ! reservoir type for levelpool - real :: levelpool_assimilated_value ! levelpool assimilated sentinel value + real*8 :: levelpool_assimilated_value ! levelpool assimilated sentinel value character(len=256) :: levelpool_assimilated_source_file ! none for levelpool assimilated source type (levelpool), pointer :: levelpool_ptr ! pointer to levelpool object diff --git a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F index 7aa568260..08ea7ad9d 100644 --- a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F +++ b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts.F @@ -213,7 +213,7 @@ subroutine run_rfc_forecasts_reservoir(this, previous_timestep_inflow, inflow, & real*8, intent(out) :: outflow ! cubic meters per second (cms) real, intent(in) :: routing_period ! seconds integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files - real, intent(out) :: assimilated_value ! value assimilated from observation or forecast + real*8, intent(out) :: assimilated_value ! value assimilated from observation or forecast character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value real*8 :: levelpool_outflow ! cubic meters per second (cms) integer :: missing_outflow_index diff --git a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F index 779dc1a04..43062c341 100644 --- a/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F +++ b/src/Routing/Reservoirs/RFC_Forecasts/module_rfc_forecasts_state.F @@ -27,7 +27,7 @@ module module_rfc_forecasts_state real :: assimilated_value character(len=256) :: assimilated_source_file integer :: levelpool_reservoir_type ! reservoir type for levelpool - real :: levelpool_assimilated_value ! levelpool assimilated sentinel value + real*8 :: levelpool_assimilated_value ! levelpool assimilated sentinel value character(len=256) :: levelpool_assimilated_source_file ! empty string for levelpool assimilated source type (levelpool), pointer :: levelpool_ptr ! pointer to levelpool object diff --git a/src/Routing/Reservoirs/module_reservoir.F b/src/Routing/Reservoirs/module_reservoir.F index d0517e192..6e7b51fbf 100644 --- a/src/Routing/Reservoirs/module_reservoir.F +++ b/src/Routing/Reservoirs/module_reservoir.F @@ -103,7 +103,7 @@ subroutine run_reservoir_interface(this, previous_timestep_inflow, inflow, & real*8, intent(out) :: outflow ! cubic meters per second (cms) real, intent(in) :: routing_period ! seconds integer, intent(out):: dynamic_reservoir_type ! dynamic reservoir type sent to lake out files - real, intent(out) :: assimilated_value ! value assimilated from observation or forecast + real*8, intent(out) :: assimilated_value ! value assimilated from observation or forecast character(len=256), intent(out) :: assimilated_source_file ! source file of assimilated value end subroutine run_reservoir_interface diff --git a/src/Routing/module_channel_routing.F b/src/Routing/module_channel_routing.F index ff2d1b299..91b65cb24 100644 --- a/src/Routing/module_channel_routing.F +++ b/src/Routing/module_channel_routing.F @@ -1733,7 +1733,7 @@ subroutine drive_CHANNEL_RSL(did, UDMP_OPT,KT, IXRT,JXRT, & real*8, allocatable,dimension(:) :: tmpQLAKEO, tmpQLAKEI real*8, allocatable,dimension(:) :: tmpRESHT integer, allocatable, dimension(:) :: tmpFinalResType - real, allocatable,dimension(:) :: tmpAssimilatedValue + real*8, allocatable,dimension(:) :: tmpAssimilatedValue character(len=256), allocatable,dimension(:) :: tmpAssimilatedSourceFile #ifdef MPP_LAND From c6b8a46fc170312ad4166bfd2c949a78bdaad33f Mon Sep 17 00:00:00 2001 From: Ryan Cabell Date: Wed, 26 Jan 2022 16:06:18 -0700 Subject: [PATCH 5/5] Add missing mpp_land_bcast_real8() and read_rst_crt_nc8() subroutines --- src/MPP/mpp_land.F | 11 ++++++++ src/Routing/module_HYDRO_io.F | 52 ++++++++++++++++++++++++++++++----- 2 files changed, 56 insertions(+), 7 deletions(-) diff --git a/src/MPP/mpp_land.F b/src/MPP/mpp_land.F index cbcd3f61e..bb95ff676 100644 --- a/src/MPP/mpp_land.F +++ b/src/MPP/mpp_land.F @@ -1096,6 +1096,17 @@ subroutine mpp_land_bcast_real(size1,inout) return end subroutine mpp_land_bcast_real + subroutine mpp_land_bcast_real8(size1,inout) + integer size1 + ! real inout(size1) + real*8 , dimension(:) :: inout + integer ierr, len + call mpi_bcast(inout,size1,MPI_REAL8, & + IO_id,HYDRO_COMM_WORLD,ierr) + call mpp_land_sync() + return + end subroutine mpp_land_bcast_real8 + subroutine mpp_land_bcast_int2d(inout) integer length1, k,length2 integer inout(:,:) diff --git a/src/Routing/module_HYDRO_io.F b/src/Routing/module_HYDRO_io.F index a03b1f8f4..6be9b5639 100644 --- a/src/Routing/module_HYDRO_io.F +++ b/src/Routing/module_HYDRO_io.F @@ -6700,7 +6700,7 @@ subroutine read_rst_gwbucket_real(ncid,outV,numbasns,& character(len=*) :: vName integer i, j,k real, dimension(gnumbasns) :: buf - call read_rst_crt_nc(ncid,dble(buf),gnumbasns,vName) + call read_rst_crt_nc(ncid,buf,gnumbasns,vName) do k = 1, numbasns outV(k) = buf(basnsInd(k)) end do @@ -6875,9 +6875,9 @@ subroutine RESTART_IN_NC(inFile,did) endif if(rt_domain(did)%NLAKES .gt. 0) then - call read_rst_crt_nc(ncid,rt_domain(did)%RESHT,rt_domain(did)%NLAKES,"resht") - call read_rst_crt_nc(ncid,rt_domain(did)%QLAKEO,rt_domain(did)%NLAKES,"qlakeo") - call read_rst_crt_nc(ncid,rt_domain(did)%QLAKEI,rt_domain(did)%NLAKES,"qlakei") + call read_rst_crt_nc8(ncid,rt_domain(did)%RESHT,rt_domain(did)%NLAKES,"resht") + call read_rst_crt_nc8(ncid,rt_domain(did)%QLAKEO,rt_domain(did)%NLAKES,"qlakeo") + call read_rst_crt_nc8(ncid,rt_domain(did)%QLAKEI,rt_domain(did)%NLAKES,"qlakei") endif if( nlst(did)%channel_only .eq. 0 .and. & @@ -7152,7 +7152,7 @@ subroutine read_rst_crt_nc(ncid,var,n,varStr) implicit none integer :: ireg, ncid, varid, n, iret !real,dimension(n) :: var - real*8,dimension(n) :: var + real,dimension(n) :: var character(len=*) :: varStr if( n .le. 0) return @@ -7179,14 +7179,52 @@ subroutine read_rst_crt_nc(ncid,var,n,varStr) #ifdef MPP_LAND endif if(n .gt. 0) then - call mpp_land_bcast(var) - !call mpp_land_bcast_real(n,var) + call mpp_land_bcast_real(n,var) endif #endif return end subroutine read_rst_crt_nc + subroutine read_rst_crt_nc8(ncid,var,n,varStr) + implicit none + integer :: ireg, ncid, varid, n, iret + !real,dimension(n) :: var + real*8,dimension(n) :: var + character(len=*) :: varStr + + if( n .le. 0) return +#ifdef MPP_LAND + if(my_id .eq. IO_id) & +#endif + iret = nf90_inq_varid(ncid, trim(varStr), varid) +#ifdef MPP_LAND + call mpp_land_bcast_int1(iret) +#endif + if (iret /= 0) then +#ifdef HYDRO_D + print*, 'variable not found: name = "', trim(varStr)//'"' +#endif + return + endif +#ifdef HYDRO_D + print*, "read restart variable ", varStr +#endif +#ifdef MPP_LAND + if(my_id .eq. IO_id) then +#endif + iret = nf90_get_var(ncid, varid, var) +#ifdef MPP_LAND + endif + if(n .gt. 0) then + call mpp_land_bcast_real8(n,var) + + endif +#endif + return + end subroutine read_rst_crt_nc8 + + subroutine read_rst_crt_stream_nc(ncid,var_out,n,varStr,gnlinks,map_l2g) implicit none integer :: ncid, varid, n, iret, gnlinks