Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DOM flow in Mosart #64

Draft
wants to merge 38 commits into
base: master
Choose a base branch
from
Draft
Changes from 1 commit
Commits
Show all changes
38 commits
Select commit Hold shift + click to select a range
9a1474f
add DOM mass balance code module
devarajun May 7, 2021
02ebfba
call DOM module
devarajun May 8, 2021
bdb1e45
Declare a type for DOM variables
devarajun May 8, 2021
f0e77af
Declare dom variables
devarajun May 8, 2021
3b0d992
add DOC/DON variables to output
devarajun May 16, 2021
dc541fb
add DOC, DON as tracer index
devarajun May 21, 2021
bd2272d
add DOM mass balance subroutines
devarajun May 24, 2021
bfd5fc2
include if condition to avoid division by zero
devarajun May 24, 2021
1089202
correction to include Type for wh
devarajun May 24, 2021
22196b0
update domRUp variable
devarajun Jun 3, 2021
22fa84c
comma addition b/w 1 and cnt
devarajun Jun 3, 2021
2233bbf
correction to eroutUp
devarajun Jun 7, 2021
691552d
add file pycache
devarajun Jun 7, 2021
f11e770
Merge branch 'dev/dom' of https://github.com/MetOs-UiO/MOSART into mo…
Oct 11, 2022
fbf614d
add doc tracer to mosart in new array first step
Nov 18, 2022
d4d9c2b
changes to mosart for DOC
Jan 5, 2023
ae92015
removed domRout and domRin
Jan 5, 2023
68e31f0
corrections to the code, typos and mistakes
Jan 6, 2023
3c18a69
changes to units , kg instead of g and corrections
Jan 10, 2023
9b245f7
comments added
Jan 10, 2023
4d1e1c4
add output concentrations of DOM, corrections of units
Jan 12, 2023
06827d4
sign corrections and fixes
Jan 16, 2023
0123953
various fixes
Jan 20, 2023
bf84861
testing
Jan 20, 2023
1b8ca3b
add surface and subsurface DOM
Jan 23, 2023
ba0506b
PRINT
Jan 31, 2023
03c7f89
unit changes
Feb 1, 2023
6439c7b
correct mass balance equations
Feb 2, 2023
9ce0d68
various changes
Feb 15, 2023
f7653d9
still too much going out of subn
Mar 13, 2023
d1fc844
various changes but one error
Mar 18, 2023
c38038c
changes
Mar 20, 2023
475c98e
changes made it worse
Mar 29, 2023
bdf66e2
corrections and description
Apr 19, 2023
7fa46ae
commented out a bunch of boundaries
May 2, 2023
dae8a26
remove the last boundary statement of DOC mass balance
May 3, 2023
bbd2aa7
added text
Jun 12, 2023
b2debee
starting to remove equations from dommas in hill and subn
Jun 20, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
correct mass balance equations
  • Loading branch information
Marius Lambert committed Feb 2, 2023
commit 6439c7b2ccde48d5c2baeaec7e50029392ce3aab
18 changes: 9 additions & 9 deletions src/riverroute/DommasbMod.F90
Original file line number Diff line number Diff line change
@@ -20,32 +20,32 @@ MODULE DommasbMod
contains

!----------------------------------------------------------------------
subroutine hillslopeRoutingDOM(iunit,nt,ntdom,theDeltaT,Darea,Dfrac)
subroutine hillslopeRoutingDOM(iunit,nt,ntdom,theDeltaT,Pwh)
! ! DESCRIPTION: solve the ODEs with Euler algorithm for hillslope routing
implicit none
integer, intent(in) :: iunit, nt, ntdom
real(r8), intent(in) :: theDeltaT, Darea, Dfrac
real(r8), intent(in) :: theDeltaT,Pwh
! assume no chemical reaction in the water hence sink term is zero implies domsur = domR*flow
! ehout is negative
Tdom%domH(iunit,ntdom) = Tdom%domH(iunit,ntdom) + (TRunoff%ehout(iunit,nt) * Darea * Dfrac * Tdom%domH(iunit,ntdom) + Tdom%domsur(iunit,ntdom)) * theDeltaT/(TRunoff%wh(iunit,nt)*Darea*Dfrac)
Tdom%domH(iunit,ntdom) = (Tdom%domH(iunit,ntdom)*Pwh + TRunoff%ehout(iunit,nt) * Tdom%domH(iunit,ntdom) + Tdom%domsur(iunit,ntdom)) * theDeltaT/TRunoff%wh(iunit,nt)
end subroutine hillslopeRoutingDOM

