From 24ad8a2ff3f36b11d48ae2194d1a8523e908606c Mon Sep 17 00:00:00 2001
From: Victor Stepanenko <vstepanenkomeister@gmail.com>
Date: Mon, 22 Aug 2022 09:40:42 +0300
Subject: [PATCH] An option tribheat=4 is added to take into account
 tributaries in large scale models (LSMs, ESMs)

---
 source/driver/driver.f90            |  30 ++--
 source/model/bathym_mod.f90         |   6 +-
 source/model/driving_params_mod.f90 |   9 +-
 source/model/lake.f90               |  23 +--
 source/model/trib.f90               | 249 +++++++++++++++++++++-------
 5 files changed, 229 insertions(+), 88 deletions(-)

diff --git a/source/driver/driver.f90 b/source/driver/driver.f90
index 1e9ba30..6144880 100644
--- a/source/driver/driver.f90
+++ b/source/driver/driver.f90
@@ -41,6 +41,7 @@
   integer(kind=iintegers), parameter :: n_select_call = 20
   integer(kind=iintegers), parameter :: ndatamax = 100
   integer(kind=iintegers), parameter :: n_netcdf_surf_files_max = 10
+  integer(kind=iintegers), parameter :: nvartribext = 4 !number of tributary variables
        
   real(kind=ireals), parameter :: day_sec = 24.*60.*60.
   real(kind=ireals), parameter :: hour_sec = 60.*60.
@@ -134,7 +135,7 @@
   real(kind=ireals) :: h_ML0
   real(kind=ireals) :: extwat, extice
   real(kind=ireals) :: kor
-  real(kind=ireals) :: trib_inflow, effl_outflow(1:ndatamax)
+  real(kind=ireals) :: trib_inflow(1:nvartribext), effl_outflow(1:ndatamax)
   real(kind=ireals) :: widthchan, lengthchan
   real(kind=ireals) :: hbot(1:2), hbotchan(1:2) ! assuming two interacting lakes
   real(kind=ireals) :: Sals0, Salb0
