New river inflow scheme implemented. To be tested

setup/Mojai2016_setup.dat

@@ -144,7 +144,7 @@ ifrad    1
 ifbubble 1
 sedim    0
 salsoil  0
-dyn_pgrad 0
+dyn_pgrad 4
 zero_model 0
 thermokarst_meth_prod 0.
 soil_meth_prod 0.
@@ -208,12 +208,14 @@ T_profile 20
 #-----------------------------------------------------------------------------------------
 #
 #
-tribheat 2
+tribheat 3
 N_tribin 1
 1
 N_triblev 24
-fileinflow  'Mojai2016_inflows.dat'
-fileoutflow 'Mojai2016_outflow.dat'
+#fileinflow  'Mojai2016_inflows.dat'
+#fileoutflow 'Mojai2016_outflow.dat'
+fileinflow  'Mojai2016_inflows_tribheat=3.dat'
+fileoutflow 'Mojai2016_outflow_tribheat=3.dat'
 iefflloc 2
 dttribupdate 1.
 #

source/model/trib.f90

@@ -115,9 +115,11 @@ SUBROUTINE TRIBTEMP(time,dt,h1,z_full,area_int,gsp,Tw,Sal,u,v,qwater,gas, &
  ! If tribheat == 1, tributaries' and effluent's data are read once, and kept constant throughout model integration time
  if (tribheat%par == 1 .and. firstcall) call TRIBUPDATE(spinup_done,gsp)
 
+ !print*, 'time', time, dt
  if (tribheat%par == 2 .or. tribheat%par == 3) then 
    work = mod(time-dt,dttribupdate*day_sec) ! time is for next time step, thus lateral inflow
                                             ! values are updated to the current time step
+   !print*, 'work',work,firstcall
    if ((work >= 0.d0 .and. work < dt) .or. firstcall .or. spinup_done) then
      ! Updating tributary data
      call TRIBUPDATE(spinup_done,gsp)
@@ -263,6 +265,7 @@ SUBROUTINE TRIBTEMP(time,dt,h1,z_full,area_int,gsp,Tw,Sal,u,v,qwater,gas, &
      outflow = outflow + U_tribout_(k,i)*width_tribout_(k,i)
    enddo
    outflow = outflow/area_int(i)
+   !print*, 'outflow', outflow, U_tribin_(1,i),T_tribin_(1,i),width_tribin_(1,i), area_int(i), invdt
 
    if (outflpar%par == 1) then
      meanTin = 0.
@@ -295,18 +298,40 @@ SUBROUTINE TRIBTEMP(time,dt,h1,z_full,area_int,gsp,Tw,Sal,u,v,qwater,gas, &
    endif
 
    !Implicit scheme for the outflow term
+   !Check: if for any of tributaries the data on a variable is missing,
+   !this variable is not updated with inlet/outlet sources
    work = invdt + outflow*(1. + outflpar%par)
-   Tw(i) = (invdt*Tw(i) + heatinflow + outflpar%par*outflow*meanTin)/work
-   Sal(i) = (invdt*Sal(i) + salinflow + outflpar%par*outflow*meansalin)/work
-   qwater(i) = (invdt*qwater(i) + methinflow + outflpar%par*outflow*meanmethin)/work
+
+   !print*, 'Tw1', Tw
+   !print*, heatinflow, outflow, meanTin
+   if (.not. any(T_tribin_(:,i) == missing_value)) &
+   & Tw(i) = (invdt*Tw(i) + heatinflow + outflpar%par*outflow*meanTin)/work
+   !print*, 'Tw2', Tw
+   !print*, heatinflow, outflow, meanTin
+   !read*
+
+   if (.not. any(Sal_tribin_(:,i) == missing_value)) &
+   & Sal(i) = (invdt*Sal(i) + salinflow + outflpar%par*outflow*meansalin)/work
+
+   if (.not. any(meth_tribin_(:,i) == missing_value)) &
+   & qwater(i) = (invdt*qwater(i) + methinflow + outflpar%par*outflow*meanmethin)/work
+
    if (carbon_model%par == 2) then
-     gas%DOC(i) = (invdt*gas%DOC(i) + DOCinflow + outflpar%par*outflow*meanDOCin)/work
-     gas%POCL(i) = (invdt*gas%POCL(i) + POCinflow + outflpar%par*outflow*meanPOCin)/work
+     if (.not. any(DOC_tribin_(:,i) == missing_value)) &
+     & gas%DOC(i) = (invdt*gas%DOC(i) + DOCinflow + outflpar%par*outflow*meanDOCin)/work
+     if (.not. any(POC_tribin_(:,i) == missing_value)) &
+     & gas%POCL(i) = (invdt*gas%POCL(i) + POCinflow + outflpar%par*outflow*meanPOCin)/work
    endif
-   gas%carbdi(i,1) = (invdt*gas%carbdi(i,1) + DICinflow + outflpar%par*outflow*meanDICin)/work
-   Sal(i) = (invdt*Sal(i) + salinflow + outflpar%par*outflow*meansalin)/work
-   u(i) = (invdt*u(i) + uinflow + outflpar%par*outflow*meanuin)/work
-   v(i) = (invdt*v(i) + vinflow + outflpar%par*outflow*meanvin)/work
+
+   if (.not. any(DIC_tribin_(:,i) == missing_value)) &
+   & gas%carbdi(i,1) = (invdt*gas%carbdi(i,1) + DICinflow + outflpar%par*outflow*meanDICin)/work
+
+   !Sal(i) = (invdt*Sal(i) + salinflow + outflpar%par*outflow*meansalin)/work
+   if (.not. any(Ux_tribin_(:,i) == missing_value)) &
+   & u(i) = (invdt*u(i) + uinflow + outflpar%par*outflow*meanuin)/work
+
+   if (.not. any(Uy_tribin_(:,i) == missing_value)) &
+   & v(i) = (invdt*v(i) + vinflow + outflpar%par*outflow*meanvin)/work
    ! qwater(i) = (invdt*qwater(i) + 0.d0)/(invdt + outflow) ! Implying zero concentration of methane in the inflow
 
  enddo
@@ -358,6 +383,8 @@ use DRIVING_PARAMS, only : &
 
 use INOUT, only : CHECK_UNIT
 
+use ARRAYS_GRID, only : gridspacing_type
+
 implicit none
 
 ! Input variables
@@ -369,7 +396,7 @@ type(gridspacing_type), intent(in) :: gsp
 ! Local variables
 real(kind=ireals), allocatable :: work(:,:)
 real(kind=ireals) :: disch, depth, width, temp, sal, &
-& dirprojX, dirprojY, DOC, POC, DIC, CH_4
+& dirprojX, dirprojY, DOC, POC, DIC, CH_4, xx
 integer(kind=iintegers), save :: nunit1 = lake_misc_unit_min, &
 & nunit2 = lake_misc_unit_min + 1
 integer(kind=iintegers) :: ncol, i, idep
@@ -377,7 +404,7 @@ integer(kind=iintegers) :: ncol, i, idep
 character(len=20) :: workchar
 
 logical, save :: firstcall = .true.
-
+!print*, 'entered', tribheat%par
 if (firstcall) then
   call CHECK_UNIT(lake_misc_unit_min,lake_misc_unit_max,nunit1)
   call CHECK_UNIT(lake_misc_unit_min,lake_misc_unit_max,nunit2)
@@ -422,12 +449,15 @@ else if (tribheat%par == 3) then
   ! timestring, discharge, depth, width, temperature, salinity, riverflow direction
   ! projection on X axis (-1,+1), riverflow direction projection on Y axis (-1,+1),
   ! DOC, POC, DIC, CH_4
+  !print*, 'entered2'
   do i = 1, N_tribin%par
 
-    read(nunit1,*) workchar, disch, depth, width, temp, sal, &
+    read(nunit1,*) xx, disch, depth, width, temp, sal, &
     & dirprojX, dirprojY, DOC, POC, DIC, CH_4
+    !print*, xx, disch, depth, width, temp, sal, &
+    !& dirprojX, dirprojY, DOC, POC, DIC, CH_4
     idep = 1
-    do while (gsp%z_half(i)/depth < 1.) 
+    do while (gsp%z_half(idep)/depth < 1.) 
       idep = idep + 1
     enddo
     idep = idep - 1
@@ -443,12 +473,18 @@ else if (tribheat%par == 3) then
     POC_tribin  (i,1:idep) = POC                           ; POC_tribin  (i,idep+1:M+1) = 0.
     DIC_tribin  (i,1:idep) = DIC                           ; DIC_tribin  (i,idep+1:M+1) = 0.
     meth_tribin (i,1:idep) = CH_4                          ; meth_tribin (i,idep+1:M+1) = 0.
+    !if (idep > 10) then
+    !print*, idep
+    !read*
+    !endif
   enddo
+  !print*, 'in', width_tribin, U_tribin, T_tribin
+  !read*
 
   ! In current version there is only ONE outflow 
   read(nunit2,*) workchar, disch, depth, width
   idep = 1
-  do while (gsp%z_half(i)/depth < 1.) 
+  do while (gsp%z_half(idep)/depth < 1.) 
     idep = idep + 1
   enddo
   idep = idep - 1