Skip to content
Snippets Groups Projects
Commit 24ad8a2f authored by Victor Stepanenko's avatar Victor Stepanenko
Browse files

An option tribheat=4 is added to take into account tributaries in large scale models (LSMs, ESMs)

parent 6548bc92
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
source/driver/driver.f90
+ 18
− 12
View file @ 24ad8a2f
......@@ -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) - &
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), &
......
source/model/bathym_mod.f90
+ 3
− 3
View file @ 24ad8a2f
......@@ -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
......
source/model/driving_params_mod.f90
+ 5
− 4
View file @ 24ad8a2f
......@@ -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.
......@@ -198,7 +199,6 @@ character(len=60) :: setupfile ! ; logical :: ok_setupfile = .false.
character(len=40) :: fileinflow = 'filenamenotgiven', &
& fileoutflow = 'filenamenotgiven'
SAVE
contains
......@@ -216,6 +216,7 @@ 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
......
source/model/lake.f90
+ 12
− 11
View file @ 24ad8a2f
......@@ -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
......@@ -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
......
source/model/trib.f90
+ 191
− 58
View file @ 24ad8a2f
......@@ -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
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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment