From fd2141473c714132972d1d8f381f6bb81168806d Mon Sep 17 00:00:00 2001 From: rniemeyer07 Date: Wed, 19 Jul 2017 13:32:42 -0700 Subject: [PATCH 01/27] inserted a comment at top of Begin.f90 to verify PR --- src/Begin.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Begin.f90 b/src/Begin.f90 index bf115d5..8c8c257 100755 --- a/src/Begin.f90 +++ b/src/Begin.f90 @@ -1,5 +1,7 @@ Subroutine BEGIN(param_file,spatial_file) ! +! Ryan make a comment to verify PR (7/19/2017) +! use Block_Energy use Block_Hydro use Block_Network From 2253a6abb8bf89e20d552ba52096cfd98bb4a1d5 Mon Sep 17 00:00:00 2001 From: rniemeyer07 Date: Tue, 25 Jul 2017 13:32:34 -0700 Subject: [PATCH 02/27] small changes to read_forcing.f90 to read in energy and flow files --- src/Read_Forcing.f90 | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index 46bdcb3..d9e96e1 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -8,7 +8,7 @@ SUBROUTINE Read_Forcing ! integer :: nc,ncell,nnd,no_flow,no_heat,nr,nrec_flow,nrec_heat real :: Q_avg,Q_dmmy - +character(40) :: line no_flow=0 no_heat=0 @@ -19,17 +19,27 @@ SUBROUTINE Read_Forcing ! nrec_flow=flow_cells*(ndays-1)+no_flow nrec_heat=heat_cells*(ndays-1)+no_heat + ! print *, 'no_heat', no_heat,'nrec_flow',nrec_flow, 'nrec_heat', nrec_heat ! - read(35,'(2i5,2f10.1,2f6.1,f7.1,f6.2)' & - ,rec=nrec_flow) nnd,ncell & - ,Q_out(no_heat),Q_dmmy,Q_diff(no_heat) & - ,depth(no_heat),width(no_heat),u(no_heat) +! read(35,'(i1,1x,i1,1x,f5.1,1x,f4.1,1x,f2.1,1x,f2.1,1x,f3.1,1x,f1.2)' & + read(35,*) & + ! , rec=nrec_flow) & + nnd,ncell & + ,Q_out(no_heat),Q_dmmy, Q_diff(no_heat) & + ,depth(no_heat) & + ,width(no_heat) & + ,u(no_heat) ! + +!print *, 'nnd', nnd, 'ncell',ncell,'Q_out',Q_out(no_heat),'Q_dmy',Q_dmmy & +!,'Qdiff',Q_diff(no_heat),'dep',depth(no_heat),'width',width(no_heat),'u',u(no_heat) if(u(no_heat).lt.0.01) u(no_heat)=0.01 if(ncell.ne.no_heat) write(*,*) 'Flow file error',ncell,no_heat ! - read(36,'(i5,2f6.1,2f7.4,f6.3,f7.1,f5.1)' & - ,rec=nrec_heat) ncell & +! read(36,'(i5,2f6.1,2f7.4,f6.3,f7.1,f5.1)' & + read(36, *) & + ! ,rec=nrec_heat) + ncell & ,dbt(no_heat),ea(no_heat) & ,Q_ns(no_heat),Q_na(no_heat),rho & ,press(no_heat),wind(no_heat) @@ -68,8 +78,9 @@ SUBROUTINE Read_Forcing ! Q_trib(nr)=Q_out(no_heat) nrec_heat=heat_cells*(ndays-1)+no_heat - read(36,'(i5,2f6.1,2f7.4,f6.3,f7.1,f5.1)' & - ,rec=nrec_heat) ncell & +! read(36,'(i5,2f6.1,2f7.4,f6.3,f7.1,f5.1)' & +! ,rec=nrec_heat) + read(36, *) ncell & ,dbt(no_heat),ea(no_heat) & ,Q_ns(no_heat),Q_na(no_heat),rho & ,press(no_heat),wind(no_heat) From f16833b888000a7c19e69a75fc9c8b3928fe7923 Mon Sep 17 00:00:00 2001 From: rniemeyer07 Date: Mon, 31 Jul 2017 07:53:27 -0700 Subject: [PATCH 03/27] removed ^M carriage returns --- src/Begin.f90 | 46 +++++----- src/Block_Energy.f90 | 2 +- src/Energy.f90 | 52 +++++------ src/Systmm.f90 | 64 +++++++------- src/Write.f90 | 2 +- src/rbm10_VIC.f90 | 200 +++++++++++++++++++++---------------------- src/tntrp.f90 | 82 +++++++++--------- 7 files changed, 224 insertions(+), 224 deletions(-) diff --git a/src/Begin.f90 b/src/Begin.f90 index 8c8c257..0da3662 100755 --- a/src/Begin.f90 +++ b/src/Begin.f90 @@ -14,7 +14,7 @@ Subroutine BEGIN(param_file,spatial_file) character (len=200):: param_file,source_file,spatial_file ! integer:: Julian - integer:: head_name,trib_cell + integer:: head_name,trib_cell integer:: jul_start,main_stem,nyear1,nyear2,nc,ncell,nseg integer:: ns_max_test,node,ncol,nrow,nr,cum_sgmnt ! @@ -88,7 +88,7 @@ Subroutine BEGIN(param_file,spatial_file) ! Read the number of cells in this reach, the headwater #, ! the number of the cell where it enters the next higher order stream, ! the headwater number of the next higher order stream it enters, and -! the river mile of the headwaters. +! the river mile of the headwaters. ! read(90,'(i5,11x,i4,10x,i5,15x,i5,15x,f10.0,i5)') no_cells(nr) & ,head_name,trib_cell,main_stem,rmile0 @@ -104,7 +104,7 @@ Subroutine BEGIN(param_file,spatial_file) if (trib_cell.gt.0) then no_tribs(trib_cell) = no_tribs(trib_cell)+1 trib(trib_cell,no_tribs(trib_cell)) = nr - end if + end if ! ! Reading Mohseni parameters for each headwaters (UW_JRY_2011/06/18) ! @@ -152,9 +152,9 @@ Subroutine BEGIN(param_file,spatial_file) ! Added variable ndelta (UW_JRY_2011/03/15) ! dx(ncell)=miles_to_ft*(rmile0-rmile1)/ndelta(ncell) - rmile0=rmile1 - nndlta=0 -200 continue + rmile0=rmile1 + nndlta=0 +200 continue nndlta=nndlta+1 nseg=nseg+1 segment_cell(nr,nseg)=ncell @@ -171,12 +171,12 @@ Subroutine BEGIN(param_file,spatial_file) ! Added variable ndelta (UW_JRY_2011/03/15) ! if(nndlta.lt.ndelta(ncell)) go to 200 - no_celm(nr)=nseg + no_celm(nr)=nseg segment_cell(nr,nseg)=ncell - x_dist(nr,nseg)=miles_to_ft*rmile1 + x_dist(nr,nseg)=miles_to_ft*rmile1 ! ! End of cell and segment loop -! +! end do ! ! If this is a reach that is tributary to another, set the confluence cell to the previous @@ -188,26 +188,26 @@ Subroutine BEGIN(param_file,spatial_file) conflnce(trib_cell,no_tribs(trib_cell)) = ncell end if -if(ns_max_test.lt.nseg) ns_max_test=nseg -! +if(ns_max_test.lt.nseg) ns_max_test=nseg +! ! End of reach loop -! +! end do if(ns_max_test.gt.ns_max) then write(*,*) 'RBM is terminating because' write(*,*) 'NS_MAX exceeded. Change NS_MAX in Block_Network to: ',ns_max_test stop -end if -! -nwpd=1 -xwpd=nwpd -dt_comp=86400./xwpd -! -! ****************************************************** -! Return to RMAIN -! ****************************************************** -! -900 continue +end if +! +nwpd=1 +xwpd=nwpd +dt_comp=86400./xwpd +! +! ****************************************************** +! Return to RMAIN +! ****************************************************** +! +900 continue ! ! end subroutine BEGIN diff --git a/src/Block_Energy.f90 b/src/Block_Energy.f90 index 3adcaa9..5f8b290 100755 --- a/src/Block_Energy.f90 +++ b/src/Block_Energy.f90 @@ -34,6 +34,6 @@ module Block_Energy real :: lvp,rb,rho real,parameter :: evap_coeff=1.5e-9 !Lake Hefner coefficient, 1/meters real,parameter :: pf=0.640,pi=3.14159 - real,parameter :: rfac=304.8 !rho/Cp kg/meter**3/Kilocalories/kg/Deg K + real,parameter :: rfac=304.8 !rho/Cp kg/meter**3/Kilocalories/kg/Deg K ! end module Block_Energy diff --git a/src/Energy.f90 b/src/Energy.f90 index fa75c23..3f9995a 100755 --- a/src/Energy.f90 +++ b/src/Energy.f90 @@ -1,37 +1,37 @@ SUBROUTINE Energy(T_surf,q_surf,ncell) - use Block_Energy + use Block_Energy implicit none integer::i,ncell,nd real::A,B,e0,q_surf,q_conv,q_evap,q_ws,td,T_surf real, dimension(2):: q_fit, T_fit ! - td=nd - T_fit(1)=T_surf-1.0 - T_fit(2)=T_surf+1.0 - do i=1,2 - e0=2.1718E8*EXP(-4157.0/(T_fit(i)+239.09)) - rb=pf*(dbt(ncell)-T_fit(i)) - lvp=597.0-0.57*T_fit(i) - q_evap=1000.*lvp*evap_coeff*wind(ncell) - if(q_evap.lt.0.0) q_evap=0.0 - q_conv=rb*q_evap - q_evap=q_evap*(e0-ea(ncell)) - q_ws=6.693E-2+1.471E-3*T_fit(i) - q_fit(i)=q_ns(ncell)+q_na(ncell)-q_ws-q_evap+q_conv - end do -! -! q=AT+B + td=nd + T_fit(1)=T_surf-1.0 + T_fit(2)=T_surf+1.0 + do i=1,2 + e0=2.1718E8*EXP(-4157.0/(T_fit(i)+239.09)) + rb=pf*(dbt(ncell)-T_fit(i)) + lvp=597.0-0.57*T_fit(i) + q_evap=1000.*lvp*evap_coeff*wind(ncell) + if(q_evap.lt.0.0) q_evap=0.0 + q_conv=rb*q_evap + q_evap=q_evap*(e0-ea(ncell)) + q_ws=6.693E-2+1.471E-3*T_fit(i) + q_fit(i)=q_ns(ncell)+q_na(ncell)-q_ws-q_evap+q_conv + end do +! +! q=AT+B ! ! Linear fit over the range of 2.0 deg C. ! These results can be used to estimate the "equilibrium" ! temperature and linear rate constant. -! - A=(q_fit(1)-q_fit(2))/(T_fit(1)-T_fit(2)) - q_surf=0.5*(q_fit(1)+q_fit(2)) - B=(q_surf/A)-(T_fit(1)+T_fit(2))/2. -! -! ****************************************************** -! Return to Subroutine RIVMOD -! ****************************************************** -! +! + A=(q_fit(1)-q_fit(2))/(T_fit(1)-T_fit(2)) + q_surf=0.5*(q_fit(1)+q_fit(2)) + B=(q_surf/A)-(T_fit(1)+T_fit(2))/2. +! +! ****************************************************** +! Return to Subroutine RIVMOD +! ****************************************************** +! END Subroutine Energy diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 1dcd026..f2aa551 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -7,7 +7,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) Implicit None ! ! -character (len=200):: temp_file +character (len=200):: temp_file character (len=200):: param_file ! integer :: ncell,nncell,ncell0,nc_head,no_flow,no_heat @@ -20,8 +20,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Indices for lagrangian interpolation ! integer :: njb,npndx,ntrp -integer, dimension(2):: ndltp=(/-1,-2/) -integer, dimension(2):: nterp=(/2,3/) +integer, dimension(2):: ndltp=(/-1,-2/) +integer, dimension(2):: nterp=(/2,3/) ! real :: dt_calc,dt_total,hpd,q_dot,q_surf,z @@ -36,7 +36,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) real,dimension(4):: ta,xa ! real,dimension(:),allocatable :: T_head,T_smth,T_trib - + logical:: DONE ! ! @@ -80,8 +80,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! ! open the output file ! - -open(20,file=TRIM(temp_file),status='unknown') + +open(20,file=TRIM(temp_file),status='unknown') ! ! ! Initialize dummy counters that facilitate updating simulated values @@ -109,13 +109,13 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Start the numbers of days-to-date counter ! ndays=ndays+1 -! +! ! Daily period loop starts -! +! DO ndd=1,nwpd - xdd = ndd + xdd = ndd time=year+(xd+(xdd-0.5)*hpd)/xd_year - + ! ! Read the hydrologic and meteorologic forcings ! @@ -137,7 +137,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! ! Variable Mohseni parameters (UW_JRY_2011/06/16) ! - T_head(nr)=mu(nr)+(alphaMu(nr) & + T_head(nr)=mu(nr)+(alphaMu(nr) & /(1.+exp(gmma(nr)*(beta(nr)-T_smth(nr))))) ! temp(nr,0,n1)=T_head(nr) @@ -156,10 +156,10 @@ SUBROUTINE SYSTMM(temp_file,param_file) call Particle_Track(nr,ns,nx_s,nx_head) ! ncell=segment_cell(nr,ns) -! -! Now do the third-order interpolation to -! establish the starting temperature values -! for each parcel +! +! Now do the third-order interpolation to +! establish the starting temperature values +! for each parcel ! nseg=nstrt_elm(ns) ! @@ -184,7 +184,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! do ntrp=nterp(npndx),1,-1 npart=nseg+ntrp+ndltp(npndx) - xa(ntrp)=x_dist(nr,npart) + xa(ntrp)=x_dist(nr,npart) ta(ntrp)=temp(nr,npart,n1) end do ! @@ -286,44 +286,44 @@ SUBROUTINE SYSTMM(temp_file,param_file) if (T_0.lt.0.5) T_0=0.5 temp(nr,ns,n2)=T_0 T_trib(nr)=T_0 -! +! ! Write all temperature output UW_JRY_11/08/2013 ! The temperature is output at the beginning of the ! reach. It is, of course, possible to get output at ! other points by some additional code that keys on the -! value of ndelta (now a vector)(UW_JRY_11/08/2013) -! +! value of ndelta (now a vector)(UW_JRY_11/08/2013) +! call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell),Q_inflow,Q_outflow) -! -! End of computational element loop ! - end do -! End of reach loop -! +! End of computational element loop +! + end do +! End of reach loop +! end do ntmp=n1 n1=n2 n2=ntmp -! +! ! End of weather period loop (NDD=1,NWPD) -! +! 4650 format(16x,12(6x,f6.0,6x)) 4700 format(f10.4,f6.0,15(f6.1,f8.3)) 4750 format(f10.4,10(i4,f8.0)) - end do -! + end do +! ! End of main loop (ND=1,365/366) -! +! end do -! +! ! End of year loop ! end do ! ! ! ****************************************************** -! return to rmain +! return to rmain ! ****************************************************** -! +! 950 return end SUBROUTINE SYSTMM diff --git a/src/Write.f90 b/src/Write.f90 index a824864..b3a8b41 100755 --- a/src/Write.f90 +++ b/src/Write.f90 @@ -9,6 +9,6 @@ SUBROUTINE WRITE(time,nd,nr,ncell,ns,T_0,T_head,dbt,Q_inflow,Q_outflow) real :: T_head real :: dbt ! -write (20,'(f12.4,4i6,3f8.2,4f25.1,f8.1,f8.4)') & +write (20,'(f12.4,4i6,3f8.2,4f25.1,f8.1,f8.4)') & time,nd,nr,ncell,ns,T_0,T_head,dbt,Q_inflow,Q_outflow end SUBROUTINE WRITE diff --git a/src/rbm10_VIC.f90 b/src/rbm10_VIC.f90 index 4d6c68e..5682adc 100755 --- a/src/rbm10_VIC.f90 +++ b/src/rbm10_VIC.f90 @@ -1,10 +1,10 @@ -! - PROGRAM RBM10_VIC -! -! Dynamic river basin model for simulating water quality in -! branching river systems with freely-flowing river segments. -! -! This version uses Reverse Particle Tracking in the Lagrangian +! + PROGRAM RBM10_VIC +! +! Dynamic river basin model for simulating water quality in +! branching river systems with freely-flowing river segments. +! +! This version uses Reverse Particle Tracking in the Lagrangian ! mode and Lagrangian interpolation in the Eulerian mode. ! ! This version of the model software has the following limited features: @@ -17,113 +17,113 @@ PROGRAM RBM10_VIC ! 2. Output is written to Unit 20 (specified in the command line for executing ! RBM) and includes the simulated temperatures for the first computational cell ! in each segment in the same order as that of the *.network file and at the -! same time step. -! -! Topology and routing is set up to be consistent with output -! from the Variable Infiltration Capacity (VIC) model developed by the +! same time step. +! +! Topology and routing is set up to be consistent with output +! from the Variable Infiltration Capacity (VIC) model developed by the ! Land Surface Hydrology Group at the University of Washington. ! Model details are described in ! ! Yearsley, J. (2012), A grid-based approach for simulating stream temperature, -! Water Resour. Res., 48, W03506, doi:10.1029/2011WR011515 -! -! For additional information contact: -! -! John Yearsley -! Land Surface Hydrology Group -! Dept. of Civil and Environmental Engineering -! Box 352700 -! University of Washington -! Seattle, Washington +! Water Resour. Res., 48, W03506, doi:10.1029/2011WR011515 +! +! For additional information contact: +! +! John Yearsley +! Land Surface Hydrology Group +! Dept. of Civil and Environmental Engineering +! Box 352700 +! University of Washington +! Seattle, Washington ! 98195-2700 -! +! !use BGIN !use SYSTM ! -implicit none -! -! -character (len=200 ):: inPrefix -character (len=200 ):: outPrefix +implicit none +! +! +character (len=200 ):: inPrefix +character (len=200 ):: outPrefix character (len=200 ):: flow_file character (len=200 ):: heat_file -character (len=200 ):: net_file -character (len=200 ):: param_file -character (len=200 ):: temp_file +character (len=200 ):: net_file +character (len=200 ):: param_file +character (len=200 ):: temp_file character (len=200 ):: spatial_file -character (len=8) :: start_data,end_data -integer iargc -integer numarg - -! -! Command line input -! -numarg = iargc ( ) -if (numarg .lt. 2) then - write (*,*) 'Too few arguments were given' - write (*,*) ' ' - write (*,*) 'First: Location and prefix of input files' - write (*,*) ' (networkfile and parameterfile)' - write (*,*) 'Second: Location and prefix of output files' - write (*,*) ' ' - write (*,*) 'eg: $ ./input/Salmon_0.50 ./output/Salmon_Test' - write (*,*) ' ' - stop -end if -call getarg ( 1, inPrefix ) -call getarg ( 2, outPrefix ) -! +character (len=8) :: start_data,end_data +integer iargc +integer numarg + +! +! Command line input +! +numarg = iargc ( ) +if (numarg .lt. 2) then + write (*,*) 'Too few arguments were given' + write (*,*) ' ' + write (*,*) 'First: Location and prefix of input files' + write (*,*) ' (networkfile and parameterfile)' + write (*,*) 'Second: Location and prefix of output files' + write (*,*) ' ' + write (*,*) 'eg: $ ./input/Salmon_0.50 ./output/Salmon_Test' + write (*,*) ' ' + stop +end if +call getarg ( 1, inPrefix ) +call getarg ( 2, outPrefix ) +! ! Identify input/output files -! -net_file = TRIM(inPrefix)//'_Network' -param_file = TRIM(inPrefix)//'_Parameters' -spatial_file = TRIM(outPrefix)//'.Spat' -temp_file = TRIM(outPrefix)//'.Temp' -! -write(*,*) 'Spatial file: ',spatial_file -write(*,*) 'Network file : ',net_file -write(*,*) 'Parameter file : ',param_file! -write(*,*) 'Temperature file: ',temp_file -! -OPEN(UNIT=90,FILE=TRIM(net_file),STATUS='OLD') -! -! Read header information from control file -! -read(90,*) -read(90,'(A)') flow_file -! -! Open file with hydrologic data -! -open(unit=35,FILE=TRIM(flow_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=60,STATUS='old') -! -! -read(90,'(A)') heat_file -! -! Open file with meteorologic data -! -open(unit=36,FILE=TRIM(heat_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=50,STATUS='old') -! -! Call systems programs to get started -! -! SUBROUTINE BEGIN reads control file, sets up topology and -! important properties of reaches -! +! +net_file = TRIM(inPrefix)//'_Network' +param_file = TRIM(inPrefix)//'_Parameters' +spatial_file = TRIM(outPrefix)//'.Spat' +temp_file = TRIM(outPrefix)//'.Temp' +! +write(*,*) 'Spatial file: ',spatial_file +write(*,*) 'Network file : ',net_file +write(*,*) 'Parameter file : ',param_file! +write(*,*) 'Temperature file: ',temp_file +! +OPEN(UNIT=90,FILE=TRIM(net_file),STATUS='OLD') +! +! Read header information from control file +! +read(90,*) +read(90,'(A)') flow_file +! +! Open file with hydrologic data +! +open(unit=35,FILE=TRIM(flow_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=60,STATUS='old') +! +! +read(90,'(A)') heat_file +! +! Open file with meteorologic data +! +open(unit=36,FILE=TRIM(heat_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=50,STATUS='old') +! +! Call systems programs to get started +! +! SUBROUTINE BEGIN reads control file, sets up topology and +! important properties of reaches +! write(*,*) 'Calling BEGIN' ! ! SUBROUTINE BEGIN reads in river system information from the NETWORK file -! -CALL BEGIN(param_file, spatial_file) -! -! SUBROUTINE SYSTMM performs the simulations -! -CALL SYSTMM(temp_file,param_file) ! (WUR_WF_MvV_2011/01/05) -! -! Close files after simulation is complete -! +! +CALL BEGIN(param_file, spatial_file) +! +! SUBROUTINE SYSTMM performs the simulations +! +CALL SYSTMM(temp_file,param_file) ! (WUR_WF_MvV_2011/01/05) +! +! Close files after simulation is complete +! write(*,*) ' Closing files after simulation' - + CLOSE(35) -CLOSE(36) -CLOSE(90) -STOP -END PROGRAM RBM10_VIC +CLOSE(36) +CLOSE(90) +STOP +END PROGRAM RBM10_VIC diff --git a/src/tntrp.f90 b/src/tntrp.f90 index 26755c7..f4a1489 100755 --- a/src/tntrp.f90 +++ b/src/tntrp.f90 @@ -1,47 +1,47 @@ -! -! Third-order polynomial interpolation using Lagrange -! polynomials. FUNCTION is SUBROUTINE POLINT from -! Numerial Recipes -! -real FUNCTION tntrp(XA,YA,X,N) +! +! Third-order polynomial interpolation using Lagrange +! polynomials. FUNCTION is SUBROUTINE POLINT from +! Numerial Recipes +! +real FUNCTION tntrp(XA,YA,X,N) ! implicit NONE ! integer:: i,m,n,ns -real :: den,dif,dift,dy,hO,hp,w,x,y - real, DIMENSION(4):: XA,YA,C,D - NS=1 - DIF=ABS(X-XA(1)) - DO 11 I=1,N - DIFT=ABS(X-XA(I)) - IF (DIFT.LT.DIF) THEN - NS=I - DIF=DIFT - ENDIF - C(I)=YA(I) - D(I)=YA(I) -11 CONTINUE - Y=YA(NS) - NS=NS-1 - DO 13 M=1,N-1 - DO 12 I=1,N-M - HO=XA(I)-X - HP=XA(I+M)-X - W=C(I+1)-D(I) - DEN=HO-HP - IF(DEN.EQ.0.) DEN=0.001 - DEN=W/DEN - D(I)=HP*DEN - C(I)=HO*DEN -12 CONTINUE - IF (2*NS.LT.N-M)THEN - DY=C(NS+1) - ELSE - DY=D(NS) - NS=NS-1 - ENDIF - Y=Y+DY -13 CONTINUE - tntrp=y +real :: den,dif,dift,dy,hO,hp,w,x,y + real, DIMENSION(4):: XA,YA,C,D + NS=1 + DIF=ABS(X-XA(1)) + DO 11 I=1,N + DIFT=ABS(X-XA(I)) + IF (DIFT.LT.DIF) THEN + NS=I + DIF=DIFT + ENDIF + C(I)=YA(I) + D(I)=YA(I) +11 CONTINUE + Y=YA(NS) + NS=NS-1 + DO 13 M=1,N-1 + DO 12 I=1,N-M + HO=XA(I)-X + HP=XA(I+M)-X + W=C(I+1)-D(I) + DEN=HO-HP + IF(DEN.EQ.0.) DEN=0.001 + DEN=W/DEN + D(I)=HP*DEN + C(I)=HO*DEN +12 CONTINUE + IF (2*NS.LT.N-M)THEN + DY=C(NS+1) + ELSE + DY=D(NS) + NS=NS-1 + ENDIF + Y=Y+DY +13 CONTINUE + tntrp=y RETURN end From f5f2430fb1fa7ef6af0862cfd5f01ddc44fed339 Mon Sep 17 00:00:00 2001 From: rniemeyer07 Date: Mon, 31 Jul 2017 10:14:00 -0700 Subject: [PATCH 04/27] updated rbm10, added Block_Reservoir --- src/Block_Hydro.f90 | 3 ++ src/Block_Reservoir.f90 | 66 +++++++++++++++++++++++++++++++++++++++++ src/Read_Forcing.f90 | 1 + src/Systmm.f90 | 2 +- src/rbm10_VIC.f90 | 33 +++++++++++++++++++-- 5 files changed, 101 insertions(+), 4 deletions(-) create mode 100644 src/Block_Reservoir.f90 diff --git a/src/Block_Hydro.f90 b/src/Block_Hydro.f90 index c8b980b..7e96580 100755 --- a/src/Block_Hydro.f90 +++ b/src/Block_Hydro.f90 @@ -19,5 +19,8 @@ module Block_Hydro real, dimension(:,:), allocatable :: temp_nps,thermal real, dimension(:,:), allocatable :: x_dist real, dimension(:,:,:), allocatable :: temp + real, dimension(:), allocatable :: Q_trib_tot + real, dimension(:), allocatable :: T_trib, T_head + real, dimension(:), allocatable :: flow_source, source_num_cell end module Block_Hydro diff --git a/src/Block_Reservoir.f90 b/src/Block_Reservoir.f90 new file mode 100644 index 0000000..eddae67 --- /dev/null +++ b/src/Block_Reservoir.f90 @@ -0,0 +1,66 @@ +module Block_Reservoir + + implicit none + + ! -------------------------------------------------------------------- + ! reservoir to RBM variables + !---------------------------------------------------------------------- + + !------------------------read in reservoir info-------------------------- + real, dimension (:), allocatable:: dam_lat, dam_lon, res_grid_lat, res_grid_lon + real, dimension (:), allocatable:: res_depth_feet, res_width_feet, res_length_feet + integer, dimension (:), allocatable:: dam_number, start_operating_year + real, dimension (:), allocatable:: res_top_vol, res_bot_vol, res_max_flow, res_min_flow + integer, dimension (:), allocatable:: res_start_node, res_end_node + ! integer, dimension (:,:), allocatable:: nodes_x !for each reach, what are all the nodes + real , dimension(:), allocatable :: rmile_node + integer, dimension(:,:), allocatable :: res_num, nodes_x, x_dist_res + logical, dimension(:,:), allocatable :: res_pres, res_upstream + logical, dimension(:), allocatable :: flag_turnover + logical :: reservoir, res_upstreamx ! the first is TRUE or FALSE in fifth line of _Network file whether reserovirs are present + integer, dimension (:), allocatable :: xres, resx + integer :: xres2, nres, nm_start, ncell0res + real, dimension(:), allocatable :: dx_res, dt_res, reservoir_storage,reservoir_storage_prev + real :: advec_tot,q_surf_tot + real, dimension (:), allocatable::diffusion_tot,advec_hyp_tot,advec_epi_tot, qsurf_tot + ! real :: Q_res_in_1, T_res_in_1 + + !---------------------------------------------------------------------- + ! variables from reservoir model + !---------------------------------------------------------------------- + + ! ----------------------- constants -------------------------- + real, parameter :: prcnt_flow_epil = 1, prcnt_flow_hypo = 0, gravity = 9.8 + real, parameter :: density = 1000, heat_c = 4180, ftsec_to_msec = 0.0283168 + real, parameter :: J_to_kcal = 0.00023885, kPa_to_mb = 10, ft_to_m = 0.3048 + real, parameter :: acrefeet_to_m3 = 1233.48 + real, parameter :: heat_c_kcal = 1 ! heat capacity in kcal/kg*C + real, dimension (:), allocatable :: K_z !diffusion coefficient (m^2/sec) + + ! -------------------- temperature and meterological variables ------ + real, dimension (:), allocatable :: temp_change_ep, temp_change_hyp, temp_out, temp_out_i ! energy + real, dimension (:), allocatable :: T_epil,T_hypo, stream_T_in + real, dimension (:), allocatable :: density_epil, density_hypo, density_in + + ! -------------------- reservoinr information -------------- + real, dimension (:), allocatable :: surface_area, depth_total, depth_e, depth_h + real, dimension (:), allocatable :: delta_vol_e_x, delta_vol_h_x + real, dimension (:), allocatable :: delta_vol_e_T_x, delta_vol_h_T_x, dV_dt_epi,dV_dt_hyp + real, parameter :: depth_e_frac=0.4, depth_h_frac=0.6 + real, dimension (:), allocatable :: Q_tot, Q_pen, Q_spill + real, dimension (:), allocatable :: depth_e_inital, volume_e_initial, depth_h_inital, volume_h_initial + real, dimension (:), allocatable :: volume_e_x,volume_h_x, T_res, T_res_in, T_trib_tot, Q_res_in + logical, dimension (:), allocatable :: res_run, res_start, trib_res ! logical to only get start of reservoir and model entire reservoir once each loop + real :: outflow_x, volume_tot, T_res_in_x, Q_trib_tot_x, T_trib_in_x + ! -------------------- energy terms ----------- + real, dimension (:), allocatable :: area + ! real :: flow_in_hyp_x, flow_in_epi_x, flow_out_epi_x, flow_out_hyp_x + ! real :: flow_epi_hyp_x + real :: epix, hypox, dif_epi_x, dif_hyp_x, energy_x + real :: advec_in_epix, advec_out_epix, advec_in_hypx, advec_out_hypx + real :: advec_epi_hyp,tntrp_x + + ! --------------- path and directories of input and output files ------- + character(len=300 ) :: reservoir_file + character(len=300 ) :: releases_file +end module Block_Reservoir diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index d9e96e1..d5fec54 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -95,6 +95,7 @@ SUBROUTINE Read_Forcing depth(no_heat)=depth(no_heat-1) width(no_heat)=width(no_heat-1) dt(no_heat)=dx(ncell)/u(no_heat) +if(nnd .gt. 400) stop ! ends run after certain day end do ! ! Call the water balance subroutine diff --git a/src/Systmm.f90 b/src/Systmm.f90 index f2aa551..8be894e 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -35,7 +35,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) real :: dt_ttotal real,dimension(4):: ta,xa ! -real,dimension(:),allocatable :: T_head,T_smth,T_trib +real,dimension(:),allocatable :: T_smth logical:: DONE ! diff --git a/src/rbm10_VIC.f90 b/src/rbm10_VIC.f90 index 5682adc..839b73e 100755 --- a/src/rbm10_VIC.f90 +++ b/src/rbm10_VIC.f90 @@ -51,6 +51,9 @@ PROGRAM RBM10_VIC character (len=200 ):: param_file character (len=200 ):: temp_file character (len=200 ):: spatial_file +character (len=200 ):: reservoir_output_file +character (len=200 ):: reservoir_info_file +character (len=200 ):: reservoir_storage_file character (len=8) :: start_data,end_data integer iargc integer numarg @@ -79,11 +82,13 @@ PROGRAM RBM10_VIC param_file = TRIM(inPrefix)//'_Parameters' spatial_file = TRIM(outPrefix)//'.Spat' temp_file = TRIM(outPrefix)//'.Temp' +reservoir_output_file = TRIM(outPrefix)//'.Reservoir' ! write(*,*) 'Spatial file: ',spatial_file write(*,*) 'Network file : ',net_file -write(*,*) 'Parameter file : ',param_file! +write(*,*) 'Parameter file : ',param_file write(*,*) 'Temperature file: ',temp_file +write(*,*) 'Reservoir Output file: ',reservoir_output_file ! OPEN(UNIT=90,FILE=TRIM(net_file),STATUS='OLD') ! @@ -94,14 +99,34 @@ PROGRAM RBM10_VIC ! ! Open file with hydrologic data ! -open(unit=35,FILE=TRIM(flow_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=60,STATUS='old') +!open(unit=35,FILE=TRIM(flow_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=60,STATUS='old') +open(unit=35,FILE=TRIM(flow_file) ,FORM='FORMATTED',ACCESS='SEQUENTIAL' ,STATUS='old') ! ! read(90,'(A)') heat_file ! ! Open file with meteorologic data ! -open(unit=36,FILE=TRIM(heat_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=50,STATUS='old') +!open(unit=36,FILE=TRIM(heat_file) ,FORM='FORMATTED',ACCESS='DIRECT' ,RECL=50,STATUS='old') +open(unit=36,FILE=TRIM(heat_file) ,FORM='FORMATTED',ACCESS='SEQUENTIAL' ,STATUS='old') +! +! +! +read(90,'(A)') reservoir_info_file +! +! open reservior information file +! +open(unit=37,FILE=TRIM(reservoir_info_file),ACCESS='SEQUENTIAL',FORM='FORMATTED', STATUS='old') +! +! +! +read(90,'(A)') reservoir_storage_file +! +! open reservior storage file +! +open(unit=38,FILE=TRIM(reservoir_storage_file),ACCESS='SEQUENTIAL',FORM='FORMATTED',STATUS='old') + + ! ! Call systems programs to get started ! @@ -124,6 +149,8 @@ PROGRAM RBM10_VIC CLOSE(35) CLOSE(36) +CLOSE(37) +CLOSE(38) CLOSE(90) STOP END PROGRAM RBM10_VIC From 3352fb08e3af3920a9e2d53d5e90cb524123aa5b Mon Sep 17 00:00:00 2001 From: rniemeyer07 Date: Tue, 1 Aug 2017 07:40:39 -0700 Subject: [PATCH 05/27] added Block_Reservoir, added variables to read in reservoir data, added if loop in Systmm --- src/Begin.f90 | 91 ++++++++++++++++++-- src/Block_Hydro.f90 | 2 + src/Block_Network.f90 | 7 ++ src/Block_Reservoir.f90 | 47 ++++++----- src/Read_Forcing.f90 | 6 +- src/Systmm.f90 | 181 +++++++++++++++++++++++++++++++--------- src/makefile | 10 ++- 7 files changed, 274 insertions(+), 70 deletions(-) mode change 100644 => 100755 src/Block_Reservoir.f90 diff --git a/src/Begin.f90 b/src/Begin.f90 index 0da3662..3982267 100755 --- a/src/Begin.f90 +++ b/src/Begin.f90 @@ -1,10 +1,10 @@ Subroutine BEGIN(param_file,spatial_file) ! -! Ryan make a comment to verify PR (7/19/2017) ! use Block_Energy use Block_Hydro use Block_Network +use Block_Reservoir ! implicit none ! @@ -12,13 +12,15 @@ Subroutine BEGIN(param_file,spatial_file) character (len=8) :: lat character (len=10):: long character (len=200):: param_file,source_file,spatial_file + character :: dam_name ! integer:: Julian integer:: head_name,trib_cell integer:: jul_start,main_stem,nyear1,nyear2,nc,ncell,nseg integer:: ns_max_test,node,ncol,nrow,nr,cum_sgmnt + integer:: nreservoir, node_res ! - logical:: first_cell,source + logical:: first_cell, res_presx ! real :: nndlta real :: rmile0,rmile1,xwpd @@ -43,7 +45,7 @@ Subroutine BEGIN(param_file,spatial_file) ! jul_start = Julian(start_year,start_month,start_day) ! -read(90,*) nreach,flow_cells,heat_cells,source +read(90,*) nreach, flow_cells, heat_cells, source, nsource, reservoir, nres ! ! Allocate dynamic arrays ! @@ -68,6 +70,50 @@ Subroutine BEGIN(param_file,spatial_file) allocate(head_cell(nreach)) allocate(segment_cell(nreach,ns_max)) allocate(x_dist(nreach,0:ns_max)) +! Allocate reservoir info + allocate(dam_lat(nres)) + allocate(dam_lon(nres)) + allocate(res_grid_lat(nres)) + allocate(res_grid_lon(nres)) + allocate(res_depth_feet(nres)) + allocate(res_width_feet(nres)) + allocate(res_length_feet(nres)) + allocate(dam_number(nres)) + allocate(start_operating_year(nres)) + allocate(res_top_vol(nres)) + allocate(res_bot_vol(nres)) + allocate(res_max_flow(nres)) + allocate(res_min_flow(nres)) + allocate(res_start_node(nres)) + allocate(res_end_node(nres)) + allocate(nodes_x(nreach,0:ns_max)) + allocate(res_num(nreach,0:ns_max)) + allocate(res_pres(nreach,0:ns_max)) + allocate(res_upstream(nreach,0:ns_max)) + allocate(rmile_node(0:ns_max)) + allocate(xres(0:ns_max)) + allocate(dx_res(0:heat_cells)) + + +! +!--------------- read in reservoir info (if reservoirs present) -------------- +! +if(reservoir) then + read(37,*) + do nreservoir = 1,nres + read(37,*) dam_number(nreservoir) & + , res_grid_lat(nreservoir), res_grid_lon(nreservoir) & + , res_top_vol(nreservoir) & + , res_bot_vol(nreservoir),res_depth_feet(nreservoir) & + , res_width_feet(nreservoir), res_length_feet(nreservoir) & + ,res_start_node(nreservoir), res_end_node(nreservoir) & + , res_min_flow(nreservoir) + + print *, 'nreservoir', nreservoir, 'dam_num', dam_number(nreservoir) & + , 'end', res_end_node(nreservoir) , 'start', res_start_node(nreservoir) + end do +end if + ! ! Start reading the reach date and initialize the reach index, NR ! and the cell index, NCELL @@ -137,8 +183,27 @@ Subroutine BEGIN(param_file,spatial_file) ! Variable ndelta read in here. At present, number of elements ! is entered manually into the network file (UW_JRY_2011/03/15) ! - read(90,'(5x,i5,5x,i5,8x,i5,6x,a8,6x,a10,7x,f10.0,f5.0)') & + + if(reservoir) then + read(90,'(5x,i5,5x,i5,8x,i5,6x,a8,6x,a10,7x,f10.0,f5.0,i6)') & !nodes for each reach + node,nrow,ncol,lat,long,rmile1,ndelta(ncell),node_res + if(node_res .gt.0) then + res_presx = .true. + else + res_presx = .false. + end if + + ! X,Y matrix with nodes in each reach + nodes_x(nr,ncell) = node ! matrix of nodes for each cell + res_num(nr,ncell) = node_res ! matrix of reservoir index values + res_pres(nr,ncell) = res_presx ! matrix for if reservoir is present for each cell + rmile_node(ncell) = rmile1 + + else ! IF no reservoir information in _Network file: + read(90,'(5x,i5,5x,i5,8x,i5,6x,a8,6x,a10,7x,f10.0,f5.0)') & node,nrow,ncol,lat,long,rmile1,ndelta(ncell) + end if + ! ! Set the number of segments of the default, if not specified ! @@ -158,8 +223,24 @@ Subroutine BEGIN(param_file,spatial_file) nndlta=nndlta+1 nseg=nseg+1 segment_cell(nr,nseg)=ncell - write(*,*) 'nndlta -- ',nr,nndlta,nseg,ncell,segment_cell(nr,nseg) x_dist(nr,nseg)=x_dist(nr,nseg-1)-dx(ncell) + + +! +! calcualte distance to next upstream reservoir +! + if(reservoir) then + if(res_presx) then !if cell in reservoir + res_upstream(nr,ncell) = .false. + + else if(any(res_pres(nr,:)) ) then ! any reservoirs upstream of this cell in this reach + xres(1:size(res_end_node)) = ncell - res_end_node + res_upstream(nr,ncell) = .true. + + end if + + end if + ! ! Write Segment List for mapping to temperature output (UW_JRY_2008/11/19) ! diff --git a/src/Block_Hydro.f90 b/src/Block_Hydro.f90 index 7e96580..586d3a5 100755 --- a/src/Block_Hydro.f90 +++ b/src/Block_Hydro.f90 @@ -23,4 +23,6 @@ module Block_Hydro real, dimension(:), allocatable :: T_trib, T_head real, dimension(:), allocatable :: flow_source, source_num_cell + integer :: nsource + end module Block_Hydro diff --git a/src/Block_Network.f90 b/src/Block_Network.f90 index 9c7b9ea..f96ef53 100755 --- a/src/Block_Network.f90 +++ b/src/Block_Network.f90 @@ -21,4 +21,11 @@ Module Block_Network real :: delta_n,n_default=2 real :: dt_comp real, dimension(:), allocatable :: ndelta +! +! Logical variables +! + logical :: source ! presence of thermal plant + + + end module Block_Network diff --git a/src/Block_Reservoir.f90 b/src/Block_Reservoir.f90 old mode 100644 new mode 100755 index eddae67..8f4c91b --- a/src/Block_Reservoir.f90 +++ b/src/Block_Reservoir.f90 @@ -1,28 +1,32 @@ -module Block_Reservoir +Module Block_Reservoir + + implicit none - + ! -------------------------------------------------------------------- ! reservoir to RBM variables !---------------------------------------------------------------------- - !------------------------read in reservoir info-------------------------- - real, dimension (:), allocatable:: dam_lat, dam_lon, res_grid_lat, res_grid_lon - real, dimension (:), allocatable:: res_depth_feet, res_width_feet, res_length_feet + !------------------------read in reservoir + !info-------------------------- + real, dimension (:), allocatable:: dam_lat, dam_lon, res_grid_lat,res_grid_lon + real, dimension (:), allocatable:: res_depth_feet, res_width_feet,res_length_feet integer, dimension (:), allocatable:: dam_number, start_operating_year - real, dimension (:), allocatable:: res_top_vol, res_bot_vol, res_max_flow, res_min_flow + real, dimension (:), allocatable:: res_top_vol, res_bot_vol,res_max_flow, res_min_flow integer, dimension (:), allocatable:: res_start_node, res_end_node - ! integer, dimension (:,:), allocatable:: nodes_x !for each reach, what are all the nodes + ! integer, dimension (:,:), allocatable:: nodes_x !for each reach, + ! what are all the nodes real , dimension(:), allocatable :: rmile_node integer, dimension(:,:), allocatable :: res_num, nodes_x, x_dist_res logical, dimension(:,:), allocatable :: res_pres, res_upstream - logical, dimension(:), allocatable :: flag_turnover + logical, dimension(:), allocatable :: flag_turnover logical :: reservoir, res_upstreamx ! the first is TRUE or FALSE in fifth line of _Network file whether reserovirs are present integer, dimension (:), allocatable :: xres, resx integer :: xres2, nres, nm_start, ncell0res - real, dimension(:), allocatable :: dx_res, dt_res, reservoir_storage,reservoir_storage_prev + real, dimension(:), allocatable :: dx_res, dt_res,reservoir_storage,reservoir_storage_prev real :: advec_tot,q_surf_tot - real, dimension (:), allocatable::diffusion_tot,advec_hyp_tot,advec_epi_tot, qsurf_tot + real, dimension (:),allocatable::diffusion_tot,advec_hyp_tot,advec_epi_tot, qsurf_tot ! real :: Q_res_in_1, T_res_in_1 !---------------------------------------------------------------------- @@ -31,26 +35,26 @@ module Block_Reservoir ! ----------------------- constants -------------------------- real, parameter :: prcnt_flow_epil = 1, prcnt_flow_hypo = 0, gravity = 9.8 - real, parameter :: density = 1000, heat_c = 4180, ftsec_to_msec = 0.0283168 - real, parameter :: J_to_kcal = 0.00023885, kPa_to_mb = 10, ft_to_m = 0.3048 - real, parameter :: acrefeet_to_m3 = 1233.48 + real, parameter :: density = 1000, heat_c = 4180, ftsec_to_msec =0.0283168 + real, parameter :: J_to_kcal = 0.00023885, kPa_to_mb = 10, ft_to_m =0.3048 + real, parameter :: acrefeet_to_m3 = 1233.48 real, parameter :: heat_c_kcal = 1 ! heat capacity in kcal/kg*C real, dimension (:), allocatable :: K_z !diffusion coefficient (m^2/sec) ! -------------------- temperature and meterological variables ------ - real, dimension (:), allocatable :: temp_change_ep, temp_change_hyp, temp_out, temp_out_i ! energy + real, dimension (:), allocatable :: temp_change_ep, temp_change_hyp,temp_out, temp_out_i ! energy real, dimension (:), allocatable :: T_epil,T_hypo, stream_T_in - real, dimension (:), allocatable :: density_epil, density_hypo, density_in + real, dimension (:), allocatable :: density_epil, density_hypo,density_in ! -------------------- reservoinr information -------------- - real, dimension (:), allocatable :: surface_area, depth_total, depth_e, depth_h + real, dimension (:), allocatable :: surface_area, depth_total, depth_e,depth_h real, dimension (:), allocatable :: delta_vol_e_x, delta_vol_h_x - real, dimension (:), allocatable :: delta_vol_e_T_x, delta_vol_h_T_x, dV_dt_epi,dV_dt_hyp + real, dimension (:), allocatable :: delta_vol_e_T_x, delta_vol_h_T_x,dV_dt_epi,dV_dt_hyp real, parameter :: depth_e_frac=0.4, depth_h_frac=0.6 real, dimension (:), allocatable :: Q_tot, Q_pen, Q_spill - real, dimension (:), allocatable :: depth_e_inital, volume_e_initial, depth_h_inital, volume_h_initial - real, dimension (:), allocatable :: volume_e_x,volume_h_x, T_res, T_res_in, T_trib_tot, Q_res_in - logical, dimension (:), allocatable :: res_run, res_start, trib_res ! logical to only get start of reservoir and model entire reservoir once each loop + real, dimension (:), allocatable :: depth_e_inital, volume_e_initial,depth_h_inital, volume_h_initial + real, dimension (:), allocatable :: volume_e_x,volume_h_x, T_res,T_res_in, T_trib_tot, Q_res_in + logical, dimension (:), allocatable :: res_run, res_start, trib_res !logical to only get start of reservoir and model entire reservoir once each loop real :: outflow_x, volume_tot, T_res_in_x, Q_trib_tot_x, T_trib_in_x ! -------------------- energy terms ----------- real, dimension (:), allocatable :: area @@ -58,9 +62,10 @@ module Block_Reservoir ! real :: flow_epi_hyp_x real :: epix, hypox, dif_epi_x, dif_hyp_x, energy_x real :: advec_in_epix, advec_out_epix, advec_in_hypx, advec_out_hypx - real :: advec_epi_hyp,tntrp_x + real :: advec_epi_hyp,tntrp_x ! --------------- path and directories of input and output files ------- character(len=300 ) :: reservoir_file character(len=300 ) :: releases_file + end module Block_Reservoir diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index d5fec54..4dd9ab4 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -77,7 +77,7 @@ SUBROUTINE Read_Forcing ! Tributary flow from this reach equals Q_out for this cell ! Q_trib(nr)=Q_out(no_heat) - nrec_heat=heat_cells*(ndays-1)+no_heat + nrec_heat=heat_cells*(ndays-1)+no_heat ! read(36,'(i5,2f6.1,2f7.4,f6.3,f7.1,f5.1)' & ! ,rec=nrec_heat) read(36, *) ncell & @@ -91,11 +91,11 @@ SUBROUTINE Read_Forcing ! Q_in(no_heat)=Q_out(no_heat-1) - u(no_heat)=u(no_heat-1) + u(no_heat)=u(no_heat-1) depth(no_heat)=depth(no_heat-1) width(no_heat)=width(no_heat-1) dt(no_heat)=dx(ncell)/u(no_heat) -if(nnd .gt. 400) stop ! ends run after certain day +!if(nnd .gt. 400) stop ! ends run after certain day end do ! ! Call the water balance subroutine diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 8be894e..3175001 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -3,6 +3,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) use Block_Energy use Block_Hydro use Block_Network +use Block_Reservoir ! Implicit None ! @@ -60,6 +61,60 @@ SUBROUTINE SYSTMM(temp_file,param_file) allocate (Q_na(heat_cells)) allocate (press(heat_cells)) allocate (wind(heat_cells)) +allocate (dt_res(2*heat_cells)) +allocate (resx(2*heat_cells)) +! +! reservoir allocatable +! +allocate (depth_e(nres)) +allocate (depth_h(nres)) +allocate (surface_area(nres)) +allocate (T_epil(nres)) +T_epil = 10 +allocate (T_hypo(nres)) +T_hypo = 10 +allocate (stream_T_in(nres)) +stream_T_in = 10 +allocate (density_epil(nres)) +density_epil = 0 +allocate (density_hypo(nres)) +density_hypo = 0 +allocate (density_in(nres)) +density_in = 0 +allocate (volume_e_x(nres)) +allocate (volume_h_x(nres)) +allocate (dV_dt_epi(nres)) +allocate (dV_dt_hyp(nres)) +allocate (K_z(nres)) +K_z = 0.1 +allocate (temp_change_ep(nres)) +temp_change_ep = 0 +allocate (temp_change_hyp(nres)) +temp_change_hyp = 0 +allocate (res_run(nres)) +res_run = .false. +allocate (res_start(nres)) +res_start = .false. +allocate (T_res(nres)) +T_res = 10 +allocate (T_res_in(nres)) +T_res_in = 10 +allocate (Q_trib_tot(heat_cells)) +allocate (T_trib_tot(heat_cells)) +allocate (Q_res_in(nres)) +allocate(temp_out(nres)) +temp_out = 10 +allocate(temp_out_i(nres)) +temp_out_i = 10 +allocate (trib_res(heat_cells)) +allocate(diffusion_tot(nres)) +allocate(advec_hyp_tot(nres)) +allocate(advec_epi_tot(nres)) +allocate(qsurf_tot(nres)) +allocate(reservoir_storage(nres)) +allocate(reservoir_storage_prev(nres)) + + ! ! Initialize some arrays ! @@ -77,6 +132,16 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! T_smth=4.0 +! +! initialize reservoir geometery variables +! +depth_e = res_depth_feet(:) * depth_e_frac * ft_to_m ! ft_to_m converst from feet to m +depth_h = res_depth_feet(:) * depth_h_frac * ft_to_m ! ft_to_m converst from feet to m +surface_area = res_width_feet(:) * res_length_feet(:) * ft_to_m * ft_to_m !ft_to_m converst from feet to m +volume_e_x = surface_area(:) * depth_e(:) +volume_h_x = surface_area(:) * depth_h(:) + + ! ! open the output file ! @@ -115,6 +180,13 @@ SUBROUTINE SYSTMM(temp_file,param_file) DO ndd=1,nwpd xdd = ndd time=year+(xd+(xdd-0.5)*hpd)/xd_year + ! ------------ reservoir arrays ----------- + res_run = .false. ! re-initialize reservoir fun T/F array + res_start = .false. ! re-initialize T/F for reservoir start + Q_trib_tot = 0 ! re-set the tributary flow to 0 + T_trib_tot = 0 ! re-set the tributary flow to 0 + trib_res = .false. + ! ! Read the hydrologic and meteorologic forcings @@ -147,46 +219,58 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! do ns=1,no_celm(nr) ! - DONE=.FALSE. + ! ----------------------------------------------------------------------- + ! + ! River loop + ! + ! ----------------------------------------------------------------------- + + + ! if segment in river reach, or the first segment of reservoir + if((res_pres(nr,segment_cell(nr,ns)) .eqv. .false.) .or. & + (any(segment_cell(nr,ns) == res_start_node(:)) .eqv. .false. .and. & + res_pres(nr,segment_cell(nr,ns-1)) .eqv. .false.) ) then + + DONE=.FALSE. -! Testing new code 8/8/2016 -! -! Establish particle tracks -! - call Particle_Track(nr,ns,nx_s,nx_head) -! - ncell=segment_cell(nr,ns) -! -! Now do the third-order interpolation to -! establish the starting temperature values -! for each parcel -! - nseg=nstrt_elm(ns) -! -! Perform polynomial interpolation -! -! -! Interpolation inside the domain -! - npndx=2 -! -! Interpolation at the upstream boundary if the -! parcel has reached that boundary -! - if(nx_head.eq.0) then - T_0 = T_head(nr) - else -! -! -! Interpolation at the upstream or downstream boundary -! - if(nseg .eq. 1 .or. nseg .eq. no_celm(nr)) npndx=1 -! - do ntrp=nterp(npndx),1,-1 + ! Testing new code 8/8/2016 + ! + ! Establish particle tracks + ! + call Particle_Track(nr,ns,nx_s,nx_head) + ! + ncell=segment_cell(nr,ns) + ! + ! Now do the third-order interpolation to + ! establish the starting temperature values + ! for each parcel + ! + nseg=nstrt_elm(ns) + ! + ! Perform polynomial interpolation + ! + ! + ! Interpolation inside the domain + ! + npndx=2 + ! + ! Interpolation at the upstream boundary if the + ! parcel has reached that boundary + ! + if(nx_head.eq.0) then + T_0 = T_head(nr) + else + ! + ! + ! Interpolation at the upstream or downstream boundary + ! + if(nseg .eq. 1 .or. nseg .eq. no_celm(nr)) npndx=1 + ! + do ntrp=nterp(npndx),1,-1 npart=nseg+ntrp+ndltp(npndx) xa(ntrp)=x_dist(nr,npart) ta(ntrp)=temp(nr,npart,n1) - end do + end do ! ! Start the cell counter for nx_s ! @@ -283,9 +367,26 @@ SUBROUTINE SYSTMM(temp_file,param_file) end if dt_total=dt_total+dt_calc end do - if (T_0.lt.0.5) T_0=0.5 + + end if ! end river loop + + ! ----------------------------------------------------------------------- + ! + ! Reservoir Loop + ! + ! ----------------------------------------------------------------------- + + ! ------------- if cell is in reservoir ---------------- + if(reservoir .and. res_pres(nr,segment_cell(nr,ns))) then + + + + end if ! end reservoir loop + + + if (T_0.lt.0.5) T_0=0.5 temp(nr,ns,n2)=T_0 - T_trib(nr)=T_0 + T_trib(nr)=T_0 ! ! Write all temperature output UW_JRY_11/08/2013 ! The temperature is output at the beginning of the @@ -294,6 +395,10 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! value of ndelta (now a vector)(UW_JRY_11/08/2013) ! call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell),Q_inflow,Q_outflow) + + +if(ncell .eq. 82) write(67,*) nyear,',', ndays,',',nd,',',ncell,',',ns,',',T_0 + ! ! End of computational element loop ! diff --git a/src/makefile b/src/makefile index 11a9706..bb0709f 100644 --- a/src/makefile +++ b/src/makefile @@ -8,7 +8,7 @@ objects = rbm10_VIC.o Begin.o Systmm.o Particle_Track.o \ Energy.o Julian.o tntrp.o Read_Forcing.o \ Block_Energy.o Block_Hydro.o Block_Network.o\ - Water_Balance.o Write.o + Block_Reservoir.o Water_Balance.o Write.o f90comp = gfortran # Makefile rbm10_VIC: $(objects) @@ -25,7 +25,11 @@ Block_Network.o: Block_Network.f90 $(f90comp) -c Block_Network.f90 block_network.mod: Block_Network.o Block_Network.f90 $(f90comp) -c Block_Network.f90 -Begin.o: block_energy.mod block_network.mod block_hydro.mod Begin.f90 +Block_Reservoir.o: Block_Reservoir.f90 + $(f90comp) -c Block_Reservoir.f90 +block_reservoir.mod: Block_Reservoir.o Block_Reservoir.f90 + $(f90comp) -c Block_Reservoir.f90 +Begin.o: block_energy.mod block_network.mod block_hydro.mod block_reservoir.mod Begin.f90 $(f90comp) -c Begin.f90 Read_Forcing.o: block_energy.mod block_hydro.mod block_network.mod Read_Forcing.f90 $(f90comp) -c Read_Forcing.f90 @@ -49,5 +53,5 @@ rbm10_VIC.o: rbm10_VIC.f90 # Cleaning everything clean: rm block_energy.mod block_hydro.mod block_network.mod\ - rbm10_VIC + rbm10_VIC block_reservoir.mod rm $(objects) From 965ade8d3fe347d03ff7932d5c4123733aadbab3 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Tue, 7 Nov 2017 10:04:48 -0800 Subject: [PATCH 06/27] fix small bugs in RBM in particle tracking --- src/Begin.f90 | 4 +-- src/Block_Network.f90 | 3 +-- src/Block_Reservoir.f90 | 1 + src/Particle_Track.f90 | 8 +++--- src/Read_Forcing.f90 | 9 ++++--- src/Systmm.f90 | 37 ++++++++++++++++++--------- src/flow_subroutine.f90 | 28 +++++++++++++++++--- src/makefile | 3 +++ src/reservoir_subroutine_implicit.f90 | 8 ++++-- 9 files changed, 72 insertions(+), 29 deletions(-) diff --git a/src/Begin.f90 b/src/Begin.f90 index 9441be4..b82809f 100755 --- a/src/Begin.f90 +++ b/src/Begin.f90 @@ -175,14 +175,14 @@ Subroutine BEGIN(param_file,spatial_file) if (reservoir) then read(90,'(5x,i5,5x,i5,8x,i5,6x,a8,6x,a10,7x,f10.0,f5.0,i6)') & node,nrow,ncol,lat,long,rmile1,ndelta(ncell),res_num(ncell) - write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell),res_num(ncell) + !write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell),res_num(ncell) if(res_num(ncell) .gt. 0) then res_pres(ncell) = .TRUE. end if else read(90,'(5x,i5,5x,i5,8x,i5,6x,a8,6x,a10,7x,f10.0,f5.0)') & node,nrow,ncol,lat,long,rmile1,ndelta(ncell) - write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell) + !write(*,*) node,nrow,ncol,lat,long,rmile1,ndelta(ncell) end if ! ! Set the number of segments of the default, if not specified diff --git a/src/Block_Network.f90 b/src/Block_Network.f90 index f96ef53..546bb8b 100755 --- a/src/Block_Network.f90 +++ b/src/Block_Network.f90 @@ -12,7 +12,7 @@ Module Block_Network ! integer:: flow_cells,heat_cells integer:: ndays,nreach,ntrb,nwpd - integer,parameter::ns_max=200 + integer,parameter::ns_max=400 integer:: start_year,start_month,start_day integer:: end_year,end_month,end_day ! @@ -24,7 +24,6 @@ Module Block_Network ! ! Logical variables ! - logical :: source ! presence of thermal plant diff --git a/src/Block_Reservoir.f90 b/src/Block_Reservoir.f90 index f32eefa..8a3588d 100755 --- a/src/Block_Reservoir.f90 +++ b/src/Block_Reservoir.f90 @@ -54,6 +54,7 @@ module Block_Reservoir real, dimension(:), allocatable :: qsurf_tot real, dimension(:), allocatable :: res_capacity_mcm real, dimension(:), allocatable :: volume_h_min + real, dimension(:), allocatable :: volume_e_min real, dimension(:,:), allocatable :: res_storage ! ! diff --git a/src/Particle_Track.f90 b/src/Particle_Track.f90 index 1b46f3e..b9707bb 100755 --- a/src/Particle_Track.f90 +++ b/src/Particle_Track.f90 @@ -26,7 +26,7 @@ SUBROUTINE Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) .and. .NOT. ns_res_pres) nx_s=nx_s+1 x_part(nx_s)=x_dist(nr,nx_part_next) - ncell=segment_cell(nr,nx_part_next) + ncell=segment_cell(nr,nx_part) dt_part(nx_s)=dt(ncell) dt_total=dt_total+dt_part(nx_s) ! @@ -35,7 +35,9 @@ SUBROUTINE Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) ! if (nx_s.gt.1) nx_part=nx_part-1 nx_part_next=nx_part_next-1 - ns_res_pres = res_pres(segment_cell(nr,nx_part_next)) + if (nx_part_next .gt. 0) then + ns_res_pres = res_pres(segment_cell(nr,nx_part_next)) + end if end do ! If water particle travels back to a reservoir, ! ns_res_num denotes the number of the grid cell. @@ -52,7 +54,7 @@ SUBROUTINE Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) dt_part(nx_s)=dt_part(nx_s)-dt_xcess x_part(nx_s)=x_dist(nr,nx_part+1)+u(ncell)*dt_part(nx_s) end if - nx_head=nx_part + nx_head=nx_part_next nx_part=max(1,nx_part) nstrt_elm(ns)=nx_part no_dt(ns)=nx_s diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index 1a654e2..3bbb359 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -11,6 +11,7 @@ SUBROUTINE Read_Forcing integer :: nreservoir,n1,n2 real :: Q_avg,Q_dmmy real :: z_temp,w_temp + real :: min_flow = 5.0 ! n1=1 n2=2 @@ -31,10 +32,10 @@ SUBROUTINE Read_Forcing ! If the streamflow is 0, recalculate water velocity based on ! local streamflow ! - if (Q_out(no_heat) .eq. 0) then - z_temp = a_z * (Q_local(no_heat)**b_z) - w_temp = a_w * (Q_local(no_heat)**b_w) - u(no_heat) = Q_local(no_heat)/(z_temp*w_temp) + if (Q_out(no_heat) .lt. min_flow) then + z_temp = a_z * (min_flow**b_z) + w_temp = a_w * (min_flow**b_w) + u(no_heat) = min_flow/(z_temp*w_temp) depth(no_heat) = z_temp width(no_heat) = w_temp end if diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 8e32da0..b75f691 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -114,6 +114,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) allocate (depth_e(nres)) allocate (depth_h(nres)) allocate (volume_h_min(nres)) + allocate (volume_e_min(nres)) allocate (res_stratif_start(nres)) allocate (res_turnover(nres)) ! @@ -169,11 +170,11 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! At the start of each year, reset the depth of epilimnion ! and hypolimnion. ! - !res_depth_meter = depth_e(:) + depth_h(:) - !depth_e = res_depth_meter(:) * depth_e_frac - !depth_h = res_depth_meter(:) * depth_h_frac - !volume_e_x = surface_area(:) * depth_e(:) - !volume_h_x = surface_area(:) * depth_h(:) + res_depth_meter = depth_e(:) + depth_h(:) + depth_e = res_depth_meter(:) * depth_e_frac + depth_h = res_depth_meter(:) * depth_h_frac + volume_e_x = surface_area(:) * depth_e(:) + volume_h_x = surface_area(:) * depth_h(:) ! ! Day loop starts ! @@ -239,6 +240,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Establish particle tracks ! call Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) + !if (ncell .gt. 860) write(*,*) nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num ! ! Now do the third-order interpolation to ! establish the starting temperature values @@ -251,6 +253,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! if(ns_res_pres) then T_0 = temp(nr,nseg-1,n1) + !if(ncell.eq.911) write(*,*) nd, ns, nseg, T_0, temp(nr,nseg-1,n1), n1, temp(300,50,:) else ! ! Perform polynomial interpolation @@ -276,7 +279,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) npart=nseg+ntrp+ndltp(npndx) xa(ntrp)=x_dist(nr,npart) ta(ntrp)=temp(nr,npart,n1) - !if(nr.eq.2 .and.ns.eq.2 .and. nd.lt.30) write(*,*) & + !if(ncell .eq. 3022) write(*,*) & ! 'nd', nd, 'ntrp', ntrp, 'npart',npart,xa(ntrp), ta(ntrp) end do ! @@ -287,6 +290,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Call the interpolation function ! T_0=tntrp(xa,ta,x,nterp(npndx)) + !if (ncell .eq. 3022) write(*,*) & + ! nd, ns, 'xa', xa, 'ta', ta, 'x', x, T_0,T_head(nr) end if end if ! @@ -311,12 +316,16 @@ SUBROUTINE SYSTMM(temp_file,param_file) !!end if ! do nm=no_dt(ns),1,-1 + !if (ncell .eq. 912) write(*,*) & + !nd, ns, nm, 'begin', T_0, dbt(nncell), nncell, nseg dt_calc=dt_part(nm) z=depth(nncell) call energy(T_0,q_surf,nncell) ! q_dot=(q_surf/(z*rfac)) T_0=T_0+q_dot*dt_calc + !if (ncell .eq. 912) write(*,*) & + ! nd, ns, nm, 'end', T_0, q_surf, z, q_dot*dt_calc if(T_0.lt.0.0) T_0=0.0 ! ! Add distributed flows @@ -400,7 +409,6 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! if(ncell .eq. res_start_node(res_no) .and. & .NOT. res_start(res_no)) then - if(res_no.eq.3 .and. nd.lt.5) write(*,*) res_start_node(res_no) ns_res_start(res_no) = ns Q_res_in(res_no) = Q_in(ncell) ! Inflow for reservoir T_res_in(res_no) = temp(nr,ns-1,n1) ! Inflow temperature @@ -464,7 +472,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! If it is the last segment in the reservoir ! if(ncell .eq. res_end_node(res_no) .and. & - .NOT. res_pres(segment_cell(nr,ns+1))) then + (.NOT. res_pres(segment_cell(nr,ns+1)) & + .or. res_num(segment_cell(nr,ns+1)) .ne. res_no)) then ns_res_end(res_no) = ns Q_res_outflow(res_no) = Q_out(ncell) !!!!TESTINFLOW(need to uncomment) ! @@ -487,7 +496,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! ! Calculate flow subroutine ! - call flow_subroutine(res_no) + call flow_subroutine(res_no, nyear, nd) ! ! Energy balance ! @@ -496,11 +505,15 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Calculate reservoir temperature ! !call reservoir_subroutine(res_no,q_surf) - call reservoir_subroutine_implicit(res_no,q_surf,nd,dbt(ncell)) + !if (res_storage(res_no,1) .gt.res_capacity_mcm(res_no)*(10**6)*0.2) then + call reservoir_subroutine_implicit(res_no,q_surf,nd,dbt(ncell)) + !else + ! call reservoir_single_subroutine(res_no,q_surf,nd) + !end if ! T_0 = T_hypo(res_no) ! In reservoir, water is released from hypolimnion if (T_0.lt.0.5) T_0=0.5 - temp(nr,ns,n2)=T_0 + temp(nr,ns_res_start(res_no):ns_res_end(res_no),n2)=T_0 T_trib(nr)=T_0 ! ! Write output @@ -543,7 +556,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) if (reservoir) then write(75,*) nyear, nd, T_epil(:), T_hypo(:), T_res_inflow(:), K_z(:), depth_e(:), depth_h(:),& Q_res_inflow(:), Q_res_outflow(:) - write(76,'(2i6, 37f15.10, 37f15.10)') nyear, nd, depth_e(:), depth_h(:), volume_e_x,volume_h_x + write(76,'(i8, i6, 240f15.10, 240f15.10)') nyear, nd, depth_e(:), depth_h(:), volume_e_x,volume_h_x write(77,*) nyear, nd, diffusion_tot, advec_hyp_tot,advec_epi_tot,qsurf_tot end if end do diff --git a/src/flow_subroutine.f90 b/src/flow_subroutine.f90 index dd9d2ca..1d913a5 100755 --- a/src/flow_subroutine.f90 +++ b/src/flow_subroutine.f90 @@ -1,4 +1,4 @@ -SUBROUTINE flow_subroutine (res_no) +SUBROUTINE flow_subroutine (res_no, nyear, nd) use Block_Hydro use Block_Network @@ -10,7 +10,7 @@ SUBROUTINE flow_subroutine (res_no) real :: ratio_sp, ratio_pen real :: Q1, Q2 !real :: res_vol_delta_x, vol_change_hyp_x, vol_change_epi_x - integer :: res_no, nnd + integer :: res_no, nnd, nyear, nd !************************************************************************* @@ -55,18 +55,36 @@ SUBROUTINE flow_subroutine (res_no) flow_out_hyp_x = Q2 ! * ftsec_to_msec * dt_comp flow_out_epi_x = 0 flow_epi_hyp_x = flow_in_epi_x + ! + ! if res_storage_post > res_storage_pre, when reservoir is storing water + ! + if (Q1 - Q2 .gt. 0 .and. flow_in_epi_x .gt. 0 & + .and. volume_e_x(res_no) .lt. res_capacity_mcm(res_no)*(10**6)*0.3) then + flow_epi_hyp_x = flow_out_hyp_x + end if + ! ! based on inflow and outflow + ! vol_change_epi_x = flow_in_epi_x - flow_out_epi_x - flow_epi_hyp_x vol_change_hyp_x = flow_in_hyp_x + flow_epi_hyp_x - flow_out_hyp_x ! ! Check whether hypolimnion volume is smaller than minimum hypolimnion volume ! - volume_h_min(res_no) = 0 !res_capacity_mcm(res_no) * (10**6) * 0.1 + !volume_h_min(res_no) = 0 !res_capacity_mcm(res_no) * (10**6) * 0.1 + volume_h_min(res_no) = res_capacity_mcm(res_no) * (10**6) * 0.05 if ((volume_h_x(res_no) + vol_change_hyp_x) .lt. volume_h_min(res_no)) then - vol_change_epi_x = vol_change_epi_x + vol_change_hyp_x + & + vol_change_epi_x = vol_change_epi_x + vol_change_hyp_x + & (volume_h_x(res_no) - volume_h_min(res_no)) vol_change_hyp_x = - (volume_h_x(res_no) - volume_h_min(res_no)) end if + + volume_e_min(res_no) = res_capacity_mcm(res_no) * (10**6) * 0.05 + if (res_no .eq. 19) write(*,*) volume_h_min(res_no), volume_e_min(res_no) + if ((volume_e_x(res_no) + vol_change_epi_x) .lt. volume_e_min(res_no)) then + vol_change_hyp_x = vol_change_epi_x + vol_change_hyp_x + & + (volume_e_x(res_no) - volume_e_min(res_no)) + vol_change_epi_x = - (volume_e_x(res_no) - volume_e_min(res_no)) + end if !----- update epilimnion and hypolimnion volume ------- volume_e_x(res_no) = volume_e_x(res_no) + vol_change_epi_x volume_h_x(res_no) = volume_h_x(res_no) + vol_change_hyp_x @@ -77,4 +95,6 @@ SUBROUTINE flow_subroutine (res_no) depth_e(res_no) = volume_e_x(res_no) / surface_area(res_no) depth_h(res_no) = volume_h_x(res_no) / surface_area(res_no) + !if(res_no .eq. 11) write(*,*) vol_change_epi_x, & + ! vol_change_hyp_x, depth_e(res_no), depth_h(res_no) END SUBROUTINE flow_subroutine diff --git a/src/makefile b/src/makefile index bba850d..d31ec82 100644 --- a/src/makefile +++ b/src/makefile @@ -11,6 +11,7 @@ objects = rbm10_VIC.o Begin.o Systmm.o Particle_Track.o \ Block_Reservoir.o \ density_subroutine.o flow_subroutine.o \ reservoir_subroutine_implicit.o \ + reservoir_single_subroutine.o \ Water_Balance.o Write.o f90comp = gfortran # Makefile @@ -50,6 +51,8 @@ flow_subroutine.o: block_hydro.mod block_network.mod block_reservoir.mod flow_su $(f90comp) -c flow_subroutine.f90 reservoir_subroutine_implicit.o: block_network.mod block_reservoir.mod reservoir_subroutine_implicit.f90 $(f90comp) -c reservoir_subroutine_implicit.f90 +reservoir_single_subroutine.o: block_network.mod block_reservoir.mod reservoir_single_subroutine.f90 + $(f90comp) -c reservoir_single_subroutine.f90 Julian.o: Julian.f90 $(f90comp) -c Julian.f90 tntrp.o: tntrp.f90 diff --git a/src/reservoir_subroutine_implicit.f90 b/src/reservoir_subroutine_implicit.f90 index 7eede87..cf5a377 100755 --- a/src/reservoir_subroutine_implicit.f90 +++ b/src/reservoir_subroutine_implicit.f90 @@ -52,7 +52,7 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) if (temp_epi - temp_hyp .gt. 0.2 .and. .NOT. res_stratif_start(res_no)) then res_stratif_start(res_no) = .TRUE. end if - + if (temp_epi - temp_hyp .lt. 1 .and. nd.gt.180 .and. .NOT. res_turnover(res_no)) then res_turnover(res_no) = .TRUE. end if @@ -93,7 +93,11 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) ! depth_e(res_no) = volume_e_x(res_no) / surface_area(res_no) depth_h(res_no) = volume_h_x(res_no) / surface_area(res_no) - ! + T_res(res_no) = (T_epil(res_no) * (volume_e_x(res_no)/volume_tot)) + & (T_hypo(res_no)*(volume_h_x(res_no)/volume_tot) ) ! weighted averge temp + ! + if(res_no .eq. 19) write(*,*) nd, T_epil(res_no), T_hypo(res_no),& +volume_e_x(res_no)/surface_area(res_no), & +volume_h_x(res_no)/surface_area(res_no) end subroutine reservoir_subroutine_implicit From 8e22071b3379ab4b56361314e34fbfb7b12bb4fb Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Mon, 29 Jan 2018 15:18:30 -0800 Subject: [PATCH 07/27] fix minor error in particle tracking branch --- src/Particle_Track.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Particle_Track.f90 b/src/Particle_Track.f90 index b9707bb..079afec 100755 --- a/src/Particle_Track.f90 +++ b/src/Particle_Track.f90 @@ -29,12 +29,13 @@ SUBROUTINE Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) ncell=segment_cell(nr,nx_part) dt_part(nx_s)=dt(ncell) dt_total=dt_total+dt_part(nx_s) + if (dt_total .le. dt_comp) nx_part_next=nx_part_next-1 ! ! Increment the segment counter if the total time is less than the ! computational interval ! if (nx_s.gt.1) nx_part=nx_part-1 - nx_part_next=nx_part_next-1 + !nx_part_next=nx_part_next-1 if (nx_part_next .gt. 0) then ns_res_pres = res_pres(segment_cell(nr,nx_part_next)) end if From 0ad3424984be76a2c6e9f209dcdf2e445a7857b5 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Tue, 30 Jan 2018 14:58:30 -0800 Subject: [PATCH 08/27] test for instability in RBM add some output sentence to check the output --- src/Particle_Track.f90 | 5 ++++- src/Read_Forcing.f90 | 4 ++++ src/Systmm.f90 | 16 ++++++++-------- src/flow_subroutine.f90 | 2 +- src/reservoir_subroutine_implicit.f90 | 6 +++--- 5 files changed, 20 insertions(+), 13 deletions(-) diff --git a/src/Particle_Track.f90 b/src/Particle_Track.f90 index b9707bb..f7745ce 100755 --- a/src/Particle_Track.f90 +++ b/src/Particle_Track.f90 @@ -29,12 +29,13 @@ SUBROUTINE Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) ncell=segment_cell(nr,nx_part) dt_part(nx_s)=dt(ncell) dt_total=dt_total+dt_part(nx_s) + if (dt_total.le.dt_comp) nx_part_next=nx_part_next-1 ! ! Increment the segment counter if the total time is less than the ! computational interval ! if (nx_s.gt.1) nx_part=nx_part-1 - nx_part_next=nx_part_next-1 + !nx_part_next=nx_part_next-1 if (nx_part_next .gt. 0) then ns_res_pres = res_pres(segment_cell(nr,nx_part_next)) end if @@ -58,4 +59,6 @@ SUBROUTINE Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) nx_part=max(1,nx_part) nstrt_elm(ns)=nx_part no_dt(ns)=nx_s + if (nr.eq.1078 .and. ns.le.2) write(*,*) 'particle track', nx_head, & + nx_part, nstrt_elm(ns), no_dt(ns), nx_s, x_part(nx_s) END SUBROUTINE Particle_Track diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index 3bbb359..b2dc20f 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -45,6 +45,10 @@ SUBROUTINE Read_Forcing if(u(no_heat).lt.0.01) u(no_heat)=0.01 if(ncell.ne.no_heat) write(*,*) 'Flow file error',ncell,no_heat ! + !#############just for test + if(no_heat .eq. 3442) write(*,*) 'output',Q_out(no_heat), depth(no_heat), & + width(no_heat), u(no_heat) + !#############just for test read(36,*) ncell & ,dbt(no_heat),ea(no_heat) & ,Q_ns(no_heat),Q_na(no_heat),rho & diff --git a/src/Systmm.f90 b/src/Systmm.f90 index b75f691..3f5e26f 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -279,8 +279,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) npart=nseg+ntrp+ndltp(npndx) xa(ntrp)=x_dist(nr,npart) ta(ntrp)=temp(nr,npart,n1) - !if(ncell .eq. 3022) write(*,*) & - ! 'nd', nd, 'ntrp', ntrp, 'npart',npart,xa(ntrp), ta(ntrp) + if(ncell .eq. 3442) write(*,*) & + 'nd', nd, 'ntrp', ntrp, 'npart',npart,xa(ntrp), ta(ntrp) end do ! ! Start the cell counter for nx_s @@ -290,8 +290,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Call the interpolation function ! T_0=tntrp(xa,ta,x,nterp(npndx)) - !if (ncell .eq. 3022) write(*,*) & - ! nd, ns, 'xa', xa, 'ta', ta, 'x', x, T_0,T_head(nr) + if (ncell .eq. 3442) write(*,*) & + nd, ns, 'xa', xa, 'ta', ta, 'x', x, T_0,T_head(nr) end if end if ! @@ -316,16 +316,16 @@ SUBROUTINE SYSTMM(temp_file,param_file) !!end if ! do nm=no_dt(ns),1,-1 - !if (ncell .eq. 912) write(*,*) & - !nd, ns, nm, 'begin', T_0, dbt(nncell), nncell, nseg + if (ncell .eq. 3442) write(*,*) & + nd, ns, nm, 'begin', T_0, T_head(nr), dbt(nncell), nncell, nseg dt_calc=dt_part(nm) z=depth(nncell) call energy(T_0,q_surf,nncell) ! q_dot=(q_surf/(z*rfac)) T_0=T_0+q_dot*dt_calc - !if (ncell .eq. 912) write(*,*) & - ! nd, ns, nm, 'end', T_0, q_surf, z, q_dot*dt_calc + if (ncell .eq. 3442) write(*,*) & + nd, ns, nm, 'end', T_0, q_surf, z, q_dot*dt_calc, dt_calc if(T_0.lt.0.0) T_0=0.0 ! ! Add distributed flows diff --git a/src/flow_subroutine.f90 b/src/flow_subroutine.f90 index 1d913a5..f96f95f 100755 --- a/src/flow_subroutine.f90 +++ b/src/flow_subroutine.f90 @@ -79,7 +79,7 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) end if volume_e_min(res_no) = res_capacity_mcm(res_no) * (10**6) * 0.05 - if (res_no .eq. 19) write(*,*) volume_h_min(res_no), volume_e_min(res_no) + !if (res_no .eq. 19) write(*,*) volume_h_min(res_no), volume_e_min(res_no) if ((volume_e_x(res_no) + vol_change_epi_x) .lt. volume_e_min(res_no)) then vol_change_hyp_x = vol_change_epi_x + vol_change_hyp_x + & (volume_e_x(res_no) - volume_e_min(res_no)) diff --git a/src/reservoir_subroutine_implicit.f90 b/src/reservoir_subroutine_implicit.f90 index cf5a377..bf31608 100755 --- a/src/reservoir_subroutine_implicit.f90 +++ b/src/reservoir_subroutine_implicit.f90 @@ -97,7 +97,7 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) T_res(res_no) = (T_epil(res_no) * (volume_e_x(res_no)/volume_tot)) + & (T_hypo(res_no)*(volume_h_x(res_no)/volume_tot) ) ! weighted averge temp ! - if(res_no .eq. 19) write(*,*) nd, T_epil(res_no), T_hypo(res_no),& -volume_e_x(res_no)/surface_area(res_no), & -volume_h_x(res_no)/surface_area(res_no) + !if(res_no .eq. 19) write(*,*) nd, T_epil(res_no), T_hypo(res_no),& + ! volume_e_x(res_no)/surface_area(res_no), & + ! volume_h_x(res_no)/surface_area(res_no) end subroutine reservoir_subroutine_implicit From 2a9817497d6ee28742958acd8b010a285e4028a3 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 1 Feb 2018 10:17:23 -0800 Subject: [PATCH 09/27] calculate the error aroused because of numerical schemes implemented --- src/Block_Energy.f90 | 2 ++ src/Energy.f90 | 17 ++++++++++++++++- src/Systmm.f90 | 36 +++++++++++++++++++++++++++++++----- 3 files changed, 49 insertions(+), 6 deletions(-) diff --git a/src/Block_Energy.f90 b/src/Block_Energy.f90 index c284301..15a2e3a 100755 --- a/src/Block_Energy.f90 +++ b/src/Block_Energy.f90 @@ -32,6 +32,8 @@ module Block_Energy ! Some important constants ! real :: lvp,rb,rho + real :: deriv_2nd,error_EE + real :: deriv_conv,deriv_evap,deriv_ws real,parameter :: evap_coeff=1.5e-9 !Lake Hefner coefficient, 1/meters real,parameter :: pf=0.640,pi=3.14159 real,parameter :: rfac=304.8 !rho/Cp kg/meter**3/Kilocalories/kg/Deg K diff --git a/src/Energy.f90 b/src/Energy.f90 index 3f9995a..908176b 100755 --- a/src/Energy.f90 +++ b/src/Energy.f90 @@ -1,8 +1,10 @@ -SUBROUTINE Energy(T_surf,q_surf,ncell) +SUBROUTINE Energy(T_surf,q_surf,ncell,z) use Block_Energy implicit none integer::i,ncell,nd real::A,B,e0,q_surf,q_conv,q_evap,q_ws,td,T_surf + real::const1,const2,z + real::deriv_1st !,deriv_conv,deriv_evap,deriv_ws real, dimension(2):: q_fit, T_fit ! td=nd @@ -29,9 +31,22 @@ SUBROUTINE Energy(T_surf,q_surf,ncell) A=(q_fit(1)-q_fit(2))/(T_fit(1)-T_fit(2)) q_surf=0.5*(q_fit(1)+q_fit(2)) B=(q_surf/A)-(T_fit(1)+T_fit(2))/2. + + deriv_1st=(q_surf/(z*rfac)) ! ! ****************************************************** ! Return to Subroutine RIVMOD ! ****************************************************** ! +! Calculate the second derivative of Temperature +! + const1=1000.*evap_coeff*wind(ncell)*pf + const2=1000.*evap_coeff*wind(ncell) + e0=2.1718E8*EXP(-4157.0/(T_surf+239.09)) + deriv_conv=const1*(-(597.0-0.57*T_surf) -0.57*(dbt(ncell)-T_surf)) + deriv_evap=const2*(-0.57*(e0-ea(ncell))+(597.0-0.57*T_surf)*e0*(4157.0/((T_surf+239.09)**2))) + deriv_ws=1.471E-3 + + deriv_2nd=deriv_1st*(deriv_conv-deriv_evap-deriv_ws)/(z*rfac) + END Subroutine Energy diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 3f5e26f..3783c82 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -25,7 +25,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) integer, dimension(2):: ndltp=(/-1,-2/) integer, dimension(2):: nterp=(/2,3/) ! - real :: dt_calc,dt_total,hpd,q_dot,q_surf,z + real :: dt_calc,dt_total,hpd,q_dot,q_surf,z,q_dot_pre real :: Q_dstrb,Q_inflow,Q_outflow,Q_ratio,Q_inflow_origin real :: Q_inflow_out, Q_outflow_out, Q_tot real :: sto_pre, sto_post @@ -320,13 +320,33 @@ SUBROUTINE SYSTMM(temp_file,param_file) nd, ns, nm, 'begin', T_0, T_head(nr), dbt(nncell), nncell, nseg dt_calc=dt_part(nm) z=depth(nncell) - call energy(T_0,q_surf,nncell) + call energy(T_0,q_surf,nncell,z) + ! + ! apply a different numerical method to solve + ! energy balance for river temperature ! q_dot=(q_surf/(z*rfac)) + ! + ! Initialize dH/dt (temperature increase because + ! of energy balance) + ! + if (nyear.eq.start_year .and. nd.eq.1) q_dot_pre = q_dot + ! + ! caluclate stream temperature based on Heun + ! mehtod + ! T_0=T_0+q_dot*dt_calc + ! + ! + ! + error_EE=deriv_2nd*dt_calc**2/2 if (ncell .eq. 3442) write(*,*) & - nd, ns, nm, 'end', T_0, q_surf, z, q_dot*dt_calc, dt_calc + nd, ns, nm, 'end', T_0, z, q_dot*dt_calc, dt_calc, error_EE, & + deriv_conv*q_dot/deriv_2nd/(z*rfac), & + deriv_evap*q_dot/deriv_2nd/(z*rfac), & + deriv_ws*q_dot/deriv_2nd/(z*rfac) if(T_0.lt.0.0) T_0=0.0 + q_dot_pre = q_dot ! ! Add distributed flows ! @@ -355,8 +375,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) end if end do end if - Q_trb_sum = Q_trb_sum/2 - T_trb_load = T_trb_load/2 + Q_trb_sum = Q_trb_sum/ndelta(nncell) + T_trb_load = T_trb_load/ndelta(nncell) ! ! Do the mass/energy balance ! @@ -395,6 +415,12 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! ! if the segment is located in reservoir ! + ! + ! test output for a specific grid cell + ! + if (ncell.eq.3442) write(34,*) nyear, nd, & + ncell,ns,T_0,T_head(nr),dbt(ncell) + ! call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell), & Q_inflow_out, Q_outflow_out) end if From 3d32e8f4ffcd821e35d330e7608f170bb5287703 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 8 Feb 2018 09:27:34 -0800 Subject: [PATCH 10/27] add estimation of second order error --- src/Systmm.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 3783c82..1012f38 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -419,7 +419,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! test output for a specific grid cell ! if (ncell.eq.3442) write(34,*) nyear, nd, & - ncell,ns,T_0,T_head(nr),dbt(ncell) + ncell,ns,T_0,T_head(nr),dbt(ncell),error_EE ! call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell), & Q_inflow_out, Q_outflow_out) From ea88e91a0b926e283156176262654b25e1aae87b Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 8 Feb 2018 09:30:56 -0800 Subject: [PATCH 11/27] comment out printing statement (for test) --- src/Systmm.f90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 1012f38..d10d2f8 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -279,8 +279,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) npart=nseg+ntrp+ndltp(npndx) xa(ntrp)=x_dist(nr,npart) ta(ntrp)=temp(nr,npart,n1) - if(ncell .eq. 3442) write(*,*) & - 'nd', nd, 'ntrp', ntrp, 'npart',npart,xa(ntrp), ta(ntrp) + !if(ncell .eq. 3442) write(*,*) & + ! 'nd', nd, 'ntrp', ntrp, 'npart',npart,xa(ntrp), ta(ntrp) end do ! ! Start the cell counter for nx_s @@ -290,8 +290,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Call the interpolation function ! T_0=tntrp(xa,ta,x,nterp(npndx)) - if (ncell .eq. 3442) write(*,*) & - nd, ns, 'xa', xa, 'ta', ta, 'x', x, T_0,T_head(nr) + !if (ncell .eq. 3442) write(*,*) & + ! nd, ns, 'xa', xa, 'ta', ta, 'x', x, T_0,T_head(nr) end if end if ! @@ -316,8 +316,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) !!end if ! do nm=no_dt(ns),1,-1 - if (ncell .eq. 3442) write(*,*) & - nd, ns, nm, 'begin', T_0, T_head(nr), dbt(nncell), nncell, nseg + !if (ncell .eq. 3442) write(*,*) & + ! nd, ns, nm, 'begin', T_0, T_head(nr), dbt(nncell), nncell, nseg dt_calc=dt_part(nm) z=depth(nncell) call energy(T_0,q_surf,nncell,z) @@ -340,11 +340,11 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! ! error_EE=deriv_2nd*dt_calc**2/2 - if (ncell .eq. 3442) write(*,*) & - nd, ns, nm, 'end', T_0, z, q_dot*dt_calc, dt_calc, error_EE, & - deriv_conv*q_dot/deriv_2nd/(z*rfac), & - deriv_evap*q_dot/deriv_2nd/(z*rfac), & - deriv_ws*q_dot/deriv_2nd/(z*rfac) + !if (ncell .eq. 3442) write(*,*) & + ! nd, ns, nm, 'end', T_0, z, q_dot*dt_calc, dt_calc, error_EE, & + ! deriv_conv*q_dot/deriv_2nd/(z*rfac), & + ! deriv_evap*q_dot/deriv_2nd/(z*rfac), & + ! deriv_ws*q_dot/deriv_2nd/(z*rfac) if(T_0.lt.0.0) T_0=0.0 q_dot_pre = q_dot ! @@ -418,8 +418,8 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! ! test output for a specific grid cell ! - if (ncell.eq.3442) write(34,*) nyear, nd, & - ncell,ns,T_0,T_head(nr),dbt(ncell),error_EE + !if (ncell.eq.3442) write(34,*) nyear, nd, & + ! ncell,ns,T_0,T_head(nr),dbt(ncell),error_EE ! call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell), & Q_inflow_out, Q_outflow_out) From 06acd9a97e605073c53436c23947f5bddb8d4917 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 8 Feb 2018 10:40:22 -0800 Subject: [PATCH 12/27] add output filtering --- src/Begin.f90 | 11 ++++++++++- src/Block_Network.f90 | 4 +++- src/Systmm.f90 | 9 +++++++-- 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/src/Begin.f90 b/src/Begin.f90 index b82809f..325f673 100755 --- a/src/Begin.f90 +++ b/src/Begin.f90 @@ -17,7 +17,7 @@ Subroutine BEGIN(param_file,spatial_file) integer:: jul_start,main_stem,nyear1,nyear2,nc,ncell,nseg integer:: ns_max_test,node,ncol,nrow,nr,cum_sgmnt ! - integer:: nreservoir + integer:: nreservoir,nseg_temp ! logical:: first_cell,source ! @@ -83,6 +83,7 @@ Subroutine BEGIN(param_file,spatial_file) allocate(res_start_node(nres)) allocate(res_end_node(nres)) allocate(res_capacity_mcm(nres)) + allocate(nseg_out(nres,400,nseg_out_num)) ! ! Start reading the reach date and initialize the reach index, NR ! and the cell index, NCELL @@ -203,6 +204,14 @@ Subroutine BEGIN(param_file,spatial_file) nndlta=nndlta+1 nseg=nseg+1 segment_cell(nr,nseg)=ncell + ! + ! Here we define the output segments + ! + do nseg_temp=1,nseg_out_num + if (ndelta(ncell)/nseg.eq.nseg_temp) then + nseg_out(nr,ncell,nseg_temp)=nseg + end if + end do !write(*,*) 'nndlta -- ',nr,nndlta,nseg,ncell,segment_cell(nr,nseg) x_dist(nr,nseg)=x_dist(nr,nseg-1)-dx(ncell) ! diff --git a/src/Block_Network.f90 b/src/Block_Network.f90 index 546bb8b..7024361 100755 --- a/src/Block_Network.f90 +++ b/src/Block_Network.f90 @@ -7,12 +7,14 @@ Module Block_Network ! integer, dimension(:,:), allocatable::conflnce,reach_cell,segment_cell,trib ! + integer, dimension(:,:,:), allocatable::nseg_out ! ! Integer variables ! integer:: flow_cells,heat_cells integer:: ndays,nreach,ntrb,nwpd - integer,parameter::ns_max=400 + integer,parameter::ns_max=4000 + integer,parameter::nseg_out_num=2 integer:: start_year,start_month,start_day integer:: end_year,end_month,end_day ! diff --git a/src/Systmm.f90 b/src/Systmm.f90 index d10d2f8..6fd549e 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -18,6 +18,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) integer :: n1,n2,nnd,nobs,nyear,nd_year,ntmp integer :: npart,nseg,nx_s,nx_part,nx_head integer :: ns_res_num, res_no, nsegment + integer :: nseg_temp ! ! Indices for lagrangian interpolation ! @@ -421,8 +422,12 @@ SUBROUTINE SYSTMM(temp_file,param_file) !if (ncell.eq.3442) write(34,*) nyear, nd, & ! ncell,ns,T_0,T_head(nr),dbt(ncell),error_EE ! - call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell), & - Q_inflow_out, Q_outflow_out) + do nseg_temp=1,nseg_out_num + if (nseg_out(nr,ncell,nseg_temp).eq.ns) then + call WRITE(time,nd,nr,ncell,nseg_temp,T_0,T_head(nr),dbt(ncell), & + Q_inflow_out, Q_outflow_out) + end if + end do end if if (res_pres(ncell)) then From b680dffe833aac3bcc7383bda4d9b4c4c7e198bb Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 8 Feb 2018 15:45:05 -0800 Subject: [PATCH 13/27] Fix the number of segments output for each grid cell. --- src/Begin.f90 | 21 +++++++++++---------- src/Block_Network.f90 | 2 +- src/Systmm.f90 | 14 +++++++++----- 3 files changed, 21 insertions(+), 16 deletions(-) diff --git a/src/Begin.f90 b/src/Begin.f90 index 325f673..6a34e35 100755 --- a/src/Begin.f90 +++ b/src/Begin.f90 @@ -17,7 +17,7 @@ Subroutine BEGIN(param_file,spatial_file) integer:: jul_start,main_stem,nyear1,nyear2,nc,ncell,nseg integer:: ns_max_test,node,ncol,nrow,nr,cum_sgmnt ! - integer:: nreservoir,nseg_temp + integer:: nreservoir,nseg_temp,nseg_cum ! logical:: first_cell,source ! @@ -83,7 +83,7 @@ Subroutine BEGIN(param_file,spatial_file) allocate(res_start_node(nres)) allocate(res_end_node(nres)) allocate(res_capacity_mcm(nres)) - allocate(nseg_out(nres,400,nseg_out_num)) + allocate(nseg_out(nreach,heat_cells,nseg_out_num)) ! ! Start reading the reach date and initialize the reach index, NR ! and the cell index, NCELL @@ -118,6 +118,7 @@ Subroutine BEGIN(param_file,spatial_file) ! Initialize NSEG, the total number of segments in this reach ! nseg=0 + nseg_cum=0 write(*,*) ' Starting to read reach ',nr ! ! Read the number of cells in this reach, the headwater #, @@ -199,19 +200,18 @@ Subroutine BEGIN(param_file,spatial_file) ! dx(ncell)=miles_to_ft*(rmile0-rmile1)/ndelta(ncell) rmile0=rmile1 + ! + ! Here we define the output segments + ! + do nseg_temp=1,nseg_out_num + nseg_out(nr,ncell,nseg_temp)=nseg_cum+ndelta(ncell)*nseg_temp/(nseg_out_num) + end do + nseg_cum = nseg_cum+ndelta(ncell) nndlta=0 200 continue nndlta=nndlta+1 nseg=nseg+1 segment_cell(nr,nseg)=ncell - ! - ! Here we define the output segments - ! - do nseg_temp=1,nseg_out_num - if (ndelta(ncell)/nseg.eq.nseg_temp) then - nseg_out(nr,ncell,nseg_temp)=nseg - end if - end do !write(*,*) 'nndlta -- ',nr,nndlta,nseg,ncell,segment_cell(nr,nseg) x_dist(nr,nseg)=x_dist(nr,nseg-1)-dx(ncell) ! @@ -228,6 +228,7 @@ Subroutine BEGIN(param_file,spatial_file) no_celm(nr)=nseg segment_cell(nr,nseg)=ncell x_dist(nr,nseg)=miles_to_ft*rmile1 + write(*,*) 'number of segment in reach', nr, nseg ! ! End of cell and segment loop ! diff --git a/src/Block_Network.f90 b/src/Block_Network.f90 index 7024361..db369cf 100755 --- a/src/Block_Network.f90 +++ b/src/Block_Network.f90 @@ -13,7 +13,7 @@ Module Block_Network ! integer:: flow_cells,heat_cells integer:: ndays,nreach,ntrb,nwpd - integer,parameter::ns_max=4000 + integer,parameter::ns_max=1000 integer,parameter::nseg_out_num=2 integer:: start_year,start_month,start_day integer:: end_year,end_month,end_day diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 6fd549e..2ea6d7d 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -424,7 +424,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! do nseg_temp=1,nseg_out_num if (nseg_out(nr,ncell,nseg_temp).eq.ns) then - call WRITE(time,nd,nr,ncell,nseg_temp,T_0,T_head(nr),dbt(ncell), & + call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell), & Q_inflow_out, Q_outflow_out) end if end do @@ -550,10 +550,14 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! Write output ! do nsegment=ns_res_start(res_no),ns_res_end(res_no) - call WRITE(time,nd,nr,ncell,nsegment,T_0, & - T_head(nr),dbt(segment_cell(nr,nsegment)), & - Q_res_inflow(res_no), Q_res_outflow(res_no), & - res_storage_post, T_res(res_no)) + do nseg_temp=1,nseg_out_num + if (nseg_out(nr,ncell,nseg_temp).eq.nsegment) then + call WRITE(time,nd,nr,ncell,nsegment,T_0, & + T_head(nr),dbt(segment_cell(nr,nsegment)), & + Q_res_inflow(res_no), Q_res_outflow(res_no), & + res_storage_post, T_res(res_no)) + end if + end do end do end if end if From 6b552fa49730d4538e4fcb7f3a9a8ce4c6da4fc2 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 8 Feb 2018 22:58:50 -0800 Subject: [PATCH 14/27] fix output statement --- src/Systmm.f90 | 7 +++++-- src/rbm10_VIC.f90 | 7 +++++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 2ea6d7d..6137ff5 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -1,4 +1,4 @@ -SUBROUTINE SYSTMM(temp_file,param_file) +SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! use Block_Energy use Block_Hydro @@ -10,6 +10,7 @@ SUBROUTINE SYSTMM(temp_file,param_file) ! character (len=200):: temp_file character (len=200):: param_file + character (len=200):: res_file ! integer :: ncell,nncell,ncell0,nc_head,no_flow,no_heat integer :: nc,nd,ndd,nm,nr,ns @@ -589,7 +590,9 @@ SUBROUTINE SYSTMM(temp_file,param_file) 4700 format(f10.4,f6.0,15(f6.1,f8.3)) 4750 format(f10.4,10(i4,f8.0)) if (reservoir) then - write(75,*) nyear, nd, T_epil(:), T_hypo(:), T_res_inflow(:), K_z(:), depth_e(:), depth_h(:),& + open(75,file=TRIM(res_file),status='unknown') + write(75, '(i8, i6, 1680f20.10)') nyear, nd, T_epil(:), T_hypo(:), T_res_inflow(:),& + depth_e(:), depth_h(:),& Q_res_inflow(:), Q_res_outflow(:) write(76,'(i8, i6, 240f15.10, 240f15.10)') nyear, nd, depth_e(:), depth_h(:), volume_e_x,volume_h_x write(77,*) nyear, nd, diffusion_tot, advec_hyp_tot,advec_epi_tot,qsurf_tot diff --git a/src/rbm10_VIC.f90 b/src/rbm10_VIC.f90 index 0044a30..eaf3bb5 100755 --- a/src/rbm10_VIC.f90 +++ b/src/rbm10_VIC.f90 @@ -50,6 +50,7 @@ PROGRAM RBM10_VIC character (len=200 ):: net_file character (len=200 ):: param_file character (len=200 ):: temp_file + character (len=200 ):: res_file character (len=200 ):: spatial_file character (len=200 ):: reservoir_file character (len=200 ):: reservoir_storage_file @@ -80,11 +81,13 @@ PROGRAM RBM10_VIC param_file = TRIM(inPrefix)//'_Parameters' spatial_file = TRIM(outPrefix)//'.Spat' temp_file = TRIM(outPrefix)//'.Temp' + res_file = TRIM(outPrefix)//'.Resv' ! - write(*,*) 'Spatial file: ',spatial_file write(*,*) 'Network file : ',net_file write(*,*) 'Parameter file : ',param_file! write(*,*) 'Temperature file: ',temp_file + write(*,*) 'Spatial file: ',spatial_file + write(*,*) 'Reservoir file: ',res_file ! OPEN(UNIT=90,FILE=TRIM(net_file),STATUS='OLD') ! @@ -130,7 +133,7 @@ PROGRAM RBM10_VIC ! ! SUBROUTINE SYSTMM performs the simulations ! - CALL SYSTMM(temp_file,param_file) ! (WUR_WF_MvV_2011/01/05) + CALL SYSTMM(temp_file,res_file,param_file) ! (WUR_WF_MvV_2011/01/05) ! ! Close files after simulation is complete ! From f7e0c062bfb4465ad91b9bc4b103cd92cf7c5787 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Mon, 12 Feb 2018 14:43:22 -0800 Subject: [PATCH 15/27] debug - allocated space is not sufficient for variables in Block_Hydro.f90 --- src/Begin.f90 | 8 ++++++-- src/Block_Hydro.f90 | 4 ++-- src/Block_Network.f90 | 2 +- src/Particle_Track.f90 | 4 ++-- src/Systmm.f90 | 25 ++----------------------- src/makefile | 2 +- src/reservoir_subroutine_implicit.f90 | 3 +++ 7 files changed, 17 insertions(+), 31 deletions(-) diff --git a/src/Begin.f90 b/src/Begin.f90 index 6a34e35..94ac68f 100755 --- a/src/Begin.f90 +++ b/src/Begin.f90 @@ -217,8 +217,12 @@ Subroutine BEGIN(param_file,spatial_file) ! ! Write Segment List for mapping to temperature output (UW_JRY_2008/11/19) ! - open(22,file=TRIM(spatial_file),status='unknown') ! (changed by WUR_WF_MvV_2011/01/05) - write(22,'(4i6,1x,a8,1x,a10,f5.0)') nr,ncell,nrow,ncol,lat,long,nndlta + do nseg_temp=1,nseg_out_num + if (nseg_out(nr,ncell,nseg_temp).eq.nseg) then + open(22,file=TRIM(spatial_file),status='unknown') ! (changed by WUR_WF_MvV_2011/01/05) + write(22,'(4i6,1x,a8,1x,a10,i5)') nr,ncell,nrow,ncol,lat,long,nseg_temp + end if + end do ! ! ! diff --git a/src/Block_Hydro.f90 b/src/Block_Hydro.f90 index 6434bd5..f40a8f2 100755 --- a/src/Block_Hydro.f90 +++ b/src/Block_Hydro.f90 @@ -2,8 +2,8 @@ ! Module for hydraulic characteristics and water quality constituents of the basin ! module Block_Hydro - integer, dimension(2000):: no_dt,nstrt_elm - real, dimension(2000) :: dt_part,x_part + integer, dimension(2500):: no_dt,nstrt_elm + real, dimension(2500) :: dt_part,x_part ! real, dimension(:), allocatable :: depth real, dimension(:), allocatable :: width diff --git a/src/Block_Network.f90 b/src/Block_Network.f90 index db369cf..a2a5b87 100755 --- a/src/Block_Network.f90 +++ b/src/Block_Network.f90 @@ -13,7 +13,7 @@ Module Block_Network ! integer:: flow_cells,heat_cells integer:: ndays,nreach,ntrb,nwpd - integer,parameter::ns_max=1000 + integer,parameter::ns_max=3000 integer,parameter::nseg_out_num=2 integer:: start_year,start_month,start_day integer:: end_year,end_month,end_day diff --git a/src/Particle_Track.f90 b/src/Particle_Track.f90 index f7745ce..da47773 100755 --- a/src/Particle_Track.f90 +++ b/src/Particle_Track.f90 @@ -59,6 +59,6 @@ SUBROUTINE Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) nx_part=max(1,nx_part) nstrt_elm(ns)=nx_part no_dt(ns)=nx_s - if (nr.eq.1078 .and. ns.le.2) write(*,*) 'particle track', nx_head, & - nx_part, nstrt_elm(ns), no_dt(ns), nx_s, x_part(nx_s) + !if (nr.eq.1247 .and. ns.eq.2017) write(*,*) 'particle track', nx_head, & + ! nx_part, nstrt_elm(ns), no_dt(ns), nx_s, x_part(nx_s) END SUBROUTINE Particle_Track diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 6137ff5..1bca8ab 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -227,8 +227,6 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) DO ns=1,no_celm(nr) ! ncell=segment_cell(nr,ns) - !if (res_pres(ncell)) then - !write(*,*) ncell !end if ! ! If this segment is not located in the reservoir @@ -242,7 +240,6 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! Establish particle tracks ! call Particle_Track(nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num) - !if (ncell .gt. 860) write(*,*) nr,ns,nx_s,nx_head,ns_res_pres,ns_res_num ! ! Now do the third-order interpolation to ! establish the starting temperature values @@ -255,7 +252,6 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! if(ns_res_pres) then T_0 = temp(nr,nseg-1,n1) - !if(ncell.eq.911) write(*,*) nd, ns, nseg, T_0, temp(nr,nseg-1,n1), n1, temp(300,50,:) else ! ! Perform polynomial interpolation @@ -281,8 +277,6 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) npart=nseg+ntrp+ndltp(npndx) xa(ntrp)=x_dist(nr,npart) ta(ntrp)=temp(nr,npart,n1) - !if(ncell .eq. 3442) write(*,*) & - ! 'nd', nd, 'ntrp', ntrp, 'npart',npart,xa(ntrp), ta(ntrp) end do ! ! Start the cell counter for nx_s @@ -292,8 +286,6 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! Call the interpolation function ! T_0=tntrp(xa,ta,x,nterp(npndx)) - !if (ncell .eq. 3442) write(*,*) & - ! nd, ns, 'xa', xa, 'ta', ta, 'x', x, T_0,T_head(nr) end if end if ! @@ -312,14 +304,7 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! ! Set initial river storage (Initial storage for first grid cell is 0) ! - !!if(nyear.eq.start_year .and. nd.eq.1 .and. ndd.eq.1 .and. ns.ge.3) then - !sto(nr,ns,n1) = width(ncell) * depth(ncell) * dx(ncell)/2 - !! temp_sto(nr,ns,n1) = T_head(nr) - !!end if - ! do nm=no_dt(ns),1,-1 - !if (ncell .eq. 3442) write(*,*) & - ! nd, ns, nm, 'begin', T_0, T_head(nr), dbt(nncell), nncell, nseg dt_calc=dt_part(nm) z=depth(nncell) call energy(T_0,q_surf,nncell,z) @@ -339,14 +324,10 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! T_0=T_0+q_dot*dt_calc ! - ! + ! estimate the error by numerical method ! error_EE=deriv_2nd*dt_calc**2/2 - !if (ncell .eq. 3442) write(*,*) & - ! nd, ns, nm, 'end', T_0, z, q_dot*dt_calc, dt_calc, error_EE, & - ! deriv_conv*q_dot/deriv_2nd/(z*rfac), & - ! deriv_evap*q_dot/deriv_2nd/(z*rfac), & - ! deriv_ws*q_dot/deriv_2nd/(z*rfac) + ! if(T_0.lt.0.0) T_0=0.0 q_dot_pre = q_dot ! @@ -420,8 +401,6 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! ! test output for a specific grid cell ! - !if (ncell.eq.3442) write(34,*) nyear, nd, & - ! ncell,ns,T_0,T_head(nr),dbt(ncell),error_EE ! do nseg_temp=1,nseg_out_num if (nseg_out(nr,ncell,nseg_temp).eq.ns) then diff --git a/src/makefile b/src/makefile index d31ec82..dae9c2f 100644 --- a/src/makefile +++ b/src/makefile @@ -13,7 +13,7 @@ objects = rbm10_VIC.o Begin.o Systmm.o Particle_Track.o \ reservoir_subroutine_implicit.o \ reservoir_single_subroutine.o \ Water_Balance.o Write.o -f90comp = gfortran +f90comp = gfortran -O0 -g # Makefile rbm10_VIC: $(objects) $(f90comp) -o rbm10_VIC $(objects) diff --git a/src/reservoir_subroutine_implicit.f90 b/src/reservoir_subroutine_implicit.f90 index bf31608..dd019a8 100755 --- a/src/reservoir_subroutine_implicit.f90 +++ b/src/reservoir_subroutine_implicit.f90 @@ -96,6 +96,9 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) T_res(res_no) = (T_epil(res_no) * (volume_e_x(res_no)/volume_tot)) + & (T_hypo(res_no)*(volume_h_x(res_no)/volume_tot) ) ! weighted averge temp + + if(T_epil(res_no).lt.0.5) T_epil(res_no) = 0.5 + if(T_hypo(res_no).lt.0.5) T_hypo(res_no) = 0.5 ! !if(res_no .eq. 19) write(*,*) nd, T_epil(res_no), T_hypo(res_no),& ! volume_e_x(res_no)/surface_area(res_no), & From 0ff5c18b2408578374a9e033e21f34719ad5b370 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 22 Feb 2018 09:48:32 -0800 Subject: [PATCH 16/27] git rid of the speeding up in make file --- src/Read_Forcing.f90 | 4 ---- src/makefile | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index b2dc20f..3bbb359 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -45,10 +45,6 @@ SUBROUTINE Read_Forcing if(u(no_heat).lt.0.01) u(no_heat)=0.01 if(ncell.ne.no_heat) write(*,*) 'Flow file error',ncell,no_heat ! - !#############just for test - if(no_heat .eq. 3442) write(*,*) 'output',Q_out(no_heat), depth(no_heat), & - width(no_heat), u(no_heat) - !#############just for test read(36,*) ncell & ,dbt(no_heat),ea(no_heat) & ,Q_ns(no_heat),Q_na(no_heat),rho & diff --git a/src/makefile b/src/makefile index dae9c2f..490c53a 100644 --- a/src/makefile +++ b/src/makefile @@ -13,7 +13,7 @@ objects = rbm10_VIC.o Begin.o Systmm.o Particle_Track.o \ reservoir_subroutine_implicit.o \ reservoir_single_subroutine.o \ Water_Balance.o Write.o -f90comp = gfortran -O0 -g +f90comp = gfortran #-O2#-O0 -g # Makefile rbm10_VIC: $(objects) $(f90comp) -o rbm10_VIC $(objects) From d4be9467b35080cea0033a546a63eda8d50a4899 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Fri, 23 Feb 2018 14:39:10 -0800 Subject: [PATCH 17/27] experiment for temperature spikes in model simulation --- src/Energy.f90 | 6 ++---- src/Systmm.f90 | 2 +- src/makefile | 2 +- 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/Energy.f90 b/src/Energy.f90 index 908176b..84290b1 100755 --- a/src/Energy.f90 +++ b/src/Energy.f90 @@ -1,13 +1,12 @@ -SUBROUTINE Energy(T_surf,q_surf,ncell,z) +SUBROUTINE Energy(T_surf,q_surf,ncell,z,nd) use Block_Energy implicit none integer::i,ncell,nd - real::A,B,e0,q_surf,q_conv,q_evap,q_ws,td,T_surf + real::A,B,e0,q_surf,q_conv,q_evap,q_ws,T_surf real::const1,const2,z real::deriv_1st !,deriv_conv,deriv_evap,deriv_ws real, dimension(2):: q_fit, T_fit ! - td=nd T_fit(1)=T_surf-1.0 T_fit(2)=T_surf+1.0 do i=1,2 @@ -31,7 +30,6 @@ SUBROUTINE Energy(T_surf,q_surf,ncell,z) A=(q_fit(1)-q_fit(2))/(T_fit(1)-T_fit(2)) q_surf=0.5*(q_fit(1)+q_fit(2)) B=(q_surf/A)-(T_fit(1)+T_fit(2))/2. - deriv_1st=(q_surf/(z*rfac)) ! ! ****************************************************** diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 1bca8ab..d5e8bce 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -307,7 +307,7 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) do nm=no_dt(ns),1,-1 dt_calc=dt_part(nm) z=depth(nncell) - call energy(T_0,q_surf,nncell,z) + call energy(T_0,q_surf,nncell,z,nd) ! ! apply a different numerical method to solve ! energy balance for river temperature diff --git a/src/makefile b/src/makefile index 490c53a..713b3d0 100644 --- a/src/makefile +++ b/src/makefile @@ -13,7 +13,7 @@ objects = rbm10_VIC.o Begin.o Systmm.o Particle_Track.o \ reservoir_subroutine_implicit.o \ reservoir_single_subroutine.o \ Water_Balance.o Write.o -f90comp = gfortran #-O2#-O0 -g +f90comp = gfortran #-O0 -g #-O2 # Makefile rbm10_VIC: $(objects) $(f90comp) -o rbm10_VIC $(objects) From 10e5b6edbd3ef45e4f0796e382428004d88bde19 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Fri, 23 Feb 2018 14:44:54 -0800 Subject: [PATCH 18/27] set the minimum wind speed to be 1 m/s --- src/Read_Forcing.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index 3bbb359..2208059 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -50,6 +50,10 @@ SUBROUTINE Read_Forcing ,Q_ns(no_heat),Q_na(no_heat),rho & ,press(no_heat),wind(no_heat) ! + ! Set the minimum wind speed to be 1 [m s-1] + ! + if(wind(no_heat).lt.1.0) wind(no_heat)=1.0 + ! if(ncell.ne.no_heat) write(*,*) 'Heat file error',ncell,no_heat ! ! Added variable ndelta (UW_JRY_2011/03/15 From 71ecae41f7b26400efa1575abdfd6870a2202f96 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Sun, 25 Feb 2018 16:46:11 -0800 Subject: [PATCH 19/27] fix small bug for wind speed threshold --- src/Read_Forcing.f90 | 1 + src/Systmm.f90 | 3 +++ 2 files changed, 4 insertions(+) diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index 2208059..9c28b8a 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -85,6 +85,7 @@ SUBROUTINE Read_Forcing ,dbt(no_heat),ea(no_heat) & ,Q_ns(no_heat),Q_na(no_heat),rho & ,press(no_heat),wind(no_heat) + if(wind(no_heat).lt.1.0) wind(no_heat)=1.0 if (nr .eq. nreach .and. nc .eq. no_cells(nr)) then read(35,*) nnd,ncell & ,Q_out(no_heat),Q_dmmy,Q_diff(no_heat) & diff --git a/src/Systmm.f90 b/src/Systmm.f90 index d5e8bce..cbcb0fb 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -328,6 +328,8 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! error_EE=deriv_2nd*dt_calc**2/2 ! + !if (ncell.eq.4089 .and. nd.gt.220.and.nd.lt.230) write(*,*) & + ! nyear,nd,ns,nm,wind(nncell),q_ns(nncell)+q_na(nncell),T_0,error_EE if(T_0.lt.0.0) T_0=0.0 q_dot_pre = q_dot ! @@ -406,6 +408,7 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) if (nseg_out(nr,ncell,nseg_temp).eq.ns) then call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell), & Q_inflow_out, Q_outflow_out) + !if (ncell.eq.4089) write(*,*)nyear,nd,nr,ncell,ns,T_0 end if end do end if From f358f92bad0b1b8d200371be57396f1a141809b6 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Mon, 26 Feb 2018 12:44:40 -0800 Subject: [PATCH 20/27] check reservoir temperature instability --- src/reservoir_subroutine_implicit.f90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/reservoir_subroutine_implicit.f90 b/src/reservoir_subroutine_implicit.f90 index dd019a8..dacd6c9 100755 --- a/src/reservoir_subroutine_implicit.f90 +++ b/src/reservoir_subroutine_implicit.f90 @@ -68,7 +68,11 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) (volume_h_x(res_no) + volume_e_x(res_no) + flow_in_epi_x + flow_in_hyp_x) temp_hyp = temp_epi end if - + if (res_no.eq.66) write(99,'(8f30.6)') & + volume_h_x(res_no),volume_e_x(res_no),flow_in_epi_x,flow_in_hyp_x, & + volume_h_x(res_no)*T_hypo(res_no), volume_e_x(res_no)*T_epil(res_no),& + (flow_in_epi_x + flow_in_hyp_x) * T_res_inflow(res_no), & + q_surf * dt_comp * surface_area(res_no) T_epil(res_no) = temp_epi T_hypo(res_no) = temp_hyp !if (res_no .eq. 3) write(49,*) flow_in_hyp_x+flow_in_epi_x, T_res_inflow(res_no) From 518aabbb90577ae755ea4979e34e218817b44c8a Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Tue, 27 Feb 2018 17:54:43 -0800 Subject: [PATCH 21/27] add Error_estimate subroutine based on numerical error to adjust the timestep to calculate stream temperature --- src/Block_Network.f90 | 3 ++ src/Block_Reservoir.f90 | 9 ++++ src/Energy.f90 | 2 + src/Error_estimate.f90 | 55 +++++++++++++++++++++++ src/Systmm.f90 | 64 ++++++++++++++++++++++----- src/flow_subroutine.f90 | 14 +++--- src/makefile | 3 ++ src/reservoir_subroutine_implicit.f90 | 47 +++++++++++++++++--- 8 files changed, 172 insertions(+), 25 deletions(-) create mode 100755 src/Error_estimate.f90 diff --git a/src/Block_Network.f90 b/src/Block_Network.f90 index a2a5b87..d97714b 100755 --- a/src/Block_Network.f90 +++ b/src/Block_Network.f90 @@ -17,11 +17,14 @@ Module Block_Network integer,parameter::nseg_out_num=2 integer:: start_year,start_month,start_day integer:: end_year,end_month,end_day + integer:: numsub !number of subdaily timestep + integer:: nsub ! ! Real variables ! real :: delta_n,n_default=2 real :: dt_comp + real :: dt_res real, dimension(:), allocatable :: ndelta ! ! Logical variables diff --git a/src/Block_Reservoir.f90 b/src/Block_Reservoir.f90 index 8a3588d..bc9e2a2 100755 --- a/src/Block_Reservoir.f90 +++ b/src/Block_Reservoir.f90 @@ -12,6 +12,8 @@ module Block_Reservoir logical, dimension(:), allocatable :: res_trib_calc ! logical, dimension(:), allocatable :: res_stratif_start ! logical, dimension(:), allocatable :: res_turnover ! + logical :: exceed_error_bound + logical :: adjust_timestep ! ! real, parameter :: depth_e_frac=0.4, depth_h_frac=0.6 @@ -23,6 +25,13 @@ module Block_Reservoir real :: flow_epi_hyp_x, outflow_x real :: res_vol_delta_x, vol_change_hyp_x, vol_change_epi_x real :: res_storage_pre, res_storage_post + ! + ! Error estimation + ! + real :: m11,m12,m21,m22 + real :: const1,const2 + real :: error_e,error_h + real, parameter :: error_threshold = 0.1 real, dimension(:), allocatable :: K_z real, dimension(:), allocatable :: depth_e, depth_h real, dimension(:), allocatable :: density_epil, density_hypo diff --git a/src/Energy.f90 b/src/Energy.f90 index 84290b1..ca5ad11 100755 --- a/src/Energy.f90 +++ b/src/Energy.f90 @@ -19,6 +19,8 @@ SUBROUTINE Energy(T_surf,q_surf,ncell,z,nd) q_evap=q_evap*(e0-ea(ncell)) q_ws=6.693E-2+1.471E-3*T_fit(i) q_fit(i)=q_ns(ncell)+q_na(ncell)-q_ws-q_evap+q_conv + if (ncell.eq.2034) write(88,*) & + q_ns(ncell),q_na(ncell),q_ws,q_evap,q_conv end do ! ! q=AT+B diff --git a/src/Error_estimate.f90 b/src/Error_estimate.f90 new file mode 100755 index 0000000..8b6c4da --- /dev/null +++ b/src/Error_estimate.f90 @@ -0,0 +1,55 @@ +SUBROUTINE Error_estimate(nd,res_no) + use Block_Network + use Block_Reservoir + ! + implicit none + ! + real :: dayx, q_surf, log_K_z, n_stability, density_dif + real :: temp_epi_pre, temp_hyp_pre + integer :: nd,res_no + + temp_epi_pre = T_epil(res_no) + temp_hyp_pre = T_hypo(res_no) + ! ---------------- turnover loop driven only by T_epil and T_hyp + ! ---------- + density_dif = density_epil(res_no) - density_hypo(res_no) + n_stability = (-1) * (gravity/density_epil(res_no)) * & + ((density_dif)/((depth_e(res_no) + depth_h(res_no))/2)) + if(n_stability .le. 0) then + n_stability = 1e-10 + end if + log_K_z = log10(n_stability) * (-1) - 5.699 ! high scenario - 2 w/ + K_z(res_no) = 10**log_K_z + ! ONLY for no stratification run: + write(57, *) res_no,K_z(res_no) + ! + ! Estimate the numerical error in reservoir temperature + ! + if (temp_epi_pre-temp_hyp_pre>0.01) then + m11 = -(flow_in_epi_x+K_z(res_no)*surface_area(res_no))/volume_e_x(res_no)/dt_res + m12 = K_z(res_no)*surface_area(res_no)/volume_e_x(res_no)/dt_res + m21 = (flow_epi_hyp_x+K_z(res_no)*surface_area(res_no))/volume_h_x(res_no)/dt_res + m22 = -(flow_in_hyp_x+flow_epi_hyp_x+K_z(res_no)*surface_area(res_no))/volume_h_x(res_no)/dt_res + const1 = ((q_surf*dt_res*surface_area(res_no))/(density*heat_c_kcal) + & + flow_in_epi_x*T_res_inflow(res_no))/volume_e_x(res_no)/dt_res + const2 = flow_in_hyp_x*T_res_inflow(res_no)/volume_h_x(res_no)/dt_res + error_e = -0.5*((m11**2+m12*m21)*temp_epi_pre + & + (m11*m12+m12*m22)*temp_hyp_pre + & + m11*const1 + m12*const2)*dt_res**2 + error_h = -0.5*((m21*m11+m22*m21)*temp_epi_pre + & + (m21*m12+m22**2)*temp_hyp_pre + & + m21*const1 + m22*const2)*dt_res**2 + else + m11 = -(flow_in_epi_x+flow_in_hyp_x)/(volume_e_x(res_no)+volume_h_x(res_no))/dt_res + const1 = ((flow_in_epi_x+flow_in_hyp_x)*T_res_inflow(res_no) + & + (q_surf*dt_res*surface_area(res_no))/(density*heat_c_kcal))& + /(volume_e_x(res_no)+volume_h_x(res_no))/dt_res + error_e = -0.5*(m11**2*temp_epi_pre+m11*const1) + error_h = error_e + end if + if (res_no.eq.66) write(*,*) res_no,nd,dt_res,error_e,error_h + if (abs(error_e) > error_threshold.or. abs(error_h)>error_threshold) then + exceed_error_bound=.TRUE. + if (res_no.eq.66) write(*,*) 'new',res_no,nd,exceed_error_bound + end if +end subroutine Error_estimate diff --git a/src/Systmm.f90 b/src/Systmm.f90 index cbcb0fb..d14ddd2 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -488,6 +488,8 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) if(ncell .eq. res_end_node(res_no) .and. & (.NOT. res_pres(segment_cell(nr,ns+1)) & .or. res_num(segment_cell(nr,ns+1)) .ne. res_no)) then + exceed_error_bound=.FALSE. + adjust_timestep=.FALSE. ns_res_end(res_no) = ns Q_res_outflow(res_no) = Q_out(ncell) !!!!TESTINFLOW(need to uncomment) ! @@ -505,22 +507,60 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! Calculate the density of inflow, compared with epilimnion/hypolimnion ! call stream_density(T_res_inflow(res_no), density_in(res_no)) - call stream_density(T_epil(res_no), density_epil(res_no)) - call stream_density(T_hypo(res_no), density_hypo(res_no)) + !call stream_density(T_epil(res_no),density_epil(res_no)) + !call stream_density(T_hypo(res_no),density_hypo(res_no)) ! - ! Calculate flow subroutine + ! Start the subdaily calculation ! - call flow_subroutine(res_no, nyear, nd) +500 continue + numsub=1 ! - ! Energy balance + ! calculate the number of timestep based + ! on error bound ! - call energy(T_epil(res_no), q_surf, ncell) - ! - ! Calculate reservoir temperature - ! - !call reservoir_subroutine(res_no,q_surf) - !if (res_storage(res_no,1) .gt.res_capacity_mcm(res_no)*(10**6)*0.2) then - call reservoir_subroutine_implicit(res_no,q_surf,nd,dbt(ncell)) + if (exceed_error_bound) then + numsub=MAX(ceiling(sqrt(abs(error_e)/error_threshold)), & + ceiling(sqrt(abs(error_h)/error_threshold))) + numsub=MIN(20,numsub) + end if + dt_res = dt_comp/numsub + DO nsub=1,numsub + ! + ! Calculate stream density + ! + call stream_density(T_epil(res_no), density_epil(res_no)) + call stream_density(T_hypo(res_no), density_hypo(res_no)) + ! + ! Calculate flow subroutine + ! + call flow_subroutine(res_no, nyear, nd) + ! + ! Energy balance + ! + call energy(T_epil(res_no), q_surf, ncell) + ! + ! Error estimation + ! + if (.not.exceed_error_bound) then + call Error_estimate(nd,res_no) + end if + ! + ! Adjust the timestep based on error + ! + if (exceed_error_bound .and..not.adjust_timestep) then + adjust_timestep=.TRUE. + go to 500 + end if + ! + ! Calculate reservoir temperature + ! + !call reservoir_subroutine(res_no,q_surf) + !if (res_storage(res_no,1) .gt.res_capacity_mcm(res_no)*(10**6)*0.2) then + call reservoir_subroutine_implicit(res_no,q_surf,nd,dbt(ncell)) + !if(res_no.eq.63) write(*,*) & + ! nyear,nd,nsub,dt_res,T_epil(res_no),T_hypo(res_no) + end do + !else ! call reservoir_single_subroutine(res_no,q_surf,nd) !end if diff --git a/src/flow_subroutine.f90 b/src/flow_subroutine.f90 index f96f95f..3c729bf 100755 --- a/src/flow_subroutine.f90 +++ b/src/flow_subroutine.f90 @@ -16,12 +16,14 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) !************************************************************************* ! read inflow from vic/rvic simulations !************************************************************************* - res_storage_pre = res_storage(res_no,2) - res_storage_post = res_storage(res_no,1) + if (.not.adjust_timestep) then + res_storage_pre = res_storage(res_no,2) + res_storage_post = res_storage(res_no,1) + end if !res_storage_pre = res_capacity_mcm(res_no) * (10**6) * 0.95 ! TESTINFLOW !res_storage_post = res_capacity_mcm(res_no) * (10**6) * 0.95 ! TESTINFLOW - Q1 = Q_res_inflow(res_no) * ftsec_to_msec * dt_comp + Q1 = Q_res_inflow(res_no) * ftsec_to_msec * dt_res ! converts ft^3/sec to m^3/sec, multiplys by seconds per time step if ( density_in(res_no) .le. density_hypo(res_no) ) then @@ -34,7 +36,7 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) ! ! calculate the reservoir storage ! - Q2 = Q_res_outflow(res_no) * ftsec_to_msec * dt_comp + Q2 = Q_res_outflow(res_no) * ftsec_to_msec * dt_res ! ! Initialization of reservoir storage ! @@ -51,8 +53,8 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) ! ! calculate water withdrawal based on inflow/outflow and storage change ! - water_withdrawal(res_no) = Q1 - Q2 - (res_storage_post - res_storage_pre) - flow_out_hyp_x = Q2 ! * ftsec_to_msec * dt_comp + if (nsub.eq.1) water_withdrawal(res_no) = Q1 - Q2 - (res_storage_post - res_storage_pre)/numsub + flow_out_hyp_x = Q2 ! * ftsec_to_msec * dt_res flow_out_epi_x = 0 flow_epi_hyp_x = flow_in_epi_x ! diff --git a/src/makefile b/src/makefile index 713b3d0..7da8270 100644 --- a/src/makefile +++ b/src/makefile @@ -10,6 +10,7 @@ objects = rbm10_VIC.o Begin.o Systmm.o Particle_Track.o \ Block_Energy.o Block_Hydro.o Block_Network.o \ Block_Reservoir.o \ density_subroutine.o flow_subroutine.o \ + Error_estimate.o \ reservoir_subroutine_implicit.o \ reservoir_single_subroutine.o \ Water_Balance.o Write.o @@ -49,6 +50,8 @@ density_subroutine.o: density_subroutine.f90 $(f90comp) -c density_subroutine.f90 flow_subroutine.o: block_hydro.mod block_network.mod block_reservoir.mod flow_subroutine.f90 $(f90comp) -c flow_subroutine.f90 +Error_estimate.o: block_network.mod block_reservoir.mod Error_estimate.f90 + $(f90comp) -c Error_estimate.f90 reservoir_subroutine_implicit.o: block_network.mod block_reservoir.mod reservoir_subroutine_implicit.f90 $(f90comp) -c reservoir_subroutine_implicit.f90 reservoir_single_subroutine.o: block_network.mod block_reservoir.mod reservoir_single_subroutine.f90 diff --git a/src/reservoir_subroutine_implicit.f90 b/src/reservoir_subroutine_implicit.f90 index dacd6c9..7b6f9b3 100755 --- a/src/reservoir_subroutine_implicit.f90 +++ b/src/reservoir_subroutine_implicit.f90 @@ -33,13 +33,43 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) K_z(res_no) = 10**log_K_z ! ONLY for no stratification run: write(57, *) res_no,K_z(res_no) + ! + ! Estimate the numerical error in reservoir temperature + ! +! if (temp_epi_pre-temp_hyp_pre>0.01) then +! m11 = -(flow_in_epi_x+K_z(res_no)*surface_area(res_no))/volume_e_x(res_no)/dt_res +! m12 = K_z(res_no)*surface_area(res_no)/volume_e_x(res_no)/dt_res +! m21 = (flow_epi_hyp_x+K_z(res_no)*surface_area(res_no))/volume_h_x(res_no)/dt_res +! m22 = -(flow_in_hyp_x+flow_epi_hyp_x+K_z(res_no)*surface_area(res_no))/volume_h_x(res_no)/dt_res +! const1 = ((q_surf*dt_res*surface_area(res_no))/(density*heat_c_kcal) + & +! flow_in_epi_x*T_res_inflow(res_no))/volume_e_x(res_no)/dt_res +! const2 = flow_in_hyp_x*T_res_inflow(res_no)/volume_h_x(res_no)/dt_res +! error_e = -0.5*((m11**2+m12*m21)*temp_epi_pre + & +! (m11*m12+m12*m22)*temp_hyp_pre + & +! m11*const1 + m12*const2)*dt_res**2 +! error_h = -0.5*((m21*m11+m22*m21)*temp_epi_pre + & +! (m21*m12+m22**2)*temp_hyp_pre + & +! m21*const1 + m22*const2)*dt_res**2 +! else +! m11 = -(flow_in_epi_x+flow_in_hyp_x)/(volume_e_x(res_no)+volume_h_x(res_no))/dt_res +! const1 = ((flow_in_epi_x+flow_in_hyp_x)*T_res_inflow(res_no) + & +! (q_surf*dt_res*surface_area(res_no))/(density*heat_c_kcal))& +! /(volume_e_x(res_no)+volume_h_x(res_no))/dt_res +! error_e = -0.5*(m11**2*temp_epi_pre+m11*const1) +! error_h = error_e +! end if +! if (res_no.eq.66) write(*,*) res_no,nd,dt_res,error_e,error_h +! if (error_e.gt.error_threshold.or.error_h.gt.error_threshold) then +! exceed_error_bound=.TRUE. +! go to 500 +! end if ! ! Implicit method to solve the energy balance equation ! coeff_e1 = volume_e_x(res_no) + flow_in_epi_x + K_z(res_no) * surface_area(res_no) coeff_h1 = - K_z(res_no) * surface_area(res_no) const_1 = volume_e_x(res_no)*temp_epi_pre + flow_in_epi_x*T_res_inflow(res_no) + & - (q_surf * dt_comp * surface_area(res_no)) / (density * heat_c_kcal) + (q_surf * dt_res * surface_area(res_no)) / (density * heat_c_kcal) coeff_e2 = -(flow_epi_hyp_x + K_z(res_no) * surface_area(res_no)) coeff_h2 = volume_h_x(res_no) + flow_in_hyp_x + flow_epi_hyp_x & + K_z(res_no) * surface_area(res_no) @@ -64,15 +94,15 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) if (temp_epi - temp_hyp .le. 0 .or. res_turnover(res_no)) then temp_epi = (volume_h_x(res_no)*T_hypo(res_no) + volume_e_x(res_no)*T_epil(res_no) + & (flow_in_epi_x + flow_in_hyp_x) * T_res_inflow(res_no) + & - (q_surf * dt_comp * surface_area(res_no)) / (density * heat_c_kcal)) / & + (q_surf * dt_res * surface_area(res_no)) / (density * heat_c_kcal)) / & (volume_h_x(res_no) + volume_e_x(res_no) + flow_in_epi_x + flow_in_hyp_x) temp_hyp = temp_epi end if - if (res_no.eq.66) write(99,'(8f30.6)') & - volume_h_x(res_no),volume_e_x(res_no),flow_in_epi_x,flow_in_hyp_x, & - volume_h_x(res_no)*T_hypo(res_no), volume_e_x(res_no)*T_epil(res_no),& - (flow_in_epi_x + flow_in_hyp_x) * T_res_inflow(res_no), & - q_surf * dt_comp * surface_area(res_no) + !if (res_no.eq.66) write(99,'(8f30.6)') & + ! volume_h_x(res_no),volume_e_x(res_no),flow_in_epi_x,flow_in_hyp_x, & + ! volume_h_x(res_no)*T_hypo(res_no), volume_e_x(res_no)*T_epil(res_no),& + ! (flow_in_epi_x + flow_in_hyp_x) * T_res_inflow(res_no), & + ! q_surf * dt_res * surface_area(res_no) T_epil(res_no) = temp_epi T_hypo(res_no) = temp_hyp !if (res_no .eq. 3) write(49,*) flow_in_hyp_x+flow_in_epi_x, T_res_inflow(res_no) @@ -103,6 +133,9 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) if(T_epil(res_no).lt.0.5) T_epil(res_no) = 0.5 if(T_hypo(res_no).lt.0.5) T_hypo(res_no) = 0.5 + !if(res_no.eq.1) write(*,*) & + ! 'volume',nsub,volume_e_x(res_no),volume_h_x(res_no),water_withdrawal(res_no),& + ! depth_e(res_no),depth_h(res_no) ! !if(res_no .eq. 19) write(*,*) nd, T_epil(res_no), T_hypo(res_no),& ! volume_e_x(res_no)/surface_area(res_no), & From 25d01919379998967ccfd02e4aa7178e3d0440d2 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Tue, 6 Mar 2018 13:16:58 -0800 Subject: [PATCH 22/27] fix bugs in dynamic time step to calculate reservoir temperature this bug is resulted from reservoir storage calculation. --- src/Block_Reservoir.f90 | 1 + src/Energy.f90 | 5 ++-- src/Error_estimate.f90 | 4 +-- src/Systmm.f90 | 8 ++--- src/flow_subroutine.f90 | 14 ++++++--- src/reservoir_subroutine_implicit.f90 | 43 --------------------------- 6 files changed, 18 insertions(+), 57 deletions(-) diff --git a/src/Block_Reservoir.f90 b/src/Block_Reservoir.f90 index bc9e2a2..ef262d1 100755 --- a/src/Block_Reservoir.f90 +++ b/src/Block_Reservoir.f90 @@ -14,6 +14,7 @@ module Block_Reservoir logical, dimension(:), allocatable :: res_turnover ! logical :: exceed_error_bound logical :: adjust_timestep + logical :: recalculate_volume ! ! real, parameter :: depth_e_frac=0.4, depth_h_frac=0.6 diff --git a/src/Energy.f90 b/src/Energy.f90 index ca5ad11..8d4b881 100755 --- a/src/Energy.f90 +++ b/src/Energy.f90 @@ -19,8 +19,9 @@ SUBROUTINE Energy(T_surf,q_surf,ncell,z,nd) q_evap=q_evap*(e0-ea(ncell)) q_ws=6.693E-2+1.471E-3*T_fit(i) q_fit(i)=q_ns(ncell)+q_na(ncell)-q_ws-q_evap+q_conv - if (ncell.eq.2034) write(88,*) & - q_ns(ncell),q_na(ncell),q_ws,q_evap,q_conv + !if (ncell.eq.4480) write(*,*) & + ! 'energy',q_ns(ncell),q_na(ncell),-q_ws,-q_evap,q_conv, & + ! 'air',dbt(ncell),wind(ncell) end do ! ! q=AT+B diff --git a/src/Error_estimate.f90 b/src/Error_estimate.f90 index 8b6c4da..a9909a4 100755 --- a/src/Error_estimate.f90 +++ b/src/Error_estimate.f90 @@ -44,12 +44,10 @@ SUBROUTINE Error_estimate(nd,res_no) const1 = ((flow_in_epi_x+flow_in_hyp_x)*T_res_inflow(res_no) + & (q_surf*dt_res*surface_area(res_no))/(density*heat_c_kcal))& /(volume_e_x(res_no)+volume_h_x(res_no))/dt_res - error_e = -0.5*(m11**2*temp_epi_pre+m11*const1) + error_e = -0.5*(m11**2*temp_epi_pre+m11*const1)*dt_res**2 error_h = error_e end if - if (res_no.eq.66) write(*,*) res_no,nd,dt_res,error_e,error_h if (abs(error_e) > error_threshold.or. abs(error_h)>error_threshold) then exceed_error_bound=.TRUE. - if (res_no.eq.66) write(*,*) 'new',res_no,nd,exceed_error_bound end if end subroutine Error_estimate diff --git a/src/Systmm.f90 b/src/Systmm.f90 index d14ddd2..914fcac 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -490,6 +490,7 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) .or. res_num(segment_cell(nr,ns+1)) .ne. res_no)) then exceed_error_bound=.FALSE. adjust_timestep=.FALSE. + recalculate_volume=.FALSE. ns_res_end(res_no) = ns Q_res_outflow(res_no) = Q_out(ncell) !!!!TESTINFLOW(need to uncomment) ! @@ -549,16 +550,13 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! if (exceed_error_bound .and..not.adjust_timestep) then adjust_timestep=.TRUE. + recalculate_volume=.TRUE. go to 500 end if ! ! Calculate reservoir temperature ! - !call reservoir_subroutine(res_no,q_surf) - !if (res_storage(res_no,1) .gt.res_capacity_mcm(res_no)*(10**6)*0.2) then - call reservoir_subroutine_implicit(res_no,q_surf,nd,dbt(ncell)) - !if(res_no.eq.63) write(*,*) & - ! nyear,nd,nsub,dt_res,T_epil(res_no),T_hypo(res_no) + call reservoir_subroutine_implicit(res_no,q_surf,nd,dbt(ncell)) end do !else diff --git a/src/flow_subroutine.f90 b/src/flow_subroutine.f90 index 3c729bf..a7f0448 100755 --- a/src/flow_subroutine.f90 +++ b/src/flow_subroutine.f90 @@ -20,6 +20,16 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) res_storage_pre = res_storage(res_no,2) res_storage_post = res_storage(res_no,1) end if + ! + ! if we need to adjust the timestep, we need to make the volume + ! go back to the beginning of the timestep. + ! + if (recalculate_volume) then + volume_e_x(res_no) = volume_e_x(res_no) - vol_change_epi_x + volume_h_x(res_no) = volume_h_x(res_no) - vol_change_hyp_x + recalculate_volume=.FALSE. ! This only needs to be done once + end if + !res_storage_pre = res_capacity_mcm(res_no) * (10**6) * 0.95 ! TESTINFLOW !res_storage_post = res_capacity_mcm(res_no) * (10**6) * 0.95 ! TESTINFLOW @@ -72,7 +82,6 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) ! ! Check whether hypolimnion volume is smaller than minimum hypolimnion volume ! - !volume_h_min(res_no) = 0 !res_capacity_mcm(res_no) * (10**6) * 0.1 volume_h_min(res_no) = res_capacity_mcm(res_no) * (10**6) * 0.05 if ((volume_h_x(res_no) + vol_change_hyp_x) .lt. volume_h_min(res_no)) then vol_change_epi_x = vol_change_epi_x + vol_change_hyp_x + & @@ -81,7 +90,6 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) end if volume_e_min(res_no) = res_capacity_mcm(res_no) * (10**6) * 0.05 - !if (res_no .eq. 19) write(*,*) volume_h_min(res_no), volume_e_min(res_no) if ((volume_e_x(res_no) + vol_change_epi_x) .lt. volume_e_min(res_no)) then vol_change_hyp_x = vol_change_epi_x + vol_change_hyp_x + & (volume_e_x(res_no) - volume_e_min(res_no)) @@ -97,6 +105,4 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) depth_e(res_no) = volume_e_x(res_no) / surface_area(res_no) depth_h(res_no) = volume_h_x(res_no) / surface_area(res_no) - !if(res_no .eq. 11) write(*,*) vol_change_epi_x, & - ! vol_change_hyp_x, depth_e(res_no), depth_h(res_no) END SUBROUTINE flow_subroutine diff --git a/src/reservoir_subroutine_implicit.f90 b/src/reservoir_subroutine_implicit.f90 index 7b6f9b3..49796fe 100755 --- a/src/reservoir_subroutine_implicit.f90 +++ b/src/reservoir_subroutine_implicit.f90 @@ -33,36 +33,6 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) K_z(res_no) = 10**log_K_z ! ONLY for no stratification run: write(57, *) res_no,K_z(res_no) - ! - ! Estimate the numerical error in reservoir temperature - ! -! if (temp_epi_pre-temp_hyp_pre>0.01) then -! m11 = -(flow_in_epi_x+K_z(res_no)*surface_area(res_no))/volume_e_x(res_no)/dt_res -! m12 = K_z(res_no)*surface_area(res_no)/volume_e_x(res_no)/dt_res -! m21 = (flow_epi_hyp_x+K_z(res_no)*surface_area(res_no))/volume_h_x(res_no)/dt_res -! m22 = -(flow_in_hyp_x+flow_epi_hyp_x+K_z(res_no)*surface_area(res_no))/volume_h_x(res_no)/dt_res -! const1 = ((q_surf*dt_res*surface_area(res_no))/(density*heat_c_kcal) + & -! flow_in_epi_x*T_res_inflow(res_no))/volume_e_x(res_no)/dt_res -! const2 = flow_in_hyp_x*T_res_inflow(res_no)/volume_h_x(res_no)/dt_res -! error_e = -0.5*((m11**2+m12*m21)*temp_epi_pre + & -! (m11*m12+m12*m22)*temp_hyp_pre + & -! m11*const1 + m12*const2)*dt_res**2 -! error_h = -0.5*((m21*m11+m22*m21)*temp_epi_pre + & -! (m21*m12+m22**2)*temp_hyp_pre + & -! m21*const1 + m22*const2)*dt_res**2 -! else -! m11 = -(flow_in_epi_x+flow_in_hyp_x)/(volume_e_x(res_no)+volume_h_x(res_no))/dt_res -! const1 = ((flow_in_epi_x+flow_in_hyp_x)*T_res_inflow(res_no) + & -! (q_surf*dt_res*surface_area(res_no))/(density*heat_c_kcal))& -! /(volume_e_x(res_no)+volume_h_x(res_no))/dt_res -! error_e = -0.5*(m11**2*temp_epi_pre+m11*const1) -! error_h = error_e -! end if -! if (res_no.eq.66) write(*,*) res_no,nd,dt_res,error_e,error_h -! if (error_e.gt.error_threshold.or.error_h.gt.error_threshold) then -! exceed_error_bound=.TRUE. -! go to 500 -! end if ! ! Implicit method to solve the energy balance equation ! @@ -98,14 +68,8 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) (volume_h_x(res_no) + volume_e_x(res_no) + flow_in_epi_x + flow_in_hyp_x) temp_hyp = temp_epi end if - !if (res_no.eq.66) write(99,'(8f30.6)') & - ! volume_h_x(res_no),volume_e_x(res_no),flow_in_epi_x,flow_in_hyp_x, & - ! volume_h_x(res_no)*T_hypo(res_no), volume_e_x(res_no)*T_epil(res_no),& - ! (flow_in_epi_x + flow_in_hyp_x) * T_res_inflow(res_no), & - ! q_surf * dt_res * surface_area(res_no) T_epil(res_no) = temp_epi T_hypo(res_no) = temp_hyp - !if (res_no .eq. 3) write(49,*) flow_in_hyp_x+flow_in_epi_x, T_res_inflow(res_no) ! ! Update reservoir storage based on water withdrawal ! @@ -133,11 +97,4 @@ SUBROUTINE reservoir_subroutine_implicit(res_no,q_surf,nd,tair) if(T_epil(res_no).lt.0.5) T_epil(res_no) = 0.5 if(T_hypo(res_no).lt.0.5) T_hypo(res_no) = 0.5 - !if(res_no.eq.1) write(*,*) & - ! 'volume',nsub,volume_e_x(res_no),volume_h_x(res_no),water_withdrawal(res_no),& - ! depth_e(res_no),depth_h(res_no) - ! - !if(res_no .eq. 19) write(*,*) nd, T_epil(res_no), T_hypo(res_no),& - ! volume_e_x(res_no)/surface_area(res_no), & - ! volume_h_x(res_no)/surface_area(res_no) end subroutine reservoir_subroutine_implicit From d80444af3fd4f5dc6b2f3fc89160a86a6abe8382 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 8 Mar 2018 14:29:40 -0800 Subject: [PATCH 23/27] Set the initial storage to be 95% of reservoir capacity --- src/Systmm.f90 | 7 +++---- src/flow_subroutine.f90 | 6 +++--- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 914fcac..7d0146d 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -551,6 +551,9 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) if (exceed_error_bound .and..not.adjust_timestep) then adjust_timestep=.TRUE. recalculate_volume=.TRUE. + if (nyear.eq.start_year.and.nd.eq.1) then + initial_storage(res_no)=.TRUE. + end if go to 500 end if ! @@ -558,10 +561,6 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! call reservoir_subroutine_implicit(res_no,q_surf,nd,dbt(ncell)) end do - - !else - ! call reservoir_single_subroutine(res_no,q_surf,nd) - !end if ! T_0 = T_hypo(res_no) ! In reservoir, water is released from hypolimnion if (T_0.lt.0.5) T_0=0.5 diff --git a/src/flow_subroutine.f90 b/src/flow_subroutine.f90 index a7f0448..731d98b 100755 --- a/src/flow_subroutine.f90 +++ b/src/flow_subroutine.f90 @@ -23,13 +23,12 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) ! ! if we need to adjust the timestep, we need to make the volume ! go back to the beginning of the timestep. - ! + ! if (recalculate_volume) then volume_e_x(res_no) = volume_e_x(res_no) - vol_change_epi_x volume_h_x(res_no) = volume_h_x(res_no) - vol_change_hyp_x recalculate_volume=.FALSE. ! This only needs to be done once end if - !res_storage_pre = res_capacity_mcm(res_no) * (10**6) * 0.95 ! TESTINFLOW !res_storage_post = res_capacity_mcm(res_no) * (10**6) * 0.95 ! TESTINFLOW @@ -52,7 +51,8 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) ! if (initial_storage(res_no)) then water_withdrawal(res_no) = 0 - res_depth_meter(res_no) = res_storage(res_no,1) / surface_area(res_no) + !res_depth_meter(res_no) = res_storage(res_no,1) / surface_area(res_no) + res_depth_meter(res_no) = res_capacity_mcm(res_no) * 0.95 * 1e6 / surface_area(res_no) depth_e(res_no) = res_depth_meter(res_no) * depth_e_frac depth_h(res_no) = res_depth_meter(res_no) * depth_h_frac volume_e_x(res_no) = surface_area(res_no) * depth_e(res_no) From 5dcc50d873c190019cf773dfb56fe71b263fb0d6 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Wed, 14 Mar 2018 14:40:08 -0700 Subject: [PATCH 24/27] 1) calculate equilibrium temperature 2) estimate how long it will take to reach equilibrium --- src/Block_Energy.f90 | 1 + src/Energy.f90 | 75 ++++++++++++++++++++++++++++++++++++++++++-- src/Systmm.f90 | 5 ++- 3 files changed, 77 insertions(+), 4 deletions(-) diff --git a/src/Block_Energy.f90 b/src/Block_Energy.f90 index 15a2e3a..903d599 100755 --- a/src/Block_Energy.f90 +++ b/src/Block_Energy.f90 @@ -34,6 +34,7 @@ module Block_Energy real :: lvp,rb,rho real :: deriv_2nd,error_EE real :: deriv_conv,deriv_evap,deriv_ws + real :: temp_equil,time_equil real,parameter :: evap_coeff=1.5e-9 !Lake Hefner coefficient, 1/meters real,parameter :: pf=0.640,pi=3.14159 real,parameter :: rfac=304.8 !rho/Cp kg/meter**3/Kilocalories/kg/Deg K diff --git a/src/Energy.f90 b/src/Energy.f90 index 8d4b881..b6c44a3 100755 --- a/src/Energy.f90 +++ b/src/Energy.f90 @@ -1,11 +1,15 @@ SUBROUTINE Energy(T_surf,q_surf,ncell,z,nd) use Block_Energy + use Block_Network implicit none integer::i,ncell,nd + integer::iter,num_timestep,tot_timestep real::A,B,e0,q_surf,q_conv,q_evap,q_ws,T_surf real::const1,const2,z real::deriv_1st !,deriv_conv,deriv_evap,deriv_ws real, dimension(2):: q_fit, T_fit + real::q_tot,q_deriv,temp_calc + real::rate_temp,error_estimate,timestep ! T_fit(1)=T_surf-1.0 T_fit(2)=T_surf+1.0 @@ -19,9 +23,9 @@ SUBROUTINE Energy(T_surf,q_surf,ncell,z,nd) q_evap=q_evap*(e0-ea(ncell)) q_ws=6.693E-2+1.471E-3*T_fit(i) q_fit(i)=q_ns(ncell)+q_na(ncell)-q_ws-q_evap+q_conv - !if (ncell.eq.4480) write(*,*) & - ! 'energy',q_ns(ncell),q_na(ncell),-q_ws,-q_evap,q_conv, & - ! 'air',dbt(ncell),wind(ncell) +! if (ncell.eq.1827.and.nd.eq.197) write(*,*) & +! 'energy',q_ns(ncell),q_na(ncell),-q_ws,-q_evap,q_conv, & +! 'air',dbt(ncell),wind(ncell) end do ! ! q=AT+B @@ -34,6 +38,7 @@ SUBROUTINE Energy(T_surf,q_surf,ncell,z,nd) q_surf=0.5*(q_fit(1)+q_fit(2)) B=(q_surf/A)-(T_fit(1)+T_fit(2))/2. deriv_1st=(q_surf/(z*rfac)) +! if (ncell.eq.1827) write(*,*) nd,A,B+T_surf,B,dbt(ncell) ! ! ****************************************************** ! Return to Subroutine RIVMOD @@ -49,5 +54,69 @@ SUBROUTINE Energy(T_surf,q_surf,ncell,z,nd) deriv_ws=1.471E-3 deriv_2nd=deriv_1st*(deriv_conv-deriv_evap-deriv_ws)/(z*rfac) +! +! ****************************************************** +! Calculate equilibrium temperature +! ****************************************************** +! +if (ncell.eq.1827) then + iter=1 + temp_equil=T_surf + q_tot=q_surf + q_deriv=deriv_conv-deriv_evap-deriv_ws + temp_equil=0 + do while (iter.lt.50 .and. abs(q_tot/q_deriv).gt.0.001) + e0=2.1718E8*EXP(-4157.0/(temp_equil+239.09)) + rb=pf*(dbt(ncell)-temp_equil) + lvp=597.0-0.57*temp_equil + q_evap=1000.*lvp*evap_coeff*wind(ncell) + if(q_evap.lt.0.0) q_evap=0.0 + q_conv=rb*q_evap + q_evap=q_evap*(e0-ea(ncell)) + q_ws=6.693E-2+1.471E-3*temp_equil + q_tot=(q_ns(ncell)+q_na(ncell)-q_ws-q_evap+q_conv) + + rate_temp=q_tot/(z*rfac) + + const1=1000.*evap_coeff*wind(ncell)*pf + const2=1000.*evap_coeff*wind(ncell) + e0=2.1718E8*EXP(-4157.0/(temp_equil+239.09)) + deriv_conv=const1*(-(597.0-0.57*temp_equil) -0.57*(dbt(ncell)-temp_equil)) + deriv_evap=const2*(-0.57*(e0-ea(ncell))+(597.0-0.57*temp_equil)*e0*(4157.0/((temp_equil+239.09)**2))) + deriv_ws=1.471E-3 + q_deriv=deriv_conv-deriv_evap-deriv_ws + temp_equil=temp_equil-(q_tot/q_deriv) + iter=iter+1 + enddo + ! + ! Estimate the time to reach equilibrium tmeperature + ! + error_estimate=deriv_2nd*dt_comp**2/2 + num_timestep=max(1,ceiling(abs(error_estimate/0.1))) + timestep=dt_comp/num_timestep + tot_timestep=0 + temp_calc=T_surf + do while (tot_timestep.le.(50*num_timestep).and.abs(temp_equil-temp_calc).ge.0.01) + e0=2.1718E8*EXP(-4157.0/(temp_calc+239.09)) + rb=pf*(dbt(ncell)-temp_calc) + lvp=597.0-0.57*temp_calc + q_evap=1000.*lvp*evap_coeff*wind(ncell) + if(q_evap.lt.0.0) q_evap=0.0 + q_conv=rb*q_evap + q_evap=q_evap*(e0-ea(ncell)) + q_ws=6.693E-2+1.471E-3*temp_calc + q_tot=(q_ns(ncell)+q_na(ncell)-q_ws-q_evap+q_conv) + rate_temp=q_tot/(z*rfac) + temp_calc=temp_calc+rate_temp*timestep + tot_timestep=tot_timestep+1 + if (nd.eq.197) then + write(88,*) 'sub',tot_timestep,temp_equil,temp_calc,rate_temp,timestep + end if + enddo + time_equil=tot_timestep*timestep + if (nd.eq.197) then + write(88,*) temp_equil,temp_calc,num_timestep,timestep,tot_timestep,error_estimate + end if +end if END Subroutine Energy diff --git a/src/Systmm.f90 b/src/Systmm.f90 index 7d0146d..b9e1441 100755 --- a/src/Systmm.f90 +++ b/src/Systmm.f90 @@ -305,6 +305,8 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! Set initial river storage (Initial storage for first grid cell is 0) ! do nm=no_dt(ns),1,-1 + if (ncell.eq.1827.and.nd.eq.197) write(88,*) & + nyear,nd,ns,nm dt_calc=dt_part(nm) z=depth(nncell) call energy(T_0,q_surf,nncell,z,nd) @@ -404,11 +406,12 @@ SUBROUTINE SYSTMM(temp_file,res_file,param_file) ! test output for a specific grid cell ! ! + if (ncell.eq.1827.and.nd.eq.197) write(*,*) nyear,nd,nr,ncell,ns,T_0 do nseg_temp=1,nseg_out_num if (nseg_out(nr,ncell,nseg_temp).eq.ns) then call WRITE(time,nd,nr,ncell,ns,T_0,T_head(nr),dbt(ncell), & Q_inflow_out, Q_outflow_out) - !if (ncell.eq.4089) write(*,*)nyear,nd,nr,ncell,ns,T_0 + if (ncell.eq.1827) write(89,*)nyear,nd,nr,ncell,ns,T_0,temp_equil,time_equil,deriv_2nd*dt_comp**2/2 end if end do end if From b7c132ed6af33daec3e1e3ba31a87906b035dfd3 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Thu, 5 Apr 2018 15:15:12 -0700 Subject: [PATCH 25/27] set the minimum river depth to be 5 feet --- src/Read_Forcing.f90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/Read_Forcing.f90 b/src/Read_Forcing.f90 index 9c28b8a..236d2fa 100755 --- a/src/Read_Forcing.f90 +++ b/src/Read_Forcing.f90 @@ -12,6 +12,7 @@ SUBROUTINE Read_Forcing real :: Q_avg,Q_dmmy real :: z_temp,w_temp real :: min_flow = 5.0 + real :: min_depth = 5.0 ! n1=1 n2=2 @@ -38,6 +39,15 @@ SUBROUTINE Read_Forcing u(no_heat) = min_flow/(z_temp*w_temp) depth(no_heat) = z_temp width(no_heat) = w_temp + Q_out(no_heat) = min_flow + end if + ! + ! If the river depth is smaller than minimum river depth, + ! recalculate water width based on river depth + ! + if (depth(no_heat).lt.min_depth) then + depth(no_heat) = min_depth + width(no_heat) = Q_out(no_heat)/(u(no_heat)*depth(no_heat)) end if ! if (Q_diff(no_heat) .gt. 0) write(*,*) 'Q_diff is not equal to 0', no_heat, Q_diff(no_heat) From 963319e1ee09ce9ed5d830327ae73ce74dda5b90 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Mon, 16 Apr 2018 14:17:01 -0700 Subject: [PATCH 26/27] bug fix for reservoir storage data --- src/flow_subroutine.f90 | 4 +++- src/reservoir_single_subroutine.f90 | 26 ++++++++++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) create mode 100755 src/reservoir_single_subroutine.f90 diff --git a/src/flow_subroutine.f90 b/src/flow_subroutine.f90 index 731d98b..37465cf 100755 --- a/src/flow_subroutine.f90 +++ b/src/flow_subroutine.f90 @@ -57,7 +57,9 @@ SUBROUTINE flow_subroutine (res_no, nyear, nd) depth_h(res_no) = res_depth_meter(res_no) * depth_h_frac volume_e_x(res_no) = surface_area(res_no) * depth_e(res_no) volume_h_x(res_no) = surface_area(res_no) * depth_h(res_no) - res_storage_pre = res_storage(res_no,1) + !res_storage_pre = res_storage(res_no,1) + res_storage_post = res_storage(res_no,1) + res_storage_pre = volume_e_x(res_no) + volume_h_x(res_no) initial_storage(res_no)=.FALSE. end if ! diff --git a/src/reservoir_single_subroutine.f90 b/src/reservoir_single_subroutine.f90 new file mode 100755 index 0000000..c45148d --- /dev/null +++ b/src/reservoir_single_subroutine.f90 @@ -0,0 +1,26 @@ +SUBROUTINE reservoir_single_subroutine(res_no,q_surf,nd) + use Block_Network + use Block_Reservoir + ! + implicit none + ! + real :: q_surf, volume_tot + real :: temp_epi, temp_hyp + integer :: nd, res_no + +temp_epi = (volume_h_x(res_no)*T_hypo(res_no) + volume_e_x(res_no)*T_epil(res_no) + & + (flow_in_epi_x + flow_in_hyp_x) * T_res_inflow(res_no) + & + (q_surf * dt_comp * surface_area(res_no)) / (density * heat_c_kcal)) / & + (volume_h_x(res_no) + volume_e_x(res_no) + flow_in_epi_x + flow_in_hyp_x) +temp_hyp = temp_epi + +T_epil(res_no) = temp_epi +T_hypo(res_no) = temp_hyp + +volume_tot = volume_e_x(res_no) + volume_h_x(res_no) - water_withdrawal(res_no) + +volume_e_x(res_no) = volume_tot / 2.0 +volume_h_x(res_no) = volume_tot / 2.0 + +end subroutine reservoir_single_subroutine + From 7c385a3e2733a5ff19987b2cadd6d8590ff4e936 Mon Sep 17 00:00:00 2001 From: "Y. Cheng" Date: Wed, 25 Apr 2018 16:27:43 -0700 Subject: [PATCH 27/27] clean up makefile --- src/makefile | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/makefile b/src/makefile index 7da8270..162c93b 100644 --- a/src/makefile +++ b/src/makefile @@ -12,7 +12,6 @@ objects = rbm10_VIC.o Begin.o Systmm.o Particle_Track.o \ density_subroutine.o flow_subroutine.o \ Error_estimate.o \ reservoir_subroutine_implicit.o \ - reservoir_single_subroutine.o \ Water_Balance.o Write.o f90comp = gfortran #-O0 -g #-O2 # Makefile @@ -54,8 +53,6 @@ Error_estimate.o: block_network.mod block_reservoir.mod Error_estimate.f90 $(f90comp) -c Error_estimate.f90 reservoir_subroutine_implicit.o: block_network.mod block_reservoir.mod reservoir_subroutine_implicit.f90 $(f90comp) -c reservoir_subroutine_implicit.f90 -reservoir_single_subroutine.o: block_network.mod block_reservoir.mod reservoir_single_subroutine.f90 - $(f90comp) -c reservoir_single_subroutine.f90 Julian.o: Julian.f90 $(f90comp) -c Julian.f90 tntrp.o: tntrp.f90