@@ -272,7 +273,7 @@
  & (1,hour0, tinteg, dt, height_T_q, height_u, interval, &
  &  h10, select_h10, l10, ls10, hs10, Ts0, Tb0, Tbb0, h_ML0, &
  &  extwat, select_extwat, extice, &
- &  kor, trib_inflow, effl_outflow, &
+ &  kor, trib_inflow(1), effl_outflow, &
  &  widthchan, lengthchan, hbot, hbotchan, &
  &  Sals0, Salb0, fetch, phi, lam, us0, vs0, &
  &  Tm, alphax, alphay, a_veg, c_veg, h_veg, area_lake, cellipt, depth_area, &
@@ -295,7 +296,7 @@
     ny = 1
   endif
 
-  vartrib_inflow = (trib_inflow == -999.) ! Tributary inflow varying with time
+  vartrib_inflow = (trib_inflow(1) == -999.) ! Tributary inflow varying with time
  
   allocate (ta(1:npoints) )
   allocate (qa(1:npoints) )
@@ -337,7 +338,7 @@
   allocate (extwat_2d(1:npoints,1:ny)) ; extwat_2d(:,1) = extwat
   allocate (extice_2d(1:npoints,1:ny)) ; extice_2d(:,1) = extice
   allocate (kor_2d(1:npoints,1:ny))    ; kor_2d(:,1) = kor
-  allocate (trib_inflow_2d(1:npoints,1:ny)) ; trib_inflow_2d(:,1) = trib_inflow
+  allocate (trib_inflow_2d(1:npoints,1:ny)) ; trib_inflow_2d(:,1) = trib_inflow(1)
   allocate (effl_outflow_2d(1:ndatamax,1:npoints,1:ny)) ;
   do i = 1, npoints 
     effl_outflow_2d(1:ndatamax,i,1) = effl_outflow(1:ndatamax)
@@ -408,7 +409,7 @@
       & (j,hour0, tinteg, dt, height_T_q, height_u, interval, &
       &  h10, select_h10, l10, ls10, hs10, Ts0, Tb0, Tbb0, h_ML0, &
       &  extwat, select_extwat, extice, &
-      &  kor, trib_inflow, effl_outflow, &
+      &  kor, trib_inflow(1), effl_outflow, &
       &  widthchan, lengthchan, hbot, hbotchan, &
       &  Sals0, Salb0, fetch, phi, lam, us0, vs0, &
       &  Tm, alphax, alphay, a_veg, c_veg, &
@@ -434,7 +435,7 @@
       extwat_2d(:,j) = extwat
       extice_2d(:,j) = extice
       kor_2d(:,j) = kor
-      trib_inflow_2d(:,j) = trib_inflow
+      trib_inflow_2d(:,j) = trib_inflow(1)
       do i = 1, npoints
         effl_outflow_2d(1:ndatamax,i,j) = effl_outflow(1:ndatamax)
       enddo
@@ -588,8 +589,8 @@
    time = time + dt
 
    if (vartrib_inflow) then ! Inflow variable with time
-     call TRIB_INFLOW_UPDATE(time,interval,dataname,trib_inflow)
-     trib_inflow_2d(:,:) = trib_inflow
+     call TRIB_INFLOW_UPDATE(time,interval,dataname,trib_inflow(1))
+     trib_inflow_2d(:,:) = trib_inflow(1)
    endif
 
    if (forc_format == 0) then ! Reading input file in ASCII format
@@ -745,16 +746,21 @@
 
      iycyc : do iy = 1, ny
 
-       tribinfl = trib_inflow_2d(ix,iy) - &
-       & OUTFLOW_DISCHARGE(ndatamax,effl_outflow_2d(1,ix,iy),h_2d(ix,iy)) + &
-       & dhdt(iy) ! Adding interaction of two lakes
+       tribinfl = trib_inflow_2d(ix,iy) 
+       if (lakinterac > 1) then
+         tribinfl = tribinfl - &
+         & OUTFLOW_DISCHARGE(ndatamax,effl_outflow_2d(1,ix,iy),h_2d(ix,iy)) + &
+         & dhdt(iy) ! Adding interaction of two lakes
+       endif
+       trib_inflow(1)             = tribinfl
+       trib_inflow(2:nvartribext) = missing_value
 
        call LAKE(ta(ix),qa(ix),pa(ix),ua(ix),va(ix),atm_rad(ix),Rad_bal(ix), &
        &  solar_rad(ix),precip(ix),Sflux(ix), &
        &  ch4_pres_atm0, co2_pres_atm0, o2_pres_atm0, zref, dt, &
        &  h10_2d(ix,iy), l10_2d(ix,iy), ls10_2d(ix,iy), hs10_2d(ix,iy), Ts0_2d(ix,iy), &
        &  Tb0_2d(ix,iy), Tbb0_2d(ix,iy), h_ML0_2d(ix,iy), extwat_2d(ix,iy), &
-       &  extice_2d(ix,iy), kor_2d(ix,iy), tribinfl, Sals0_2d(ix,iy), &
+       &  extice_2d(ix,iy), kor_2d(ix,iy), trib_inflow, Sals0_2d(ix,iy), &
        &  Salb0_2d(ix,iy), fetch_2d(ix,iy), phi_2d(ix,iy), lam_2d(ix,iy), &
        &  us0_2d(ix,iy), vs0_2d(ix,iy), Tm_2d(ix,iy), alphax_2d(ix,iy), &
        &  alphay_2d(ix,iy), a_veg_2d(ix,iy), c_veg_2d(ix,iy), h_veg_2d(ix,iy), &
diff --git a/source/model/bathym_mod.f90 b/source/model/bathym_mod.f90
index a25a071..e1a70f1 100644
--- a/source/model/bathym_mod.f90
+++ b/source/model/bathym_mod.f90
@@ -13,7 +13,7 @@ SUBROUTINE BATHYMETRY(M,Mice,ns,ix,iy,ndatamax,month,day,lakeform,hour,dt, &
 use LAKE_DATATYPES, only : &
 & ireals, iintegers
 
-use DRIVING_PARAMS, only : soilcolconjtype
+use DRIVING_PARAMS, only : soilcolconjtype, nvartribext
 
 use PHYS_CONSTANTS, only : g, row0
 
@@ -65,7 +65,7 @@ real(kind=ireals),       intent(in) :: area_lake, cellipt
 
 logical, intent(in) :: multsoil
 
-real(kind=ireals), intent(in) :: trib_inflow
+real(kind=ireals), intent(in) :: trib_inflow(1:nvartribext)
 real(kind=ireals), intent(inout) :: dhwtrib
 
 real(kind=ireals), intent(out) :: vol, botar
@@ -349,7 +349,7 @@ if (water_abstraction .and. month == 12 .and. day == 31 &
 & .and. (24. - hour)*hour_sec <= dt) then
 ! Water abstraction for the next year, m**3/s
   watabstr = ABSTR(vol)
-  if (trib_inflow /= -9999.) then
+  if (trib_inflow(1) /= -9999.) then
     dhwtrib = dhwtrib - watabstr*dt/area_int(1)
   endif
 endif
diff --git a/source/model/driving_params_mod.f90 b/source/model/driving_params_mod.f90
index 060e546..f7499b0 100644
--- a/source/model/driving_params_mod.f90
+++ b/source/model/driving_params_mod.f90
@@ -73,6 +73,7 @@ real(kind=ireals) :: dttribupdate = missing_value
 
 !     The group of tributaries characteristics
 type(intpar), target :: N_tribin
+integer(kind=iintegers), parameter :: nvartribext = 4 !Driving parameter, not provided in user's files
 integer(kind=iintegers), allocatable :: itribloc(:)
 type(intpar), target :: N_triblev 
 integer(kind=iintegers) :: N_tribout  = 0!; logical :: ok_N_tribout = .false.
@@ -195,9 +196,8 @@ logical :: all_par_present = .true.
 
 character(len=40) :: path ; logical :: ok_path = .false.
 character(len=60) :: setupfile ! ; logical :: ok_setupfile = .false.
-character(len=40) :: fileinflow='filenamenotgiven', &
-& fileoutflow='filenamenotgiven'
-
+character(len=40) :: fileinflow  = 'filenamenotgiven', &
+&                    fileoutflow = 'filenamenotgiven'
 
 SAVE
 
@@ -211,11 +211,12 @@ implicit none
 integer(kind=iintegers) :: i ! Loop index
 character(len=200) :: line
 logical :: firstcall
-data firstcall, line/.true., 'begin file'/
+data firstcall, line /.true., 'begin file'/
 
 
 if (firstcall) then
 
+  !Driving parameters, not provided in the user's input files
   soilcolconjtype = 2
 
   dertypepar_real(1)%p => dt_out
diff --git a/source/model/lake.f90 b/source/model/lake.f90
index c86a384..16f371b 100644
--- a/source/model/lake.f90
+++ b/source/model/lake.f90
@@ -313,7 +313,7 @@ real(kind=ireals), intent(in) :: extwat, extice
 real(kind=ireals), intent(in) :: kor_in
 
 
-real(kind=ireals), intent(in) :: trib_inflow !> Tributary inflow discharge, m**3/s
+real(kind=ireals), intent(in) :: trib_inflow(1:nvartribext) !> The group of tributary data: discharge, ...
 real(kind=ireals), intent(in) :: Sals0, Salb0
 real(kind=ireals), intent(in) :: fetch
 real(kind=ireals), intent(in) :: phi, lam
@@ -343,8 +343,8 @@ integer(kind=iintegers), intent(in) :: init_T
 !> MPI parameters
 integer(kind=iintegers), intent(in) :: comm3d, rank_comm3d, coords(1:3) 
 !> Control point (CP) switch: 0 - do nothing
-!!                       1 - write CP                             
-!!                       2 - read CP as initial conditions                                           
+!!                            1 - write CP                             
+!!                            2 - read CP as initial conditions                                           
 integer(kind=iintegers), intent(in) :: icp 
 
                                            
@@ -753,10 +753,9 @@ endif
 call BATHYMETRY(M,Mice,ns,ix,iy,ndatamax,month,day,lakeform,hour,dt, &
                & depth_area,area_lake,cellipt, &
                & multsoil,trib_inflow,dhwtrib,vol,botar)
-! Tributary inflow
-if (trib_inflow /= -9999.) then
-  dhwtrib = trib_inflow*dt/area_int(1)
-  ! Will be corrected for water abstraction in BATHYMETRY 
+! Tributary inflow (may be overwritten later in BATHYMETRY and/or in TRIBTEMP)
+if (trib_inflow(1) /= -9999.) then
+  dhwtrib = trib_inflow(1)*dt/area_int(1)
 else
   dhwtrib = (h10 - h1)/10. ! 10 - arbitrary value, the "time" of depth damping towards initial depth
 endif
@@ -835,8 +834,9 @@ turb_density_flux = TURB_DENS_FLUX(eflux0_kinem,0.e0_ireals,Tw1(1),0.e0_ireals)
 Buoyancy0 = g/row0*turb_density_flux
 
 ! Updating tributaries data and adding inflow of all substances
-if (tribheat%par > 0) then
-  call TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spinup_done)
+if (tribheat%par > 0 .and. .not. (tribheat%par == 4 .and. trib_inflow(1) == -9999.) ) then
+  call TRIBTEMP(time,dt,h1,h10,dhwtrib,trib_inflow, &
+  &             z_full,area_int,area_half,gsp,gas,wst,spinup_done)
   ! Adding tendency due to mean vertical velocity
   call VERTADV_UPDATE(M,dt,h1,gsp,bathymwater,wst,gas)
 endif
@@ -1180,8 +1180,9 @@ PEMF = 0.e0_ireals
 pt_down_f = 0.e0_ireals
 
 ! Updating tributaries data and adding inflow of all substances
-if (tribheat%par > 0) then
-  call TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spinup_done)
+if (tribheat%par > 0 .and. .not. (tribheat%par == 4 .and. trib_inflow(1) == -9999.)) then
+  call TRIBTEMP(time,dt,h1,h10,dhwtrib,trib_inflow, &
+  &             z_full,area_int,area_half,gsp,gas,wst,spinup_done)
   ! Adding tendency due to mean vertical velocity
   call VERTADV_UPDATE(M,dt,h1,gsp,bathymwater,wst,gas)
 endif
diff --git a/source/model/trib.f90 b/source/model/trib.f90
index f5c2734..f72a5f1 100644
--- a/source/model/trib.f90
+++ b/source/model/trib.f90
@@ -24,7 +24,8 @@ contains
 !!Tw --- the temperature profile in lake, C
 !!Sal --- salinity, kg/kg
 !!(u,v) --- velocity components, m/s
-SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spinup_done)
+SUBROUTINE TRIBTEMP(time,dt,h1,h10,dhwtrib,trib_inflow, &
+&                   z_full,area_int,area_half,gsp,gas,wst,spinup_done)
 
  use DRIVING_PARAMS , only : &
  & N_tribin, & 
@@ -50,7 +51,8 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
  & missing_value, &
  & outflpar, &
  & tribheat, &
- & carbon_model
+ & carbon_model, &
+ & nvartribext
 
  use ARRAYS, only : gas_type
 
@@ -58,8 +60,9 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
 
  implicit none
 
- real(kind=ireals), intent(in)    :: time, dt, h1 
+ real(kind=ireals), intent(in)    :: time, dt, h1, h10
  real(kind=ireals), intent(inout) :: dhwtrib
+ real(kind=ireals), intent(in)    :: trib_inflow(1:nvartribext)
  real(kind=ireals), intent(in)    :: z_full(1:M+1)
  real(kind=ireals), intent(in)    :: area_int(1:M+1), area_half(1:M)
 
@@ -98,7 +101,7 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
    ! Allocation of inflow and outflow profiles
    if (tribheat%par == 1 .or. tribheat%par == 2) then
      itriblev = N_triblev%par
-   elseif (tribheat%par == 3) then
+   elseif (tribheat%par == 3 .or. tribheat%par == 4) then
      itriblev = M+1
    endif
    allocate (width_tribin  (1:N_tribin%par ,1:itriblev) )
@@ -118,20 +121,23 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
    allocate (U_tribout     (1:N_tribout    ,1:itriblev) ) 
  endif 
 
- ! 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,wst,area_int,h1,dhwtrib,dt,level)
+ ! 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,wst,area_int,h1,dhwtrib,trib_inflow,dt,level)
 
