subroutine gwflow_simulate !! ~ ~ ~ PURPOSE ~ ~ ~ !! this subroutine calculates new groundwater storage and solute mass for each gwflow grid cell; !! also, computes and write out daily/annual/average annual fluxes and mass balance error use gwflow_module use hydrograph_module use hru_module use sd_channel_module use time_module use soil_module use water_body_module, only: res_wat_d use maximum_data_module, only : db_mx use calibration_data_module, only : lsu_out implicit none !counters and general information integer :: i = 0 ! |counter integer :: j = 0 ! |counter integer :: k = 0 ! |counter integer :: n = 0 ! |counter integer :: s = 0 ! |solute counter integer :: dum = 0 ! |dummy variable integer :: num_ts = 0 ! |number of flow time steps during the daily time step integer :: cell_id = 0 ! |id of gwflow cell real :: sum = 0. ! |general summation !lateral flow calculations real :: area1 = 0. !m2 |spatial area of first connected cell real :: area2 = 0. !m2 |spatial area of second connected cell real :: area = 0. !m2 |smaller of the two areas real :: conn_length = 0. !m |length of connection between two adjacent cells real :: dist_x = 0. !m |x distance between centroids of adjacent cells real :: dist_y = 0. !m |y distance between centroids of adjacent cells real :: grad_distance = 0. !m |distance between centroids of two connected cells real :: Q_cell = 0. !m3 |groundater flow between two cells real :: face_K = 0. !m/day |hydraulic conductivity at cell interface real :: sat_thick1 = 0. !m |saturated thickness of connected cell real :: sat_thick2 = 0. !m |saturated thickness of cell real :: face_sat = 0. !m |saturated thickness at cell interface !storage calculations real :: stor_change = 0. !m3 |daily change in groundwater storage for a cell real :: sat_change = 0. !m |daily change in saturated thickness for a cell real :: flow_area = 0. !m2 |groundwater flow area between cells real :: Q = 0. !m3 |total flow in/out of cell real :: gradient = 0. !m/m |groundwater gradient between cells !groundwater conditions real :: frac_sat = 0. ! |fraction of cells that are saturated (i.e., water table at ground surface) real :: depth_wt_avg = 0. !m |average water table depth from all grid cells !water balance analysis real :: mass_error = 0. ! |mass error in groundwater balance and solute mass balance !tile drainage outflow real :: sum_tile(50) = 0. !m3 |summation of flow from tile cell groups real :: sum_mass(50,100) = 0. !g |total solute mass in tile cell groups real :: c_tile(50,100) = 0. !g/m3 |average concentration of solute mass in tile cell groups !grid totals of fluxes for day (m3, mm) real*8 :: vbef_grid = 0.d0 real*8 :: vaft_grid = 0.d0 real*8 :: rech_grid = 0.d0 real*8 :: gwet_grid = 0.d0 real*8 :: gwsw_grid = 0.d0 real*8 :: swgw_grid = 0.d0 real*8 :: satx_grid = 0.d0 real*8 :: soil_grid = 0.d0 real*8 :: latl_grid = 0.d0 real*8 :: bndr_grid = 0.d0 real*8 :: ppag_grid = 0.d0 real*8 :: ppdf_grid = 0.d0 real*8 :: ppex_grid = 0.d0 real*8 :: tile_grid = 0.d0 real*8 :: resv_grid = 0.d0 real*8 :: wetl_grid = 0.d0 real*8 :: canl_grid = 0.d0 real*8 :: fpln_grid = 0.d0 !groundwater solutes integer :: t = 0 ! |counter for solute transport time steps real :: gw_trans_time_step = 0. !days |length of solute transport time steps real :: time_fraction = 0. ! |fraction of flow time step real :: gw_volume_old = 0. !m3 |cell groundwater volume from previous flow time step real :: gw_volume_new = 0. !m3 |cell groundwater volume from current flow time step real :: gw_volume_inter = 0. !m3 |interpolated cell groundwater volume for the current transport time step real :: mass_adv(100) = 0. !g |solute mass advected into/out of cell real :: mass_dsp(100) = 0. !g |solute mass dispersed into/out of cell real :: m_change(100) = 0. !g |change in cell's solute mass, for one transport time step real :: del_no_sorp = 0. !g |change in cell's solute mass, without sorption real :: mass_sorb(100) = 0. !g |solute mass sorbed !total grid values for solute mass (kg) real :: sol_grid_mbef = 0. real :: sol_grid_maft = 0. real :: sol_grid_rech = 0. real :: sol_grid_gwsw = 0. real :: sol_grid_swgw = 0. real :: sol_grid_satx = 0. real :: sol_grid_advn = 0. real :: sol_grid_disp = 0. real :: sol_grid_rcti = 0. real :: sol_grid_rcto = 0. real :: sol_grid_minl = 0. real :: sol_grid_sorb = 0. real :: sol_grid_ppag = 0. real :: sol_grid_ppex = 0. real :: sol_grid_tile = 0. real :: sol_grid_soil = 0. real :: sol_grid_resv = 0. real :: sol_grid_wetl = 0. real :: sol_grid_canl = 0. real :: sol_grid_fpln = 0. !usgs observation wells real :: head_sum = 0. real :: head_avg = 0. real :: head_mae = 0. real :: error_sum = 0. real :: head_residual = 0. real :: error_sum_well = 0. real :: head_mae_well(1000,2) = 0. real :: sat_div_well(1000,2) = 0. integer :: write_yr = 0 integer :: usgs_yr = 0 integer :: val_count = 0 integer :: val_count_well = 0 integer :: num_yrs_calb = 0 integer :: num_yrs_test = 0 real :: head_mae_calb = 0. real :: head_mae_test = 0. real :: sat_sum = 0. real :: sat_avg = 0. real :: sum_sat = 0. real :: sum_sat_well = 0. real :: sat_div = 0. real :: sat_div_calb = 0. real :: sat_div_test = 0. integer :: num_gw_meas = 0 integer :: num_gw_meas_calb = 0 integer :: num_gw_meas_test = 0 integer :: num_gw_meas_well(500,2) = 0 !stream observations integer :: chan_count = 0 integer :: month_count = 0 integer :: month_count_calb = 0 integer :: month_count_test = 0 real :: sum_obs_flow = 0. real :: sum_sim_flow = 0. real :: mean_obs_flow = 0. real :: mean_sim_flow = 0. real :: sum_resi_nse = 0. real :: sum_diff_nse = 0. real :: sum_resi_nse1 = 0. real :: sum_diff_nse1 = 0. real :: sum_num = 0. real :: sum_den = 0. real :: sum_den1 = 0. real :: sum_den2 = 0. real :: cc = 0. real :: sum_sim_sd = 0. real :: sum_obs_sd = 0. real :: sim_sd = 0. real :: obs_sd = 0. integer :: num_months_calb = 0 integer :: num_months_test = 0 !record file write(out_gw,*) 'gwflow subroutine called:',time%yrc,time%day !1. calculate the available volume of groundwater (m3) in each cell --------------------------------------------- do i=1,ncell if(gw_state(i)%stat.gt.0) then if(gw_state(i)%head > gw_state(i)%botm) then gw_state(i)%stor = ((gw_state(i)%head - gw_state(i)%botm) * gw_state(i)%area) * gw_state(i)%spyd else gw_state(i)%stor = 0. endif endif enddo !2. calculate groundwater sources and sinks (m3) for each cell -------------------------------------------------- !recharge --------------------------------------------------------------- call gwflow_rech !groundwater ET --------------------------------------------------------- call gwflow_gwet !groundwater-channel exchange ------------------------------------------- !gwflow_gwsw called in sd_channel_control !groundwater saturation excess flow ------------------------------------- !gwflow_satx called in sd_channel_control !groundwater-->soil exchange -------------------------------------------- !gwflow_soil called in swr_percmain !groundwater pumping (irrigation) --------------------------------------- !gwflow_ppag called in wallo_withdraw !track and print out pumping for HRUs do i=1,sp_ob%hru hru_pump_mo(i) = hru_pump_mo(i) + hru_pump(i) hru_pump_yr(i) = hru_pump_yr(i) + hru_pump(i) enddo if (hru_pump_flag == 1) then !pumping output for specified HRUs do i=1,num_hru_pump_obs hru_pump_obs(i) = hru_pump(hru_pump_ids(i)) enddo write(out_hru_pump_obs,119) time%yrc,time%day, & (hru_pump_obs(i),i=1,num_hru_pump_obs) endif hru_pump = 0. !groundwater pumping (external use) ------------------------------------- call gwflow_ppex !discharge from groundwater to tile drains ------------------------------ !gwflow_tile called in sd_channel_control !retrieve information for tile cell groups if (gw_tile_flag == 1) then !computer flow rate and solute concentration for the tile cell groups if(gw_tile_group_flag == 1) then do i=1,gw_tile_num_group sum_tile(i) = 0. do j=1,num_tile_cells(i) sum_tile(i) = sum_tile(i) + gw_ss(gw_tile_groups(i,j))%tile !m3 enddo sum_tile(i) = (sum_tile(i)*(-1)) / 86400. !m3 --> m3/sec if (gw_solute_flag == 1) then do s=1,gw_nsolute !loop through the solutes sum_mass(i,s) = 0. enddo do j=1,num_tile_cells(i) do s=1,gw_nsolute !loop through the solutes sum_mass(i,s) = sum_mass(i,s) + & gwsol_ss(gw_tile_groups(i,j))%solute(s)%tile !g enddo enddo if(sum_tile(i).ne.0) then do s=1,gw_nsolute !loop through the solutes c_tile(i,s) = sum_mass(i,s) / sum_tile(i) !g/m3 enddo else do s=1,gw_nsolute !loop through the solutes c_tile(i,s) = 0. !g/m3 enddo endif endif enddo if (gw_solute_flag == 1) then write(out_tile_cells,102) time%day,time%yrc, & (sum_tile(i),i=1,gw_tile_num_group), & (c_tile(i,1),i=1,gw_tile_num_group), & (c_tile(i,2),i=1,gw_tile_num_group) !only no3 and p else write(out_tile_cells,102) time%day,time%yrc, & (sum_tile(i),i=1,gw_tile_num_group) endif endif endif !groundwater exchange with reservoirs ----------------------------------- !gwflow_resv called in res_control !groundwater-wetland exchange ------------------------------------------- !gwflow_wetl called in wetland_control !groundwater-channel exchange from floodplain cells --------------------- !gwflow_fpln called in sd_channel_control !seepage from irrigation canals to groundwater -------------------------- !gwflow_canl called in sd_channel_control !seepage from irrigation canals to groundwater -------------------------- !(for canals that originate outside the model boundary) call gwflow_canl_out !3. sum sources/sinks for each grid cell ------------------------------------------------------------------------ !m3 for water; g for solutes do i=1,ncell if(gw_state(i)%stat == 1) then gw_ss(i)%totl = gw_ss(i)%rech + gw_ss(i)%gwet + gw_ss(i)%gwsw + gw_ss(i)%swgw + & gw_ss(i)%satx + gw_ss(i)%soil + & gw_ss(i)%ppag + gw_ss(i)%ppex + gw_ss(i)%tile + & gw_ss(i)%resv + gw_ss(i)%wetl + gw_ss(i)%canl + gw_ss(i)%fpln endif enddo if (gw_solute_flag == 1) then do i=1,ncell if(gw_state(i)%stat == 1) then do s=1,gw_nsolute !loop through the solutes gwsol_ss(i)%solute(s)%totl = gwsol_ss(i)%solute(s)%rech + & !recharge gwsol_ss(i)%solute(s)%gwsw + & !gw-->channel gwsol_ss(i)%solute(s)%swgw + & !channel-->gw gwsol_ss(i)%solute(s)%satx + & !gw-->channel gwsol_ss(i)%solute(s)%soil + & !gw-->soil gwsol_ss(i)%solute(s)%ppag + & !pumping (irrigation) gwsol_ss(i)%solute(s)%ppex + & !pumping (external) gwsol_ss(i)%solute(s)%tile + & !tile outflow gwsol_ss(i)%solute(s)%resv + & !reservoir exchange gwsol_ss(i)%solute(s)%wetl + & !wetland gwsol_ss(i)%solute(s)%canl + & !canal exchange gwsol_ss(i)%solute(s)%fpln !floodplain exchange enddo endif enddo endif ! 4. calculate new groundwater storage and head for each grid cell ---------------------------------------------- !compute groundwater volume and solute mass at beginning of day do i=1,ncell if(gw_state(i)%stat == 1) then gw_state(i)%vbef = (gw_state(i)%head-gw_state(i)%botm) * gw_state(i)%area * gw_state(i)%spyd !m3 endif enddo !compute cell solute mass at the beginning of the day if (gw_solute_flag == 1) then do i=1,ncell if(gw_state(i)%stat == 1) then do s=1,gw_nsolute !loop through the solutes gwsol_state(i)%solute(s)%mbef = gwsol_state(i)%solute(s)%mass enddo endif enddo endif !determine number of flow time steps; determine size of transport time step num_ts = int(1./gw_time_step) if (gw_solute_flag == 1) then gw_trans_time_step = gw_time_step / num_ts_transport endif !prepare arrays do i=1,ncell gw_state(i)%hnew = 0. gw_state(i)%hold = 0. enddo if (gw_solute_flag == 1) then do i=1,ncell do s=1,gw_nsolute gwsol_state(i)%solute(s)%cnew = 0. enddo enddo endif !calculate new storage and head for each cell do n=1,num_ts do i=1,ncell !only proceed if the cell is active if(gw_state(i)%stat > 0) then !if the cell is interior (not a boundry cell) if(gw_state(i)%stat == 1) then !loop through the cells connected to the current cell Q = 0. do k=1,gw_state(i)%ncon !id of the connected cell cell_id = cell_con(i)%cell_id(k) Q_cell = 0. !calculate groundwater flow between the cells, using Darcy's Law if(gw_state(cell_id)%stat == 0) then Q_cell = 0. elseif(gw_state(cell_id)%stat == 2 .and. bc_type == 2) then !boundary cell Q_cell = 0. else !length of connection between the two cells area1 = gw_state(cell_id)%area !area of connected cell area2 = gw_state(i)%area !area of current cell area = min(area1,area2) !smaller of the two conn_length = sqrt(area) !K along the interface (harmonic mean) face_K = conn_length / (((conn_length/2.)/gw_state(cell_id)%hydc) + ((conn_length/2.)/gw_state(i)%hydc)) !saturated thickness of connected cell if(gw_state(cell_id)%head > gw_state(cell_id)%botm) then sat_thick1 = gw_state(cell_id)%head - gw_state(cell_id)%botm else sat_thick1 = 0. endif !saturated thickness of current cell if(gw_state(i)%head > gw_state(i)%botm) then sat_thick2 = gw_state(i)%head - gw_state(i)%botm else sat_thick2 = 0. endif !saturated thickness at the interface (m) face_sat = (sat_thick1 + sat_thick2) / 2. !groundwater hydraulic gradient (m/m) dist_x = gw_state(i)%xcrd - gw_state(cell_id)%xcrd !m dist_y = gw_state(i)%ycrd - gw_state(cell_id)%ycrd !m grad_distance = sqrt((dist_x)**2 + (dist_y)**2) gradient = (gw_state(cell_id)%head - gw_state(i)%head) / grad_distance !groundwater flow area (m2) flow_area = face_sat * conn_length !groundwater flow rate (m3/day) Q_cell = face_K * gradient * flow_area !Darcy's Law !store for solute transport mass balance calculations cell_con(i)%latl(k) = Q_cell cell_con(i)%sat(k) = face_sat endif !store for cell water balance if(gw_state(cell_id)%stat == 2) then !boundary flow gw_ss(i)%bndr = gw_ss(i)%bndr + (Q_cell*gw_time_step) else gw_ss(i)%latl = gw_ss(i)%latl + (Q_cell*gw_time_step) gw_ss_sum(i)%latl = gw_ss_sum(i)%latl + (Q_cell*gw_time_step) endif !sum total flow to/from current cell Q = Q + Q_cell enddo !go to next connected cell !update storage and head for the cell stor_change = (Q + gw_ss(i)%totl) * gw_time_step !change in storage (m3) gw_state(i)%stor = gw_state(i)%stor + stor_change !new storage (m3) sat_change = stor_change / (gw_state(i)%spyd * gw_state(i)%area) !change in saturated thickness (m3) gw_state(i)%hnew = gw_state(i)%head + sat_change !new groundwater head (m) elseif(gw_state(i)%stat == 2) then !constant head cell gw_state(i)%hnew = gw_state(i)%init gw_state(i)%stor = ((gw_state(i)%hnew - gw_state(i)%botm) * gw_state(i)%area) * gw_state(i)%spyd endif endif !check if cell is active enddo !go to next cell !store new head values into regular head array do i=1,ncell gw_state(i)%hold = gw_state(i)%head gw_state(i)%head = gw_state(i)%hnew enddo !simulate fate and transport of solutes - calculate new concentrations if (gw_solute_flag == 1) then do t=1,num_ts_transport do i=1,ncell if(gw_state(i)%stat == 1) then !interior cell !calculate old and new groundwater volume in the cell (m3) if(gw_state(i)%hold > gw_state(i)%botm) then gw_volume_old = gw_state(i)%area * (gw_state(i)%hold - gw_state(i)%botm) * gw_state(i)%spyd else gw_volume_old = 0. endif if(gw_state(i)%head > gw_state(i)%botm) then gw_volume_new = gw_state(i)%area * (gw_state(i)%head - gw_state(i)%botm) * gw_state(i)%spyd else gw_volume_new = 0. endif !calculate groundwater volume for the current transport time step (via interpolation) time_fraction = real(t)/real(num_ts_transport) gw_volume_inter = gw_volume_old + ((gw_volume_new-gw_volume_old)*time_fraction) !advection transport mass_adv = 0. do k=1,gw_state(i)%ncon !loop through the connected cells cell_id = cell_con(i)%cell_id(k) !id of the connected cell Q_cell = cell_con(i)%latl(k) !m3/day flow between cell and connected cell if(Q_cell > 0) then !mass entering cell do s=1,gw_nsolute mass_adv(s) = mass_adv(s) + (Q_cell * gwsol_state(cell_id)%solute(s)%conc) !g enddo else !mass leaving cell do s=1,gw_nsolute mass_adv(s) = mass_adv(s) + (Q_cell * gwsol_state(i)%solute(s)%conc) !g enddo endif enddo !go to next connected cell !dispersion transport mass_dsp = 0. do k=1,gw_state(i)%ncon !loop through the connected cells cell_id = cell_con(i)%cell_id(k) !id of the connected cell face_sat = cell_con(i)%sat(k) !m saturated thickness at interface between cells area1 = gw_state(cell_id)%area !area of connected cell area2 = gw_state(i)%area !area of current cell area = min(area1,area2) !smaller of the two conn_length = sqrt(area) do s=1,gw_nsolute !loop through the solutes mass_dsp(s) = mass_dsp(s) + (gw_long_disp * ((gwsol_state(cell_id)%solute(s)%conc - & gwsol_state(i)%solute(s)%conc)/conn_length) * face_sat) !g enddo enddo !go to next connected cell !chemical reactions --> fill in mass_rct and mass_min cell_id = i mass_rct = 0. mass_min = 0. call gwflow_chem(cell_id,gw_volume_inter) !calculate change in mass (g) do s=1,gw_nsolute !loop through the solutes m_change(s) = (mass_adv(s) + mass_dsp(s) + mass_rct(s) + mass_min(s) + gwsol_ss(i)%solute(s)%totl) * & (gw_trans_time_step/gwsol_sorb(s)) enddo !calculate mass removed due to sorption (g) do s=1,gw_nsolute !loop through the solutes del_no_sorp = (mass_adv(s) + mass_dsp(s) + mass_rct(s) + mass_min(s) + gwsol_ss(i)%solute(s)%totl) * & gw_trans_time_step mass_sorb(s) = del_no_sorp - m_change(s) enddo !calculate new mass in the cell (g) do s=1,gw_nsolute !loop through the solutes gwsol_state(i)%solute(s)%mass = gwsol_state(i)%solute(s)%mass + m_change(s) if(gwsol_state(i)%solute(s)%mass < 0) then gwsol_state(i)%solute(s)%mass = 0. endif enddo !calculate new concentration (g/m3) if(gw_volume_inter > 0) then do s=1,gw_nsolute !loop through the solutes gwsol_state(i)%solute(s)%cnew = gwsol_state(i)%solute(s)%mass / gw_volume_inter enddo else do s=1,gw_nsolute !loop through the solutes gwsol_state(i)%solute(s)%cnew = 0. gwsol_state(i)%solute(s)%mass = 0. enddo endif !store values for mass budget analysis do s=1,gw_nsolute !loop through the solutes gwsol_ss(i)%solute(s)%advn = gwsol_ss(i)%solute(s)%advn + (mass_adv(s)*(gw_trans_time_step/gwsol_sorb(s))) gwsol_ss(i)%solute(s)%disp = gwsol_ss(i)%solute(s)%disp + (mass_dsp(s)*(gw_trans_time_step/gwsol_sorb(s))) if(mass_rct(s) > 0) then gwsol_ss(i)%solute(s)%rcti = gwsol_ss(i)%solute(s)%rcti + (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !produced else gwsol_ss(i)%solute(s)%rcto = gwsol_ss(i)%solute(s)%rcto + (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !consumed endif gwsol_ss(i)%solute(s)%minl = gwsol_ss(i)%solute(s)%minl + (mass_min(s)*(gw_trans_time_step/gwsol_sorb(s))) gwsol_ss(i)%solute(s)%sorb = gwsol_ss(i)%solute(s)%sorb + mass_sorb(s) if(mass_sorb(3) > 0) then dum = 10 endif enddo !track for annual write-out do s=1,gw_nsolute !loop through the solutes if(mass_rct(s) > 0) then gwsol_ss_sum(i)%solute(s)%rcti = gwsol_ss_sum(i)%solute(s)%rcti + & (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !produced else gwsol_ss_sum(i)%solute(s)%rcto = gwsol_ss_sum(i)%solute(s)%rcto + & (mass_rct(s)*(gw_trans_time_step/gwsol_sorb(s))) !consumed endif gwsol_ss_sum(i)%solute(s)%minl = gwsol_ss_sum(i)%solute(s)%minl + & (mass_min(s)*(gw_trans_time_step/gwsol_sorb(s))) gwsol_ss_sum(i)%solute(s)%sorb = gwsol_ss_sum(i)%solute(s)%sorb + mass_sorb(s) enddo elseif(gw_state(i)%stat == 2) then !constant concentration cell do s=1,gw_nsolute gwsol_state(i)%solute(s)%cnew = 0. enddo endif enddo !go to next cell !store new concentration values into regular array do i=1,ncell do s=1,gw_nsolute !loop through the solutes gwsol_state(i)%solute(s)%conc = gwsol_state(i)%solute(s)%cnew enddo enddo enddo !go to next transport time step endif !check if solute transport is being simulated enddo !next flow time step -------------------------------------------------------------------------------------- !5. save/write out head and solute concentrations --------------------------------------------------------------- !save head values and solute concentration values for monthly and annual averages do i=1,ncell gw_state(i)%hdmo = gw_state(i)%hdmo + gw_state(i)%head gw_state(i)%hdyr = gw_state(i)%hdyr + gw_state(i)%head enddo if (gw_solute_flag == 1) then do i=1,ncell do s=1,gw_nsolute !loop through the solutes gwsol_state(i)%solute(s)%cnmo = gwsol_state(i)%solute(s)%cnmo + gwsol_state(i)%solute(s)%conc gwsol_state(i)%solute(s)%cnyr = gwsol_state(i)%solute(s)%cnyr + gwsol_state(i)%solute(s)%conc enddo enddo endif !open(96622,file='gwflow_head_daily') !if(grid_type == "structured") then ! grid_val = 0. ! do i=1,grid_nrow ! do j=1,grid_ncol ! if(cell_id_usg(i,j) > 0) then ! grid_val(i,j) = gw_state(cell_id_usg(i,j))%head ! endif ! enddo ! enddo ! do i=1,grid_nrow ! write(96622,100) (grid_val(i,j),j=1,grid_ncol) ! enddo ! endif !write(96622,*) !print out new head values and solute concentration values, if requested if(gw_output_index.le.gw_num_output) then if(gw_output_yr(gw_output_index).eq.time%yrc .and. gw_output_day(gw_output_index).eq.time%day) then write(out_gwheads,*) 'Groundwater Head for:',time%yrc,time%day if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_state(cell_id_usg(i,j))%head endif enddo enddo do i=1,grid_nrow write(out_gwheads,100) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gwheads,120) (gw_state(i)%head,i=1,ncell) endif write(out_gwheads,*) if (gw_solute_flag == 1) then do s=1,gw_nsolute !loop through the solutes write(out_gwconc,*) gwsol_nm(s),'concentration for:',time%yrc,time%day if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_state(cell_id_usg(i,j))%solute(s)%conc endif enddo enddo do i=1,grid_nrow write(out_gwconc,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gwconc,120) (gwsol_state(i)%solute(s)%conc,i=1,ncell) endif write(out_gwconc,*) enddo endif gw_output_index = gw_output_index + 1 endif endif !calculate the average depth to water table sum = 0. do i=1,ncell if(gw_state(i)%stat == 1) then sum = sum + (gw_state(i)%elev - gw_state(i)%head) endif enddo depth_wt_avg = sum / num_active !print out head values and solute concentration values for observation cells (each time step) do k=1,gw_num_obs_wells gw_obs_head(k) = gw_state(gw_obs_cells(k))%head if (gw_solute_flag == 1) then do s=1,gw_nsolute !loop through the solutes gw_obs_solute(k,s) = gwsol_state(gw_obs_cells(k))%solute(s)%conc enddo endif !is usgs wells, store for end-of-year analysis if(usgs_obs == 1) then gw_obs_head_annual(k,time%day) = gw_state(gw_obs_cells(k))%head gw_obs_sat_annual(k,time%day) = gw_state(gw_obs_cells(k))%head - & gw_state(gw_obs_cells(k))%botm endif enddo write(out_gwobs,119) time%yrc,time%day,(gw_obs_head(k),k=1,gw_num_obs_wells) if (gw_solute_flag == 1) then write(out_gwobs_sol,119) time%yrc,time%day,(gw_obs_solute(k,1),k=1,gw_num_obs_wells), & (gw_obs_solute(k,2),k=1,gw_num_obs_wells) !need to continue if there are more solutes... endif !if the end of the month has been reached, then store flow rates for specified channels if (stream_obs == 1) then if(gw_num_obs_chan.gt.0) then if (time%end_mo == 1) then do chan_count=1,gw_num_obs_chan sim_flow_vals(chan_count,sim_month) = ch_out_m(obs_channels(chan_count))%flo / time%day_mo enddo sim_month = sim_month + 1 endif endif endif !compute groundwater volumes at the end of the day do i=1,ncell if(gw_state(i)%stat == 1) then gw_state(i)%vaft = ((gw_state(i)%head - gw_state(i)%botm) * gw_state(i)%area) * gw_state(i)%spyd endif enddo !compute solute mass at the end of the day if (gw_solute_flag == 1) then do i=1,ncell if(gw_state(i)%stat == 1) then do s=1,gw_nsolute !loop through the solutes gwsol_state(i)%solute(s)%maft = gwsol_state(i)%solute(s)%mass !g enddo endif enddo endif !print out source/sink information for a specified cell !gw_cell_obs_ss_vals(1) = gw_state(gw_cell_obs_ss)%head !gw_cell_obs_ss_vals(2) = gw_state(gw_cell_obs_ss)%vbef !gw_cell_obs_ss_vals(3) = gw_state(gw_cell_obs_ss)%vaft !gw_cell_obs_ss_vals(4) = gw_ss(gw_cell_obs_ss)%rech !gw_cell_obs_ss_vals(5) = gw_ss(gw_cell_obs_ss)%gwet !gw_cell_obs_ss_vals(6) = gw_ss(gw_cell_obs_ss)%gwsw !gw_cell_obs_ss_vals(7) = gw_ss(gw_cell_obs_ss)%swgw !gw_cell_obs_ss_vals(8) = gw_ss(gw_cell_obs_ss)%satx !gw_cell_obs_ss_vals(9) = gw_ss(gw_cell_obs_ss)%soil !gw_cell_obs_ss_vals(10) = gw_ss(gw_cell_obs_ss)%latl !gw_cell_obs_ss_vals(11) = gw_ss(gw_cell_obs_ss)%ppag !gw_cell_obs_ss_vals(12) = gw_ss(gw_cell_obs_ss)%ppex !gw_cell_obs_ss_vals(13) = gw_ss(gw_cell_obs_ss)%tile !gw_cell_obs_ss_vals(14) = gw_ss(gw_cell_obs_ss)%resv !gw_cell_obs_ss_vals(15) = gw_ss(gw_cell_obs_ss)%wetl !gw_cell_obs_ss_vals(16) = gw_ss(gw_cell_obs_ss)%canl !gw_cell_obs_ss_vals(17) = gw_ss(gw_cell_obs_ss)%fpln !write(out_gwobs_ss,102) time%yrc,time%day,(gw_cell_obs_ss_vals(i),i=1,17) !sum groundwater budget terms for each HUC12 (if national model mode) if (nat_model == 1) then do n=1,sp_ob%outlet !loop through the HUC12 catchments do k=1,huc12_ncell(n) !loop through the cells within each HUC12 catchment cell_id = huc12_cells(n,k) if(gw_state(cell_id)%stat.eq.1) then gw_huc12_wb(1,n) = gw_huc12_wb(1,n) + gw_ss(cell_id)%rech gw_huc12_wb(2,n) = gw_huc12_wb(2,n) + gw_ss(cell_id)%gwet gw_huc12_wb(3,n) = gw_huc12_wb(3,n) + gw_ss(cell_id)%gwsw gw_huc12_wb(4,n) = gw_huc12_wb(4,n) + gw_ss(cell_id)%swgw gw_huc12_wb(5,n) = gw_huc12_wb(5,n) + gw_ss(cell_id)%satx gw_huc12_wb(6,n) = gw_huc12_wb(6,n) + gw_ss(cell_id)%soil gw_huc12_wb(7,n) = gw_huc12_wb(7,n) + gw_ss(cell_id)%latl gw_huc12_wb(8,n) = gw_huc12_wb(8,n) + gw_ss(cell_id)%ppag gw_huc12_wb(9,n) = gw_huc12_wb(9,n) + gw_ss(cell_id)%ppex gw_huc12_wb(10,n) = gw_huc12_wb(10,n) + gw_ss(cell_id)%tile gw_huc12_wb(11,n) = gw_huc12_wb(11,n) + gw_ss(cell_id)%resv gw_huc12_wb(12,n) = gw_huc12_wb(12,n) + gw_ss(cell_id)%wetl gw_huc12_wb(13,n) = gw_huc12_wb(13,n) + gw_ss(cell_id)%canl gw_huc12_wb(14,n) = gw_huc12_wb(14,n) + gw_ss(cell_id)%fpln gw_huc12_wb(15,n) = gw_huc12_wb(15,n) + gw_ss(cell_id)%ppdf gw_huc12_wb_mo(1,n) = gw_huc12_wb_mo(1,n) + gw_ss(cell_id)%rech gw_huc12_wb_mo(2,n) = gw_huc12_wb_mo(2,n) + gw_ss(cell_id)%gwet gw_huc12_wb_mo(3,n) = gw_huc12_wb_mo(3,n) + gw_ss(cell_id)%gwsw gw_huc12_wb_mo(4,n) = gw_huc12_wb_mo(4,n) + gw_ss(cell_id)%swgw gw_huc12_wb_mo(5,n) = gw_huc12_wb_mo(5,n) + gw_ss(cell_id)%satx gw_huc12_wb_mo(6,n) = gw_huc12_wb_mo(6,n) + gw_ss(cell_id)%soil gw_huc12_wb_mo(7,n) = gw_huc12_wb_mo(7,n) + gw_ss(cell_id)%latl gw_huc12_wb_mo(8,n) = gw_huc12_wb_mo(8,n) + gw_ss(cell_id)%ppag gw_huc12_wb_mo(9,n) = gw_huc12_wb_mo(9,n) + gw_ss(cell_id)%ppex gw_huc12_wb_mo(10,n) = gw_huc12_wb_mo(10,n) + gw_ss(cell_id)%tile gw_huc12_wb_mo(11,n) = gw_huc12_wb_mo(11,n) + gw_ss(cell_id)%resv gw_huc12_wb_mo(12,n) = gw_huc12_wb_mo(12,n) + gw_ss(cell_id)%wetl gw_huc12_wb_mo(13,n) = gw_huc12_wb_mo(13,n) + gw_ss(cell_id)%canl gw_huc12_wb_mo(14,n) = gw_huc12_wb_mo(14,n) + gw_ss(cell_id)%fpln gw_huc12_wb_mo(15,n) = gw_huc12_wb_mo(15,n) + gw_ss(cell_id)%ppdf endif enddo enddo endif !6. compute and write out daily fluxes and mass balance error !calculate values for entire grid (all cells) vbef_grid = 0. vaft_grid = 0. rech_grid = 0. gwet_grid = 0. gwsw_grid = 0. swgw_grid = 0. satx_grid = 0. soil_grid = 0. latl_grid = 0. bndr_grid = 0. ppag_grid = 0. ppdf_grid = 0. ppex_grid = 0. tile_grid = 0. resv_grid = 0. wetl_grid = 0. canl_grid = 0. fpln_grid = 0. !open(1356,file='gwflow_daily_cell_balance') do i=1,ncell if(gw_state(i)%stat == 1) then vbef_grid = vbef_grid + gw_state(i)%vbef vaft_grid = vaft_grid + gw_state(i)%vaft rech_grid = rech_grid + gw_ss(i)%rech gwet_grid = gwet_grid + gw_ss(i)%gwet gwsw_grid = gwsw_grid + gw_ss(i)%gwsw swgw_grid = swgw_grid + gw_ss(i)%swgw satx_grid = satx_grid + gw_ss(i)%satx soil_grid = soil_grid + gw_ss(i)%soil latl_grid = latl_grid + gw_ss(i)%latl bndr_grid = bndr_grid + gw_ss(i)%bndr ppag_grid = ppag_grid + gw_ss(i)%ppag ppdf_grid = ppdf_grid + gw_ss(i)%ppdf ppex_grid = ppex_grid + gw_ss(i)%ppex tile_grid = tile_grid + gw_ss(i)%tile resv_grid = resv_grid + gw_ss(i)%resv wetl_grid = wetl_grid + gw_ss(i)%wetl canl_grid = canl_grid + gw_ss(i)%canl fpln_grid = fpln_grid + gw_ss(i)%fpln !write(1356,100) gw_state(i)%vbef,gw_state(i)%vaft,gw_ss(i)%rech,gw_ss(i)%gwet,gw_ss(i)%gwsw, & ! gw_ss(i)%swgw,gw_ss(i)%satx,gw_ss(i)%soil, & ! gw_ss(i)%latl,gw_ss(i)%bndr,gw_ss(i)%ppag, & ! gw_ss(i)%ppdf,gw_ss(i)%ppex,gw_ss(i)%tile, & ! gw_ss(i)%resv,gw_ss(i)%wetl,gw_ss(i)%canl, & ! gw_ss(i)%fpln endif enddo !write(1356,*) mass_error = (1-((vbef_grid + rech_grid + gwet_grid + gwsw_grid + swgw_grid - satx_grid - soil_grid + & latl_grid + bndr_grid + ppag_grid + ppex_grid + tile_grid + resv_grid + wetl_grid + & canl_grid + fpln_grid) & /vaft_grid)) * 100 !print out daily information !first, normalize volumes to the watershed area (m3 --> mm) vbef_grid = (vbef_grid / (bsn%area_tot_ha*10000.)) * 1000. vaft_grid = (vaft_grid / (bsn%area_tot_ha*10000.)) * 1000. rech_grid = (rech_grid / (bsn%area_tot_ha*10000.)) * 1000. gwet_grid = (gwet_grid / (bsn%area_tot_ha*10000.)) * 1000. gwsw_grid = (gwsw_grid / (bsn%area_tot_ha*10000.)) * 1000. swgw_grid = (swgw_grid / (bsn%area_tot_ha*10000.)) * 1000. satx_grid = (satx_grid / (bsn%area_tot_ha*10000.)) * 1000. soil_grid = (soil_grid / (bsn%area_tot_ha*10000.)) * 1000. latl_grid = (latl_grid / (bsn%area_tot_ha*10000.)) * 1000. bndr_grid = (bndr_grid / (bsn%area_tot_ha*10000.)) * 1000. ppag_grid = (ppag_grid / (bsn%area_tot_ha*10000.)) * 1000. ppdf_grid = (ppdf_grid / (bsn%area_tot_ha*10000.)) * 1000. ppex_grid = (ppex_grid / (bsn%area_tot_ha*10000.)) * 1000. tile_grid = (tile_grid / (bsn%area_tot_ha*10000.)) * 1000. resv_grid = (resv_grid / (bsn%area_tot_ha*10000.)) * 1000. wetl_grid = (wetl_grid / (bsn%area_tot_ha*10000.)) * 1000. canl_grid = (canl_grid / (bsn%area_tot_ha*10000.)) * 1000. fpln_grid = (fpln_grid / (bsn%area_tot_ha*10000.)) * 1000. frac_sat = real(satx_count) / real(num_active) !calculate the fraction of active grid cells that are fully saturated if(gwflag_day.eq.1) then write(out_gwbal,102) time%yrc,time%day,gw_time_step,vbef_grid,vaft_grid,rech_grid,gwet_grid,gwsw_grid,swgw_grid, & satx_grid,soil_grid,latl_grid,bndr_grid,ppag_grid,ppex_grid, & tile_grid,resv_grid,wetl_grid,canl_grid,fpln_grid, & mass_error,frac_sat,depth_wt_avg,ppdf_grid endif !add daily water balance volumes to yearly values ss_grid_yr%chng = ss_grid_yr%chng + (vaft_grid-vbef_grid) ss_grid_yr%rech = ss_grid_yr%rech + rech_grid ss_grid_yr%gwet = ss_grid_yr%gwet + gwet_grid ss_grid_yr%gwsw = ss_grid_yr%gwsw + gwsw_grid ss_grid_yr%swgw = ss_grid_yr%swgw + swgw_grid ss_grid_yr%satx = ss_grid_yr%satx + satx_grid ss_grid_yr%soil = ss_grid_yr%soil + soil_grid ss_grid_yr%latl = ss_grid_yr%latl + latl_grid ss_grid_yr%bndr = ss_grid_yr%bndr + bndr_grid ss_grid_yr%ppag = ss_grid_yr%ppag + ppag_grid ss_grid_yr%ppdf = ss_grid_yr%ppdf + ppdf_grid ss_grid_yr%ppex = ss_grid_yr%ppex + ppex_grid ss_grid_yr%tile = ss_grid_yr%tile + tile_grid ss_grid_yr%resv = ss_grid_yr%resv + resv_grid ss_grid_yr%wetl = ss_grid_yr%wetl + wetl_grid ss_grid_yr%canl = ss_grid_yr%canl + canl_grid ss_grid_yr%fpln = ss_grid_yr%fpln + fpln_grid !add daily water balance volumes to total values ss_grid_tt%chng = ss_grid_tt%chng + (vaft_grid-vbef_grid) ss_grid_tt%rech = ss_grid_tt%rech + rech_grid ss_grid_tt%gwet = ss_grid_tt%gwet + gwet_grid ss_grid_tt%gwsw = ss_grid_tt%gwsw + gwsw_grid ss_grid_tt%swgw = ss_grid_tt%swgw + swgw_grid ss_grid_tt%satx = ss_grid_tt%satx + satx_grid ss_grid_tt%soil = ss_grid_tt%soil + soil_grid ss_grid_tt%latl = ss_grid_tt%latl + latl_grid ss_grid_tt%bndr = ss_grid_tt%bndr + bndr_grid ss_grid_tt%ppag = ss_grid_tt%ppag + ppag_grid ss_grid_tt%ppdf = ss_grid_tt%ppdf + ppdf_grid ss_grid_tt%ppex = ss_grid_tt%ppex + ppex_grid ss_grid_tt%tile = ss_grid_tt%tile + tile_grid ss_grid_tt%resv = ss_grid_tt%resv + resv_grid ss_grid_tt%wetl = ss_grid_tt%wetl + wetl_grid ss_grid_tt%canl = ss_grid_tt%canl + canl_grid ss_grid_tt%fpln = ss_grid_tt%fpln + fpln_grid if (gw_solute_flag == 1) then !if solutes are simulated !loop through the solutes do s=1,gw_nsolute sol_grid_mbef = 0. sol_grid_maft = 0. sol_grid_rech = 0. sol_grid_gwsw = 0. sol_grid_swgw = 0. sol_grid_satx = 0. sol_grid_soil = 0. sol_grid_advn = 0. sol_grid_disp = 0. sol_grid_rcti = 0. sol_grid_rcto = 0. sol_grid_minl = 0. sol_grid_sorb = 0. sol_grid_ppag = 0. sol_grid_ppex = 0. sol_grid_tile = 0. sol_grid_soil = 0. sol_grid_resv = 0. sol_grid_wetl = 0. sol_grid_canl = 0. sol_grid_fpln = 0. !add up mass for the grid (convert g-->kg) do i=1,ncell if(gw_state(i)%stat == 1) then sol_grid_mbef = sol_grid_mbef + (gwsol_state(i)%solute(s)%mbef / 1000.) sol_grid_maft = sol_grid_maft + (gwsol_state(i)%solute(s)%maft / 1000.) sol_grid_rech = sol_grid_rech + (gwsol_ss(i)%solute(s)%rech / 1000.) sol_grid_gwsw = sol_grid_gwsw + (gwsol_ss(i)%solute(s)%gwsw / 1000.) sol_grid_swgw = sol_grid_swgw + (gwsol_ss(i)%solute(s)%swgw / 1000.) sol_grid_satx = sol_grid_satx + (gwsol_ss(i)%solute(s)%satx / 1000.) sol_grid_advn = sol_grid_advn + (gwsol_ss(i)%solute(s)%advn / 1000.) sol_grid_disp = sol_grid_disp + (gwsol_ss(i)%solute(s)%disp / 1000.) sol_grid_rcti = sol_grid_rcti + (gwsol_ss(i)%solute(s)%rcti / 1000.) sol_grid_rcto = sol_grid_rcto + (gwsol_ss(i)%solute(s)%rcto / 1000.) sol_grid_minl = sol_grid_minl + (gwsol_ss(i)%solute(s)%minl / 1000.) sol_grid_sorb = sol_grid_sorb + (gwsol_ss(i)%solute(s)%sorb / 1000.) sol_grid_ppag = sol_grid_ppag + (gwsol_ss(i)%solute(s)%ppag / 1000.) sol_grid_ppex = sol_grid_ppex + (gwsol_ss(i)%solute(s)%ppex / 1000.) sol_grid_tile = sol_grid_tile + (gwsol_ss(i)%solute(s)%tile / 1000.) sol_grid_soil = sol_grid_soil + (gwsol_ss(i)%solute(s)%soil / 1000.) sol_grid_resv = sol_grid_resv + (gwsol_ss(i)%solute(s)%resv / 1000.) sol_grid_wetl = sol_grid_wetl + (gwsol_ss(i)%solute(s)%wetl / 1000.) sol_grid_canl = sol_grid_canl + (gwsol_ss(i)%solute(s)%canl / 1000.) sol_grid_fpln = sol_grid_fpln + (gwsol_ss(i)%solute(s)%fpln / 1000.) endif enddo sol_grid_sorb = sol_grid_sorb * (-1) !leaving groundwater (sorbing to aquifer material) !calculate mass error if(sol_grid_maft > 0) then mass_error = (1- ((sol_grid_mbef + sol_grid_rech + sol_grid_gwsw + sol_grid_swgw + & sol_grid_satx + sol_grid_advn + sol_grid_disp + & sol_grid_rcti + sol_grid_rcto + sol_grid_minl + & sol_grid_ppag + sol_grid_ppex + sol_grid_tile + sol_grid_soil + & sol_grid_resv + sol_grid_wetl + sol_grid_canl + sol_grid_fpln) / sol_grid_maft)) * 100 endif !print out daily values for the solute if(gwflag_day == 1) then write(out_solbal_dy+s,102) time%yrc,time%day,gw_time_step, & sol_grid_mbef,sol_grid_maft,sol_grid_rech,sol_grid_gwsw,sol_grid_swgw, & sol_grid_satx,sol_grid_soil,sol_grid_advn,sol_grid_disp, & sol_grid_rcti,sol_grid_rcto,sol_grid_minl, & sol_grid_sorb,sol_grid_ppag,sol_grid_ppex,sol_grid_tile,sol_grid_resv, & sol_grid_wetl,sol_grid_canl,sol_grid_fpln, & mass_error endif !add grid values to yearly and total mass values !yearly (kg) sol_grid_chng_yr(s) = sol_grid_chng_yr(s) + (sol_grid_maft-sol_grid_mbef) sol_grid_rech_yr(s) = sol_grid_rech_yr(s) + sol_grid_rech sol_grid_gwsw_yr(s) = sol_grid_gwsw_yr(s) + sol_grid_gwsw sol_grid_swgw_yr(s) = sol_grid_swgw_yr(s) + sol_grid_swgw sol_grid_satx_yr(s) = sol_grid_satx_yr(s) + sol_grid_satx sol_grid_advn_yr(s) = sol_grid_advn_yr(s) + sol_grid_advn sol_grid_disp_yr(s) = sol_grid_disp_yr(s) + sol_grid_disp sol_grid_rcti_yr(s) = sol_grid_rcti_yr(s) + sol_grid_rcti sol_grid_rcto_yr(s) = sol_grid_rcto_yr(s) + sol_grid_rcto sol_grid_minl_yr(s) = sol_grid_minl_yr(s) + sol_grid_minl sol_grid_sorb_yr(s) = sol_grid_sorb_yr(s) + sol_grid_sorb sol_grid_ppag_yr(s) = sol_grid_ppag_yr(s) + sol_grid_ppag sol_grid_ppex_yr(s) = sol_grid_ppex_yr(s) + sol_grid_ppex sol_grid_tile_yr(s) = sol_grid_tile_yr(s) + sol_grid_tile sol_grid_soil_yr(s) = sol_grid_soil_yr(s) + sol_grid_soil sol_grid_resv_yr(s) = sol_grid_resv_yr(s) + sol_grid_resv sol_grid_wetl_yr(s) = sol_grid_wetl_yr(s) + sol_grid_wetl sol_grid_canl_yr(s) = sol_grid_canl_yr(s) + sol_grid_canl sol_grid_fpln_yr(s) = sol_grid_fpln_yr(s) + sol_grid_fpln !total (kg) sol_grid_chng_tt(s) = sol_grid_chng_tt(s) + (sol_grid_maft-sol_grid_mbef) sol_grid_rech_tt(s) = sol_grid_rech_tt(s) + sol_grid_rech sol_grid_gwsw_tt(s) = sol_grid_gwsw_tt(s) + sol_grid_gwsw sol_grid_swgw_tt(s) = sol_grid_swgw_tt(s) + sol_grid_swgw sol_grid_satx_tt(s) = sol_grid_satx_tt(s) + sol_grid_satx sol_grid_advn_tt(s) = sol_grid_advn_tt(s) + sol_grid_advn sol_grid_disp_tt(s) = sol_grid_disp_tt(s) + sol_grid_disp sol_grid_rcti_tt(s) = sol_grid_rcti_tt(s) + sol_grid_rcti sol_grid_rcto_tt(s) = sol_grid_rcto_tt(s) + sol_grid_rcto sol_grid_minl_tt(s) = sol_grid_minl_tt(s) + sol_grid_minl sol_grid_sorb_tt(s) = sol_grid_sorb_tt(s) + sol_grid_sorb sol_grid_ppag_tt(s) = sol_grid_ppag_tt(s) + sol_grid_ppag sol_grid_ppex_tt(s) = sol_grid_ppex_tt(s) + sol_grid_ppex sol_grid_tile_tt(s) = sol_grid_tile_tt(s) + sol_grid_tile sol_grid_soil_tt(s) = sol_grid_soil_tt(s) + sol_grid_soil sol_grid_resv_tt(s) = sol_grid_resv_tt(s) + sol_grid_resv sol_grid_wetl_tt(s) = sol_grid_wetl_tt(s) + sol_grid_wetl sol_grid_canl_tt(s) = sol_grid_canl_tt(s) + sol_grid_canl sol_grid_fpln_tt(s) = sol_grid_fpln_tt(s) + sol_grid_fpln enddo !go to next solute endif !check for solute transport !7. end of month: average head, average solute concentration, sum pumped groundwater (m3) for HRUs if (time%end_mo == 1) then !huc12 values, if national model is in use if(nat_model == 1) then do n=1,sp_ob%outlet !loop through the huc12 catchments write(out_huc12wb_mo,125) time%yrc,time%mo,huc12(n),(gw_huc12_wb_mo(i,n),i=1,15) enddo gw_huc12_wb_mo = 0. endif !monthly average groundwater head do i=1,ncell gw_state(i)%hdmo = gw_state(i)%hdmo / time%day_mo enddo write(out_head_mo,*) time%yrc,time%mo if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_state(cell_id_usg(i,j))%hdmo endif enddo enddo do i=1,grid_nrow write(out_head_mo,100) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_head_mo,121) (gw_state(i)%hdmo,i=1,ncell) endif write(out_head_mo,*) !zero out for next month do i=1,ncell gw_state(i)%hdmo = 0. enddo !monthly average solute concentration if(gw_solute_flag == 1) then write(out_conc_mo,*) time%yrc,time%mo do s=1,gw_nsolute !calculate average concentration do i=1,ncell gwsol_state(i)%solute(s)%cnmo = gwsol_state(i)%solute(s)%cnmo / time%day_mo enddo !write out write(out_conc_mo,*) gwsol_nm(s) !solute name if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_state(cell_id_usg(i,j))%solute(s)%cnmo endif enddo enddo do i=1,grid_nrow write(out_conc_mo,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_conc_mo,121) (gwsol_state(i)%solute(s)%cnmo,i=1,ncell) endif !zero out for next month do i=1,ncell gwsol_state(i)%solute(s)%cnmo = 0. enddo enddo !next solute write(out_conc_mo,*) endif !pumping (irrigation) (for HRUs) do i=1,sp_ob%hru hru_pump_mo_all(i,((time%yrs-1)*12)+time%mo) = hru_pump_mo(i) enddo hru_pump_mo = 0. endif !8. end of year: write out annual flux values if(time%end_yr == 1) then !annual average groundwater head do i=1,ncell gw_state(i)%hdyr = gw_state(i)%hdyr / time%day enddo write(out_head_yr,*) time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_state(cell_id_usg(i,j))%hdyr endif enddo enddo do i=1,grid_nrow write(out_head_yr,100) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_head_yr,121) (gw_state(i)%hdyr,i=1,ncell) endif write(out_head_yr,*) !zero out for next year do i=1,ncell gw_state(i)%hdyr = 0. enddo !annual average solute concentration if (gw_solute_flag == 1) then write(out_conc_yr,*) time%yrc do s=1,gw_nsolute !calculate average concentration do i=1,ncell gwsol_state(i)%solute(s)%cnyr = gwsol_state(i)%solute(s)%cnyr / time%day enddo !write out write(out_conc_yr,*) gwsol_nm(s) !solute name if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_state(cell_id_usg(i,j))%solute(s)%cnyr endif enddo enddo do i=1,grid_nrow write(out_conc_yr,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_conc_yr,121) (gwsol_state(i)%solute(s)%cnyr,i=1,ncell) endif !zero out for next year do i=1,ncell gwsol_state(i)%solute(s)%cnyr = 0. enddo enddo !next solute write(out_conc_yr,*) endif !check for solutes !compute average daily groundwater fluxes (m3/day) for the year do i=1,ncell gw_ss_sum(i)%rech = gw_ss_sum(i)%rech / time%day_end_yr gw_ss_sum(i)%gwet = gw_ss_sum(i)%gwet / time%day_end_yr gw_ss_sum(i)%gwsw = gw_ss_sum(i)%gwsw / time%day_end_yr gw_ss_sum(i)%swgw = gw_ss_sum(i)%swgw / time%day_end_yr gw_ss_sum(i)%satx = gw_ss_sum(i)%satx / time%day_end_yr gw_ss_sum(i)%soil = gw_ss_sum(i)%soil / time%day_end_yr gw_ss_sum(i)%latl = gw_ss_sum(i)%latl / time%day_end_yr gw_ss_sum(i)%bndr = gw_ss_sum(i)%bndr / time%day_end_yr gw_ss_sum(i)%ppag = gw_ss_sum(i)%ppag / time%day_end_yr gw_ss_sum(i)%ppdf = gw_ss_sum(i)%ppdf / time%day_end_yr gw_ss_sum(i)%ppex = gw_ss_sum(i)%ppex / time%day_end_yr gw_ss_sum(i)%tile = gw_ss_sum(i)%tile / time%day_end_yr gw_ss_sum(i)%resv = gw_ss_sum(i)%resv / time%day_end_yr gw_ss_sum(i)%wetl = gw_ss_sum(i)%wetl / time%day_end_yr gw_ss_sum(i)%fpln = gw_ss_sum(i)%fpln / time%day_end_yr gw_ss_sum(i)%canl = gw_ss_sum(i)%canl / time%day_end_yr enddo !compute average daily solute fluxes (kg/day) for the year if (gw_solute_flag == 1) then do i=1,ncell do s=1,gw_nsolute gwsol_ss_sum(i)%solute(s)%rech = (gwsol_ss_sum(i)%solute(s)%rech/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%gwsw = (gwsol_ss_sum(i)%solute(s)%gwsw/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%swgw = (gwsol_ss_sum(i)%solute(s)%swgw/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%soil = (gwsol_ss_sum(i)%solute(s)%soil/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%satx = (gwsol_ss_sum(i)%solute(s)%satx/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%ppex = (gwsol_ss_sum(i)%solute(s)%ppex/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%tile = (gwsol_ss_sum(i)%solute(s)%tile/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%resv = (gwsol_ss_sum(i)%solute(s)%resv/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%wetl = (gwsol_ss_sum(i)%solute(s)%wetl/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%fpln = (gwsol_ss_sum(i)%solute(s)%fpln/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%canl = (gwsol_ss_sum(i)%solute(s)%canl/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%advn = (gwsol_ss_sum(i)%solute(s)%advn/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%disp = (gwsol_ss_sum(i)%solute(s)%disp/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%rcti = (gwsol_ss_sum(i)%solute(s)%rcti/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%rcto = (gwsol_ss_sum(i)%solute(s)%rcto/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%minl = (gwsol_ss_sum(i)%solute(s)%minl/1000.) / time%day_end_yr !g --> kg gwsol_ss_sum(i)%solute(s)%sorb = (gwsol_ss_sum(i)%solute(s)%sorb/1000.) / time%day_end_yr !g --> kg enddo enddo endif !recharge write(out_gw_rech,*) 'Recharge for year (m3/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%rech endif enddo enddo do i=1,grid_nrow write(out_gw_rech,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_rech,121) (gw_ss_sum(i)%rech,i=1,ncell) endif write(out_gw_rech,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_rech,*) gwsol_nm(s),'recharge flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%rech endif enddo enddo do i=1,grid_nrow write(out_sol_rech,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_rech,121) (gwsol_ss_sum(i)%solute(s)%rech,i=1,ncell) endif write(out_sol_rech,*) enddo endif !groundwater ET write(out_gw_et,*) 'Groundwater ET for year:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%gwet endif enddo enddo do i=1,grid_nrow write(out_gw_et,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_et,121) (gw_ss_sum(i)%gwet,i=1,ncell) endif write(out_gw_et,*) !gw-sw exchange rates write(out_gwsw,*) 'GW-SW Exchange Rates for year:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%gwsw endif enddo enddo do i=1,grid_nrow write(out_gwsw,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gwsw,121) (gw_ss_sum(i)%gwsw,i=1,ncell) endif write(out_gwsw,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_gwsw,*) gwsol_nm(s),'gw-channel flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%gwsw endif enddo enddo do i=1,grid_nrow write(out_sol_gwsw,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_gwsw,121) (gwsol_ss_sum(i)%solute(s)%gwsw,i=1,ncell) endif write(out_sol_gwsw,*) enddo endif !saturation excess flow if(gw_satx_flag.eq.1) then write(out_gw_satex,*) 'Saturation Excess Volumes for:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%satx endif enddo enddo do i=1,grid_nrow write(out_gw_satex,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_satex,121) (gw_ss_sum(i)%satx,i=1,ncell) endif write(out_gw_satex,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_satx,*) gwsol_nm(s),'sat. excess flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%satx endif enddo enddo do i=1,grid_nrow write(out_sol_satx,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_satx,121) (gwsol_ss_sum(i)%solute(s)%satx,i=1,ncell) endif write(out_sol_satx,*) enddo endif endif !groundwater --> soil transfer if(gw_soil_flag.eq.1) then write(out_gw_soil,*) 'Groundwater --> Soil Transfer for:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%soil endif enddo enddo do i=1,grid_nrow write(out_gw_soil,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_soil,121) (gw_ss_sum(i)%soil,i=1,ncell) endif write(out_gw_soil,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_soil,*) gwsol_nm(s),'gw-->soil flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%soil endif enddo enddo do i=1,grid_nrow write(out_sol_soil,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_soil,121) (gwsol_ss_sum(i)%solute(s)%soil,i=1,ncell) endif write(out_sol_soil,*) enddo endif endif !lateral flow write(out_lateral,*) 'Lateral flow for year:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%latl endif enddo enddo do i=1,grid_nrow write(out_lateral,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_lateral,121) (gw_ss_sum(i)%latl,i=1,ncell) endif write(out_lateral,*) !tile drain flow if(gw_tile_flag == 1) then write(out_gw_tile,*) 'Tile Drain Outflow Volumes for:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%tile endif enddo enddo do i=1,grid_nrow write(out_gw_tile,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_tile,121) (gw_ss_sum(i)%tile,i=1,ncell) endif write(out_gw_tile,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_tile,*) gwsol_nm(s),'tile drain flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%tile endif enddo enddo do i=1,grid_nrow write(out_sol_tile,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_tile,121) (gwsol_ss_sum(i)%solute(s)%tile,i=1,ncell) endif write(out_sol_tile,*) enddo endif endif !pumping (irrigation) write(out_gw_pumpag,*) 'Pumping rate for year (m3/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%ppag endif enddo enddo do i=1,grid_nrow write(out_gw_pumpag,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_pumpag,121) (gw_ss_sum(i)%ppag,i=1,ncell) endif write(out_gw_pumpag,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_ppag,*) gwsol_nm(s),'ag pumping flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%ppag endif enddo enddo do i=1,grid_nrow write(out_sol_ppag,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_ppag,121) (gwsol_ss_sum(i)%solute(s)%ppag,i=1,ncell) endif write(out_sol_ppag,*) enddo endif !pumping (irrigation) (for HRUs) do i=1,sp_ob%hru hru_pump_yr_all(i,time%yrs) = hru_pump_yr(i) enddo hru_pump_yr = 0. !pumping deficit (not satisfied) (irrigation) write(out_gw_pumpdef,*) 'Pumping rate not satisfied for year (m3/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%ppdf endif enddo enddo do i=1,grid_nrow write(out_gw_pumpdef,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_pumpdef,121) (gw_ss_sum(i)%ppdf,i=1,ncell) endif write(out_gw_pumpdef,*) !pumping (user specified) if (gw_pumpex_flag == 1) then write(out_gw_pumpex,*) 'Pumping rate for year (m3/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%ppex endif enddo enddo do i=1,grid_nrow write(out_gw_pumpex,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_pumpex,121) (gw_ss_sum(i)%ppex,i=1,ncell) endif write(out_gw_pumpex,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_ppex,*) gwsol_nm(s),'ex pumping flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%ppex endif enddo enddo do i=1,grid_nrow write(out_sol_ppex,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_ppex,121) (gwsol_ss_sum(i)%solute(s)%ppex,i=1,ncell) endif write(out_sol_ppex,*) enddo endif endif !groundwater-reservoir exchange if (gw_res_flag == 1) then write(out_gw_res,*) 'Groundwater-Reservoir Exchange Volumes for:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%resv endif enddo enddo do i=1,grid_nrow write(out_gw_res,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_res,121) (gw_ss_sum(i)%resv,i=1,ncell) endif write(out_gw_res,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_resv,*) gwsol_nm(s),'gw-reservoir flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%resv endif enddo enddo do i=1,grid_nrow write(out_sol_resv,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_resv,121) (gwsol_ss_sum(i)%solute(s)%resv,i=1,ncell) endif write(out_sol_resv,*) enddo endif endif !groundwater-wetland exchange if (gw_wet_flag == 1) then write(out_gw_wet,*) 'Groundwater outflow to wetlands for:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%wetl endif enddo enddo do i=1,grid_nrow write(out_gw_wet,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_wet,121) (gw_ss_sum(i)%wetl,i=1,ncell) endif write(out_gw_wet,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_wetl,*) gwsol_nm(s),'gw-wetland flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%wetl endif enddo enddo do i=1,grid_nrow write(out_sol_wetl,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_wetl,121) (gwsol_ss_sum(i)%solute(s)%wetl,i=1,ncell) endif write(out_sol_wetl,*) enddo endif endif !groundwater-canal exchange if (gw_canal_flag == 1) then write(out_gw_canal,*) 'Groundwater-Canal Exchange Volumes for:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%canl endif enddo enddo do i=1,grid_nrow write(out_gw_canal,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_canal,121) (gw_ss_sum(i)%canl,i=1,ncell) endif write(out_gw_canal,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_canl,*) gwsol_nm(s),'gw-canal flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%canl endif enddo enddo do i=1,grid_nrow write(out_sol_canl,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_canl,121) (gwsol_ss_sum(i)%solute(s)%canl,i=1,ncell) endif write(out_sol_canl,*) enddo endif endif !floodplain exchange if (gw_fp_flag == 1) then write(out_gw_fp,*) 'Floodplain Exchange Volumes for:',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gw_ss_sum(cell_id_usg(i,j))%fpln endif enddo enddo do i=1,grid_nrow write(out_gw_fp,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_gw_fp,121) (gw_ss_sum(i)%fpln,i=1,ncell) endif write(out_gw_fp,*) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_fpln,*) gwsol_nm(s),'gw-floodplain flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%fpln endif enddo enddo do i=1,grid_nrow write(out_sol_fpln,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_fpln,121) (gwsol_ss_sum(i)%solute(s)%fpln,i=1,ncell) endif write(out_sol_fpln,*) enddo endif endif !chemical reaction (produced = positive values) if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_rcti,*) gwsol_nm(s),'chem. reaction flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%rcti endif enddo enddo do i=1,grid_nrow write(out_sol_rcti,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_rcti,121) (gwsol_ss_sum(i)%solute(s)%rcti,i=1,ncell) endif write(out_sol_rcti,*) enddo endif !chemical reaction (consumed = negative values) if(gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_rcto,*) gwsol_nm(s),'chem. reaction flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%rcto endif enddo enddo do i=1,grid_nrow write(out_sol_rcto,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_rcto,121) (gwsol_ss_sum(i)%solute(s)%rcto,i=1,ncell) endif write(out_sol_rcto,*) enddo endif !precipitation-dissolution if(gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_minl,*) gwsol_nm(s),'mineral dissolved mass for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%minl endif enddo enddo do i=1,grid_nrow write(out_sol_minl,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_minl,121) (gwsol_ss_sum(i)%solute(s)%minl,i=1,ncell) endif write(out_sol_minl,*) enddo endif !sorption if (gw_solute_flag == 1) then !solute mass flux do s=1,gw_nsolute write(out_sol_sorb,*) gwsol_nm(s),'sorption flux for year (kg/day):',time%yrc if(grid_type == "structured") then grid_val = 0. do i=1,grid_nrow do j=1,grid_ncol if(cell_id_usg(i,j) > 0) then grid_val(i,j) = gwsol_ss_sum(cell_id_usg(i,j))%solute(s)%sorb endif enddo enddo do i=1,grid_nrow write(out_sol_sorb,101) (grid_val(i,j),j=1,grid_ncol) enddo else write(out_sol_sorb,121) (gwsol_ss_sum(i)%solute(s)%sorb,i=1,ncell) endif write(out_sol_sorb,*) enddo endif !zero out flux sums to prepare for the next year !flow do i=1,ncell gw_ss_sum(i)%rech = 0. gw_ss_sum(i)%gwet = 0. gw_ss_sum(i)%gwsw = 0. gw_ss_sum(i)%swgw = 0. gw_ss_sum(i)%satx = 0. gw_ss_sum(i)%soil = 0. gw_ss_sum(i)%latl = 0. gw_ss_sum(i)%bndr = 0. gw_ss_sum(i)%ppag = 0. gw_ss_sum(i)%ppdf = 0. gw_ss_sum(i)%ppex = 0. gw_ss_sum(i)%tile = 0. gw_ss_sum(i)%resv = 0. gw_ss_sum(i)%wetl = 0. gw_ss_sum(i)%fpln = 0. gw_ss_sum(i)%canl = 0. enddo !solute do i=1,ncell do s=1,gw_nsolute gwsol_ss_sum(i)%solute(s)%rech = 0. gwsol_ss_sum(i)%solute(s)%gwsw = 0. gwsol_ss_sum(i)%solute(s)%swgw = 0. gwsol_ss_sum(i)%solute(s)%satx = 0. gwsol_ss_sum(i)%solute(s)%soil = 0. gwsol_ss_sum(i)%solute(s)%ppag = 0. gwsol_ss_sum(i)%solute(s)%ppex = 0. gwsol_ss_sum(i)%solute(s)%tile = 0. gwsol_ss_sum(i)%solute(s)%resv = 0. gwsol_ss_sum(i)%solute(s)%wetl = 0. gwsol_ss_sum(i)%solute(s)%fpln = 0. gwsol_ss_sum(i)%solute(s)%canl = 0. gwsol_ss_sum(i)%solute(s)%advn = 0. gwsol_ss_sum(i)%solute(s)%disp = 0. gwsol_ss_sum(i)%solute(s)%rcti = 0. gwsol_ss_sum(i)%solute(s)%rcto = 0. gwsol_ss_sum(i)%solute(s)%minl = 0. gwsol_ss_sum(i)%solute(s)%sorb = 0. enddo enddo !yearly water balance if(gwflag_yr.eq.1) then write(out_gwbal_yr,105) time%yrc, & ss_grid_yr%chng,ss_grid_yr%rech,ss_grid_yr%gwet,ss_grid_yr%gwsw,ss_grid_yr%swgw, & ss_grid_yr%satx,ss_grid_yr%soil,ss_grid_yr%latl,ss_grid_yr%bndr,ss_grid_yr%ppag, & ss_grid_yr%ppex,ss_grid_yr%tile,ss_grid_yr%resv,ss_grid_yr%wetl,ss_grid_yr%canl, & ss_grid_yr%fpln,ss_grid_yr%ppdf endif !if usgs wells are used, calculate and store annual average values for usgs wells and observation cells if (usgs_obs == 1) then do k=1,gw_num_obs_wells !compute the average head value for the year, for the current observation cell !also, compute the average saturated thickness for the year head_sum = 0. sat_sum = 0. do i=1,time%day head_sum = head_sum + gw_obs_head_annual(k,i) sat_sum = sat_sum + gw_obs_sat_annual(k,i) enddo head_avg = head_sum / time%day sim_head_vals(k,time%yrs) = head_avg sat_avg = sat_sum / time%day sim_sat_vals(k,time%yrs) = sat_avg enddo endif !zero out annual arrays ss_grid_yr%chng = 0. ss_grid_yr%rech = 0. ss_grid_yr%gwet = 0. ss_grid_yr%gwsw = 0. ss_grid_yr%swgw = 0. ss_grid_yr%satx = 0. ss_grid_yr%soil = 0. ss_grid_yr%latl = 0. ss_grid_yr%bndr = 0. ss_grid_yr%ppag = 0. ss_grid_yr%ppdf = 0. ss_grid_yr%ppex = 0. ss_grid_yr%tile = 0. ss_grid_yr%resv = 0. ss_grid_yr%wetl = 0. ss_grid_yr%canl = 0. ss_grid_yr%fpln = 0. !solute mass values if (gw_solute_flag == 1) then do s=1,gw_nsolute !loop through the solutes !write out annual values if(gwflag_yr.eq.1) then write(out_solbal_yr+s,105) time%yrc, & sol_grid_chng_yr,sol_grid_rech_yr,sol_grid_gwsw_yr,sol_grid_swgw_yr,sol_grid_satx_yr, & sol_grid_soil_yr,sol_grid_advn_yr,sol_grid_disp_yr, & sol_grid_rcti_yr,sol_grid_rcto_yr,sol_grid_minl_yr,sol_grid_sorb_yr, & sol_grid_ppag_yr,sol_grid_ppex_yr,sol_grid_tile_yr,sol_grid_resv_yr,sol_grid_wetl_yr, & sol_grid_canl_yr,sol_grid_fpln_yr endif !zero out values for next year sol_grid_chng_yr = 0. sol_grid_rech_yr = 0. sol_grid_gwsw_yr = 0. sol_grid_swgw_yr = 0. sol_grid_satx_yr = 0. sol_grid_soil_yr = 0. sol_grid_advn_yr = 0. sol_grid_disp_yr = 0. sol_grid_rcti_yr = 0. sol_grid_rcto_yr = 0. sol_grid_minl_yr = 0. sol_grid_sorb_yr = 0. sol_grid_ppag_yr = 0. sol_grid_ppex_yr = 0. sol_grid_tile_yr = 0. sol_grid_resv_yr = 0. sol_grid_wetl_yr = 0. sol_grid_canl_yr = 0. sol_grid_fpln_yr = 0. enddo !go to next solute endif endif !9. last day of the simulationreached: print out average annual values; calculate head and flow metrics if(time%yrc == time%yrc_end .and. time%day == time%day_end) then !pumping for HRUs num_months = time%nbyr * 12 do i=1,sp_ob%hru write(out_hru_pump_mo,105) i,(hru_pump_mo_all(i,j),j=1,num_months) write(out_hru_pump_yr,105) i,(hru_pump_yr_all(i,j),j=1,time%nbyr) enddo !average annual water balance ss_grid_tt%chng = ss_grid_tt%chng + (vaft_grid-vbef_grid) ss_grid_tt%rech = ss_grid_tt%rech / time%nbyr ss_grid_tt%gwet = ss_grid_tt%gwet / time%nbyr ss_grid_tt%gwsw = ss_grid_tt%gwsw / time%nbyr ss_grid_tt%swgw = ss_grid_tt%swgw / time%nbyr ss_grid_tt%satx = ss_grid_tt%satx / time%nbyr ss_grid_tt%soil = ss_grid_tt%soil / time%nbyr ss_grid_tt%latl = ss_grid_tt%latl / time%nbyr ss_grid_tt%bndr = ss_grid_tt%bndr / time%nbyr ss_grid_tt%ppag = ss_grid_tt%ppag / time%nbyr ss_grid_tt%ppdf = ss_grid_tt%ppdf / time%nbyr ss_grid_tt%ppex = ss_grid_tt%ppex / time%nbyr ss_grid_tt%tile = ss_grid_tt%tile / time%nbyr ss_grid_tt%resv = ss_grid_tt%resv / time%nbyr ss_grid_tt%wetl = ss_grid_tt%wetl / time%nbyr ss_grid_tt%canl = ss_grid_tt%canl / time%nbyr ss_grid_tt%fpln = ss_grid_tt%fpln / time%nbyr if(gwflag_aa.eq.1) then write(out_gwbal_aa,105) time%yrc, & ss_grid_tt%chng,ss_grid_tt%rech,ss_grid_tt%gwet,ss_grid_tt%gwsw,ss_grid_tt%swgw, & ss_grid_tt%satx,ss_grid_tt%soil,ss_grid_tt%latl,ss_grid_tt%bndr,ss_grid_tt%ppag, & ss_grid_tt%ppex,ss_grid_tt%tile,ss_grid_tt%resv,ss_grid_tt%wetl,ss_grid_tt%canl, & ss_grid_tt%fpln,ss_grid_tt%ppdf endif !average annual solute values if (gw_solute_flag == 1) then do s=1,gw_nsolute sol_grid_chng_tt(s) = sol_grid_chng_tt(s) + (sol_grid_maft-sol_grid_mbef) sol_grid_rech_tt(s) = sol_grid_rech_tt(s) / time%nbyr sol_grid_gwsw_tt(s) = sol_grid_gwsw_tt(s) / time%nbyr sol_grid_swgw_tt(s) = sol_grid_swgw_tt(s) / time%nbyr sol_grid_satx_tt(s) = sol_grid_satx_tt(s) / time%nbyr sol_grid_advn_tt(s) = sol_grid_advn_tt(s) / time%nbyr sol_grid_disp_tt(s) = sol_grid_disp_tt(s) / time%nbyr sol_grid_rcti_tt(s) = sol_grid_rcti_tt(s) / time%nbyr sol_grid_rcto_tt(s) = sol_grid_rcto_tt(s) / time%nbyr sol_grid_minl_tt(s) = sol_grid_minl_tt(s) / time%nbyr sol_grid_sorb_tt(s) = sol_grid_sorb_tt(s) / time%nbyr sol_grid_ppag_tt(s) = sol_grid_ppag_tt(s) / time%nbyr sol_grid_ppex_tt(s) = sol_grid_ppex_tt(s) / time%nbyr sol_grid_tile_tt(s) = sol_grid_tile_tt(s) / time%nbyr sol_grid_soil_tt(s) = sol_grid_soil_tt(s) / time%nbyr sol_grid_resv_tt(s) = sol_grid_resv_tt(s) / time%nbyr sol_grid_wetl_tt(s) = sol_grid_wetl_tt(s) / time%nbyr sol_grid_canl_tt(s) = sol_grid_canl_tt(s) / time%nbyr sol_grid_fpln_tt(s) = sol_grid_fpln_tt(s) / time%nbyr if(gwflag_aa.eq.1) then write(out_solbal_aa+s,105) time%yrc, & sol_grid_chng_tt(s),sol_grid_rech_tt(s),sol_grid_gwsw_tt(s),sol_grid_swgw_tt(s), & sol_grid_satx_tt(s),sol_grid_soil_tt(s),sol_grid_advn_tt(s),sol_grid_disp_tt(s), & sol_grid_rcti_tt(s),sol_grid_rcto_tt(s),sol_grid_minl_tt(s),sol_grid_sorb_tt(s), & sol_grid_ppag_tt(s),sol_grid_ppex_tt(s),sol_grid_tile_tt(s),sol_grid_resv_tt(s), & sol_grid_wetl_tt(s),sol_grid_canl_tt(s),sol_grid_fpln_tt(s) endif enddo !next solute endif !write out observed and simulated annual averaged gw head values at USGS well locations if (usgs_obs == 1) then if(gw_flow_cal.eq.0) then !first, calculate mean absolute error (MAE) (meters) for gw head values val_count = 0 error_sum = 0. sum_sat = 0. do k=1,gw_num_obs_wells usgs_yr = time%yrc_start - 1920 + 1 !get the right year from the USGS data set val_count_well = 0 error_sum_well = 0. sum_sat_well = 0. do j=1,time%nbyr if(usgs_head_vals(k,usgs_yr).ne.(-9.99)) then head_residual = usgs_head_vals(k,usgs_yr) - sim_head_vals(k,j) sum_sat = sum_sat + sim_sat_vals(k,j) sum_sat_well = sum_sat_well + sim_sat_vals(k,j) error_sum = error_sum + abs(head_residual) error_sum_well = error_sum_well + head_residual val_count = val_count + 1 val_count_well = val_count_well + 1 endif usgs_yr = usgs_yr + 1 enddo if(val_count_well.gt.0) then head_mae_well(k,1) = error_sum_well / val_count_well sat_div_well(k,1) = abs((sum_sat_well/val_count_well) / head_mae_well(k,1)) num_gw_meas_well(k,1) = val_count_well else head_mae_well(k,1) = 0. sat_div_well(k,1) = 0. num_gw_meas_well(k,1) = 0 endif enddo if(val_count.gt.0) then head_mae = error_sum / val_count sat_div = (sum_sat/val_count) / head_mae num_gw_meas = val_count else head_mae = 0. sat_div = 0. num_gw_meas = 0 endif !write out results write(out_gwobs_usgs,*) 'MAE (m), SAT/MAE, #Meas' write(out_gwobs_usgs,*) head_mae,sat_div,num_gw_meas write(out_gwobs_usgs,*) write(out_gwobs_usgs,*) 'MAE (m), SAT/MAE, #Meas per well' write(out_gwobs_usgs,*) 'Observed values - Simulated values' do k=1,gw_num_obs_wells write(out_gwobs_usgs,111) usgs_id(k),head_mae_well(k,1),sat_div_well(k,1),num_gw_meas_well(k,1) enddo else num_yrs_calb = gw_flow_cal_yrs num_yrs_test = time%nbyr - num_yrs_calb !calibration period !first, calculate mean absolute error (MAE) (meters) for gw head values val_count = 0 error_sum = 0. sum_sat = 0. do k=1,gw_num_obs_wells usgs_yr = time%yrc_start - 1920 + 1 !get the right year from the USGS data set val_count_well = 0 error_sum_well = 0. sum_sat_well = 0. do j=1,num_yrs_calb if(usgs_head_vals(k,usgs_yr).ne.(-9.99)) then head_residual = usgs_head_vals(k,usgs_yr) - sim_head_vals(k,j) sum_sat = sum_sat + sim_sat_vals(k,j) sum_sat_well = sum_sat_well + sim_sat_vals(k,j) error_sum = error_sum + abs(head_residual) error_sum_well = error_sum_well + head_residual val_count = val_count + 1 val_count_well = val_count_well + 1 endif usgs_yr = usgs_yr + 1 enddo if(val_count_well.gt.0) then head_mae_well(k,1) = error_sum_well / val_count_well sat_div_well(k,1) = abs((sum_sat_well/val_count_well) / head_mae_well(k,1)) num_gw_meas_well(k,1) = val_count_well else head_mae_well(k,1) = 0. sat_div_well(k,1) = 0. num_gw_meas_well(k,1) = 0 endif enddo if(val_count.gt.0) then head_mae_calb = error_sum / val_count sat_div_calb = (sum_sat/val_count) / head_mae_calb num_gw_meas_calb = val_count else head_mae_calb = 0. sat_div_calb = 0. num_gw_meas_calb = 0 endif !testing period val_count = 0 error_sum = 0. sum_sat = 0. do k=1,gw_num_obs_wells usgs_yr = time%yrc_start - 1920 + 1 + num_yrs_calb !get the right year from the USGS data set val_count_well = 0 error_sum_well = 0. sum_sat_well = 0. do j=1,num_yrs_test if(usgs_head_vals(k,usgs_yr).ne.(-9.99)) then head_residual = usgs_head_vals(k,usgs_yr) - sim_head_vals(k,j+num_yrs_calb) sum_sat = sum_sat + sim_sat_vals(k,j+num_yrs_calb) sum_sat_well = sum_sat_well + sim_sat_vals(k,j+num_yrs_calb) error_sum = error_sum + abs(head_residual) error_sum_well = error_sum_well + head_residual val_count = val_count + 1 val_count_well = val_count_well + 1 endif usgs_yr = usgs_yr + 1 enddo if(val_count_well.gt.0) then head_mae_well(k,2) = error_sum_well / val_count_well sat_div_well(k,2) = abs((sum_sat_well/val_count_well) / head_mae_well(k,1)) num_gw_meas_well(k,2) = val_count_well else head_mae_well(k,2) = 0. sat_div_well(k,2) = 0. num_gw_meas_well(k,2) = 0 endif enddo if(val_count.gt.0) then head_mae_test = error_sum / val_count sat_div_test = (sum_sat/val_count) / head_mae_test num_gw_meas_test = val_count else head_mae_test = 0. sat_div_test = 0. num_gw_meas_test = 0 endif !write out results write(out_gwobs_usgs,*) 'MAE_calib,MAE_test,SAT/MAE_calib,SAT_MAE_test,#calib,#test' write(out_gwobs_usgs,115) head_mae_calb,head_mae_test,sat_div_calb,sat_div_test,num_gw_meas_calb,num_gw_meas_test write(out_gwobs_usgs,*) write(out_gwobs_usgs,*) 'MAE (m), SAT/MAE per well, #Meas per well' write(out_gwobs_usgs,*) 'Observed values - Simulated values' do k=1,gw_num_obs_wells write(out_gwobs_usgs,116) usgs_id(k),head_mae_well(k,1),head_mae_well(k,2),sat_div_well(k,1),sat_div_well(k,2), & num_gw_meas_well(k,1),num_gw_meas_well(k,2) enddo endif write(out_gwobs_usgs,*) write(out_gwobs_usgs,*) 'Year, Well, Row, Col, Obs, Sim, SatThick' do k=1,gw_num_obs_wells write_yr = time%yrc_start usgs_yr = time%yrc_start - 1920 + 1 !get the right year from the USGS data set do j=1,time%nbyr write(out_gwobs_usgs,110) write_yr,usgs_id(k),gw_obs_cells(k),usgs_head_vals(k,usgs_yr), & sim_head_vals(k,j),sim_sat_vals(k,j) write_yr = write_yr + 1 usgs_yr = usgs_yr + 1 enddo enddo endif !write out monthly flow rates for specified channels if (stream_obs == 1) then if(gw_num_obs_chan.gt.0) then !first, calculate Nash-Sutcliffe model efficiciency coefficient (NSE) for each specified channel allocate (stream_nse(gw_num_obs_chan,2), source = 0.) allocate (stream_nse1(gw_num_obs_chan,2), source = 0.) allocate (stream_nnse(gw_num_obs_chan,2), source = 0.) allocate (stream_kg(gw_num_obs_chan,2), source = 0.) allocate (stream_pbias(gw_num_obs_chan,2), source = 0.) stream_nse = -9.99 stream_nse1 = -9.99 stream_nnse = -9.99 stream_kg = -9.99 stream_pbias = -9.99 if(gw_flow_cal.eq.0) then do i=1,gw_num_obs_chan !calculate mean of observed discharges; also count the number of months with streamflow data sum_obs_flow = 0. sum_sim_flow = 0. month_count = 0 do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_obs_flow = sum_obs_flow + obs_flow_vals(i,j) sum_sim_flow = sum_sim_flow + sim_flow_vals(i,j) month_count = month_count + 1 endif enddo if(month_count.gt.30) then !only calculate if number of data points > 30 mean_obs_flow = sum_obs_flow / month_count mean_sim_flow = sum_sim_flow / month_count !calculate NSE, NNSE, and NSE1 sum_resi_nse = 0. sum_diff_nse = 0. sum_resi_nse1 = 0. sum_diff_nse1 = 0. do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_resi_nse = sum_resi_nse + (obs_flow_vals(i,j) - sim_flow_vals(i,j))**2 sum_diff_nse = sum_diff_nse + (obs_flow_vals(i,j) - mean_obs_flow)**2 sum_resi_nse1 = sum_resi_nse1 + abs(obs_flow_vals(i,j) - sim_flow_vals(i,j)) sum_diff_nse1 = sum_diff_nse1 + abs(obs_flow_vals(i,j) - mean_obs_flow) endif enddo if(sum_diff_nse.ne.0) then stream_nse(i,1) = 1 - (sum_resi_nse / sum_diff_nse) stream_nse1(i,1) = 1 - (sum_resi_nse1 / sum_diff_nse1) stream_nnse(i,1) = 1 / (2 - stream_nse(i,1)) endif !calculate PBIAS sum_num = 0. sum_den = 0. do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_num = sum_num + (obs_flow_vals(i,j) - sim_flow_vals(i,j)) sum_den = sum_den + obs_flow_vals(i,j) endif enddo if(sum_den.ne.0) then stream_pbias(i,1) = (sum_num*100.) / sum_den endif !calculate KLING-GUPTA efficiency sum_num = 0. sum_den1 = 0. sum_den2 = 0. do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_num = sum_num + ((obs_flow_vals(i,j)-mean_obs_flow)*(sim_flow_vals(i,j)-mean_sim_flow)) sum_den1 = sum_den1 + (obs_flow_vals(i,j) - mean_obs_flow)**2 sum_den2 = sum_den2 + (sim_flow_vals(i,j) - mean_sim_flow)**2 endif enddo if((sum_den1+sum_den2).ne.0) then cc = sum_num / (sqrt(sum_den1)*sqrt(sum_den2)) endif sum_sim_sd = 0. sum_obs_sd = 0. do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_sim_sd = sum_sim_sd + (sim_flow_vals(i,j)-mean_sim_flow)**2 sum_obs_sd = sum_obs_sd + (obs_flow_vals(i,j)-mean_obs_flow)**2 endif enddo sim_sd = sqrt(sum_sim_sd/month_count) obs_sd = sqrt(sum_obs_sd/month_count) stream_kg(i,1) = 1 - sqrt((cc-1)**2 + ((sim_sd/obs_sd)-1)**2 + ((mean_sim_flow/mean_obs_flow)-1)**2) endif enddo !write out results write(out_strobs,*) 'Statistics for specified channels' write(out_strobs,*) 'Channel, count, NSE, NSE1, NNSE, PBIAS, KGE' do i=1,gw_num_obs_chan write(out_strobs,114) obs_channels(i),month_count,stream_nse(i,1),stream_nse1(i,1),stream_nnse(i,1), & stream_pbias(i,1),stream_kg(i,1) enddo else !calculate statistics for calibration and testing periods num_months_calb = 12 * gw_flow_cal_yrs num_months_test = num_months - num_months_calb do i=1,gw_num_obs_chan !calibration period !calculate mean of observed discharges sum_obs_flow = 0. sum_sim_flow = 0. month_count_calb = 0 do j=1,num_months_calb if(obs_flow_vals(i,j).gt.0) then sum_obs_flow = sum_obs_flow + obs_flow_vals(i,j) sum_sim_flow = sum_sim_flow + sim_flow_vals(i,j) month_count_calb = month_count_calb + 1 endif enddo if(month_count_calb.gt.30) then mean_obs_flow = sum_obs_flow / month_count_calb mean_sim_flow = sum_sim_flow / month_count_calb !calculate NSE, NNSE, and NSE1 sum_resi_nse = 0. sum_diff_nse = 0. sum_resi_nse1 = 0. sum_diff_nse1 = 0. do j=1,num_months_calb if(obs_flow_vals(i,j).gt.0) then sum_resi_nse = sum_resi_nse + (obs_flow_vals(i,j) - sim_flow_vals(i,j))**2 sum_diff_nse = sum_diff_nse + (obs_flow_vals(i,j) - mean_obs_flow)**2 sum_resi_nse1 = sum_resi_nse1 + abs(obs_flow_vals(i,j) - sim_flow_vals(i,j)) sum_diff_nse1 = sum_diff_nse1 + abs(obs_flow_vals(i,j) - mean_obs_flow) endif enddo if(sum_diff_nse.ne.0) then stream_nse(i,1) = 1 - (sum_resi_nse / sum_diff_nse) stream_nse1(i,1) = 1 - (sum_resi_nse1 / sum_diff_nse1) stream_nnse(i,1) = 1 / (2 - stream_nse(i,1)) endif endif !calculate PBIAS sum_num = 0. sum_den = 0. do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_num = sum_num + (obs_flow_vals(i,j) - sim_flow_vals(i,j)) sum_den = sum_den + obs_flow_vals(i,j) endif enddo if(sum_den.ne.0) then stream_pbias(i,1) = (sum_num*100.) / sum_den endif !calculate KLING-GUPTA efficiency sum_num = 0. sum_den1 = 0. sum_den2 = 0. do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_num = sum_num + ((obs_flow_vals(i,j)-mean_obs_flow)*(sim_flow_vals(i,j)-mean_sim_flow)) sum_den1 = sum_den1 + (obs_flow_vals(i,j) - mean_obs_flow)**2 sum_den2 = sum_den2 + (sim_flow_vals(i,j) - mean_sim_flow)**2 endif enddo if((sum_den1+sum_den2).ne.0) then cc = sum_num / (sqrt(sum_den1)*sqrt(sum_den2)) endif sum_sim_sd = 0. sum_obs_sd = 0. do j=1,num_months if(obs_flow_vals(i,j).gt.0) then sum_sim_sd = sum_sim_sd + (sim_flow_vals(i,j)-mean_sim_flow)**2 sum_obs_sd = sum_obs_sd + (obs_flow_vals(i,j)-mean_obs_flow)**2 endif enddo sim_sd = sqrt(sum_sim_sd/month_count_calb) obs_sd = sqrt(sum_obs_sd/month_count_calb) stream_kg(i,1) = 1 - sqrt((cc-1)**2 + ((sim_sd/obs_sd)-1)**2 + ((mean_sim_flow/mean_obs_flow)-1)**2) !testing period !calculate mean of observed discharges sum_obs_flow = 0. sum_sim_flow = 0. month_count_test = 0 do j=1,num_months_test if(obs_flow_vals(i,j+num_months_calb).gt.0) then sum_obs_flow = sum_obs_flow + obs_flow_vals(i,j+num_months_calb) sum_sim_flow = sum_sim_flow + sim_flow_vals(i,j+num_months_calb) month_count_test = month_count_test + 1 endif enddo if(month_count_test.gt.30) then mean_obs_flow = sum_obs_flow / month_count_test mean_sim_flow = sum_sim_flow / month_count_test !calculate NSE, NNSE, and NSE1 sum_resi_nse = 0. sum_diff_nse = 0. sum_resi_nse1 = 0. sum_diff_nse1 = 0. do j=1,num_months_test if(obs_flow_vals(i,j+num_months_calb).gt.0) then sum_resi_nse = sum_resi_nse + (obs_flow_vals(i,j+num_months_calb) - sim_flow_vals(i,j+num_months_calb))**2 sum_diff_nse = sum_diff_nse + (obs_flow_vals(i,j+num_months_calb) - mean_obs_flow)**2 sum_resi_nse1 = sum_resi_nse1 + abs(obs_flow_vals(i,j+num_months_calb) - sim_flow_vals(i,j+num_months_calb)) sum_diff_nse1 = sum_diff_nse1 + abs(obs_flow_vals(i,j+num_months_calb) - mean_obs_flow) endif enddo if(sum_diff_nse.ne.0) then stream_nse(i,2) = 1 - (sum_resi_nse / sum_diff_nse) stream_nse1(i,2) = 1 - (sum_resi_nse1 / sum_diff_nse1) stream_nnse(i,2) = 1 / (2 - stream_nse(i,1)) endif !calculate PBIAS sum_num = 0. sum_den = 0. do j=1,num_months if(obs_flow_vals(i,j+num_months_calb).gt.0) then sum_num = sum_num + (obs_flow_vals(i,j+num_months_calb) - sim_flow_vals(i,j+num_months_calb)) sum_den = sum_den + obs_flow_vals(i,j+num_months_calb) endif enddo if(sum_den.ne.0) then stream_pbias(i,2) = (sum_num*100.) / sum_den endif !calculate KLING-GUPTA efficiency sum_num = 0. sum_den1 = 0. sum_den2 = 0. do j=1,num_months if(obs_flow_vals(i,j+num_months_calb).gt.0) then sum_num = sum_num + ((obs_flow_vals(i,j+num_months_calb)-mean_obs_flow)* & (sim_flow_vals(i,j+num_months_calb)-mean_sim_flow)) sum_den1 = sum_den1 + (obs_flow_vals(i,j+num_months_calb) - mean_obs_flow)**2 sum_den2 = sum_den2 + (sim_flow_vals(i,j+num_months_calb) - mean_sim_flow)**2 endif enddo if((sum_den1+sum_den2).ne.0) then cc = sum_num / (sqrt(sum_den1)*sqrt(sum_den2)) endif sum_sim_sd = 0. sum_obs_sd = 0. do j=1,num_months if(obs_flow_vals(i,j+num_months_calb).gt.0) then sum_sim_sd = sum_sim_sd + (sim_flow_vals(i,j+num_months_calb)-mean_sim_flow)**2 sum_obs_sd = sum_obs_sd + (obs_flow_vals(i,j+num_months_calb)-mean_obs_flow)**2 endif enddo sim_sd = sqrt(sum_sim_sd/month_count_test) obs_sd = sqrt(sum_obs_sd/month_count_test) stream_kg(i,2) = 1 - sqrt((cc-1)**2 + ((sim_sd/obs_sd)-1)**2 + ((mean_sim_flow/mean_obs_flow)-1)**2) endif enddo !write out results write(out_strobs,*) 'Statistics for specified channels (calib/testing)' write(out_strobs,*) 'Channel, count, NSE, NSE1, NNSE, PBIAS, KGE' do i=1,gw_num_obs_chan write(out_strobs,114) obs_channels(i),month_count_calb,stream_nse(i,1),stream_nse1(i,1),stream_nnse(i,1), & stream_pbias(i,1),stream_kg(i,1) write(out_strobs,114) obs_channels(i),month_count_test,stream_nse(i,2),stream_nse1(i,2),stream_nnse(i,2), & stream_pbias(i,2),stream_kg(i,2) enddo endif write(out_strobs,*) write(out_strobs,*) 'Monthly values for each specified channel' do i=1,gw_num_obs_chan write(out_strobs,*) 'Channel:',obs_channels(i) write(out_strobs,*) 'Month Observed Simulated' do j=1,num_months write(out_strobs,113) j,obs_flow_vals(i,j),sim_flow_vals(i,j) enddo write(out_strobs,*) enddo endif endif !write out average annual groundwater budgets for each huc12 catchment if(nat_model == 1) then do n=1,sp_ob%outlet !loop through the huc12 catchments write(out_huc12wb,112) huc12(n),(gw_huc12_wb(i,n),i=1,15) enddo endif !if soft calibration, prepare for next simulation sim_month = 1 if(stream_obs == 1) then deallocate(stream_nse) deallocate(stream_nse1) deallocate(stream_nnse) deallocate(stream_kg) deallocate(stream_pbias) endif endif !end of simulation check !10. zero out flux arrays for next day !flow do i=1,ncell gw_ss(i)%rech = 0. gw_ss(i)%gwet = 0. gw_ss(i)%gwsw = 0. gw_ss(i)%swgw = 0. gw_ss(i)%satx = 0. gw_ss(i)%soil = 0. gw_ss(i)%latl = 0. gw_ss(i)%bndr = 0. gw_ss(i)%ppag = 0. gw_ss(i)%ppdf = 0. gw_ss(i)%ppex = 0. gw_ss(i)%tile = 0. gw_ss(i)%resv = 0. gw_ss(i)%wetl = 0. gw_ss(i)%fpln = 0. gw_ss(i)%canl = 0. gw_ss(i)%totl = 0. enddo satx_count = 0 !solutes if (gw_solute_flag == 1) then do i=1,ncell do s=1,2 !only for no3 and p gwsol_ss(i)%solute(s)%rech = 0. gwsol_ss(i)%solute(s)%gwsw = 0. gwsol_ss(i)%solute(s)%swgw = 0. gwsol_ss(i)%solute(s)%satx = 0. gwsol_ss(i)%solute(s)%soil = 0. gwsol_ss(i)%solute(s)%ppag = 0. gwsol_ss(i)%solute(s)%ppex = 0. gwsol_ss(i)%solute(s)%tile = 0. gwsol_ss(i)%solute(s)%resv = 0. gwsol_ss(i)%solute(s)%wetl = 0. gwsol_ss(i)%solute(s)%fpln = 0. gwsol_ss(i)%solute(s)%canl = 0. gwsol_ss(i)%solute(s)%advn = 0. gwsol_ss(i)%solute(s)%disp = 0. gwsol_ss(i)%solute(s)%rcti = 0. gwsol_ss(i)%solute(s)%rcto = 0. gwsol_ss(i)%solute(s)%minl = 0. gwsol_ss(i)%solute(s)%sorb = 0. gwsol_ss(i)%solute(s)%totl = 0. enddo do s=3,gw_nsolute !for salt and cs: do not zero out for recharge, reactions, or sorption (needed in salt_balance and cs_balance) !gwsol_ss(i)%solute(s)%rech = 0. gwsol_ss(i)%solute(s)%gwsw = 0. gwsol_ss(i)%solute(s)%swgw = 0. gwsol_ss(i)%solute(s)%satx = 0. gwsol_ss(i)%solute(s)%soil = 0. gwsol_ss(i)%solute(s)%ppag = 0. gwsol_ss(i)%solute(s)%ppex = 0. gwsol_ss(i)%solute(s)%tile = 0. gwsol_ss(i)%solute(s)%resv = 0. gwsol_ss(i)%solute(s)%wetl = 0. gwsol_ss(i)%solute(s)%fpln = 0. gwsol_ss(i)%solute(s)%canl = 0. gwsol_ss(i)%solute(s)%advn = 0. gwsol_ss(i)%solute(s)%disp = 0. !gwsol_ss(i)%solute(s)%rcti = 0. !gwsol_ss(i)%solute(s)%rcto = 0. !gwsol_ss(i)%solute(s)%minl = 0. !gwsol_ss(i)%solute(s)%sorb = 0. gwsol_ss(i)%solute(s)%totl = 0. enddo enddo mass_adv = 0. mass_dsp = 0. mass_rct = 0. mass_sorb = 0. endif 100 format(10000(f12.3)) 101 format(10000(e12.3)) 102 format(i8,i8,f10.3,e16.7,e16.7,1000(e13.4)) !*** tu Wunused-label: 103 format(i8,i8,i8,i8,i8,i8,i8,50(f15.3)) !*** tu Wunused-label: 104 format(10000(f12.2)) 105 format(i8,50(e13.4)) !*** tu Wunused-label: 106 format(i8,i8,i8,50(f12.3)) !*** tu Wunused-label: 108 format(i8,2x,50(e12.4)) !*** tu Wunused-label: 109 format(i8,i8,1000(e12.3)) 110 format(i8,f20.1,i8,f12.3,f12.3,f12.3) 111 format(f20.1,f12.3,f12.3,i8) 112 format(f15.1,50(e13.4)) 113 format(i8,1000(f12.3)) 114 format(i8,i8,1000(e12.3)) 115 format(f12.3,f12.3,f12.3,f12.3,i8,i8) 116 format(f20.1,f12.3,f12.3,f12.3,f12.3,i8,i8) 119 format(i8,i8,1000(f12.3)) !120 format(<out_cols>(f12.3)) !121 format(<out_cols>(e12.3)) 120 format(f12.3) 121 format(e12.3) 125 format(3x,i8,2x,i8,7x,f15.1,50(e13.4)) return end subroutine gwflow_simulate