Changes for constituents, salt, gwflow, and heat

.gitignore

@@ -114,3 +114,4 @@ ford.md
 /data/Osu_1hru/basin_wb_yr.txt
 /data/Osu_1hru/aquifer_aa.txt
 [Tt]mp/
+/src/command.f90

CMakeLists.txt

@@ -12,7 +12,7 @@ project(swatplus VERSION 0.2 LANGUAGES Fortran)
 include(CTest)
 
 if(NOT TAG)
-    find_package(Git REQUIRED)
+    # find_package(Git REQUIRED)
     execute_process(
       COMMAND ${GIT_EXECUTABLE} describe --tags
       WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
@@ -58,9 +58,9 @@ elseif(WIN32)
         set(FFC             "ifo")
         link_libraries("-static")
     elseif(CMAKE_Fortran_COMPILER_ID STREQUAL IntelLLVM)
-        set(fdialect        "/fpp free /fpe0 /traceback")
+        set(fdialect        "/fpp /free /fpe0 /traceback")
         set(fdebug          "/warn:all")
-        set(frelease        "/O")
+        set(frelease        "/O2 /libs:static /threads")
         set(FFC             "ifx")
         # link_libraries("-static")
     elseif(CMAKE_Fortran_COMPILER_ID MATCHES GNU)

src/aqu_1d_control.f90

@@ -255,7 +255,10 @@ subroutine aqu_1d_control
       ob(icmd)%hd(1)%solp = aqu_d(iaq)%minp
 
       !! temperature of aquifer flow
-      ob(icmd)%hd(1)%temp = w_temp%gw
+      !ob(icmd)%hd(1)%temp = w_temp%gw  
+
+      !! separate gw flo !KPeters added for stream temperature routine
+      ob(icmd)%hdsep%flo_gwsw = ob(icmd)%hd(1)%flo
 
       !! compute fraction of flow to each channel in the aquifer
       !! if connected to aquifer - add flow

src/aqu_read_init_cs.f90

@@ -36,6 +36,7 @@ subroutine aqu_read_init_cs !rtb salt/cs
       real :: aqu_bd = 0.
       real :: aqu_mass = 0.
       real :: mass_sorbed = 0.
+      real :: init_conc = 0.
       eof = 0
       imax = 0
 
@@ -100,28 +101,25 @@ subroutine aqu_read_init_cs !rtb salt/cs
           enddo
 
           !! initial pathogens
-          do isp_ini = 1, imax
-            if (aqu_init_dat_c_cs(isp_ini)%name == aqudb(iaq)%aqu_ini) then
-              do ics = 1, db_mx%pathw_ini
-                if (aqu_init_dat_c_cs(isp_ini)%path == path_init_name(ics)) then
-                  !! initialize pathogens in aquifer water and benthic from input data
-                  do ipath = 1, cs_db%num_paths
-                    cs_aqu(iaq)%path(ipath) = path_soil_ini(ics)%soil(ipath)
-                  end do
-                  exit
-                end if
-              end do  
-            endif
-          enddo
+          do ics = 1, db_mx%pathw_ini
+            if (aqu_init_dat_c_cs(iaq)%path == path_init_name(ics)) then
+              !! initialize pathogens in aquifer water and benthic from input data
+              do ipath = 1, cs_db%num_paths
+                cs_aqu(iaq)%path(ipath) = path_soil_ini(ics)%soil(ipath)
+              end do
+              exit
+            end if
+          end do
 
           !! initial salts !rtb salt
           if(cs_db%num_salts > 0) then
           do ics = 1, db_mx%salt_gw_ini
             if (aqu_init_dat_c_cs(iaq)%salt == salt_aqu_ini(ics)%name) then
               !loop for salt ions
               do isalt = 1,cs_db%num_salts
-                cs_aqu(iaq)%saltc(isalt) = salt_aqu_ini(ics)%conc(isalt) !g/m3 (mg/L)
-                cs_aqu(iaq)%salt(isalt) = (salt_aqu_ini(ics)%conc(isalt)*gw_volume) / 1000. !g/m3 --> kg
+                init_conc = salt_aqu_ini(ics)%conc(isalt)
+                cs_aqu(iaq)%saltc(isalt) = init_conc !g/m3 (mg/L)
+                cs_aqu(iaq)%salt(isalt) = (init_conc*gw_volume) / 1000. !g/m3 --> kg
               enddo
               ! loop for salt mineral fractions
               do isalt = 1,5

src/ch_cs_module.f90

@@ -19,7 +19,11 @@ module ch_cs_module !rtb cs
         type (ch_cs_balance), dimension (:), allocatable :: cs         !cs hydrographs
       end type ch_cs_output
       type (ch_cs_balance) :: ch_csbz
-
+
+      real :: cs_conc_in(8) = 0.        !kg         !constituent concentration in channel water
+
+      real, allocatable :: div_conc_cs(:,:)  !g/m3  |solute concentration in diverted water
+
       type (ch_cs_output), dimension(:), allocatable, save :: chcs_d
       type (ch_cs_output), dimension(:), allocatable, save :: chcs_m
       type (ch_cs_output), dimension(:), allocatable, save :: chcs_y

src/ch_read.f90

@@ -21,6 +21,7 @@ subroutine ch_read
       integer :: ihyd = 0               !none          |counter
       integer :: ised = 0               !none          |counter
       integer :: inut = 0               !none          |counter
+      integer :: itemp = 0              !none          |counter
 
       eof = 0
       imax = 0
@@ -90,10 +91,18 @@ subroutine ch_read
            end if
          end do   
 
+         do itemp = 1, db_mx%w_temp 
+           if (w_temp(itemp)%name == ch_dat_c(ichi)%temp) then
+             ch_dat(ichi)%temp = itemp
+             exit
+           end if
+         end do 
+
        if (ch_dat(ichi)%init == 0) write (9001,*) ch_dat_c(ichi)%init, " not found (initial.cha)"
        if (ch_dat(ichi)%hyd == 0) write (9001,*) ch_dat_c(ichi)%hyd, " not found (hydrology.cha)"
        if (ch_dat(ichi)%sed == 0) write (9001,*) ch_dat_c(ichi)%sed, " not found (sediment.cha)"
        if (ch_dat(ichi)%nut == 0) write (9001,*) ch_dat_c(ichi)%nut, " not found (nutrients.cha)"      
+       if (ch_dat(ichi)%nut == 0) write (9001,*) ch_dat_c(ichi)%temp, " not found (temperature.cha)"      
 
        end do
 
@@ -104,4 +113,5 @@ subroutine ch_read
 
       return
 
-    end subroutine ch_read
+	end subroutine ch_read
+

src/ch_read_temp.f90

@@ -1,6 +1,7 @@
       subroutine ch_read_temp
 
       use basin_module
+      use time_module
       use input_file_module
       use maximum_data_module
       use channel_data_module
@@ -19,45 +20,48 @@ subroutine ch_read_temp
 
       inquire (file=in_cha%temp, exist=i_exist)
       if (.not. i_exist .or. in_cha%temp == "null") then
-        imax = 0    !allocate (w_temp(0:0))
+        allocate (w_temp(0:0))
       else   
       do
-       open (105,file=in_cha%temp)
-       read (105,*,iostat=eof) titldum
-       if (eof < 0) exit
-       read (105,*,iostat=eof) header
-       if (eof < 0) exit
-       read (105,*,iostat=eof) w_temp
-       if (w_temp%airlag_d < 1) w_temp%airlag_d = 6
-       !if (w_temp%hex_coef1 < 1) w_temp%hex_coef1 = 6
-       if (eof < 0) exit
-
-        !do while (eof == 0)
-        !  read (105,*,iostat=eof) titldum
-        !  if (eof < 0) exit
-        !  imax = imax + 1
-        !end do
+        open (105,file=in_cha%temp)
+        read (105,*,iostat=eof) titldum
+        if (eof < 0) exit
+        read (105,*,iostat=eof) header
+        !if (eof < 0) exit
+        !read (105,*,iostat=eof) wa_temp
+        !if (wa_temp%airlag_d < 1) wa_temp%airlag_d = 6
+        !if (wa_temp%hex_coef1 < 1) wa_temp%hex_coef1 = 6
 
-      db_mx%ch_temp = imax
+        if (eof < 0) exit     
+          do while (eof == 0)
+            read (105,*,iostat=eof) titldum
+            if (eof < 0) exit
+            imax = imax + 1
+          end do
+
+        db_mx%w_temp = imax
+
+        if (allocated(w_temp)) deallocate(w_temp)
 
-      !allocate (ch_temp(0:imax))
-      !rewind (105)
-      !read (105,*,iostat=eof) titldum
-      !if (eof < 0) exit
-      !read (105,*,iostat=eof) header
-      !if (eof < 0) exit
+        allocate (w_temp(0:imax))
+        rewind (105)
+        read (105,*,iostat=eof) titldum
+        if (eof < 0) exit
+        read (105,*,iostat=eof) header
+        if (eof < 0) exit
 
-       !do ich = 1, db_mx%ch_temp
-         !read (105,*,iostat=eof) titldum
-         !if (eof < 0) exit
-         !backspace (105)
-         !read (105,*,iostat=eof) w_temp
-         !if (eof < 0) exit
-       !end do
-       close (105)
-      exit
+        do ich = 1, db_mx%w_temp
+            read (105,*,iostat=eof) titldum
+            if (eof < 0) exit
+            backspace (105)
+            read (105,*,iostat=eof) w_temp
+            if (eof < 0) exit
+        end do
+
+        exit
       enddo
       endif
-
+      close (105)
       return    
-      end subroutine ch_read_temp
+	end subroutine ch_read_temp
+

src/ch_rtmusk.f90

@@ -23,6 +23,9 @@ subroutine ch_rtmusk
       use water_body_module
       use hru_module, only : hru
       use conditional_module
+      use ch_salt_module, only : salt_conc_in,chsalt_d
+      use ch_cs_module, only : cs_conc_in,chcs_d
+      use constituent_mass_module
 
       implicit none
 
@@ -35,6 +38,8 @@ subroutine ch_rtmusk
       integer :: ires = 0
       integer :: ihyd = 0
       integer :: irel = 0
+      integer :: isalt = 0
+      integer :: ics = 0
 
       real :: ch_stor_init = 0. !m3             |storage in channel at beginning of day
       real :: fp_stor_init = 0. !m3             |storage in flood plain above wetlands emergency spillway at beginning of day
@@ -72,6 +77,7 @@ subroutine ch_rtmusk
       real :: wet_evol = 0. 
       real :: bf_flow = 0.               !m3/s           |bankfull flow rate * adjustment factor
       real :: pk_rto = 0.                !ratio          |peak to mean flow rate ratio
+      real :: seep_mass                  !kg             |salt or constituent mass in channel seepage water
 
       jrch = isdch
       jhyd = sd_dat(jrch)%hyd
@@ -137,6 +143,19 @@ subroutine ch_rtmusk
           rto = Max(0., rto)
           rto = Min(1., rto)
           tot_stor(jrch) = tot_stor(jrch) + rto * ht1
+
+          !! rtb: subdaily input of salt and constituent mass
+          if(cs_db%num_salts > 0) then
+            do isalt=1,cs_db%num_salts
+              ch_water(jrch)%salt(isalt) = ch_water(jrch)%salt(isalt) + (rto * hcs1%salt(isalt))
+            enddo
+          endif
+          if(cs_db%num_cs > 0) then
+            do ics=1,cs_db%num_cs
+              ch_water(jrch)%cs(ics) = ch_water(jrch)%cs(ics) + (rto * hcs1%cs(ics))
+            enddo
+          endif
+
         end if    ! ht1%flo > 1.e-6
 
         !! interpolate rating curve using inflow rates
@@ -183,7 +202,20 @@ subroutine ch_rtmusk
           ob(icmd)%hyd_flo(1,irtstep) = ob(icmd)%hyd_flo(1,irtstep) + outflo
           !! subtract outflow from total storage
           tot_stor(jrch) = (1. - rto) * tot_stor(jrch)
-
+          !! rtb: add salt and constituent mass to outflow; then subtract from channel mass storage
+          if(cs_db%num_salts > 0) then
+            do isalt=1,cs_db%num_salts
+              hcs2%salt(isalt) = hcs2%salt(isalt) + (rto * ch_water(jrch)%salt(isalt))
+              ch_water(jrch)%salt(isalt) = ch_water(jrch)%salt(isalt) - hcs2%salt(isalt)
+            enddo
+          endif
+          if(cs_db%num_cs > 0) then
+            do ics=1,cs_db%num_cs
+              hcs2%cs(ics) = hcs2%cs(ics) + (rto * ch_water(jrch)%cs(ics))
+              ch_water(jrch)%cs(ics) = ch_water(jrch)%cs(ics) - hcs2%cs(ics)
+            enddo
+          endif
+
           !! set rating curve for next time step
           ch_rcurv(jrch)%in1 = ch_rcurv(jrch)%in2
           ch_rcurv(jrch)%out1 = ch_rcurv(jrch)%out2
@@ -216,8 +248,34 @@ subroutine ch_rtmusk
         rto = trans_loss / ch_stor(jrch)%flo
         ch_stor(jrch) = (1. - rto) * ch_stor(jrch)
       end if
-      ch_wat_d(ich)%seep = trans_loss
-
+      ch_wat_d(jrch)%seep = trans_loss
+      !rtb gwflow: if gwflow active, compute seepage in gwflow routine
+      if(bsn_cc%gwflow == 0) then
+	    ch_wat_d(jrch)%seep = trans_loss
+	  else
+        ch_wat_d(jrch)%seep = 0.
+      endif
+      !rtb: salt and constituent mass in seepage
+      do isalt=1,cs_db%num_salts
+        seep_mass = salt_conc_in(isalt) * ch_wat_d(jrch)%seep !g/m3 * m3 = g
+        seep_mass = seep_mass / 1000. !kg
+        if(seep_mass > ch_water(jrch)%salt(isalt)) then
+          seep_mass = ch_water(jrch)%salt(isalt)
+        endif
+        ch_water(jrch)%salt(isalt) = ch_water(jrch)%salt(isalt) - seep_mass !kg
+        chsalt_d(jrch)%salt(isalt)%seep = seep_mass !kg (channel salt output)
+      enddo
+      !constituent mass in seepage
+      do ics=1,cs_db%num_cs
+        seep_mass = cs_conc_in(ics) * ch_wat_d(jrch)%seep !g/m3 * m3 = g
+        seep_mass = seep_mass / 1000. !kg
+        if(seep_mass > ch_water(jrch)%cs(ics)) then
+          seep_mass = ch_water(jrch)%cs(ics)
+        endif
+        ch_water(jrch)%cs(ics) = ch_water(jrch)%cs(ics) - seep_mass !kg
+        chcs_d(jrch)%cs(ics)%seep = seep_mass !kg (channel consituent output)
+      enddo
+
       !! calculate evaporation losses
       if (ch_stor(jrch)%flo > 1.e-6) then
         !! calculate width of channel at water level - flood plain evap calculated in wetlands
@@ -233,7 +291,7 @@ subroutine ch_rtmusk
         rto = evap / ch_stor(jrch)%flo
         ch_stor(jrch)%flo = (1. - rto) * ch_stor(jrch)%flo
       end if
-      ch_wat_d(ich)%evap = evap
+      ch_wat_d(jrch)%evap = evap
 
       tot_stor(jrch) = ch_stor(jrch) + fp_stor(jrch)
 
@@ -261,5 +319,6 @@ subroutine ch_rtmusk
       ch_fp_wb(jrch)%wet_stor_init = wet_stor_init
       ch_fp_wb(jrch)%wet_stor = wet_stor(jrch)%flo
 
+
       return
       end subroutine ch_rtmusk

src/ch_salt_module.f90

@@ -24,7 +24,11 @@ module ch_salt_module
       type (ch_salt_output), dimension(:), allocatable, save :: chsalt_m
       type (ch_salt_output), dimension(:), allocatable, save :: chsalt_y
       type (ch_salt_output), dimension(:), allocatable, save :: chsalt_a
-
+
+      real :: salt_conc_in(8) = 0.         !kg         !salt concentration in channel water
+
+      real, allocatable :: div_conc_salt(:,:)  !g/m3  |salt ion concentration in diverted water
+
       type ch_salt_header
           character (len=6) :: day =        "  jday"
           character (len=6) :: mo =         "   mon"