subroutine subnetworkRoutingDOM(iunit,nt,ntdom,theDeltaT)
subroutine subnetworkRoutingDOM(iunit,nt,ntdom,theDeltaT,Pwt)
! solve the ODEs with Euler algorithm for subnetwork routing
! etin is positive and etout is negative
implicit none
integer, intent(in) :: iunit, nt, ntdom
real(r8), intent(in) :: theDeltaT
Tdom%domT(iunit,ntdom) = Tdom%domT(iunit,ntdom) + (Tdom%domsub(iunit,ntdom) + TRunoff%etin(iunit,nt) * Tdom%domH(iunit,ntdom) + TRunoff%etout(iunit,nt) * Tdom%domT(iunit,ntdom)) * theDeltaT/TRunoff%wt(iunit,nt)
real(r8), intent(in) :: theDeltaT,Pwt
Tdom%domT(iunit,ntdom) = (Tdom%domT(iunit,ntdom)*Pwt + Tdom%domsub(iunit,ntdom) + TRunoff%etin(iunit,nt) * Tdom%domH(iunit,ntdom) + TRunoff%etout(iunit,nt) * Tdom%domT(iunit,ntdom)) * theDeltaT/TRunoff%wt(iunit,nt)
end subroutine subnetworkRoutingDOM

subroutine mainchannelRoutingDOM(iunit,nt,ntdom,theDeltaT)
subroutine mainchannelRoutingDOM(iunit,nt,ntdom,theDeltaT,Pwr)
! solve the ODE with Euler algorithm for main-channel routing
! erout is negative, while erlateral and erin are positive
implicit none
integer, intent(in) :: iunit, nt, ntdom
real(r8), intent(in) :: theDeltaT
Tdom%domR(iunit,ntdom) = Tdom%domR(iunit,ntdom) + (TRunoff%erlateral(iunit,nt)*Tdom%domT(iunit,ntdom) + Tdom%domRUp(iunit,ntdom) + TRunoff%erout(iunit,nt)*Tdom%domR(iunit,ntdom))*theDeltaT/TRunoff%wr(iunit,nt)
real(r8), intent(in) :: theDeltaT,Pwr
Tdom%domR(iunit,ntdom) = (Tdom%domR(iunit,ntdom)*Pwr + TRunoff%erlateral(iunit,nt)*Tdom%domT(iunit,ntdom) + Tdom%domRUp(iunit,ntdom) + TRunoff%erout(iunit,nt)*Tdom%domR(iunit,ntdom))*theDeltaT/TRunoff%wr(iunit,nt)
end subroutine mainchannelRoutingDOM
!-------------------------------------------------------------------------
end MODULE DommasbMod
42 changes: 29 additions & 13 deletions src/riverroute/MOSART_physics_mod.F90
Original file line number Diff line number Diff line change
@@ -49,7 +49,7 @@ subroutine Euler
implicit none

integer :: iunit, m, k, unitUp, cnt, ier !local index
real(r8) :: temp_erout, localDeltaT
real(r8) :: temp_erout, localDeltaT,Pwh,Pwt,Pwr
real(r8) :: negchan

