      subroutine gwflow_output_init

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine opens all gwflow output files and writes headers
!!    (extracted from gwflow_read)

      use gwflow_module
      use hydrograph_module
      use sd_channel_module
      use time_module
      use constituent_mass_module, only : cs_db
      use basin_module, only : pco, bsn

      implicit none

      logical :: i_exist
      integer :: in_gw = 1230
      integer :: s = 0
      character*4 aString
      character(len=16) :: hydsep_hdr(10) = ""
      !general
      character(len=13) :: header = ""
      character*100 file_name(50)
      character*100 file_name_scalar
      integer :: i = 0
      integer :: j = 0
      integer :: n = 0
      integer :: max_num = 0
      integer :: wb_cell = 0
      real :: group_area = 0.

      !initialize groundwater balance ---------------------------------------------------------------------------------------------------------------
      write(out_gw,*)
      write(out_gw,*) '     initialize groundwater balance files and arrays'

      !write cell definition file (maps cell index to spatial location for all output)
      call gwflow_write_celldef
      write(out_gw,*) '          cell definition written to gwflow_cell_definition.txt'

      !--- print.prt overrides: when gwflow entries are listed in print.prt, override gwflow.input flags ---
      if(pco%gwflow_wb%already_read_in) then
        gwflag_day = 0; gwflag_mon = 0; gwflag_yr = 0; gwflag_aa = 0
        if(pco%gwflow_wb%d == 'y') gwflag_day = 1
        if(pco%gwflow_wb%m == 'y') gwflag_mon = 1
        if(pco%gwflow_wb%y == 'y') gwflag_yr  = 1
        if(pco%gwflow_wb%a == 'y') gwflag_aa  = 1
        write(out_gw,*) '          print.prt gwflow_wb overrides gwflow.input flags'
      endif
      if(pco%gwflow_obs%already_read_in) then
        gwflag_obs = 0
        if(pco%gwflow_obs%d == 'y') gwflag_obs = 1
        write(out_gw,*) '          print.prt gwflow_obs controls observation well output'
      endif
      if(pco%gwflow_pump%already_read_in) then
        gwflag_pump = 0
        if(pco%gwflow_pump%d == 'y' .or. pco%gwflow_pump%m == 'y' .or. &
           pco%gwflow_pump%y == 'y' .or. pco%gwflow_pump%a == 'y') gwflag_pump = 1
        write(out_gw,*) '          print.prt gwflow_pump controls HRU pumping output'
      endif
      if(pco%gwflow_heat%already_read_in) then
        gwflag_heat = 0
        if(pco%gwflow_heat%d == 'y' .or. pco%gwflow_heat%y == 'y' .or. &
           pco%gwflow_heat%a == 'y') gwflag_heat = 1
        write(out_gw,*) '          print.prt gwflow_heat controls heat balance output'
      endif
      if(pco%gwflow_solute%already_read_in) then
        gwflag_solute = 0
        if(pco%gwflow_solute%d == 'y' .or. pco%gwflow_solute%m == 'y' .or. &
           pco%gwflow_solute%y == 'y' .or. pco%gwflow_solute%a == 'y') gwflag_solute = 1
        write(out_gw,*) '          print.prt gwflow_solute controls solute balance output'
      endif
      if(pco%gwflow_flux%already_read_in) then
        gwflag_flux = 0
        if(pco%gwflow_flux%d == 'y' .or. pco%gwflow_flux%y == 'y') gwflag_flux = 1
        write(out_gw,*) '          print.prt gwflow_flux controls canal/pond/tile/chan diagnostic output'
      endif

      !open file to track daily groundwater water balance
      if(gwflag_day.eq.1) then
        open(out_gwbal,file='gwflow_basin_wb_day.txt')
        write(out_gwbal,*) 'Groundwater watershed-wide fluxes for each day (mm)'
        write(out_gwbal,8000) 'jday','mon','day','yr','unit','gis_id', &
          'name','ts','vbef','vaft','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral','boundary','pump_allo','pump_ext', &
          'tile','reservoir','wetland','canal','floodplain','pond', &
          'phytorem','error','sat_frac','wt_depth','pump_def'
        write(out_gwbal,8000) '','','','','','','','---','mm','mm', &
          'mm','mm','mm','mm','mm','mm','mm','mm','mm','mm', &
          'mm','mm','mm','mm','mm','mm','mm','%','frac','m','mm'
      endif

      !open file to track monthly groundwater water balance
      if(gwflag_mon.eq.1) then
        open(out_gwbal_mon,file='gwflow_basin_wb_mon.txt')
        write(out_gwbal_mon,*) 'Groundwater watershed-wide fluxes for each month (mm)'
        write(out_gwbal_mon,8000) 'jday','mon','day','yr','unit','gis_id', &
          'name','dvol','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral','boundary','pump_allo','pump_ext', &
          'tile','reservoir','wetland','canal','floodplain','pond', &
          'phytorem','pump_def'
        write(out_gwbal_mon,8000) '','','','','','','','mm','mm','mm', &
          'mm','mm','mm','mm','mm','mm','mm','mm','mm','mm', &
          'mm','mm','mm','mm','mm','mm'
      endif

      !open file to track yearly groundwater water balance
      if(gwflag_yr.eq.1) then
        open(out_gwbal_yr,file='gwflow_basin_wb_yr.txt')
        write(out_gwbal_yr,*) 'Groundwater watershed-wide fluxes for each year (mm)'
        write(out_gwbal_yr,8000) 'jday','mon','day','yr','unit','gis_id', &
          'name','dvol','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral','boundary','pump_allo','pump_ext', &
          'tile','reservoir','wetland','canal','floodplain','pond', &
          'phytorem','pump_def'
        write(out_gwbal_yr,8000) '','','','','','','','mm','mm','mm', &
          'mm','mm','mm','mm','mm','mm','mm','mm','mm','mm', &
          'mm','mm','mm','mm','mm','mm'
      endif

      !open file to write out average annual groundwater water balance
        if(gwflag_aa.eq.1) then
        open(out_gwbal_aa,file='gwflow_basin_wb_aa.txt')
        write(out_gwbal_aa,*) 'Average annual groundwater watershed-wide fluxes (mm)'
        write(out_gwbal_aa,8000) 'jday','mon','day','yr','unit','gis_id', &
          'name','dvol','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral','boundary','pump_allo','pump_ext', &
          'tile','reservoir','wetland','canal','floodplain','pond', &
          'phytorem','pump_def'
        write(out_gwbal_aa,8000) '','','','','','','','mm','mm','mm', &
          'mm','mm','mm','mm','mm','mm','mm','mm','mm','mm', &
          'mm','mm','mm','mm','mm','mm'
        endif

      !open files to track daily groundwater water balance for selected groups of cells
      inquire(file='gwflow.wbgroups',exist=i_exist)
      if(i_exist) then
        gw_group_flag = 1
        open(in_gw,file='gwflow.wbgroups')
        read(in_gw,*) header
        read(in_gw,*) gw_wb_grp_num !number of cell groups
        read(in_gw,*) max_num !maximum number of cells in a group
        allocate(gw_wb_grp_ncell(gw_wb_grp_num))
        allocate(gw_wb_grp_cells(gw_wb_grp_num,max_num))
        gw_wb_grp_ncell = 0
        !loop through the cell groups
        do i=1,gw_wb_grp_num
          !read in the cells in each group
          read(in_gw,*) header
          read(in_gw,*) gw_wb_grp_ncell(i)
          group_area = 0.
          do j=1,gw_wb_grp_ncell(i)
            read(in_gw,*) wb_cell
            if(wb_cell > 0) then
              if(grid_type == "structured") then
                wb_cell = cell_id_list(wb_cell)
                gw_wb_grp_cells(i,j) = wb_cell
              elseif(grid_type == "unstructured") then
                gw_wb_grp_cells(i,j) = wb_cell
              endif
              group_area = group_area + gw_state(wb_cell)%area !m2
            endif
          enddo
          !open the water balance output file
          write(aString,1091) i
          file_name_scalar = 'gwflow_group_wb_day_'//aString
          open(out_gwbal_grp+i,file=file_name_scalar)
          write(out_gwbal_grp+i,*) 'Groundwater fluxes for cell group (m3)',i
          write(out_gwbal_grp+i,8000) 'jday','mon','day','yr','unit', &
            'gis_id','name','ts','vbef','vaft','recharge','gw_et', &
            'gw_sw','sw_gw','sat_excess','soil','lateral','boundary', &
            'pump_allo','pump_ext','tile','reservoir','wetland','canal', &
            'floodplain','pond','phytorem','error','wt_depth','pump_def'
          write(out_gwbal_grp+i,8000) '','','','','','','','---', &
            'm3','m3','m3','m3','m3','m3','m3','m3','m3','m3', &
            'm3','m3','m3','m3','m3','m3','m3','m3','m3','%', &
            'm','m3'
        enddo !go to next cell group
        close(in_gw)
  1091  format(i2)
      endif



      !initialize groundwater heat balance ----------------------------------------------------------------------------------------------------------
      if(gw_heat_flag == 1) then
      write(out_gw,*)
      write(out_gw,*) '     initialize groundwater heat balance files and arrays'

      !open file to track daily groundwater heat balance
      if(gwflag_day.eq.1) then
        open(out_heatbal_dy,file='gwflow_basin_heat_day.txt')
        write(out_heatbal_dy,*) 'Groundwater watershed-wide heat fluxes for each day (MJ)'
        write(out_heatbal_dy,8000) 'jday','mon','day','yr','unit','gis_id', &
          'name','ts','hbef','haft','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral','dispersion','boundary', &
          'pump_allo','pump_ext','tile','reservoir','wetland','canal', &
          'floodplain','pond','error','tavg'
        write(out_heatbal_dy,8000) '','','','','','','','---', &
          'MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ', &
          'MJ','MJ','MJ','MJ','MJ','MJ','MJ','%','degC'
      endif

      !open file to track yearly groundwater heat balance
      if(gwflag_yr.eq.1) then
        open(out_heatbal_yr,file='gwflow_basin_heat_yr.txt')
        write(out_heatbal_yr,*) 'Groundwater watershed-wide heat fluxes for each year (MJ)'
        write(out_heatbal_yr,8000) 'jday','mon','day','yr','unit','gis_id', &
          'name','hdel','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral','dispersion','boundary', &
          'pump_allo','pump_ext','tile','reservoir','wetland','canal', &
          'floodplain','pond'
        write(out_heatbal_yr,8000) '','','','','','','','MJ','MJ','MJ', &
          'MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ', &
          'MJ','MJ','MJ'
      endif

      !open file to write out average annual groundwater heat balance
      if(gwflag_aa.eq.1) then
        open(out_heatbal_aa,file='gwflow_basin_heat_aa.txt')
        write(out_heatbal_aa,*) 'Groundwater watershed-wide heat fluxes average annual (MJ)'
        write(out_heatbal_aa,8000) 'jday','mon','day','yr','unit','gis_id', &
          'name','hdel','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral','dispersion','boundary', &
          'pump_allo','pump_ext','tile','reservoir','wetland','canal', &
          'floodplain','pond'
        write(out_heatbal_aa,8000) '','','','','','','','MJ','MJ','MJ', &
          'MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ','MJ', &
          'MJ','MJ','MJ'
      endif

      endif



      !initialize solute mass balance ---------------------------------------------------------------------------------------------------------------
      if(gw_solute_flag == 1) then

        !allocate monthly, yearly, and total arrays
        !monthly
        allocate(sol_grid_chng_mo(gw_nsolute))
        allocate(sol_grid_rech_mo(gw_nsolute))
        allocate(sol_grid_gwsw_mo(gw_nsolute))
        allocate(sol_grid_swgw_mo(gw_nsolute))
        allocate(sol_grid_satx_mo(gw_nsolute))
        allocate(sol_grid_advn_mo(gw_nsolute))
        allocate(sol_grid_disp_mo(gw_nsolute))
        allocate(sol_grid_rcti_mo(gw_nsolute))
        allocate(sol_grid_rcto_mo(gw_nsolute))
        allocate(sol_grid_minl_mo(gw_nsolute))
        allocate(sol_grid_sorb_mo(gw_nsolute))
        allocate(sol_grid_ppag_mo(gw_nsolute))
        allocate(sol_grid_ppex_mo(gw_nsolute))
        allocate(sol_grid_tile_mo(gw_nsolute))
        allocate(sol_grid_soil_mo(gw_nsolute))
        allocate(sol_grid_resv_mo(gw_nsolute))
        allocate(sol_grid_wetl_mo(gw_nsolute))
        allocate(sol_grid_canl_mo(gw_nsolute))
        allocate(sol_grid_fpln_mo(gw_nsolute))
        allocate(sol_grid_pond_mo(gw_nsolute))
        !yearly
        allocate(sol_grid_chng_yr(gw_nsolute))
        allocate(sol_grid_rech_yr(gw_nsolute))
        allocate(sol_grid_gwsw_yr(gw_nsolute))
        allocate(sol_grid_swgw_yr(gw_nsolute))
        allocate(sol_grid_satx_yr(gw_nsolute))
        allocate(sol_grid_advn_yr(gw_nsolute))
        allocate(sol_grid_disp_yr(gw_nsolute))
        allocate(sol_grid_rcti_yr(gw_nsolute))
        allocate(sol_grid_rcto_yr(gw_nsolute))
        allocate(sol_grid_minl_yr(gw_nsolute))
        allocate(sol_grid_sorb_yr(gw_nsolute))
        allocate(sol_grid_ppag_yr(gw_nsolute))
        allocate(sol_grid_ppex_yr(gw_nsolute))
        allocate(sol_grid_tile_yr(gw_nsolute))
        allocate(sol_grid_soil_yr(gw_nsolute))
        allocate(sol_grid_resv_yr(gw_nsolute))
        allocate(sol_grid_wetl_yr(gw_nsolute))
        allocate(sol_grid_canl_yr(gw_nsolute))
        allocate(sol_grid_fpln_yr(gw_nsolute))
        allocate(sol_grid_pond_yr(gw_nsolute))
        !total
        allocate(sol_grid_chng_tt(gw_nsolute))
        allocate(sol_grid_rech_tt(gw_nsolute))
        allocate(sol_grid_gwsw_tt(gw_nsolute))
        allocate(sol_grid_swgw_tt(gw_nsolute))
        allocate(sol_grid_satx_tt(gw_nsolute))
        allocate(sol_grid_advn_tt(gw_nsolute))
        allocate(sol_grid_disp_tt(gw_nsolute))
        allocate(sol_grid_rcti_tt(gw_nsolute))
        allocate(sol_grid_rcto_tt(gw_nsolute))
        allocate(sol_grid_minl_tt(gw_nsolute))
        allocate(sol_grid_sorb_tt(gw_nsolute))
        allocate(sol_grid_ppag_tt(gw_nsolute))
        allocate(sol_grid_ppex_tt(gw_nsolute))
        allocate(sol_grid_tile_tt(gw_nsolute))
        allocate(sol_grid_soil_tt(gw_nsolute))
        allocate(sol_grid_resv_tt(gw_nsolute))
        allocate(sol_grid_wetl_tt(gw_nsolute))
        allocate(sol_grid_canl_tt(gw_nsolute))
        allocate(sol_grid_fpln_tt(gw_nsolute))
        allocate(sol_grid_pond_tt(gw_nsolute))

        !loop through the solutes
        do n=1,gw_nsolute

          !daily solute mass balance
          if(gwflag_day.eq.1) then
            !prepare solute mass balance output files (daily output)
            file_name(1) = 'gwflow_basin_sol_no3_day.txt'
            file_name(2) = 'gwflow_basin_sol_p_day.txt'
            if(cs_db%num_salts > 0) then
              file_name(3) = 'gwflow_basin_sol_so4_day.txt'
              file_name(4) = 'gwflow_basin_sol_ca_day.txt'
              file_name(5) = 'gwflow_basin_sol_mg_day.txt'
              file_name(6) = 'gwflow_basin_sol_na_day.txt'
              file_name(7) = 'gwflow_basin_sol_k_day.txt'
              file_name(8) = 'gwflow_basin_sol_cl_day.txt'
              file_name(9) = 'gwflow_basin_sol_co3_day.txt'
              file_name(10) = 'gwflow_basin_sol_hco3_day.txt'
            endif
            if(cs_db%num_cs > 0) then
              file_name(11) = 'gwflow_basin_sol_seo4_day.txt'
              file_name(12) = 'gwflow_basin_sol_seo3_day.txt'
            endif
            open(out_solbal_dy+n,file=file_name(n))
            write(out_solbal_dy+n,*) 'Solute basin daily loads (kg):',gwsol_nm(n)
            write(out_solbal_dy+n,8000) 'jday','mon','day','yr','unit', &
              'gis_id','name','ts','mbef','maft','recharge','gw_sw', &
              'sw_gw','sat_excess','soil','advection','dispersion', &
              'react_in','react_out','mineral','sorption','pump_allo', &
              'pump_ext','tile','reservoir','wetland','canal', &
              'floodplain','pond','error'
            write(out_solbal_dy+n,8000) '','','','','','','','---', &
              'kg','kg','kg','kg','kg','kg','kg','kg','kg','kg','kg', &
              'kg','kg','kg','kg','kg','kg','kg','kg','kg','kg','%'
          endif

          !monthly solute mass balance
          if(gwflag_mon.eq.1) then
            !prepare solute mass balance output files (monthly output)
            file_name(1) = 'gwflow_basin_sol_no3_mon.txt'
            file_name(2) = 'gwflow_basin_sol_p_mon.txt'
            if(cs_db%num_salts > 0) then
              file_name(3) = 'gwflow_basin_sol_so4_mon.txt'
              file_name(4) = 'gwflow_basin_sol_ca_mon.txt'
              file_name(5) = 'gwflow_basin_sol_mg_mon.txt'
              file_name(6) = 'gwflow_basin_sol_na_mon.txt'
              file_name(7) = 'gwflow_basin_sol_k_mon.txt'
              file_name(8) = 'gwflow_basin_sol_cl_mon.txt'
              file_name(9) = 'gwflow_basin_sol_co3_mon.txt'
              file_name(10) = 'gwflow_basin_sol_hco3_mon.txt'
            endif
            if(cs_db%num_cs > 0) then
              file_name(11) = 'gwflow_basin_sol_seo4_mon.txt'
              file_name(12) = 'gwflow_basin_sol_seo3_mon.txt'
            endif
            open(out_solbal_mo+n,file=file_name(n))
            write(out_solbal_mo+n,*) 'Solute basin monthly loads (kg):',gwsol_nm(n)
            write(out_solbal_mo+n,8000) 'jday','mon','day','yr','unit', &
              'gis_id','name','delm','recharge','gw_sw','sw_gw', &
              'sat_excess','soil','advection','dispersion','react_in', &
              'react_out','mineral','sorption','pump_allo','pump_ext', &
              'tile','reservoir','wetland','canal','floodplain','pond'
            write(out_solbal_mo+n,8000) '','','','','','','','kg','kg', &
              'kg','kg','kg','kg','kg','kg','kg','kg','kg','kg','kg', &
              'kg','kg','kg','kg','kg','kg','kg'
            !zero out monthly arrays
            sol_grid_chng_mo(n) = 0.
            sol_grid_rech_mo(n) = 0.
            sol_grid_gwsw_mo(n) = 0.
            sol_grid_swgw_mo(n) = 0.
            sol_grid_satx_mo(n) = 0.
            sol_grid_advn_mo(n) = 0.
            sol_grid_disp_mo(n) = 0.
            sol_grid_rcti_mo(n) = 0.
            sol_grid_rcto_mo(n) = 0.
            sol_grid_minl_mo(n) = 0.
            sol_grid_sorb_mo(n) = 0.
            sol_grid_ppag_mo(n) = 0.
            sol_grid_ppex_mo(n) = 0.
            sol_grid_tile_mo(n) = 0.
            sol_grid_soil_mo(n) = 0.
            sol_grid_resv_mo(n) = 0.
            sol_grid_wetl_mo(n) = 0.
            sol_grid_canl_mo(n) = 0.
            sol_grid_fpln_mo(n) = 0.
            sol_grid_pond_mo(n) = 0.
          endif

          !yearly solute mass balance
          if(gwflag_yr.eq.1) then
            !prepare solute mass balance output files (yearly output)
            file_name(1) = 'gwflow_basin_sol_no3_yr.txt'
            file_name(2) = 'gwflow_basin_sol_p_yr.txt'
            if(cs_db%num_salts > 0) then
              file_name(3) = 'gwflow_basin_sol_so4_yr.txt'
              file_name(4) = 'gwflow_basin_sol_ca_yr.txt'
              file_name(5) = 'gwflow_basin_sol_mg_yr.txt'
              file_name(6) = 'gwflow_basin_sol_na_yr.txt'
              file_name(7) = 'gwflow_basin_sol_k_yr.txt'
              file_name(8) = 'gwflow_basin_sol_cl_yr.txt'
              file_name(9) = 'gwflow_basin_sol_co3_yr.txt'
              file_name(10) = 'gwflow_basin_sol_hco3_yr.txt'
            endif
            if(cs_db%num_cs > 0) then
              file_name(11) = 'gwflow_basin_sol_seo4_yr.txt'
              file_name(12) = 'gwflow_basin_sol_seo3_yr.txt'
            endif
            open(out_solbal_yr+n,file=file_name(n))
            write(out_solbal_yr+n,*) 'Solute basin annual loads (kg):',gwsol_nm(n)
            write(out_solbal_yr+n,8000) 'jday','mon','day','yr','unit', &
              'gis_id','name','delm','recharge','gw_sw','sw_gw', &
              'sat_excess','soil','advection','dispersion','react_in', &
              'react_out','mineral','sorption','pump_allo','pump_ext', &
              'tile','reservoir','wetland','canal','floodplain','pond'
            write(out_solbal_yr+n,8000) '','','','','','','','kg','kg', &
              'kg','kg','kg','kg','kg','kg','kg','kg','kg','kg','kg', &
              'kg','kg','kg','kg','kg','kg','kg'
            !zero out yearly arrays
            sol_grid_chng_yr(n) = 0.
            sol_grid_rech_yr(n) = 0.
            sol_grid_gwsw_yr(n) = 0.
            sol_grid_swgw_yr(n) = 0.
            sol_grid_satx_yr(n) = 0.
            sol_grid_advn_yr(n) = 0.
            sol_grid_disp_yr(n) = 0.
            sol_grid_rcti_yr(n) = 0.
            sol_grid_rcto_yr(n) = 0.
            sol_grid_minl_yr(n) = 0.
            sol_grid_sorb_yr(n) = 0.
            sol_grid_ppag_yr(n) = 0.
            sol_grid_ppex_yr(n) = 0.
            sol_grid_tile_yr(n) = 0.
            sol_grid_soil_yr(n) = 0.
            sol_grid_resv_yr(n) = 0.
            sol_grid_wetl_yr(n) = 0.
            sol_grid_canl_yr(n) = 0.
            sol_grid_fpln_yr(n) = 0.
            sol_grid_pond_yr(n) = 0.
          endif

          !average annual solute mass balance
          if(gwflag_aa.eq.1) then
            !prepare solute mass balance output files (average annual output)
            file_name(1) = 'gwflow_basin_sol_no3_aa.txt'
            file_name(2) = 'gwflow_basin_sol_p_aa.txt'
            if(cs_db%num_salts > 0) then
              file_name(3) = 'gwflow_basin_sol_so4_aa.txt'
              file_name(4) = 'gwflow_basin_sol_ca_aa.txt'
              file_name(5) = 'gwflow_basin_sol_mg_aa.txt'
              file_name(6) = 'gwflow_basin_sol_na_aa.txt'
              file_name(7) = 'gwflow_basin_sol_k_aa.txt'
              file_name(8) = 'gwflow_basin_sol_cl_aa.txt'
              file_name(9) = 'gwflow_basin_sol_co3_aa.txt'
              file_name(10) = 'gwflow_basin_sol_hco3_aa.txt'
            endif
            if(cs_db%num_cs > 0) then
              file_name(11) = 'gwflow_basin_sol_seo4_aa.txt'
              file_name(12) = 'gwflow_basin_sol_seo3_aa.txt'
            endif
            open(out_solbal_aa+n,file=file_name(n))
            write(out_solbal_aa+n,*) 'Solute basin average annual loads (kg):',gwsol_nm(n)
            write(out_solbal_aa+n,8000) 'jday','mon','day','yr','unit', &
              'gis_id','name','delm','recharge','gw_sw','sw_gw', &
              'sat_excess','soil','advection','dispersion','react_in', &
              'react_out','mineral','sorption','pump_allo','pump_ext', &
              'tile','reservoir','wetland','canal','floodplain','pond'
            write(out_solbal_aa+n,8000) '','','','','','','','kg','kg', &
              'kg','kg','kg','kg','kg','kg','kg','kg','kg','kg','kg', &
              'kg','kg','kg','kg','kg','kg','kg'
            !zero out yearly arrays
            sol_grid_chng_tt(n) = 0.
            sol_grid_rech_tt(n) = 0.
            sol_grid_gwsw_tt(n) = 0.
            sol_grid_swgw_tt(n) = 0.
            sol_grid_satx_tt(n) = 0.
            sol_grid_advn_tt(n) = 0.
            sol_grid_disp_tt(n) = 0.
            sol_grid_rcti_tt(n) = 0.
            sol_grid_rcto_tt(n) = 0.
            sol_grid_minl_tt(n) = 0.
            sol_grid_sorb_tt(n) = 0.
            sol_grid_ppag_tt(n) = 0.
            sol_grid_ppex_tt(n) = 0.
            sol_grid_tile_tt(n) = 0.
            sol_grid_soil_tt(n) = 0.
            sol_grid_resv_tt(n) = 0.
            sol_grid_wetl_tt(n) = 0.
            sol_grid_canl_tt(n) = 0.
            sol_grid_fpln_tt(n) = 0.
            sol_grid_pond_tt(n) = 0.
          endif

        enddo !go to next solute

      endif !check for solutes


      !open cell-level long-format output files -----------------------------------------
      if(gwflag_day.eq.1) then
        open(out_gwcell_day,file='gwflow_cell_wb_day.txt')
        write(out_gwcell_day,*) 'gwflow cell-level daily output'
        write(out_gwcell_day,*) 'fluxes in m3/day; head and wtdepth in m'
        write(out_gwcell_day,'(a8,a6,a6,a8,a8,a10,a12,20a13)') &
          'jday','mon','day','yr','unit','gis_id','name', &
          'head','wt_depth','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral', &
          'pump_allo','pump_ext','tile','reservoir','wetland', &
          'floodplain','canal','pond','phytorem'
      endif
      if(gwflag_mon.eq.1) then
        open(out_gwcell_mon,file='gwflow_cell_wb_mon.txt')
        write(out_gwcell_mon,*) 'gwflow cell-level monthly output'
        write(out_gwcell_mon,*) 'values are monthly average daily rates' &
          //' (m3/day for fluxes, m for head/wtdepth)'
        write(out_gwcell_mon,'(a8,a6,a6,a8,a8,a10,a12,20a13)') &
          'jday','mon','day','yr','unit','gis_id','name', &
          'head','wt_depth','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral', &
          'pump_allo','pump_ext','tile','reservoir','wetland', &
          'floodplain','canal','pond','phytorem'
      endif
      if(gwflag_yr.eq.1) then
        open(out_gwcell_yr,file='gwflow_cell_wb_yr.txt')
        write(out_gwcell_yr,*) 'gwflow cell-level annual output'
        write(out_gwcell_yr,*) 'values are annual average daily rates' &
          //' (m3/day for fluxes, m for head/wtdepth)'
        write(out_gwcell_yr,'(a8,a6,a6,a8,a8,a10,a12,20a13)') &
          'jday','mon','day','yr','unit','gis_id','name', &
          'head','wt_depth','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral', &
          'pump_allo','pump_ext','tile','reservoir','wetland', &
          'floodplain','canal','pond','phytorem'
      endif
      if(gwflag_aa.eq.1) then
        open(out_gwcell_aa,file='gwflow_cell_wb_aa.txt')
        write(out_gwcell_aa,*) 'gwflow cell-level average annual output'
        write(out_gwcell_aa,*) 'values are avg annual daily rates' &
          //' (m3/day for fluxes, m for head/wtdepth)'
        write(out_gwcell_aa,'(a8,a6,a6,a8,a8,a10,a12,20a13)') &
          'jday','mon','day','yr','unit','gis_id','name', &
          'head','wt_depth','recharge','gw_et','gw_sw','sw_gw', &
          'sat_excess','soil','lateral', &
          'pump_allo','pump_ext','tile','reservoir','wetland', &
          'floodplain','canal','pond','phytorem'
      endif


      !format statement for standard SWAT+ basin-level headers