- !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
+   work = mod(time - dt, dttribupdate*day_sec) ! time is for next time step, thus lateral inflow
+                                               ! values are updated to the current time step
    if ((work >= 0.d0 .and. work < dt) .or. firstcall .or. spinup_done) then
      ! Updating tributary data
-     call TRIBUPDATE(spinup_done,gsp,wst,area_int,h1,dhwtrib,dt,level)
+     call TRIBUPDATE(spinup_done,gsp,wst,area_int,h1,dhwtrib,trib_inflow,dt,level)
    endif
  endif
 
+ if (tribheat%par == 4) &
+ &  call TRIBUPDATE(spinup_done,gsp,wst,area_int,h1,dhwtrib,trib_inflow,dt,level)
+
  allocate( U_tribin_(1:N_tribin%par,1:M+1), T_tribin_(1:N_tribin%par,1:M+1), &
  & Sal_tribin_(1:N_tribin%par,1:M+1), Ux_tribin_(1:N_tribin%par,1:M+1), &
  & Uy_tribin_(1:N_tribin%par,1:M+1), U_tribout_(1:N_tribout,1:M+1), &
@@ -153,22 +159,34 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
      !call LININTERPOL(z_,U_tribout(k,1:itriblev),itriblev,z_full,U_tribout_(k,1:M+1),M+1,flag)
      !call LININTERPOL(z_,width_tribin(k,1:itriblev),itriblev,z_full,width_tribin_(k,1:M+1),M+1,flag)
      !call LININTERPOL(z_,width_tribout(k,1:itriblev),itriblev,z_full,width_tribout_(k,1:M+1),M+1,flag)
-     call PIECEWISECONSTINT(z_,U_tribin    (k,1:itriblev),itriblev,z_full,U_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,T_tribin    (k,1:itriblev),itriblev,z_full,T_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,Sal_tribin  (k,1:itriblev),itriblev,z_full,Sal_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,meth_tribin (k,1:itriblev),itriblev,z_full,meth_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,oxyg_tribin (k,1:itriblev),itriblev,z_full,oxyg_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,DIC_tribin  (k,1:itriblev),itriblev,z_full,DIC_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,DOC_tribin  (k,1:itriblev),itriblev,z_full,DOC_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,POC_tribin  (k,1:itriblev),itriblev,z_full,POC_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,Ux_tribin   (k,1:itriblev),itriblev,z_full,Ux_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,Uy_tribin   (k,1:itriblev),itriblev,z_full,Uy_tribin_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,width_tribin(k,1:itriblev),itriblev,z_full,width_tribin_(k,1:M+1),M+1,2_iintegers)
+     call PIECEWISECONSTINT(z_    , U_tribin (k,1:itriblev),itriblev, &
+     &                      z_full, U_tribin_(k,1:M+1     ),M+1     , 2_iintegers)
+     call PIECEWISECONSTINT(z_    ,T_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,T_tribin_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,Sal_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,Sal_tribin_(k,1:M+1      ),M+1    ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,meth_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,meth_tribin_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,oxyg_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,oxyg_tribin_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,DIC_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,DIC_tribin_(k,1:M+1)     ,M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,DOC_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,DOC_tribin_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,POC_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,POC_tribin_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,Ux_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,Ux_tribin_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,Uy_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,Uy_tribin_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,width_tribin (k,1:itriblev),itriblev, &
+     &                      z_full,width_tribin_(k,1:M+1)     ,M+1     ,2_iintegers)
    enddo
    do k = 1, N_tribout
-     call PIECEWISECONSTINT(z_,U_tribout(k,1:itriblev),itriblev,z_full,U_tribout_(k,1:M+1),M+1,2_iintegers)
-     call PIECEWISECONSTINT(z_,width_tribout(k,1:itriblev),itriblev,z_full,width_tribout_(k,1:M+1),M+1,2_iintegers)
-     !print*, U_tribout, U_tribout_
+     call PIECEWISECONSTINT(z_    ,U_tribout (k,1:itriblev),itriblev, &
+     &                      z_full,U_tribout_(k,1:M+1     ),M+1     ,2_iintegers)
+     call PIECEWISECONSTINT(z_    ,width_tribout (k,1:itriblev),itriblev, &
+     &                      z_full,width_tribout_(k,1:M+1     ),M+1     ,2_iintegers)
    enddo 
    deallocate(z_)
  else
@@ -277,11 +295,10 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
      do k = 1, N_tribout
        outflow = outflow + U_tribout_(k,i)*width_tribout_(k,i)
      enddo
-     wst%wArea(i-1) = wst%wArea(i) - gsp%ddz05(i-1)*h1*(inflow - outflow) !Solution of horizontally averaged continuity
-     !print*, i, wst%wArea(i-1), wst%wArea(i), gsp%ddz05(i),h1,inflow , outflow !Solution of horizontally averaged continuity
+     !Solution of horizontally averaged continuity equation
+     wst%wArea(i-1) = wst%wArea(i) - gsp%ddz05(i-1)*h1*(inflow - outflow) 
      outflow = outflow/area_int(i)
      outflow = max(outflow,0.) !Ensuring positivity of outflow discharge, used for advection only
-     !print*, 'outflow', outflow, U_tribin_(1,i),T_tribin_(1,i),width_tribin_(1,i), area_int(i), invdt
 
      meanTin = 0.
      meansalin = 0.
@@ -396,15 +413,18 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
  end select 
  foutflow = foutflow/area_int(1) !Methane flux is normalized by lake surface area
 
- !Diagnosing vertical speed
+ !Diagnosing horizontally averaged vertical speed
  wst%w(0) = wst%wArea(0)/area_int(1)
  wst%w(M+1) = wst%wArea(M+1)/area_int(M+1)
  forall(i = 1:M) wst%w(i) = wst%wArea(i)/area_half(i)
 
+ !Correction of dhwtrib (level change due to tributaries)
  if (level /= missing_value) then
    dhwtrib = (sum(disch_tribin(:,:)) - sum(disch_tribout(:,:)))* &
    & dt/area_int(1)
  endif
+ !Currently, tribheat = 4 includes nudging water level to initial value (constant)
+ if (tribheat%par == 4) dhwtrib = (h10 - h1)/10.
 
  deallocate(U_tribin_, T_tribin_, Sal_tribin_, Ux_tribin_, Uy_tribin_, U_tribout_, &
  & width_tribin_, width_tribout_, DIC_tribin_, DOC_tribin_, &
@@ -416,7 +436,7 @@ SUBROUTINE TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spi
 
 
 !> Subroutine TRIBUPDATE updates tributary data
-SUBROUTINE TRIBUPDATE(spinup_done,gsp,wst,area_int,h1,dhwtrib,dt,level)
+SUBROUTINE TRIBUPDATE(spinup_done,gsp,wst,area_int,h1,dhwtrib,trib_inflow,dt,level)
 
 use INOUT_PARAMETERS, only : &
 & lake_misc_unit_min, &
@@ -448,6 +468,7 @@ use DRIVING_PARAMS, only : &
 & fileinflow, &
 & fileoutflow, &
 & dttribupdate, &
+& nvartribext, &
 & READPROFILE
 
 use INOUT, only : CHECK_UNIT
@@ -465,9 +486,10 @@ type(waterstate_type) , intent(in) :: wst
 real(kind=ireals), intent(in) :: area_int(1:M+1)
 real(kind=ireals), intent(in) :: h1 
 real(kind=ireals), intent(inout) :: dhwtrib
+real(kind=ireals), intent(in) :: trib_inflow(1:nvartribext)
 real(kind=ireals), intent(in) :: dt
 !Output variables
-real(kind=ireals), intent(out) :: level
+real(kind=ireals), intent(inout) :: level
 
 ! Local variables
 integer(kind=iintegers), parameter :: mixtype = 0 !> Indicates the type of mixing of tributaries water at the
@@ -483,7 +505,7 @@ integer(kind=iintegers) :: ncol, i, j, idep, idep0
 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)
@@ -523,37 +545,25 @@ if_tribform : if (tribheat%par == 2) then
   width_tribout(1,:) = work (:,2)
   U_tribout    (1,:) = work (:,3)
   deallocate(work)
+
 else if (tribheat%par == 3) then
+
   ! Reading bulk characteristics of inlets and an outlet.
   ! The data for each inlet in an input file should be given in rows containing:
   ! 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'
 
   disch_inflow = 0.
   do i = 1, N_tribin%par
 
     read(nunit1,*) xx, disch, depth, width, temp, sal, &
     & dirprojX, dirprojY, DOC, POC, DIC, CH_4, O_2
-    !print*, xx, disch, depth, width, temp, sal, &
-    !& dirprojX, dirprojY, DOC, POC, DIC, CH_4
-    !disch = 0.
+
     if (mixtype == 1 .and. temp /= missing_value) then
       ! Assuming that inlet properties are advected at levels from surface to the lake depth 
       ! where the water density equals to that of an inlet
-      allocate(work(1:M+1,1))
-      xx = sal
-      if (sal == missing_value) xx = 0.
-      call DENSITY_W(0_iintegers,eos%par,lindens%par,(/temp/),(/xx/),(/pref/),work(1,1))
-      rowtrib = work(1,1)
-      work(1:M+1,1) = abs(wst%row(1:M+1) - rowtrib)
-      j = minloc(work(1:M+1,1),1)
-      depth = gsp%z_half(j)
-      !print*, work
-      !print*, 'riv depth', j
-      !read*
-      deallocate(work)
+      depth = DEPTH_TRIB_MIX()
     endif
 
     idep = 1
@@ -564,10 +574,10 @@ else if (tribheat%par == 3) then
     if (idep /= 1_iintegers) idep = idep - 1
     depth = gsp%z_half(idep)
 
-    width_tribin(i,1:idep) = width                         ; 
-    U_tribin    (i,1:idep) = disch/(depth*width)           ; 
-    T_tribin    (i,1:idep) = temp                          ; 
-    Sal_tribin  (i,1:idep) = sal                           ; 
+    width_tribin(i,1:idep) = width                          
+    U_tribin    (i,1:idep) = disch/(depth*width)            
+    T_tribin    (i,1:idep) = temp                           
+    Sal_tribin  (i,1:idep) = sal                            
     if (dirprojX /= missing_value) then
       Ux_tribin   (i,1:idep) = disch/(depth*width)*dirprojX 
     else
@@ -578,11 +588,11 @@ else if (tribheat%par == 3) then
     else
       Uy_tribin   (i,1:idep) = missing_value
     endif
-    DOC_tribin  (i,1:idep) = DOC                           ; 
-    POC_tribin  (i,1:idep) = POC                           ; 
-    DIC_tribin  (i,1:idep) = DIC                           ; 
-    meth_tribin (i,1:idep) = CH_4                          ;     
-    oxyg_tribin (i,1:idep) = O_2                           ;     
+    DOC_tribin  (i,1:idep) = DOC                            
+    POC_tribin  (i,1:idep) = POC                            
+    DIC_tribin  (i,1:idep) = DIC                            
+    meth_tribin (i,1:idep) = CH_4                               
+    oxyg_tribin (i,1:idep) = O_2                                
     if (idep < M+1) then
       width_tribin(i,idep+1:M+1) = 0.
       U_tribin    (i,idep+1:M+1) = 0.
@@ -625,7 +635,7 @@ else if (tribheat%par == 3) then
   if (idep /= 1_iintegers) idep = idep - 1
   depth = gsp%z_half(idep)
 
-  !If data on the water level is available, outlet discharge is chosen to ensure
+  !If data on the water level is available, outlet discharge is chosen to keep
   !this level (i.e. as a residual of water balance equation)
   if (level /= missing_value) then
     disch = - area_int(1)*(level-h1)/(dttribupdate*day_sec)+ disch_inflow
@@ -642,9 +652,132 @@ else if (tribheat%par == 3) then
     U_tribout    (1,idep+1:M+1) = 0.
   endif
 
+else if (tribheat%par == 4) then
+
+  ! Taking bulk characteristics of a single inlet from trib_inflow, provided by the LAKE model driver.
+  ! The data for the inlet is discharge, depth, width, temperature, ... (possibly other scalars)
+
+  if (N_tribin%par /= 1_iintegers) print*, 'N_tribin must be 1 under tribheat = 4'
+
+  disch = trib_inflow(1)
+  depth = trib_inflow(2)
+  width = trib_inflow(3)
+  temp  = trib_inflow(4)
+  dirprojX = 1. !assuming the inflow directs eastward
+  dirprojY = 0.
+  !Other scalars, which can be added:
+  !sal, dirprojX, dirprojY, DOC, POC, DIC, CH_4, O_2
+  sal = missing_value
+  DOC = missing_value
+  POC = missing_value
+  DIC = missing_value
+  CH_4 = missing_value
+  O_2 = missing_value
+
+  if (mixtype == 1 .and. temp /= missing_value) then
+    ! Assuming that inlet properties are advected at levels from surface to the lake depth 
+    ! where the water density equals to that of an inlet
+    depth = DEPTH_TRIB_MIX()
+  endif
+
+  idep = 1
+  do while (gsp%z_half(idep)/depth < 1.) 
+    idep = idep + 1
+    if (idep == M+1) exit
+  enddo
+  if (idep /= 1_iintegers) idep = idep - 1
+  depth = gsp%z_half(idep)
+
+  i = 1
+  width_tribin(i,1:idep) = width                          
+  U_tribin    (i,1:idep) = disch/(depth*width)            
+  T_tribin    (i,1:idep) = temp                           
+  Sal_tribin  (i,1:idep) = sal                            
+  if (dirprojX /= missing_value) then
+    Ux_tribin   (i,1:idep) = disch/(depth*width)*dirprojX 
+  else
+    Ux_tribin   (i,1:idep) = missing_value
+  endif
+  if (dirprojY /= missing_value) then
+    Uy_tribin   (i,1:idep) = disch/(depth*width)*dirprojY
+  else
+    Uy_tribin   (i,1:idep) = missing_value
+  endif
+  DOC_tribin  (i,1:idep) = DOC                            
+  POC_tribin  (i,1:idep) = POC                            
+  DIC_tribin  (i,1:idep) = DIC                            
+  meth_tribin (i,1:idep) = CH_4                               
+  oxyg_tribin (i,1:idep) = O_2                                
+  if (idep < M+1) then
+    width_tribin(i,idep+1:M+1) = 0.
+    U_tribin    (i,idep+1:M+1) = 0.
+    T_tribin    (i,idep+1:M+1) = 0.
+    Sal_tribin  (i,idep+1:M+1) = 0.
+    Ux_tribin   (i,idep+1:M+1) = 0.
+    Uy_tribin   (i,idep+1:M+1) = 0.
+    DOC_tribin  (i,idep+1:M+1) = 0.
+    POC_tribin  (i,idep+1:M+1) = 0.
+    DIC_tribin  (i,idep+1:M+1) = 0.
+    meth_tribin (i,idep+1:M+1) = 0.
+    oxyg_tribin (i,idep+1:M+1) = 0.
+  endif
+  disch_inflow = disch
+
+  depth0 = 0.
+  disch = disch_inflow ! The inflow and outflow are equal, ensuring constant level
+  !Depth of upper bound of outlet stream
+  idep0 = 1
+  !do while (gsp%z_half(idep0)/(depth0 + small_value) < 1.) 
+  !  idep0 = idep0 + 1
+  !  if (idep0 == M+1) exit
+  !enddo
+  !if (idep0 /= 1_iintegers) idep0 = idep0 - 1
+  !if (idep0 == 1) then
+  !  depth0 = 0.
+  !else
+  !  depth0 = gsp%z_half(idep0-1)
+  !endif
+  !Depth of lower bound of outlet stream
+  !idep = 1
+  !do while (gsp%z_half(idep)/depth < 1.) 
+  !  idep = idep + 1
+  !  if (idep == M+1) exit
+  !enddo
+  !if (idep /= 1_iintegers) idep = idep - 1
+  !depth = gsp%z_half(idep)
+
+  width_tribout(1,idep0:idep) = width  
+  U_tribout    (1,idep0:idep) = disch/((depth-depth0)*width)
+  if (idep0 > 1) then
+    width_tribout(1,1:idep0-1) = 0.
+    U_tribout    (1,1:idep0-1) = 0.
+  endif
+  if (idep < M+1) then
+    width_tribout(1,idep+1:M+1) = 0.
+    U_tribout    (1,idep+1:M+1) = 0.
+  endif
+
 endif if_tribform
 
 if (firstcall) firstcall = .false.
+contains 
+
+!> Computing the lake depth 
+!! where the water density equals to that of an inlet
+FUNCTION DEPTH_TRIB_MIX() result(depth_mix)
+implicit none
+real(kind=ireals) :: depth_mix
+allocate(work(1:M+1,1))
+xx = sal
+if (sal == missing_value) xx = 0.
+call DENSITY_W(0_iintegers,eos%par,lindens%par,(/temp/),(/xx/),(/pref/),work(1,1))
+rowtrib = work(1,1)
+work(1:M+1,1) = abs(wst%row(1:M+1) - rowtrib)
+j = minloc(work(1:M+1,1),1)
+depth_mix = gsp%z_half(j)
+deallocate(work)
+END FUNCTION DEPTH_TRIB_MIX
+
 END SUBROUTINE TRIBUPDATE 
 
 !> Calculates artificial abstraction rate (m**3/s) 
-- 
GitLab