!------------------
@@ -61,13 +61,21 @@ subroutine Euler
if (TUnit%euler_calc(nt)) then
do iunit=rtmCTL%begr,rtmCTL%endr
if(TUnit%mask(iunit) > 0) then
write(iulog,*) 'wh before',TRunoff%wh(iunit,nt)
Pwh=TRunoff%wh(iunit,nt)
call hillslopeRouting(iunit,nt,Tctl%DeltaT)
TRunoff%wh(iunit,nt) = TRunoff%wh(iunit,nt) + TRunoff%dwh(iunit,nt) * Tctl%DeltaT
call UpdateState_hillslope(iunit,nt)
TRunoff%etin(iunit,nt) = (-TRunoff%ehout(iunit,nt) + TRunoff%qsub(iunit,nt)) * TUnit%area(iunit) * TUnit%frac(iunit)
if (TRunoff%wh(iunit, nt) > 0._r8 .and. nt==1) then ! if LIQ tracer and there is water
if (nt==1) then ! if LIQ tracer
do ntdom=1,nt_rtm_dom ! loop over DOM tracers
call hillslopeRoutingDOM(iunit,nt,ntdom,Tctl%DeltaT,TUnit%area(iunit),TUnit%frac(iunit))
write(iulog,*) 'domsur',Tdom%domsur(iunit,ntdom),'wh',TRunoff%wh(iunit,nt),'ehout',TRunoff%ehout(iunit,nt),'qsur',TRunoff%qsur(iunit,nt),'domH',Tdom%domH(iunit,ntdom),'time',Tctl%DeltaT
call hillslopeRoutingDOM(iunit,nt,ntdom,Tctl%DeltaT,Pwh)
Tdom%domsub(iunit,ntdom) = Tdom%domsub(iunit,ntdom) * TUnit%area(iunit) * TUnit%frac(iunit)
write(iulog,*) 'after domH',Tdom%domH(iunit,ntdom), 'iunit',iunit
if (Tdom%domH(iunit,ntdom)>0.3 .or. TRunoff%wh(iunit,nt)<0._r8) then
write(iulog,*) 'SHIT'
end if
end do
endif
endif
@@ -104,17 +112,21 @@ subroutine Euler
if(TUnit%mask(iunit) > 0) then
localDeltaT = Tctl%DeltaT/Tctl%DLevelH2R/TUnit%numDT_t(iunit)
do k=1,TUnit%numDT_t(iunit)
Pwt=TRunoff%wt(iunit,nt)
call subnetworkRouting(iunit,nt,localDeltaT)
TRunoff%wt(iunit,nt) = TRunoff%wt(iunit,nt) + TRunoff%dwt(iunit,nt) * localDeltaT
call UpdateState_subnetwork(iunit,nt)
TRunoff%erlateral(iunit,nt) = TRunoff%erlateral(iunit,nt)-TRunoff%etout(iunit,nt)
if (nt==1) then
do ntdom=1,nt_rtm_dom ! loop over DOM tracers
call subnetworkRoutingDOM(iunit,nt,ntdom,localDeltaT,Pwt)
if (Tdom%domT(iunit,ntdom) >0.3 .or. TRunoff%wt(iunit,nt)<0._r8) then
write(iulog,*) 'SHIT',Tdom%domT(iunit,ntdom),TRunoff%wt(iunit,nt)
end if
end do
endif
end do ! numDT_t
TRunoff%erlateral(iunit,nt) = TRunoff%erlateral(iunit,nt) / TUnit%numDT_t(iunit)
if (TRunoff%wt(iunit,nt) > 0._r8 .and. nt==1) then
do ntdom=1,nt_rtm_dom ! loop over DOM tracers
call subnetworkRoutingDOM(iunit,nt,ntdom,localDeltaT)
end do
endif
endif
end do ! iunit
endif ! euler_calc
@@ -194,6 +206,7 @@ subroutine Euler
localDeltaT = Tctl%DeltaT/Tctl%DLevelH2R/TUnit%numDT_r(iunit)
temp_erout = 0._r8
do k=1,TUnit%numDT_r(iunit)
Pwr=TRunoff%wr(iunit,nt)
call mainchannelRouting(iunit,nt,localDeltaT)
TRunoff%wr(iunit,nt) = TRunoff%wr(iunit,nt) + TRunoff%dwr(iunit,nt) * localDeltaT
! check for negative channel storage
@@ -202,16 +215,19 @@ subroutine Euler
! call shr_sys_abort('mosart: negative channel storage')
! end if
call UpdateState_mainchannel(iunit,nt)
if (nt==1) then
do ntdom=1,nt_rtm_dom ! loop over DOM tracers
call mainchannelRoutingDOM(iunit,nt,ntdom,localDeltaT,Pwr)
if (Tdom%domR(iunit,ntdom) >0.3 .or. TRunoff%wr(iunit,nt)<0._r8) then
write(iulog,*) 'SHIT',Tdom%domR(iunit,ntdom),TRunoff%wr(iunit,nt)
end if
end do
end if
temp_erout = temp_erout + TRunoff%erout(iunit,nt) ! erout here might be inflow to some downstream subbasin, so treat it differently than erlateral
end do
temp_erout = temp_erout / TUnit%numDT_r(iunit)
TRunoff%erout(iunit,nt) = temp_erout
TRunoff%flow(iunit,nt) = TRunoff%flow(iunit,nt) - TRunoff%erout(iunit,nt)
if (TRunoff%wr(iunit,nt) > 0._r8 .and. nt==1) then
do ntdom=1,nt_rtm_dom ! loop over DOM tracers
call mainchannelRoutingDOM(iunit,nt,ntdom,localDeltaT)
end do
end if
endif
end do ! iunit
endif ! euler_calc
23 changes: 15 additions & 8 deletions src/riverroute/RtmMod.F90
Original file line number Diff line number Diff line change
@@ -1327,6 +1327,11 @@ subroutine Rtmini(flood_active)
TRunoff%wt = rtmCTL%wt
TRunoff%wr = rtmCTL%wr
TRunoff%erout= rtmCTL%erout
Tdom%domH = rtmCTL%domH
Tdom%domT = rtmCTL%domT
Tdom%domR = rtmCTL%domR
Tdom%domRUp = rtmCTL%domRUp
write(iulog,*) 'UPDATED MARIUS'
else
! do nt = 1,nt_rtm
! do nr = rtmCTL%begr,rtmCTL%endr
@@ -1821,14 +1826,16 @@ subroutine Rtmrun(rstwr,nlend,rdate)
call t_stopf('mosartr_budget')
! endif