8000  format(4a6,2a8,a18,50a13)

      return
      end subroutine gwflow_output_init


      subroutine gwflow_output_day

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine computes and writes daily gwflow output:
!!    observation wells, cell groups, basin water balance, heat balance,
!!    solute mass balance; accumulates to monthly/yearly/aa totals

      use gwflow_module
      use hydrograph_module
      use sd_channel_module
      use time_module
      use basin_module, only : bsn

      implicit none

      integer :: i, j, k, s, iob
      integer :: cell_id, gis_id
      real :: sum
      real :: obs_temp, obs_no3, obs_p
      character(len=16) :: obs_name
      real :: frac_sat
      real :: depth_wt_avg
      real :: depth_wt_avg_grp
      real :: mass_error
      real :: temp_avg
      !grid totals of fluxes for day (m3, mm)
      !note: vbef_grid, vaft_grid are module variables (shared with output_aa)
      real*8 :: rech_grid
      real*8 :: gwet_grid
      real*8 :: gwsw_grid
      real*8 :: swgw_grid
      real*8 :: satx_grid
      real*8 :: soil_grid
      real*8 :: latl_grid
      real*8 :: bndr_grid
      real*8 :: ppag_grid
      real*8 :: ppdf_grid
      real*8 :: ppex_grid
      real*8 :: tile_grid
      real*8 :: resv_grid
      real*8 :: wetl_grid
      real*8 :: canl_grid
      real*8 :: fpln_grid
      real*8 :: pond_grid
      real*8 :: phyt_grid
      !cell group totals of fluxes for day (m3)
      real*8 :: vbef_grp
      real*8 :: vaft_grp
      real*8 :: rech_grp
      real*8 :: gwet_grp
      real*8 :: gwsw_grp
      real*8 :: swgw_grp
      real*8 :: satx_grp
      real*8 :: soil_grp
      real*8 :: latl_grp
      real*8 :: bndr_grp
      real*8 :: ppag_grp
      real*8 :: ppdf_grp
      real*8 :: ppex_grp
      real*8 :: tile_grp
      real*8 :: resv_grp
      real*8 :: wetl_grp
      real*8 :: canl_grp
      real*8 :: fpln_grp
      real*8 :: pond_grp
      real*8 :: phyt_grp
      !heat grid totals
      !note: heat_hbef_grid, heat_haft_grid are module variables (shared with output_aa)
      real*8 :: heat_rech_grid
      real*8 :: heat_gwet_grid
      real*8 :: heat_gwsw_grid
      real*8 :: heat_swgw_grid
      real*8 :: heat_satx_grid
      real*8 :: heat_soil_grid
      real*8 :: heat_latl_grid
      real*8 :: heat_disp_grid
      real*8 :: heat_bndr_grid
      real*8 :: heat_ppag_grid
      real*8 :: heat_ppex_grid
      real*8 :: heat_tile_grid
      real*8 :: heat_resv_grid
      real*8 :: heat_wetl_grid
      real*8 :: heat_canl_grid
      real*8 :: heat_fpln_grid
      real*8 :: heat_pond_grid
      !solute grid totals (kg)
      !note: sol_grid_mbef, sol_grid_maft are module variables (shared with output_aa)
      real :: sol_grid_rech
      real :: sol_grid_gwsw
      real :: sol_grid_swgw
      real :: sol_grid_satx
      real :: sol_grid_advn
      real :: sol_grid_disp
      real :: sol_grid_rcti
      real :: sol_grid_rcto
      real :: sol_grid_minl
      real :: sol_grid_sorb
      real :: sol_grid_ppag
      real :: sol_grid_ppex
      real :: sol_grid_tile
      real :: sol_grid_soil
      real :: sol_grid_resv
      real :: sol_grid_wetl
      real :: sol_grid_canl
      real :: sol_grid_fpln
      real :: sol_grid_pond
      integer :: count


      !compute average depth to water table across watershed
      sum = 0.
      count = 0
      do i=1,ncell
        if(gw_state(i)%stat == 1) then
          sum = sum + (gw_state(i)%elev - gw_state(i)%head)
          count = count + 1
        endif
      enddo
      depth_wt_avg = sum / count

      !compute average groundwater temperature
      temp_avg = 0.
      if(gw_heat_flag == 1) then
        sum = 0.
        count = 0
        do i=1,ncell
          if(gw_state(i)%stat == 1) then
            sum = sum + gwheat_state(i)%temp
            count = count + 1
          endif
        enddo
        temp_avg = sum / count
      endif

      !print out head values and solute concentration values for observation cells
      do k=1,gw_num_obs_wells
        gw_obs_head(k) = gw_state(gw_obs_cells(k))%head
        if(gw_heat_flag == 1) then
          gw_obs_temp(k) = gwheat_state(gw_obs_cells(k))%temp
        endif
        if(gw_solute_flag == 1) then
          do s=1,gw_nsolute
            gw_obs_solute(k,s) = gwsol_state(gw_obs_cells(k))%solute(s)%conc
          enddo
        endif
      enddo
      if(gwflag_obs == 1) then
        do k=1,gw_num_obs_wells
          obs_temp = -99.
          obs_no3 = -99.
          obs_p = -99.
          if(gw_heat_flag == 1) obs_temp = gw_obs_temp(k)
          if(gw_solute_flag == 1) then
            obs_no3 = gw_obs_solute(k,1)
            obs_p = gw_obs_solute(k,2)
          endif
          write(obs_name,'(a4,i4.4)') 'obs_',k
          write(out_gwobs,8102) time%day,time%mo,time%day_mo,time%yrc, &
            k,gw_obs_cells(k),obs_name, &
            gw_obs_head(k), &
            gw_state(gw_obs_cells(k))%elev - gw_obs_head(k), &
            obs_temp, obs_no3, obs_p
        enddo
      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 groundwater heat at the end of the day
      if(gw_heat_flag == 1) then
        do i=1,ncell
          if(gw_state(i)%stat == 1) then
            gwheat_state(i)%haft = gwheat_state(i)%stor
          endif
        enddo
      endif
      !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
              gwsol_state(i)%solute(s)%maft = gwsol_state(i)%solute(s)%mass !g
            enddo
          endif
        enddo
      endif


      !--- groundwater water balance ---

      !calculate values for cell groups (if specified)
      if(gw_group_flag == 1) then
        do i=1,gw_wb_grp_num
          vbef_grp = 0.
          vaft_grp = 0.
          rech_grp = 0.
          gwet_grp = 0.
          gwsw_grp = 0.
          swgw_grp = 0.
          satx_grp = 0.
          soil_grp = 0.
          latl_grp = 0.
          bndr_grp = 0.
          ppag_grp = 0.
          ppdf_grp = 0.
          ppex_grp = 0.
          tile_grp = 0.
          resv_grp = 0.
          wetl_grp = 0.
          canl_grp = 0.
          fpln_grp = 0.
          pond_grp = 0.
          phyt_grp = 0.
          do j=1,gw_wb_grp_ncell(i)
            cell_id = gw_wb_grp_cells(i,j)
            if(cell_id > 0) then
            sum = 0.
            if(gw_state(cell_id)%stat == 1) then
              vbef_grp = vbef_grp + gw_state(cell_id)%vbef
              vaft_grp = vaft_grp + gw_state(cell_id)%vaft
              rech_grp = rech_grp + gw_hyd_ss(cell_id)%rech
              gwet_grp = gwet_grp + gw_hyd_ss(cell_id)%gwet
              gwsw_grp = gwsw_grp + gw_hyd_ss(cell_id)%gwsw
              swgw_grp = swgw_grp + gw_hyd_ss(cell_id)%swgw
              satx_grp = satx_grp + gw_hyd_ss(cell_id)%satx
              soil_grp = soil_grp + gw_hyd_ss(cell_id)%soil
              latl_grp = latl_grp + gw_hyd_ss(cell_id)%latl
              bndr_grp = bndr_grp + gw_hyd_ss(cell_id)%bndr
              ppag_grp = ppag_grp + gw_hyd_ss(cell_id)%ppag
              ppdf_grp = ppdf_grp + gw_hyd_ss(cell_id)%ppdf
              ppex_grp = ppex_grp + gw_hyd_ss(cell_id)%ppex
              tile_grp = tile_grp + gw_hyd_ss(cell_id)%tile
              resv_grp = resv_grp + gw_hyd_ss(cell_id)%resv
              wetl_grp = wetl_grp + gw_hyd_ss(cell_id)%wetl
              canl_grp = canl_grp + gw_hyd_ss(cell_id)%canl
              fpln_grp = fpln_grp + gw_hyd_ss(cell_id)%fpln
              pond_grp = pond_grp + gw_hyd_ss(cell_id)%pond
              phyt_grp = phyt_grp + gw_hyd_ss(cell_id)%phyt
              sum = sum + (gw_state(cell_id)%elev - gw_state(cell_id)%head)
            endif
            endif
          enddo
          depth_wt_avg_grp = sum / gw_wb_grp_ncell(i)
          !calculate mass error
          mass_error = 0.
          mass_error = (1-((vbef_grp + rech_grp + gwet_grp + gwsw_grp + swgw_grp + satx_grp + soil_grp + &
                            latl_grp + bndr_grp + ppag_grp + ppex_grp + tile_grp + resv_grp + wetl_grp + &
                            canl_grp + fpln_grp + pond_grp + phyt_grp) &
                           /vaft_grp)) * 100
          !print out daily information
          write(out_gwbal_grp+i,8100) time%day,time%mo,time%day_mo, &
            time%yrc,"       1","       1",bsn%name,gw_time_step, &
            vbef_grp,vaft_grp,rech_grp,gwet_grp,gwsw_grp,swgw_grp, &
            satx_grp,soil_grp,latl_grp,bndr_grp,ppag_grp,ppex_grp, &
            tile_grp,resv_grp,wetl_grp,canl_grp,fpln_grp,pond_grp,phyt_grp, &
            mass_error,depth_wt_avg_grp,ppdf_grp
        enddo !go to next cell group
      endif

      !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.
      pond_grid = 0.
      phyt_grid = 0.
      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_hyd_ss(i)%rech
          gwet_grid = gwet_grid + gw_hyd_ss(i)%gwet
          gwsw_grid = gwsw_grid + gw_hyd_ss(i)%gwsw
          swgw_grid = swgw_grid + gw_hyd_ss(i)%swgw
          satx_grid = satx_grid + gw_hyd_ss(i)%satx
          soil_grid = soil_grid + gw_hyd_ss(i)%soil
          latl_grid = latl_grid + gw_hyd_ss(i)%latl
          bndr_grid = bndr_grid + gw_hyd_ss(i)%bndr
          ppag_grid = ppag_grid + gw_hyd_ss(i)%ppag
          ppdf_grid = ppdf_grid + gw_hyd_ss(i)%ppdf
          ppex_grid = ppex_grid + gw_hyd_ss(i)%ppex
          tile_grid = tile_grid + gw_hyd_ss(i)%tile
          resv_grid = resv_grid + gw_hyd_ss(i)%resv
          wetl_grid = wetl_grid + gw_hyd_ss(i)%wetl
          canl_grid = canl_grid + gw_hyd_ss(i)%canl
          fpln_grid = fpln_grid + gw_hyd_ss(i)%fpln
          pond_grid = pond_grid + gw_hyd_ss(i)%pond
          phyt_grid = phyt_grid + gw_hyd_ss(i)%phyt
        endif
      enddo
      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 + pond_grid + phyt_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.
      pond_grid = (pond_grid / (bsn%area_tot_ha*10000.)) * 1000.
      phyt_grid = (phyt_grid / (bsn%area_tot_ha*10000.)) * 1000.
      frac_sat = real(satx_count) / real(num_active)
      if(gwflag_day == 1) then
        write(out_gwbal,8100) time%day,time%mo,time%day_mo,time%yrc, &
          "       1","       1",bsn%name,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,pond_grid, &
          phyt_grid,mass_error,frac_sat,depth_wt_avg,ppdf_grid
      endif

      !add daily water balance volumes to monthly values
      gw_hyd_grid_mo%chng = gw_hyd_grid_mo%chng + (vaft_grid-vbef_grid)
      gw_hyd_grid_mo%rech = gw_hyd_grid_mo%rech + rech_grid
      gw_hyd_grid_mo%gwet = gw_hyd_grid_mo%gwet + gwet_grid
      gw_hyd_grid_mo%gwsw = gw_hyd_grid_mo%gwsw + gwsw_grid
      gw_hyd_grid_mo%swgw = gw_hyd_grid_mo%swgw + swgw_grid
      gw_hyd_grid_mo%satx = gw_hyd_grid_mo%satx + satx_grid
      gw_hyd_grid_mo%soil = gw_hyd_grid_mo%soil + soil_grid
      gw_hyd_grid_mo%latl = gw_hyd_grid_mo%latl + latl_grid
      gw_hyd_grid_mo%bndr = gw_hyd_grid_mo%bndr + bndr_grid
      gw_hyd_grid_mo%ppag = gw_hyd_grid_mo%ppag + ppag_grid
      gw_hyd_grid_mo%ppdf = gw_hyd_grid_mo%ppdf + ppdf_grid
      gw_hyd_grid_mo%ppex = gw_hyd_grid_mo%ppex + ppex_grid
      gw_hyd_grid_mo%tile = gw_hyd_grid_mo%tile + tile_grid
      gw_hyd_grid_mo%resv = gw_hyd_grid_mo%resv + resv_grid
      gw_hyd_grid_mo%wetl = gw_hyd_grid_mo%wetl + wetl_grid
      gw_hyd_grid_mo%canl = gw_hyd_grid_mo%canl + canl_grid
      gw_hyd_grid_mo%fpln = gw_hyd_grid_mo%fpln + fpln_grid
      gw_hyd_grid_mo%pond = gw_hyd_grid_mo%pond + pond_grid
      gw_hyd_grid_mo%phyt = gw_hyd_grid_mo%phyt + phyt_grid
      !add daily water balance volumes to yearly values
      gw_hyd_grid_yr%chng = gw_hyd_grid_yr%chng + (vaft_grid-vbef_grid)
      gw_hyd_grid_yr%rech = gw_hyd_grid_yr%rech + rech_grid
      gw_hyd_grid_yr%gwet = gw_hyd_grid_yr%gwet + gwet_grid
      gw_hyd_grid_yr%gwsw = gw_hyd_grid_yr%gwsw + gwsw_grid
      gw_hyd_grid_yr%swgw = gw_hyd_grid_yr%swgw + swgw_grid
      gw_hyd_grid_yr%satx = gw_hyd_grid_yr%satx + satx_grid
      gw_hyd_grid_yr%soil = gw_hyd_grid_yr%soil + soil_grid
      gw_hyd_grid_yr%latl = gw_hyd_grid_yr%latl + latl_grid
      gw_hyd_grid_yr%bndr = gw_hyd_grid_yr%bndr + bndr_grid
      gw_hyd_grid_yr%ppag = gw_hyd_grid_yr%ppag + ppag_grid
      gw_hyd_grid_yr%ppdf = gw_hyd_grid_yr%ppdf + ppdf_grid
      gw_hyd_grid_yr%ppex = gw_hyd_grid_yr%ppex + ppex_grid
      gw_hyd_grid_yr%tile = gw_hyd_grid_yr%tile + tile_grid
      gw_hyd_grid_yr%resv = gw_hyd_grid_yr%resv + resv_grid
      gw_hyd_grid_yr%wetl = gw_hyd_grid_yr%wetl + wetl_grid
      gw_hyd_grid_yr%canl = gw_hyd_grid_yr%canl + canl_grid
      gw_hyd_grid_yr%fpln = gw_hyd_grid_yr%fpln + fpln_grid
      gw_hyd_grid_yr%pond = gw_hyd_grid_yr%pond + pond_grid
      gw_hyd_grid_yr%phyt = gw_hyd_grid_yr%phyt + phyt_grid
      !add daily water balance volumes to total values
      gw_hyd_grid_aa%chng = gw_hyd_grid_aa%chng + (vaft_grid-vbef_grid)
      gw_hyd_grid_aa%rech = gw_hyd_grid_aa%rech + rech_grid
      gw_hyd_grid_aa%gwet = gw_hyd_grid_aa%gwet + gwet_grid
      gw_hyd_grid_aa%gwsw = gw_hyd_grid_aa%gwsw + gwsw_grid
      gw_hyd_grid_aa%swgw = gw_hyd_grid_aa%swgw + swgw_grid
      gw_hyd_grid_aa%satx = gw_hyd_grid_aa%satx + satx_grid
      gw_hyd_grid_aa%soil = gw_hyd_grid_aa%soil + soil_grid
      gw_hyd_grid_aa%latl = gw_hyd_grid_aa%latl + latl_grid
      gw_hyd_grid_aa%bndr = gw_hyd_grid_aa%bndr + bndr_grid
      gw_hyd_grid_aa%ppag = gw_hyd_grid_aa%ppag + ppag_grid
      gw_hyd_grid_aa%ppdf = gw_hyd_grid_aa%ppdf + ppdf_grid
      gw_hyd_grid_aa%ppex = gw_hyd_grid_aa%ppex + ppex_grid
      gw_hyd_grid_aa%tile = gw_hyd_grid_aa%tile + tile_grid
      gw_hyd_grid_aa%resv = gw_hyd_grid_aa%resv + resv_grid
      gw_hyd_grid_aa%wetl = gw_hyd_grid_aa%wetl + wetl_grid
      gw_hyd_grid_aa%canl = gw_hyd_grid_aa%canl + canl_grid
      gw_hyd_grid_aa%fpln = gw_hyd_grid_aa%fpln + fpln_grid
      gw_hyd_grid_aa%pond = gw_hyd_grid_aa%pond + pond_grid
      gw_hyd_grid_aa%phyt = gw_hyd_grid_aa%phyt + phyt_grid

      !groundwater heat balance output
      if(gw_heat_flag == 1) then
        !calculate values for entire grid (all cells)
        heat_hbef_grid = 0.
        heat_haft_grid = 0.
        heat_rech_grid = 0.
        heat_gwet_grid = 0.
        heat_gwsw_grid = 0.
        heat_swgw_grid = 0.
        heat_satx_grid = 0.
        heat_soil_grid = 0.
        heat_latl_grid = 0.
        heat_disp_grid = 0.
        heat_bndr_grid = 0.
        heat_ppag_grid = 0.
        heat_ppex_grid = 0.
        heat_tile_grid = 0.
        heat_resv_grid = 0.
        heat_wetl_grid = 0.
        heat_canl_grid = 0.
        heat_fpln_grid = 0.
        heat_pond_grid = 0.
        do i=1,ncell
          if(gw_state(i)%stat == 1) then
            heat_hbef_grid = heat_hbef_grid + gwheat_state(i)%hbef
            heat_haft_grid = heat_haft_grid + gwheat_state(i)%haft
            heat_rech_grid = heat_rech_grid + gw_heat_ss(i)%rech
            heat_gwet_grid = heat_gwet_grid + gw_heat_ss(i)%gwet
            heat_gwsw_grid = heat_gwsw_grid + gw_heat_ss(i)%gwsw
            heat_swgw_grid = heat_swgw_grid + gw_heat_ss(i)%swgw
            heat_satx_grid = heat_satx_grid + gw_heat_ss(i)%satx
            heat_soil_grid = heat_soil_grid + gw_heat_ss(i)%soil
            heat_latl_grid = heat_latl_grid + gw_heat_ss(i)%latl
            heat_disp_grid = heat_disp_grid + gw_heat_ss(i)%disp
            heat_bndr_grid = heat_bndr_grid + gw_heat_ss(i)%bndr
            heat_ppag_grid = heat_ppag_grid + gw_heat_ss(i)%ppag
            heat_ppex_grid = heat_ppex_grid + gw_heat_ss(i)%ppex
            heat_tile_grid = heat_tile_grid + gw_heat_ss(i)%tile
            heat_resv_grid = heat_resv_grid + gw_heat_ss(i)%resv
            heat_wetl_grid = heat_wetl_grid + gw_heat_ss(i)%wetl
            heat_canl_grid = heat_canl_grid + gw_heat_ss(i)%canl
            heat_fpln_grid = heat_fpln_grid + gw_heat_ss(i)%fpln
            heat_pond_grid = heat_pond_grid + gw_heat_ss(i)%pond
          endif
        enddo
        mass_error = (1-((heat_hbef_grid + heat_rech_grid + heat_gwet_grid + heat_gwsw_grid + &
                          heat_swgw_grid + heat_satx_grid + heat_soil_grid + heat_latl_grid + &
                          heat_disp_grid + heat_bndr_grid + heat_ppag_grid + heat_ppex_grid + &
                          heat_tile_grid + heat_resv_grid + heat_wetl_grid + heat_canl_grid + &
                          heat_fpln_grid + heat_pond_grid) / heat_haft_grid)) * 100
        !print out daily information
        !first, divide by 1,000,000 to get MJ
        heat_hbef_grid = heat_hbef_grid / 1000000.
        heat_haft_grid = heat_haft_grid / 1000000.
        heat_rech_grid = heat_rech_grid / 1000000.
        heat_gwet_grid = heat_gwet_grid / 1000000.
        heat_gwsw_grid = heat_gwsw_grid / 1000000.
        heat_swgw_grid = heat_swgw_grid / 1000000.
        heat_satx_grid = heat_satx_grid / 1000000.
        heat_soil_grid = heat_soil_grid / 1000000.
        heat_latl_grid = heat_latl_grid / 1000000.
        heat_disp_grid = heat_disp_grid / 1000000.
        heat_bndr_grid = heat_bndr_grid / 1000000.
        heat_ppag_grid = heat_ppag_grid / 1000000.
        heat_ppex_grid = heat_ppex_grid / 1000000.
        heat_tile_grid = heat_tile_grid / 1000000.
        heat_resv_grid = heat_resv_grid / 1000000.
        heat_wetl_grid = heat_wetl_grid / 1000000.
        heat_canl_grid = heat_canl_grid / 1000000.
        heat_fpln_grid = heat_fpln_grid / 1000000.
        heat_pond_grid = heat_pond_grid / 1000000.
        if(gwflag_day == 1) then
          write(out_heatbal_dy,8100) time%day,time%mo,time%day_mo, &
            time%yrc,"       1","       1",bsn%name,gw_time_step, &
            heat_hbef_grid,heat_haft_grid,heat_rech_grid,heat_gwet_grid, &
            heat_gwsw_grid,heat_swgw_grid,heat_satx_grid,heat_soil_grid, &
            heat_latl_grid,heat_disp_grid,heat_bndr_grid,heat_ppag_grid, &
            heat_ppex_grid,heat_tile_grid,heat_resv_grid,heat_wetl_grid, &
            heat_canl_grid,heat_fpln_grid,heat_pond_grid, &
            mass_error,temp_avg
        endif
        !add daily heat fluxes to yearly values
        gw_heat_grid_yr%chng = gw_heat_grid_yr%chng + (heat_haft_grid-heat_hbef_grid)
        gw_heat_grid_yr%rech = gw_heat_grid_yr%rech + heat_rech_grid
        gw_heat_grid_yr%gwet = gw_heat_grid_yr%gwet + heat_gwet_grid
        gw_heat_grid_yr%gwsw = gw_heat_grid_yr%gwsw + heat_gwsw_grid
        gw_heat_grid_yr%swgw = gw_heat_grid_yr%swgw + heat_swgw_grid
        gw_heat_grid_yr%satx = gw_heat_grid_yr%satx + heat_satx_grid
        gw_heat_grid_yr%soil = gw_heat_grid_yr%soil + heat_soil_grid
        gw_heat_grid_yr%latl = gw_heat_grid_yr%latl + heat_latl_grid
        gw_heat_grid_yr%disp = gw_heat_grid_yr%disp + heat_disp_grid
        gw_heat_grid_yr%bndr = gw_heat_grid_yr%bndr + heat_bndr_grid
        gw_heat_grid_yr%ppag = gw_heat_grid_yr%ppag + heat_ppag_grid
        gw_heat_grid_yr%ppex = gw_heat_grid_yr%ppex + heat_ppex_grid
        gw_heat_grid_yr%tile = gw_heat_grid_yr%tile + heat_tile_grid
        gw_heat_grid_yr%resv = gw_heat_grid_yr%resv + heat_resv_grid
        gw_heat_grid_yr%wetl = gw_heat_grid_yr%wetl + heat_wetl_grid
        gw_heat_grid_yr%canl = gw_heat_grid_yr%canl + heat_canl_grid
        gw_heat_grid_yr%fpln = gw_heat_grid_yr%fpln + heat_fpln_grid
        gw_heat_grid_yr%pond = gw_heat_grid_yr%pond + heat_pond_grid
        !add daily heat fluxes to total values
        gw_heat_grid_aa%chng = gw_heat_grid_aa%chng + (heat_haft_grid-heat_hbef_grid)
        gw_heat_grid_aa%rech = gw_heat_grid_aa%rech + heat_rech_grid
        gw_heat_grid_aa%gwet = gw_heat_grid_aa%gwet + heat_gwet_grid
        gw_heat_grid_aa%gwsw = gw_heat_grid_aa%gwsw + heat_gwsw_grid
        gw_heat_grid_aa%swgw = gw_heat_grid_aa%swgw + heat_swgw_grid
        gw_heat_grid_aa%satx = gw_heat_grid_aa%satx + heat_satx_grid
        gw_heat_grid_aa%soil = gw_heat_grid_aa%soil + heat_soil_grid
        gw_heat_grid_aa%latl = gw_heat_grid_aa%latl + heat_latl_grid
        gw_heat_grid_aa%disp = gw_heat_grid_aa%disp + heat_disp_grid
        gw_heat_grid_aa%bndr = gw_heat_grid_aa%bndr + heat_bndr_grid
        gw_heat_grid_aa%ppag = gw_heat_grid_aa%ppag + heat_ppag_grid
        gw_heat_grid_aa%ppex = gw_heat_grid_aa%ppex + heat_ppex_grid
        gw_heat_grid_aa%tile = gw_heat_grid_aa%tile + heat_tile_grid
        gw_heat_grid_aa%resv = gw_heat_grid_aa%resv + heat_resv_grid
        gw_heat_grid_aa%wetl = gw_heat_grid_aa%wetl + heat_wetl_grid
        gw_heat_grid_aa%canl = gw_heat_grid_aa%canl + heat_canl_grid
        gw_heat_grid_aa%fpln = gw_heat_grid_aa%fpln + heat_fpln_grid
        gw_heat_grid_aa%pond = gw_heat_grid_aa%pond + heat_pond_grid
      endif

      !solute mass balance output
      if(gw_solute_flag == 1) then

        !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_resv = 0.
          sol_grid_wetl = 0.
          sol_grid_canl = 0.
          sol_grid_fpln = 0.
          sol_grid_pond = 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.)
              sol_grid_pond = sol_grid_pond + (gwsol_ss(i)%solute(s)%pond / 1000.)
            endif
          enddo
          sol_grid_sorb = sol_grid_sorb * (-1) !leaving groundwater (sorbing to aquifer material)
          !calculate mass error
          mass_error = 0.
          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_pond) / sol_grid_maft)) * 100
          endif
          !print out daily values for the solute
          if(gwflag_day == 1) then
            write(out_solbal_dy+s,8100) time%day,time%mo,time%day_mo, &
              time%yrc,"       1","       1",bsn%name,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,sol_grid_pond, &
              mass_error
          endif
          !add grid values to monthly, yearly, and total mass values
          !monthly (kg)
          sol_grid_chng_mo(s) = sol_grid_chng_mo(s) + (sol_grid_maft-sol_grid_mbef)
          sol_grid_rech_mo(s) = sol_grid_rech_mo(s) + sol_grid_rech
          sol_grid_gwsw_mo(s) = sol_grid_gwsw_mo(s) + sol_grid_gwsw
          sol_grid_swgw_mo(s) = sol_grid_swgw_mo(s) + sol_grid_swgw
          sol_grid_satx_mo(s) = sol_grid_satx_mo(s) + sol_grid_satx
          sol_grid_advn_mo(s) = sol_grid_advn_mo(s) + sol_grid_advn
          sol_grid_disp_mo(s) = sol_grid_disp_mo(s) + sol_grid_disp
          sol_grid_rcti_mo(s) = sol_grid_rcti_mo(s) + sol_grid_rcti
          sol_grid_rcto_mo(s) = sol_grid_rcto_mo(s) + sol_grid_rcto
          sol_grid_minl_mo(s) = sol_grid_minl_mo(s) + sol_grid_minl
          sol_grid_sorb_mo(s) = sol_grid_sorb_mo(s) + sol_grid_sorb
          sol_grid_ppag_mo(s) = sol_grid_ppag_mo(s) + sol_grid_ppag
          sol_grid_ppex_mo(s) = sol_grid_ppex_mo(s) + sol_grid_ppex
          sol_grid_tile_mo(s) = sol_grid_tile_mo(s) + sol_grid_tile
          sol_grid_soil_mo(s) = sol_grid_soil_mo(s) + sol_grid_soil
          sol_grid_resv_mo(s) = sol_grid_resv_mo(s) + sol_grid_resv
          sol_grid_wetl_mo(s) = sol_grid_wetl_mo(s) + sol_grid_wetl
          sol_grid_canl_mo(s) = sol_grid_canl_mo(s) + sol_grid_canl
          sol_grid_fpln_mo(s) = sol_grid_fpln_mo(s) + sol_grid_fpln
          sol_grid_pond_mo(s) = sol_grid_pond_mo(s) + sol_grid_pond
          !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
          sol_grid_pond_yr(s) = sol_grid_pond_yr(s) + sol_grid_pond
          !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
          sol_grid_pond_tt(s) = sol_grid_pond_tt(s) + sol_grid_pond
        enddo !go to next solute
      endif !check for solute transport


      !--- cell-level daily long-format output ---
      if(gwflag_day == 1) then
        do i=1,ncell
          if(gw_state(i)%stat == 1) then
            gis_id = cell_gis_id(i)
            write(out_gwcell_day,140) time%day, time%mo, time%day_mo, &
              time%yrc, i, gis_id, &
              cell_name(i), &
              gw_state(i)%head, &
              gw_state(i)%elev - gw_state(i)%head, &
              gw_hyd_ss(i)%rech, gw_hyd_ss(i)%gwet, &
              gw_hyd_ss(i)%gwsw, gw_hyd_ss(i)%swgw, &
              gw_hyd_ss(i)%satx, gw_hyd_ss(i)%soil, &
              gw_hyd_ss(i)%latl, &
              gw_hyd_ss(i)%ppag, gw_hyd_ss(i)%ppex, &
              gw_hyd_ss(i)%tile, gw_hyd_ss(i)%resv, &
              gw_hyd_ss(i)%wetl, gw_hyd_ss(i)%fpln, &
              gw_hyd_ss(i)%canl, gw_hyd_ss(i)%pond, &
              gw_hyd_ss(i)%phyt
          endif
        enddo
      endif


      !--- HRU pumping daily long-format output ---
      if(gwflag_pump == 1) then
        do i=1,sp_ob%hru
          if(hru_pump(i) > 0.) then
            iob = sp_ob1%hru + i - 1
            write(out_hru_pump_day,8101) time%day,time%mo,time%day_mo, &
              time%yrc,i,ob(iob)%gis_id,ob(iob)%name,hru_pump(i)
          endif
        enddo
      endif

      !format statements (subroutine-local)