do nt = 1,nt_rtm

do nr = rtmCTL%begr,rtmCTL%endr
do nt = 1,nt_rtm
TRunoff%qsur(nr,nt) = TRunoff%qsur(nr,nt) / rtmCTL%area(nr)
TRunoff%qsub(nr,nt) = TRunoff%qsub(nr,nt) / rtmCTL%area(nr)
TRunoff%qgwl(nr,nt) = TRunoff%qgwl(nr,nt) / rtmCTL%area(nr)
if (nt==1) then
!write(iulog,*) 'MOSART CHECK',Tdom%domsur(nr,1),TRunoff%qsur(nr,1),Tdom%domsub(nr,1),TRunoff%qsub(nr,1)
endif
enddo
do nt = 1,nt_rtm_dom
Tdom%domsur(nr,nt) = Tdom%domsur(nr,nt) / rtmCTL%area(nr)
Tdom%domsub(nr,nt) = Tdom%domsub(nr,nt) / rtmCTL%area(nr)
enddo
enddo

@@ -1907,6 +1914,10 @@ subroutine Rtmrun(rstwr,nlend,rdate)
rtmCTL%wt = TRunoff%wt
rtmCTL%wr = TRunoff%wr
rtmCTL%erout = TRunoff%erout
rtmCTL%domH = Tdom%domH
rtmCTL%domT = Tdom%domT
rtmCTL%domR = Tdom%domR
rtmCTL%domRUp = Tdom%domRUp

do nt = 1,nt_rtm
do nr = rtmCTL%begr,rtmCTL%endr
@@ -1918,10 +1929,6 @@ subroutine Rtmrun(rstwr,nlend,rdate)
rtmCTL%dommas(nr,ntdom)=(TRunoff%wh(nr,nt)*rtmCTL%area(nr)*Tdom%domH(nr,ntdom) + &
TRunoff%wt(nr,nt)*Tdom%domT(nr,ntdom) + &
TRunoff%wr(nr,nt)*Tdom%domR(nr,ntdom))
rtmCTL%domH(nr,ntdom)=Tdom%domH(nr,ntdom)
rtmCTL%domT(nr,ntdom)=Tdom%domT(nr,ntdom)
rtmCTL%domR(nr,ntdom)=Tdom%domR(nr,ntdom)
rtmCTL%domRUp(nr,ntdom)=Tdom%domRUp(nr,ntdom)
enddo
rtmCTL%erin(nr,nt)=TRunoff%erin(nr,nt)
rtmCTL%erlateral(nr,nt)=TRunoff%erlateral(nr,nt)
10 changes: 8 additions & 2 deletions src/riverroute/RtmRestFile.F90
Original file line number Diff line number Diff line change
@@ -443,7 +443,7 @@ subroutine RtmRestart(ncid, flag)
enddo


do nv = 8,10
do nv = 8,11
do ntdom = 1,nt_rtm_dom
if (nv == 8) then
vname = 'RTM_DOMH_'//trim(rtm_tracers(ntdom))
@@ -457,9 +457,14 @@ subroutine RtmRestart(ncid, flag)
dfld => rtmCTL%domT(:,ntdom)
elseif (nv == 10) then
vname = 'RTM_DOMR_'//trim(rtm_tracers(ntdom))
lname = 'DOM storage in main channel in cell'
lname = 'DOM storage in main channel in cell'
uname = 'kg/m3'
dfld => rtmCTL%domR(:,ntdom)
elseif (nv == 11) then
vname = 'RTM_DOMRUP_'//trim(rtm_tracers(ntdom))
lname = 'DOM storage in upstream main channels'
uname = 'kg/m3'
dfld => rtmCTL%domRUp(:,ntdom)
else
write(iulog,*) 'Rtm ERROR: illegal nv value a ',nv
call shr_sys_abort()
@@ -500,6 +505,7 @@ subroutine RtmRestart(ncid, flag)
if (abs(rtmCTL%domH(n,ntdom)) > 1.e30) rtmCTL%domH(n,ntdom) = 0.
if (abs(rtmCTL%domT(n,ntdom)) > 1.e30) rtmCTL%domT(n,ntdom) = 0.
if (abs(rtmCTL%domR(n,ntdom)) > 1.e30) rtmCTL%domR(n,ntdom) = 0.
if (abs(rtmCTL%domRUp(n,ntdom)) > 1.e30) rtmCTL%domR(n,ntdom) = 0.
end do
endif
end do