119   format(i8,i8,i8,1000(f12.3))
130   format(i8,i8,1000(e13.4))
140   format(i8,i6,i6,i8,i8,i10,4x,a12,2f13.3,16e13.4)
8100  format(4i6,2a,2x,a16,f10.3,50e13.4)
8101  format(4i6,2i8,a18,e13.4)
8102  format(4i6,2i8,a18,5e13.4)

      return
      end subroutine gwflow_output_day


      subroutine gwflow_output_mon

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine writes monthly gwflow output in SWAT+ long format:
!!    one row per active cell with average head, wtdepth, and average
!!    daily flow rates; basin-level water/heat/solute balance; HRU pumping

      use gwflow_module
      use hydrograph_module
      use sd_channel_module
      use time_module
      use basin_module, only : bsn

      implicit none

      integer :: i, j, k, s, iob
      integer :: gis_id
      real :: wtdepth
      real :: day_mo_r
      real :: obs_temp, obs_no3, obs_p
      character(len=16) :: obs_name

      if (time%end_mo == 1) then

      day_mo_r = real(time%day_mo)

        !monthly average groundwater head -----------------------------------------------
        do i=1,ncell
          gw_state(i)%hdmo = gw_state(i)%hdmo / day_mo_r
        enddo
        !monthly average groundwater temperature ----------------------------------------
        if(gw_heat_flag == 1) then
          do i=1,ncell
            gwheat_state(i)%tpmo = gwheat_state(i)%tpmo / day_mo_r
          enddo
          !zero out for next month
          do i=1,ncell
            gwheat_state(i)%tpmo = 0.
          enddo
        endif

        !monthly average solute concentration -------------------------------------------
        if(gw_solute_flag == 1) then
          do s=1,gw_nsolute
            do i=1,ncell
              gwsol_state(i)%solute(s)%cnmo = gwsol_state(i)%solute(s)%cnmo / day_mo_r
            enddo
            !zero out for next month
            do i=1,ncell
              gwsol_state(i)%solute(s)%cnmo = 0.
            enddo
          enddo !next solute
        endif

        !pumping (irrigation) (for HRUs) -- monthly long-format output ----------------
        if(gwflag_pump == 1) then
          do i=1,sp_ob%hru
            if(hru_pump_mo(i) > 0.) then
              iob = sp_ob1%hru + i - 1
              write(out_hru_pump_mo,8101) time%day,time%mo,time%day_mo, &
                time%yrc,i,ob(iob)%gis_id,ob(iob)%name,hru_pump_mo(i)
            endif
          enddo
        endif
        hru_pump_mo = 0.

        !compute average daily groundwater fluxes (m3/day) for the month ----------------
        do i=1,ncell
          gw_hyd_ss_mo(i)%rech = gw_hyd_ss_mo(i)%rech / day_mo_r
          gw_hyd_ss_mo(i)%gwet = gw_hyd_ss_mo(i)%gwet / day_mo_r
          gw_hyd_ss_mo(i)%gwsw = gw_hyd_ss_mo(i)%gwsw / day_mo_r
          gw_hyd_ss_mo(i)%swgw = gw_hyd_ss_mo(i)%swgw / day_mo_r
          gw_hyd_ss_mo(i)%satx = gw_hyd_ss_mo(i)%satx / day_mo_r
          gw_hyd_ss_mo(i)%soil = gw_hyd_ss_mo(i)%soil / day_mo_r
          gw_hyd_ss_mo(i)%latl = gw_hyd_ss_mo(i)%latl / day_mo_r
          gw_hyd_ss_mo(i)%bndr = gw_hyd_ss_mo(i)%bndr / day_mo_r
          gw_hyd_ss_mo(i)%ppag = gw_hyd_ss_mo(i)%ppag / day_mo_r
          gw_hyd_ss_mo(i)%ppdf = gw_hyd_ss_mo(i)%ppdf / day_mo_r
          gw_hyd_ss_mo(i)%ppex = gw_hyd_ss_mo(i)%ppex / day_mo_r
          gw_hyd_ss_mo(i)%tile = gw_hyd_ss_mo(i)%tile / day_mo_r
          gw_hyd_ss_mo(i)%resv = gw_hyd_ss_mo(i)%resv / day_mo_r
          gw_hyd_ss_mo(i)%wetl = gw_hyd_ss_mo(i)%wetl / day_mo_r
          gw_hyd_ss_mo(i)%canl = gw_hyd_ss_mo(i)%canl / day_mo_r
          gw_hyd_ss_mo(i)%fpln = gw_hyd_ss_mo(i)%fpln / day_mo_r
          gw_hyd_ss_mo(i)%pond = gw_hyd_ss_mo(i)%pond / day_mo_r
          gw_hyd_ss_mo(i)%phyt = gw_hyd_ss_mo(i)%phyt / day_mo_r
        enddo

        !--- cell-level monthly long-format output (one row per active cell) ---
        if(gwflag_mon == 1) then
          do i=1,ncell
            if(gw_state(i)%stat == 1) then
              gis_id = cell_gis_id(i)
              wtdepth = gw_state(i)%elev - gw_state(i)%hdmo
              write(out_gwcell_mon,140) time%day, time%mo, &
                time%day_mo, time%yrc, i, gis_id, &
                cell_name(i), &
                gw_state(i)%hdmo, wtdepth, &
                gw_hyd_ss_mo(i)%rech, gw_hyd_ss_mo(i)%gwet, &
                gw_hyd_ss_mo(i)%gwsw, gw_hyd_ss_mo(i)%swgw, &
                gw_hyd_ss_mo(i)%satx, gw_hyd_ss_mo(i)%soil, &
                gw_hyd_ss_mo(i)%latl, &
                gw_hyd_ss_mo(i)%ppag, gw_hyd_ss_mo(i)%ppex, &
                gw_hyd_ss_mo(i)%tile, gw_hyd_ss_mo(i)%resv, &
                gw_hyd_ss_mo(i)%wetl, gw_hyd_ss_mo(i)%fpln, &
                gw_hyd_ss_mo(i)%canl, gw_hyd_ss_mo(i)%pond, &
                gw_hyd_ss_mo(i)%phyt
            endif
          enddo
        endif

        !--- obs well monthly output ---
        if(gwflag_obs == 1 .and. gw_num_obs_wells > 0) then
          do k=1,gw_num_obs_wells
            i = gw_obs_cells(k)
            obs_temp = -99.; obs_no3 = -99.; obs_p = -99.
            if(gw_heat_flag == 1) obs_temp = gwheat_state(i)%tpmo
            if(gw_solute_flag == 1) then
              obs_no3 = gwsol_state(i)%solute(1)%cnmo
              obs_p = gwsol_state(i)%solute(2)%cnmo
            endif
            write(obs_name,'(a4,i4.4)') 'obs_',k
            write(out_gwobs_mon,8102) time%day,time%mo,time%day_mo, &
              time%yrc,k,gw_obs_cells(k),obs_name, &
              gw_state(i)%hdmo, gw_state(i)%elev - gw_state(i)%hdmo, &
              obs_temp, obs_no3, obs_p
          enddo
        endif

        !zero out head average and flux sums for next month
        do i=1,ncell
          gw_state(i)%hdmo = 0.
          gw_hyd_ss_mo(i)%rech = 0.
          gw_hyd_ss_mo(i)%gwet = 0.
          gw_hyd_ss_mo(i)%gwsw = 0.
          gw_hyd_ss_mo(i)%swgw = 0.
          gw_hyd_ss_mo(i)%satx = 0.
          gw_hyd_ss_mo(i)%soil = 0.
          gw_hyd_ss_mo(i)%latl = 0.
          gw_hyd_ss_mo(i)%bndr = 0.
          gw_hyd_ss_mo(i)%ppag = 0.
          gw_hyd_ss_mo(i)%ppdf = 0.
          gw_hyd_ss_mo(i)%ppex = 0.
          gw_hyd_ss_mo(i)%tile = 0.
          gw_hyd_ss_mo(i)%resv = 0.
          gw_hyd_ss_mo(i)%wetl = 0.
          gw_hyd_ss_mo(i)%canl = 0.
          gw_hyd_ss_mo(i)%fpln = 0.
          gw_hyd_ss_mo(i)%pond = 0.
          gw_hyd_ss_mo(i)%phyt = 0.
        enddo

        !solute monthly zeroing
        if(gw_solute_flag == 1) then
          do i=1,ncell
            do s=1,gw_nsolute
              gwsol_ss_sum_mo(i)%solute(s)%rech = 0.
              gwsol_ss_sum_mo(i)%solute(s)%gwsw = 0.
              gwsol_ss_sum_mo(i)%solute(s)%swgw = 0.
              gwsol_ss_sum_mo(i)%solute(s)%satx = 0.
              gwsol_ss_sum_mo(i)%solute(s)%soil = 0.
              gwsol_ss_sum_mo(i)%solute(s)%ppag = 0.
              gwsol_ss_sum_mo(i)%solute(s)%ppex = 0.
              gwsol_ss_sum_mo(i)%solute(s)%tile = 0.
              gwsol_ss_sum_mo(i)%solute(s)%resv = 0.
              gwsol_ss_sum_mo(i)%solute(s)%wetl = 0.
              gwsol_ss_sum_mo(i)%solute(s)%canl = 0.
              gwsol_ss_sum_mo(i)%solute(s)%fpln = 0.
              gwsol_ss_sum_mo(i)%solute(s)%pond = 0.
              gwsol_ss_sum_mo(i)%solute(s)%advn = 0.
              gwsol_ss_sum_mo(i)%solute(s)%disp = 0.
              gwsol_ss_sum_mo(i)%solute(s)%rcti = 0.
              gwsol_ss_sum_mo(i)%solute(s)%rcto = 0.
              gwsol_ss_sum_mo(i)%solute(s)%minl = 0.
              gwsol_ss_sum_mo(i)%solute(s)%sorb = 0.
            enddo
          enddo
        endif

        !monthly groundwater balance (basin) --------------------------------------------
        if(gwflag_mon == 1) then
          write(out_gwbal_mon,8100) time%day,time%mo,time%day_mo, &
            time%yrc,"       1","       1",bsn%name, &
            gw_hyd_grid_mo%chng,gw_hyd_grid_mo%rech,gw_hyd_grid_mo%gwet, &
            gw_hyd_grid_mo%gwsw,gw_hyd_grid_mo%swgw, &
            gw_hyd_grid_mo%satx,gw_hyd_grid_mo%soil,gw_hyd_grid_mo%latl, &
            gw_hyd_grid_mo%bndr,gw_hyd_grid_mo%ppag, &
            gw_hyd_grid_mo%ppex,gw_hyd_grid_mo%tile,gw_hyd_grid_mo%resv, &
            gw_hyd_grid_mo%wetl,gw_hyd_grid_mo%canl, &
            gw_hyd_grid_mo%fpln,gw_hyd_grid_mo%pond,gw_hyd_grid_mo%phyt, &
            gw_hyd_grid_mo%ppdf
        endif
        !zero out basin monthly accumulators
        gw_hyd_grid_mo%chng = 0.
        gw_hyd_grid_mo%rech = 0.
        gw_hyd_grid_mo%gwet = 0.
        gw_hyd_grid_mo%gwsw = 0.
        gw_hyd_grid_mo%swgw = 0.
        gw_hyd_grid_mo%satx = 0.
        gw_hyd_grid_mo%soil = 0.
        gw_hyd_grid_mo%latl = 0.
        gw_hyd_grid_mo%bndr = 0.
        gw_hyd_grid_mo%ppag = 0.
        gw_hyd_grid_mo%ppdf = 0.
        gw_hyd_grid_mo%ppex = 0.
        gw_hyd_grid_mo%tile = 0.
        gw_hyd_grid_mo%resv = 0.
        gw_hyd_grid_mo%wetl = 0.
        gw_hyd_grid_mo%canl = 0.
        gw_hyd_grid_mo%fpln = 0.
        gw_hyd_grid_mo%pond = 0.
        gw_hyd_grid_mo%phyt = 0.

        !solute mass values (basin) -----------------------------------------------------
        if(gw_solute_flag == 1) then
          do s=1,gw_nsolute
            if(gwflag_mon == 1) then
              write(out_solbal_mo+s,8100) time%day,time%mo,time%day_mo, &
                time%yrc,"       1","       1",bsn%name, &
                sol_grid_chng_mo(s),sol_grid_rech_mo(s), &
                sol_grid_gwsw_mo(s),sol_grid_swgw_mo(s), &
                sol_grid_satx_mo(s), &
                sol_grid_soil_mo(s),sol_grid_advn_mo(s), &
                sol_grid_disp_mo(s), &
                sol_grid_rcti_mo(s),sol_grid_rcto_mo(s), &
                sol_grid_minl_mo(s),sol_grid_sorb_mo(s), &
                sol_grid_ppag_mo(s),sol_grid_ppex_mo(s), &
                sol_grid_tile_mo(s),sol_grid_resv_mo(s), &
                sol_grid_wetl_mo(s), &
                sol_grid_canl_mo(s),sol_grid_fpln_mo(s), &
                sol_grid_pond_mo(s)
            endif
            sol_grid_chng_mo(s) = 0.
            sol_grid_rech_mo(s) = 0.
            sol_grid_gwsw_mo(s) = 0.
            sol_grid_swgw_mo(s) = 0.
            sol_grid_satx_mo(s) = 0.
            sol_grid_soil_mo(s) = 0.
            sol_grid_advn_mo(s) = 0.
            sol_grid_disp_mo(s) = 0.
            sol_grid_rcti_mo(s) = 0.
            sol_grid_rcto_mo(s) = 0.
            sol_grid_minl_mo(s) = 0.
            sol_grid_sorb_mo(s) = 0.
            sol_grid_ppag_mo(s) = 0.
            sol_grid_ppex_mo(s) = 0.
            sol_grid_tile_mo(s) = 0.
            sol_grid_resv_mo(s) = 0.
            sol_grid_wetl_mo(s) = 0.
            sol_grid_canl_mo(s) = 0.
            sol_grid_fpln_mo(s) = 0.
            sol_grid_pond_mo(s) = 0.
          enddo !go to next solute
        endif

      endif !end_mo

      !format statements (subroutine-local)
140   format(i8,i6,i6,i8,i8,i10,4x,a12,2f13.3,16e13.4)
8100  format(4i6,2a,2x,a16,50e13.4)
8101  format(4i6,2i8,a18,e13.4)
8102  format(4i6,2i8,a18,5e13.4)

      return
      end subroutine gwflow_output_mon


      subroutine gwflow_output_yr

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine writes annual gwflow output in SWAT+ long format:
!!    one row per active cell with average head, wtdepth, and average
!!    daily flow rates; basin-level water/heat/solute balance; HRU pumping

      use gwflow_module
      use hydrograph_module
      use sd_channel_module
      use time_module
      use basin_module, only : bsn

      implicit none

      integer :: i, j, k, s, iob
      integer :: gis_id
      real :: wtdepth
      real :: day_yr_r
      real :: obs_temp, obs_no3, obs_p
      character(len=16) :: obs_name

      !--- only execute at end of year ---
      if(time%end_yr /= 1) return

      day_yr_r = real(time%day_end_yr)

      !annual average groundwater head ------------------------------------------------
      do i=1,ncell
        gw_state(i)%hdyr = gw_state(i)%hdyr / day_yr_r
      enddo

      !annual average groundwater temperature -----------------------------------------
      if(gw_heat_flag == 1) then
        do i=1,ncell
          gwheat_state(i)%tpyr = gwheat_state(i)%tpyr / day_yr_r
        enddo
        !zero out for next year
        do i=1,ncell
          gwheat_state(i)%tpyr = 0.
        enddo
      endif

      !annual average solute concentration --------------------------------------------
      if(gw_solute_flag == 1) then
        do s=1,gw_nsolute
          do i=1,ncell
            gwsol_state(i)%solute(s)%cnyr = gwsol_state(i)%solute(s)%cnyr / day_yr_r
          enddo
          !zero out for next year
          do i=1,ncell
            gwsol_state(i)%solute(s)%cnyr = 0.
          enddo
        enddo !next solute
      endif !check for solutes

      !compute average daily groundwater fluxes (m3/day) for the year -----------------
      do i=1,ncell
        gw_hyd_ss_yr(i)%rech = gw_hyd_ss_yr(i)%rech / day_yr_r
        gw_hyd_ss_yr(i)%gwet = gw_hyd_ss_yr(i)%gwet / day_yr_r
        gw_hyd_ss_yr(i)%gwsw = gw_hyd_ss_yr(i)%gwsw / day_yr_r
        gw_hyd_ss_yr(i)%swgw = gw_hyd_ss_yr(i)%swgw / day_yr_r
        gw_hyd_ss_yr(i)%satx = gw_hyd_ss_yr(i)%satx / day_yr_r
        gw_hyd_ss_yr(i)%soil = gw_hyd_ss_yr(i)%soil / day_yr_r
        gw_hyd_ss_yr(i)%latl = gw_hyd_ss_yr(i)%latl / day_yr_r
        gw_hyd_ss_yr(i)%bndr = gw_hyd_ss_yr(i)%bndr / day_yr_r
        gw_hyd_ss_yr(i)%ppag = gw_hyd_ss_yr(i)%ppag / day_yr_r
        gw_hyd_ss_yr(i)%ppdf = gw_hyd_ss_yr(i)%ppdf / day_yr_r
        gw_hyd_ss_yr(i)%ppex = gw_hyd_ss_yr(i)%ppex / day_yr_r
        gw_hyd_ss_yr(i)%tile = gw_hyd_ss_yr(i)%tile / day_yr_r
        gw_hyd_ss_yr(i)%resv = gw_hyd_ss_yr(i)%resv / day_yr_r
        gw_hyd_ss_yr(i)%wetl = gw_hyd_ss_yr(i)%wetl / day_yr_r
        gw_hyd_ss_yr(i)%canl = gw_hyd_ss_yr(i)%canl / day_yr_r
        gw_hyd_ss_yr(i)%fpln = gw_hyd_ss_yr(i)%fpln / day_yr_r
        gw_hyd_ss_yr(i)%pond = gw_hyd_ss_yr(i)%pond / day_yr_r
        gw_hyd_ss_yr(i)%phyt = gw_hyd_ss_yr(i)%phyt / day_yr_r
      enddo

      !--- cell-level annual long-format output (one row per active cell) ---
      if(gwflag_yr == 1) then
        do i=1,ncell
          if(gw_state(i)%stat == 1) then
            gis_id = cell_gis_id(i)
            wtdepth = gw_state(i)%elev - gw_state(i)%hdyr
            write(out_gwcell_yr,140) time%day, time%mo, &
              time%day_mo, time%yrc, i, gis_id, &
              cell_name(i), &
              gw_state(i)%hdyr, wtdepth, &
              gw_hyd_ss_yr(i)%rech, gw_hyd_ss_yr(i)%gwet, &
              gw_hyd_ss_yr(i)%gwsw, gw_hyd_ss_yr(i)%swgw, &
              gw_hyd_ss_yr(i)%satx, gw_hyd_ss_yr(i)%soil, &
              gw_hyd_ss_yr(i)%latl, &
              gw_hyd_ss_yr(i)%ppag, gw_hyd_ss_yr(i)%ppex, &
              gw_hyd_ss_yr(i)%tile, gw_hyd_ss_yr(i)%resv, &
              gw_hyd_ss_yr(i)%wetl, gw_hyd_ss_yr(i)%fpln, &
              gw_hyd_ss_yr(i)%canl, gw_hyd_ss_yr(i)%pond, &
              gw_hyd_ss_yr(i)%phyt
          endif
        enddo
      endif

      !accumulate yearly values into per-cell AA arrays (before zeroing)
      do i=1,ncell
        gw_head_sum_aa(i) = gw_head_sum_aa(i) + gw_state(i)%hdyr
        gw_hyd_ss_aa(i)%rech = gw_hyd_ss_aa(i)%rech + gw_hyd_ss_yr(i)%rech
        gw_hyd_ss_aa(i)%gwet = gw_hyd_ss_aa(i)%gwet + gw_hyd_ss_yr(i)%gwet
        gw_hyd_ss_aa(i)%gwsw = gw_hyd_ss_aa(i)%gwsw + gw_hyd_ss_yr(i)%gwsw
        gw_hyd_ss_aa(i)%swgw = gw_hyd_ss_aa(i)%swgw + gw_hyd_ss_yr(i)%swgw
        gw_hyd_ss_aa(i)%satx = gw_hyd_ss_aa(i)%satx + gw_hyd_ss_yr(i)%satx
        gw_hyd_ss_aa(i)%soil = gw_hyd_ss_aa(i)%soil + gw_hyd_ss_yr(i)%soil
        gw_hyd_ss_aa(i)%latl = gw_hyd_ss_aa(i)%latl + gw_hyd_ss_yr(i)%latl
        gw_hyd_ss_aa(i)%ppag = gw_hyd_ss_aa(i)%ppag + gw_hyd_ss_yr(i)%ppag
        gw_hyd_ss_aa(i)%ppex = gw_hyd_ss_aa(i)%ppex + gw_hyd_ss_yr(i)%ppex
        gw_hyd_ss_aa(i)%tile = gw_hyd_ss_aa(i)%tile + gw_hyd_ss_yr(i)%tile
        gw_hyd_ss_aa(i)%resv = gw_hyd_ss_aa(i)%resv + gw_hyd_ss_yr(i)%resv
        gw_hyd_ss_aa(i)%wetl = gw_hyd_ss_aa(i)%wetl + gw_hyd_ss_yr(i)%wetl
        gw_hyd_ss_aa(i)%fpln = gw_hyd_ss_aa(i)%fpln + gw_hyd_ss_yr(i)%fpln
        gw_hyd_ss_aa(i)%canl = gw_hyd_ss_aa(i)%canl + gw_hyd_ss_yr(i)%canl
        gw_hyd_ss_aa(i)%pond = gw_hyd_ss_aa(i)%pond + gw_hyd_ss_yr(i)%pond
        gw_hyd_ss_aa(i)%phyt = gw_hyd_ss_aa(i)%phyt + gw_hyd_ss_yr(i)%phyt
      enddo

      !--- obs well yearly output + AA accumulation ---
      if(gwflag_obs == 1 .and. gw_num_obs_wells > 0) then
        do k=1,gw_num_obs_wells
          i = gw_obs_cells(k)
          obs_temp = -99.; obs_no3 = -99.; obs_p = -99.
          if(gw_heat_flag == 1) then
            obs_temp = gwheat_state(i)%tpyr
            gw_obs_temp_aa(k) = gw_obs_temp_aa(k) + obs_temp
          endif
          if(gw_solute_flag == 1) then
            obs_no3 = gwsol_state(i)%solute(1)%cnyr
            obs_p = gwsol_state(i)%solute(2)%cnyr
            gw_obs_sol_aa(k,1) = gw_obs_sol_aa(k,1) + obs_no3
            gw_obs_sol_aa(k,2) = gw_obs_sol_aa(k,2) + obs_p
          endif
          write(obs_name,'(a4,i4.4)') 'obs_',k
          write(out_gwobs_yr,8102) time%day,time%mo,time%day_mo, &
            time%yrc,k,gw_obs_cells(k),obs_name, &
            gw_state(i)%hdyr, gw_state(i)%elev - gw_state(i)%hdyr, &
            obs_temp, obs_no3, obs_p
        enddo
      endif

      !zero out head average for next year
      do i=1,ncell
        gw_state(i)%hdyr = 0.
      enddo

      !compute average daily heat fluxes (MJ/day) for the year
      if(gw_heat_flag == 1) then
        do i=1,ncell
          gw_heat_ss_yr(i)%rech = (gw_heat_ss_yr(i)%rech/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%gwet = (gw_heat_ss_yr(i)%gwet/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%gwsw = (gw_heat_ss_yr(i)%gwsw/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%swgw = (gw_heat_ss_yr(i)%swgw/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%satx = (gw_heat_ss_yr(i)%satx/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%soil = (gw_heat_ss_yr(i)%soil/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%latl = (gw_heat_ss_yr(i)%latl/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%disp = (gw_heat_ss_yr(i)%disp/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%bndr = (gw_heat_ss_yr(i)%bndr/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%ppag = (gw_heat_ss_yr(i)%ppag/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%ppex = (gw_heat_ss_yr(i)%ppex/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%tile = (gw_heat_ss_yr(i)%tile/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%resv = (gw_heat_ss_yr(i)%resv/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%wetl = (gw_heat_ss_yr(i)%wetl/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%canl = (gw_heat_ss_yr(i)%canl/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%fpln = (gw_heat_ss_yr(i)%fpln/1000000.) / day_yr_r !J --> MJ
          gw_heat_ss_yr(i)%pond = (gw_heat_ss_yr(i)%pond/1000000.) / day_yr_r !J --> MJ
        enddo
      endif

      !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.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%gwsw = (gwsol_ss_sum(i)%solute(s)%gwsw/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%swgw = (gwsol_ss_sum(i)%solute(s)%swgw/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%soil = (gwsol_ss_sum(i)%solute(s)%soil/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%satx = (gwsol_ss_sum(i)%solute(s)%satx/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%ppex = (gwsol_ss_sum(i)%solute(s)%ppex/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%tile = (gwsol_ss_sum(i)%solute(s)%tile/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%resv = (gwsol_ss_sum(i)%solute(s)%resv/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%wetl = (gwsol_ss_sum(i)%solute(s)%wetl/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%canl = (gwsol_ss_sum(i)%solute(s)%canl/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%fpln = (gwsol_ss_sum(i)%solute(s)%fpln/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%pond = (gwsol_ss_sum(i)%solute(s)%pond/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%advn = (gwsol_ss_sum(i)%solute(s)%advn/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%disp = (gwsol_ss_sum(i)%solute(s)%disp/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%rcti = (gwsol_ss_sum(i)%solute(s)%rcti/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%rcto = (gwsol_ss_sum(i)%solute(s)%rcto/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%minl = (gwsol_ss_sum(i)%solute(s)%minl/1000.) / day_yr_r !g --> kg
            gwsol_ss_sum(i)%solute(s)%sorb = (gwsol_ss_sum(i)%solute(s)%sorb/1000.) / day_yr_r !g --> kg
          enddo
        enddo
      endif

      !pumping (irrigation) (for HRUs) -- yearly long-format output ----------------------
      if(gwflag_pump == 1) then
        do i=1,sp_ob%hru
          if(hru_pump_yr(i) > 0.) then
            iob = sp_ob1%hru + i - 1
            write(out_hru_pump_yr,8101) time%day,time%mo,time%day_mo, &
              time%yrc,i,ob(iob)%gis_id,ob(iob)%name,hru_pump_yr(i)
          endif
        enddo
      endif
      !accumulate into AA before zeroing
      do i=1,sp_ob%hru
        hru_pump_aa(i) = hru_pump_aa(i) + hru_pump_yr(i)
      enddo
      hru_pump_yr = 0.

      !(individual grid-file writes per flux component removed -- replaced by long-format above)

      !zero out flux sums to prepare for the next year ------------------------------------------
      !flow
      do i=1,ncell
        gw_hyd_ss_yr(i)%rech = 0.
        gw_hyd_ss_yr(i)%gwet = 0.
        gw_hyd_ss_yr(i)%gwsw = 0.
        gw_hyd_ss_yr(i)%swgw = 0.
        gw_hyd_ss_yr(i)%satx = 0.
        gw_hyd_ss_yr(i)%soil = 0.
        gw_hyd_ss_yr(i)%latl = 0.
        gw_hyd_ss_yr(i)%bndr = 0.
        gw_hyd_ss_yr(i)%ppag = 0.
        gw_hyd_ss_yr(i)%ppdf = 0.
        gw_hyd_ss_yr(i)%ppex = 0.
        gw_hyd_ss_yr(i)%tile = 0.
        gw_hyd_ss_yr(i)%resv = 0.
        gw_hyd_ss_yr(i)%wetl = 0.
        gw_hyd_ss_yr(i)%canl = 0.
        gw_hyd_ss_yr(i)%fpln = 0.
        gw_hyd_ss_yr(i)%pond = 0.
        gw_hyd_ss_yr(i)%phyt = 0.
      enddo
      !heat flux
      if(gw_heat_flag == 1) then
        do i=1,ncell
          gw_heat_ss_yr(i)%rech = 0.
          gw_heat_ss_yr(i)%gwet = 0.
          gw_heat_ss_yr(i)%gwsw = 0.
          gw_heat_ss_yr(i)%swgw = 0.
          gw_heat_ss_yr(i)%satx = 0.
          gw_heat_ss_yr(i)%soil = 0.
          gw_heat_ss_yr(i)%latl = 0.
          gw_heat_ss_yr(i)%disp = 0.
          gw_heat_ss_yr(i)%bndr = 0.
          gw_heat_ss_yr(i)%ppag = 0.
          gw_heat_ss_yr(i)%ppex = 0.
          gw_heat_ss_yr(i)%tile = 0.
          gw_heat_ss_yr(i)%resv = 0.
          gw_heat_ss_yr(i)%wetl = 0.
          gw_heat_ss_yr(i)%canl = 0.
          gw_heat_ss_yr(i)%fpln = 0.
          gw_heat_ss_yr(i)%pond = 0.
        enddo
      endif
      !solute
      if(gw_solute_flag == 1) then
        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)%canl = 0.
            gwsol_ss_sum(i)%solute(s)%fpln = 0.
            gwsol_ss_sum(i)%solute(s)%pond = 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
      endif

      !yearly water balance (basin) -------------------------------------------------------------
      if(gwflag_yr == 1) then
        write(out_gwbal_yr,8100) time%day,time%mo,time%day_mo, &
          time%yrc,"       1","       1",bsn%name, &
          gw_hyd_grid_yr%chng,gw_hyd_grid_yr%rech,gw_hyd_grid_yr%gwet, &
          gw_hyd_grid_yr%gwsw,gw_hyd_grid_yr%swgw, &
          gw_hyd_grid_yr%satx,gw_hyd_grid_yr%soil,gw_hyd_grid_yr%latl, &
          gw_hyd_grid_yr%bndr,gw_hyd_grid_yr%ppag, &
          gw_hyd_grid_yr%ppex,gw_hyd_grid_yr%tile,gw_hyd_grid_yr%resv, &
          gw_hyd_grid_yr%wetl,gw_hyd_grid_yr%canl, &
          gw_hyd_grid_yr%fpln,gw_hyd_grid_yr%pond,gw_hyd_grid_yr%phyt, &
          gw_hyd_grid_yr%ppdf
      endif

      !zero out annual basin arrays
      gw_hyd_grid_yr%chng = 0.
      gw_hyd_grid_yr%rech = 0.
      gw_hyd_grid_yr%gwet = 0.
      gw_hyd_grid_yr%gwsw = 0.
      gw_hyd_grid_yr%swgw = 0.
      gw_hyd_grid_yr%satx = 0.
      gw_hyd_grid_yr%soil = 0.
      gw_hyd_grid_yr%latl = 0.
      gw_hyd_grid_yr%bndr = 0.
      gw_hyd_grid_yr%ppag = 0.
      gw_hyd_grid_yr%ppdf = 0.
      gw_hyd_grid_yr%ppex = 0.
      gw_hyd_grid_yr%tile = 0.
      gw_hyd_grid_yr%resv = 0.
      gw_hyd_grid_yr%wetl = 0.
      gw_hyd_grid_yr%canl = 0.
      gw_hyd_grid_yr%fpln = 0.
      gw_hyd_grid_yr%pond = 0.
      gw_hyd_grid_yr%phyt = 0.

      !heat flux values (basin) -----------------------------------------------------------------
      if(gw_heat_flag == 1) then
        if(gwflag_yr == 1) then
          write(out_heatbal_yr,8100) time%day,time%mo,time%day_mo, &
            time%yrc,"       1","       1",bsn%name, &
            gw_heat_grid_yr%chng, &
            gw_heat_grid_yr%rech,gw_heat_grid_yr%gwet, &
            gw_heat_grid_yr%gwsw, &
            gw_heat_grid_yr%swgw,gw_heat_grid_yr%satx, &
            gw_heat_grid_yr%soil, &
            gw_heat_grid_yr%latl,gw_heat_grid_yr%disp, &
            gw_heat_grid_yr%bndr, &
            gw_heat_grid_yr%ppag,gw_heat_grid_yr%ppex, &
            gw_heat_grid_yr%tile, &
            gw_heat_grid_yr%resv,gw_heat_grid_yr%wetl, &
            gw_heat_grid_yr%canl, &
            gw_heat_grid_yr%fpln,gw_heat_grid_yr%pond
        endif
        gw_heat_grid_yr%chng = 0.
        gw_heat_grid_yr%rech = 0.
        gw_heat_grid_yr%gwet = 0.
        gw_heat_grid_yr%gwsw = 0.
        gw_heat_grid_yr%swgw = 0.
        gw_heat_grid_yr%satx = 0.
        gw_heat_grid_yr%soil = 0.
        gw_heat_grid_yr%latl = 0.
        gw_heat_grid_yr%disp = 0.
        gw_heat_grid_yr%bndr = 0.
        gw_heat_grid_yr%ppag = 0.
        gw_heat_grid_yr%ppex = 0.
        gw_heat_grid_yr%tile = 0.
        gw_heat_grid_yr%resv = 0.
        gw_heat_grid_yr%wetl = 0.
        gw_heat_grid_yr%canl = 0.
        gw_heat_grid_yr%fpln = 0.
        gw_heat_grid_yr%pond = 0.
      endif

      !solute mass values (basin) ---------------------------------------------------------------
      if(gw_solute_flag == 1) then
        do s=1,gw_nsolute
          if(gwflag_yr == 1) then
            write(out_solbal_yr+s,8100) time%day,time%mo,time%day_mo, &
              time%yrc,"       1","       1",bsn%name, &
              sol_grid_chng_yr(s),sol_grid_rech_yr(s), &
              sol_grid_gwsw_yr(s),sol_grid_swgw_yr(s), &
              sol_grid_satx_yr(s), &
              sol_grid_soil_yr(s),sol_grid_advn_yr(s), &
              sol_grid_disp_yr(s), &
              sol_grid_rcti_yr(s),sol_grid_rcto_yr(s), &
              sol_grid_minl_yr(s),sol_grid_sorb_yr(s), &
              sol_grid_ppag_yr(s),sol_grid_ppex_yr(s), &
              sol_grid_tile_yr(s),sol_grid_resv_yr(s), &
              sol_grid_wetl_yr(s), &
              sol_grid_canl_yr(s),sol_grid_fpln_yr(s), &
              sol_grid_pond_yr(s)
          endif
          sol_grid_chng_yr(s) = 0.
          sol_grid_rech_yr(s) = 0.
          sol_grid_gwsw_yr(s) = 0.
          sol_grid_swgw_yr(s) = 0.
          sol_grid_satx_yr(s) = 0.
          sol_grid_soil_yr(s) = 0.
          sol_grid_advn_yr(s) = 0.
          sol_grid_disp_yr(s) = 0.
          sol_grid_rcti_yr(s) = 0.
          sol_grid_rcto_yr(s) = 0.
          sol_grid_minl_yr(s) = 0.
          sol_grid_sorb_yr(s) = 0.
          sol_grid_ppag_yr(s) = 0.
          sol_grid_ppex_yr(s) = 0.
          sol_grid_tile_yr(s) = 0.
          sol_grid_resv_yr(s) = 0.
          sol_grid_wetl_yr(s) = 0.
          sol_grid_canl_yr(s) = 0.
          sol_grid_fpln_yr(s) = 0.
          sol_grid_pond_yr(s) = 0.
        enddo !go to next solute
      endif

      !format statements (subroutine-local)
140   format(i8,i6,i6,i8,i8,i10,4x,a12,2f13.3,16e13.4)
8100  format(4i6,2a,2x,a16,50e13.4)
8101  format(4i6,2i8,a18,e13.4)
8102  format(4i6,2i8,a18,5e13.4)

      return
      end subroutine gwflow_output_yr

      !NOTE: old individual grid-file writes per flux component removed.
      !      (formerly ~400 lines of per-variable grid write blocks)
      !      All cell-level data now written in long-format above.
      !

      subroutine gwflow_output_aa

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine writes average annual gwflow output in SWAT+ long format:
!!    one row per active cell with average annual head, wtdepth, and average
!!    daily flow rates; basin-level water/heat/solute balance; HRU pumping;
!!    transit time grids.

      use gwflow_module
      use hydrograph_module
      use sd_channel_module
      use time_module
      use basin_module, only : bsn

      implicit none

      integer :: i, j, k, s, iob
      integer :: num_months
      integer :: gis_id
      real :: wtdepth
      real :: nbyr_r
      real :: obs_temp, obs_no3, obs_p
      character(len=16) :: obs_name
      real, allocatable :: temp_array(:)
      external :: gwflow_write_cell_array

      !--- only execute on the last day of the simulation ---
      if(time%yrc /= time%yrc_end .or. time%day /= time%day_end) return

      nbyr_r = real(time%nbyr)
      allocate(temp_array(ncell))

      !pumping for HRUs -- average annual long-format output
      if(gwflag_pump == 1) then
        do i=1,sp_ob%hru
          if(hru_pump_aa(i) > 0.) then
            iob = sp_ob1%hru + i - 1
            write(out_hru_pump_aa,8101) time%day,time%mo,time%day_mo, &
              time%yrc,i,ob(iob)%gis_id,ob(iob)%name,hru_pump_aa(i)/nbyr_r
          endif
        enddo
      endif

      !--- obs well average annual output ---
      if(gwflag_obs == 1 .and. gw_num_obs_wells > 0) then
        do k=1,gw_num_obs_wells
          i = gw_obs_cells(k)
          obs_temp = -99.; obs_no3 = -99.; obs_p = -99.
          if(gw_heat_flag == 1) obs_temp = gw_obs_temp_aa(k) / nbyr_r
          if(gw_solute_flag == 1) then
            obs_no3 = gw_obs_sol_aa(k,1) / nbyr_r
            obs_p = gw_obs_sol_aa(k,2) / nbyr_r
          endif
          write(obs_name,'(a4,i4.4)') 'obs_',k
          write(out_gwobs_aa,8102) time%day,time%mo,time%day_mo, &
            time%yrc,k,gw_obs_cells(k),obs_name, &
            gw_head_sum_aa(i)/nbyr_r, &
            gw_state(i)%elev - gw_head_sum_aa(i)/nbyr_r, &
            obs_temp, obs_no3, obs_p
        enddo
      endif

      !--- cell-level average annual long-format output ---
      !      Per-cell yearly averages were accumulated into gw_hyd_ss_aa in _yr;
      !      divide by nbyr to get average annual values.
      if(gwflag_aa == 1) then
        do i=1,ncell
          if(gw_state(i)%stat == 1) then
            gis_id = cell_gis_id(i)
            wtdepth = gw_state(i)%elev - (gw_head_sum_aa(i) / nbyr_r)
            write(out_gwcell_aa,140) time%day, time%mo, &
              time%day_mo, time%yrc, i, gis_id, &
              cell_name(i), &
              gw_head_sum_aa(i) / nbyr_r, wtdepth, &
              gw_hyd_ss_aa(i)%rech / nbyr_r, &
              gw_hyd_ss_aa(i)%gwet / nbyr_r, &
              gw_hyd_ss_aa(i)%gwsw / nbyr_r, &
              gw_hyd_ss_aa(i)%swgw / nbyr_r, &
              gw_hyd_ss_aa(i)%satx / nbyr_r, &
              gw_hyd_ss_aa(i)%soil / nbyr_r, &
              gw_hyd_ss_aa(i)%latl / nbyr_r, &
              gw_hyd_ss_aa(i)%ppag / nbyr_r, &
              gw_hyd_ss_aa(i)%ppex / nbyr_r, &
              gw_hyd_ss_aa(i)%tile / nbyr_r, &
              gw_hyd_ss_aa(i)%resv / nbyr_r, &
              gw_hyd_ss_aa(i)%wetl / nbyr_r, &
              gw_hyd_ss_aa(i)%fpln / nbyr_r, &
              gw_hyd_ss_aa(i)%canl / nbyr_r, &
              gw_hyd_ss_aa(i)%pond / nbyr_r, &
              gw_hyd_ss_aa(i)%phyt / nbyr_r
          endif
        enddo
      endif

      !average annual water balance (basin) -----------------------------------------------------
      gw_hyd_grid_aa%chng = gw_hyd_grid_aa%chng + (vaft_grid-vbef_grid)
      gw_hyd_grid_aa%rech = gw_hyd_grid_aa%rech / nbyr_r
      gw_hyd_grid_aa%gwet = gw_hyd_grid_aa%gwet / nbyr_r
      gw_hyd_grid_aa%gwsw = gw_hyd_grid_aa%gwsw / nbyr_r
      gw_hyd_grid_aa%swgw = gw_hyd_grid_aa%swgw / nbyr_r
      gw_hyd_grid_aa%satx = gw_hyd_grid_aa%satx / nbyr_r
      gw_hyd_grid_aa%soil = gw_hyd_grid_aa%soil / nbyr_r
      gw_hyd_grid_aa%latl = gw_hyd_grid_aa%latl / nbyr_r
      gw_hyd_grid_aa%bndr = gw_hyd_grid_aa%bndr / nbyr_r
      gw_hyd_grid_aa%ppag = gw_hyd_grid_aa%ppag / nbyr_r
      gw_hyd_grid_aa%ppdf = gw_hyd_grid_aa%ppdf / nbyr_r
      gw_hyd_grid_aa%ppex = gw_hyd_grid_aa%ppex / nbyr_r
      gw_hyd_grid_aa%tile = gw_hyd_grid_aa%tile / nbyr_r
      gw_hyd_grid_aa%resv = gw_hyd_grid_aa%resv / nbyr_r
      gw_hyd_grid_aa%wetl = gw_hyd_grid_aa%wetl / nbyr_r
      gw_hyd_grid_aa%canl = gw_hyd_grid_aa%canl / nbyr_r
      gw_hyd_grid_aa%fpln = gw_hyd_grid_aa%fpln / nbyr_r
      gw_hyd_grid_aa%pond = gw_hyd_grid_aa%pond / nbyr_r
      gw_hyd_grid_aa%phyt = gw_hyd_grid_aa%phyt / nbyr_r
      if(gwflag_aa == 1) then
        write(out_gwbal_aa,8100) time%day,time%mo,time%day_mo, &
          time%yrc,"       1","       1",bsn%name, &
          gw_hyd_grid_aa%chng,gw_hyd_grid_aa%rech,gw_hyd_grid_aa%gwet, &
          gw_hyd_grid_aa%gwsw,gw_hyd_grid_aa%swgw, &
          gw_hyd_grid_aa%satx,gw_hyd_grid_aa%soil,gw_hyd_grid_aa%latl, &
          gw_hyd_grid_aa%bndr,gw_hyd_grid_aa%ppag, &
          gw_hyd_grid_aa%ppex,gw_hyd_grid_aa%tile,gw_hyd_grid_aa%resv, &
          gw_hyd_grid_aa%wetl,gw_hyd_grid_aa%canl, &
          gw_hyd_grid_aa%fpln,gw_hyd_grid_aa%pond,gw_hyd_grid_aa%phyt, &
          gw_hyd_grid_aa%ppdf
      endif

      !average annual heat fluxes (basin) -------------------------------------------------------
      gw_heat_grid_aa%chng = gw_heat_grid_aa%chng + (heat_haft_grid-heat_hbef_grid)
      gw_heat_grid_aa%rech = gw_heat_grid_aa%rech / nbyr_r
      gw_heat_grid_aa%gwet = gw_heat_grid_aa%gwet / nbyr_r
      gw_heat_grid_aa%gwsw = gw_heat_grid_aa%gwsw / nbyr_r
      gw_heat_grid_aa%swgw = gw_heat_grid_aa%swgw / nbyr_r
      gw_heat_grid_aa%satx = gw_heat_grid_aa%satx / nbyr_r
      gw_heat_grid_aa%soil = gw_heat_grid_aa%soil / nbyr_r
      gw_heat_grid_aa%latl = gw_heat_grid_aa%latl / nbyr_r
      gw_heat_grid_aa%disp = gw_heat_grid_aa%disp / nbyr_r
      gw_heat_grid_aa%bndr = gw_heat_grid_aa%bndr / nbyr_r
      gw_heat_grid_aa%ppag = gw_heat_grid_aa%ppag / nbyr_r
      gw_heat_grid_aa%ppex = gw_heat_grid_aa%ppex / nbyr_r
      gw_heat_grid_aa%tile = gw_heat_grid_aa%tile / nbyr_r
      gw_heat_grid_aa%resv = gw_heat_grid_aa%resv / nbyr_r
      gw_heat_grid_aa%wetl = gw_heat_grid_aa%wetl / nbyr_r
      gw_heat_grid_aa%canl = gw_heat_grid_aa%canl / nbyr_r
      gw_heat_grid_aa%fpln = gw_heat_grid_aa%fpln / nbyr_r
      gw_heat_grid_aa%pond = gw_heat_grid_aa%pond / nbyr_r
      if(gwflag_aa == 1) then
        write(out_heatbal_aa,8100) time%day,time%mo,time%day_mo, &
          time%yrc,"       1","       1",bsn%name, &
          gw_heat_grid_aa%chng, &
          gw_heat_grid_aa%rech,gw_heat_grid_aa%gwet,gw_heat_grid_aa%gwsw, &
          gw_heat_grid_aa%swgw,gw_heat_grid_aa%satx,gw_heat_grid_aa%soil, &
          gw_heat_grid_aa%latl,gw_heat_grid_aa%disp,gw_heat_grid_aa%bndr, &
          gw_heat_grid_aa%ppag,gw_heat_grid_aa%ppex,gw_heat_grid_aa%tile, &
          gw_heat_grid_aa%resv,gw_heat_grid_aa%wetl,gw_heat_grid_aa%canl, &
          gw_heat_grid_aa%fpln,gw_heat_grid_aa%pond
      endif

      !average annual solute values (basin) -----------------------------------------------------
      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) / nbyr_r
          sol_grid_gwsw_tt(s) = sol_grid_gwsw_tt(s) / nbyr_r
          sol_grid_swgw_tt(s) = sol_grid_swgw_tt(s) / nbyr_r
          sol_grid_satx_tt(s) = sol_grid_satx_tt(s) / nbyr_r
          sol_grid_advn_tt(s) = sol_grid_advn_tt(s) / nbyr_r
          sol_grid_disp_tt(s) = sol_grid_disp_tt(s) / nbyr_r
          sol_grid_rcti_tt(s) = sol_grid_rcti_tt(s) / nbyr_r
          sol_grid_rcto_tt(s) = sol_grid_rcto_tt(s) / nbyr_r
          sol_grid_minl_tt(s) = sol_grid_minl_tt(s) / nbyr_r
          sol_grid_sorb_tt(s) = sol_grid_sorb_tt(s) / nbyr_r
          sol_grid_ppag_tt(s) = sol_grid_ppag_tt(s) / nbyr_r
          sol_grid_ppex_tt(s) = sol_grid_ppex_tt(s) / nbyr_r
          sol_grid_tile_tt(s) = sol_grid_tile_tt(s) / nbyr_r
          sol_grid_soil_tt(s) = sol_grid_soil_tt(s) / nbyr_r
          sol_grid_resv_tt(s) = sol_grid_resv_tt(s) / nbyr_r
          sol_grid_wetl_tt(s) = sol_grid_wetl_tt(s) / nbyr_r
          sol_grid_canl_tt(s) = sol_grid_canl_tt(s) / nbyr_r
          sol_grid_fpln_tt(s) = sol_grid_fpln_tt(s) / nbyr_r
          sol_grid_pond_tt(s) = sol_grid_pond_tt(s) / nbyr_r
          if(gwflag_aa == 1) then
            write(out_solbal_aa+s,8100) time%day,time%mo,time%day_mo, &
              time%yrc,"       1","       1",bsn%name, &
              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), &
              sol_grid_pond_tt(s)
          endif
        enddo !next solute
      endif

      !if soft calibration, prepare for next simulation
      sim_month = 1

      !write out groundwater transit time to channels and tiles
      if(gw_ttime == 1) then
        do k=1,ncell
          temp_array(k) = gw_cell_chan_time(k)
        enddo
        call gwflow_write_cell_array(out_gw_transit_chan, temp_array, ncell, 2)
        if(gw_tile_flag == 1) then
          do k=1,ncell
            temp_array(k) = gw_cell_tile_time(k)
          enddo
          call gwflow_write_cell_array(out_gw_transit_tile, temp_array, ncell, 2)
        endif
      endif

      if(allocated(temp_array)) deallocate(temp_array)

      !format statements (subroutine-local)
105   format(i8,1000(e13.4))
140   format(i8,i6,i6,i8,i8,i10,4x,a12,2f13.3,16e13.4)
8100  format(4i6,2a,2x,a16,50e13.4)
8101  format(4i6,2i8,a18,e13.4)
8102  format(4i6,2i8,a18,5e13.4)

      return
      end subroutine gwflow_output_aa


!     ==========================================================================
!     Helper subroutines for output redesign
!     ==========================================================================


      subroutine gwflow_write_celldef
!!    Writes gwflow_cell_definition.txt once during initialization.
!!    Maps cell index to spatial location for all output files.
      use gwflow_module
      implicit none
      integer :: i

      open(out_gw_celldef, file='gwflow_cell_definition.txt')
      write(out_gw_celldef,'(a)') 'cell_id  row  col        x_coord        y_coord  zone  status          area'
      do i=1,ncell
        if(gw_state(i)%stat > 0) then
          if(grid_type == "structured") then
            write(out_gw_celldef,'(i8,2i6,2f15.1,2i6,e15.4)') &
              i, cell_row(i), cell_col(i), gw_state(i)%xcrd, gw_state(i)%ycrd, &
              gw_state(i)%zone, gw_state(i)%stat, gw_state(i)%area
          else
            write(out_gw_celldef,'(i8,2i6,2f15.1,2i6,e15.4)') &
              i, 0, 0, gw_state(i)%xcrd, gw_state(i)%ycrd, &
              gw_state(i)%zone, gw_state(i)%stat, gw_state(i)%area
          endif
        endif
      enddo
      close(out_gw_celldef)

      return
      end subroutine gwflow_write_celldef


      subroutine gwflow_write_cell_array(iunit, values, ncell_in, fmt_code)
!!    Writes active cell values as a single row.
!!    fmt_code: 1=f12.3 (heads), 2=e12.3 (fluxes), 3=e12.6 (high precision)
      use gwflow_module, only : gw_state
      implicit none
      integer, intent(in) :: iunit
      real, dimension(ncell_in), intent(in) :: values
      integer, intent(in) :: ncell_in
      integer, intent(in) :: fmt_code
      integer :: i

      select case(fmt_code)
      case(1) !heads
        write(iunit,101) (values(i), i=1,ncell_in)
      case(2) !fluxes
        write(iunit,102) (values(i), i=1,ncell_in)
      case(3) !high precision
        write(iunit,103) (values(i), i=1,ncell_in)
      end select

101   format(99999(f12.3))
102   format(99999(e12.3))
103   format(99999(e12.6))

      return
      end subroutine gwflow_write_cell_array


