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

1.Vertical advection added for autochtonous DOC 2.MUSCL scheme added for POCD...

1.Vertical advection added for autochtonous DOC 2.MUSCL scheme added for POCD sedimentation 3. Improvements in plotting scripts
parent 9df02a27
No related branches found
No related tags found
No related merge requests found
......@@ -24,7 +24,7 @@
PREPROCESS_KEY = -cpp
opt_keys = -fopenmp -O3 #-fp-model source
#check_keys = -g -ffpe-trap=invalid,zero,overflow,denormal -fcheck=all -fbacktrace
check_keys = #-g -ffpe-trap=invalid,zero,overflow -fcheck=all -fbacktrace
#check_keys = -g -ffpe-trap=invalid,zero,overflow -fcheck=all -fbacktrace
endif
objfiles_path = ../objfiles/
......@@ -62,13 +62,13 @@
$(objfiles_path)methane_mod.o \
$(objfiles_path)init_var_mod.o \
$(objfiles_path)soil_mod.o \
$(objfiles_path)vertadv_mod.o \
$(objfiles_path)salinity_mod.o \
$(objfiles_path)control_point_mod.o \
$(objfiles_path)carbon_dioxide.o \
$(objfiles_path)phosphorus_mod.o \
$(objfiles_path)trib.o \
$(objfiles_path)bathym_mod.o \
$(objfiles_path)vertadv_mod.o \
\
$(objfiles_path)massflux_convection_v10.o \
$(objfiles_path)convectpar.o \
......
......@@ -40,6 +40,8 @@ use DRIVING_PARAMS, only : carbon_model
use ARRAYS_OXYGEN, only : lampart
use VERTADV, only : MUSCL_FLUXES
implicit none
!> Grid spacing group
......@@ -88,7 +90,7 @@ real(kind=ireals), intent(out) :: Flux_atm_out
! Local variables
real(kind=ireals), parameter :: small_value = 1.d-20
real(kind=ireals), allocatable :: a(:), b(:), c(:), f(:), y(:)
real(kind=ireals), allocatable :: a(:), b(:), c(:), f(:), y(:), work(:), fluxes(:)
real(kind=ireals), allocatable, target :: conc(:)
real(kind=ireals), pointer :: lam(:)
......@@ -215,16 +217,36 @@ do i = 1, nspec
enddo
if (carbon_model%par == 2) then
!Settlement of POCD, explicit (Euler), central-differences scheme
!> @todo: do different sedimentation speed depending on turbulence regime
!> @todo: check scheme for sedimentation term in MyLake model
x = W_SEDIM(dpart,1_iintegers) ! the same laminar-limit sedimentation speed at all depths
conc(1) = gas%POCD(1) - dt*x*(gas%POCD(2) - gas%POCD(1))/(ls%h1*gsp%ddz(1))
do i = 2, gs%M
conc(i) = gas%POCD(i) - dt*x*(gas%POCD(i+1) - gas%POCD(i-1))/(ls%h1*gsp%ddz2(i-1))
!Sedimentation of POCD by MUSCL scheme
allocate(work(0:gs%M+1),fluxes(0:gs%M+1))
work(1:gs%M) = x
work(0) = 0.
work(gs%M+1) = 0.
call MUSCL_FLUXES(gs%M,gas%POCD,work,fluxes)
!print*, fluxes
!read*
do i = 1, gs%M+1
gas%POCD(i) = gas%POCD(i) - dt * (fluxes(i) - fluxes(i-1)) / (ls%h1*gsp%ddz05(i-1))
enddo
conc(gs%M+1) = gas%POCD(gs%M+1) - dt*x*(gas%POCD(gs%M+1) - gas%POCD(gs%M))/(ls%h1*gsp%ddz(gs%M))
gas%POCD(1:gs%M+1) = conc(1:gs%M+1)
if (maxval(gas%POCD) > 100.) then
print*, gas%POCD
print*, fluxes
read*
endif
deallocate(work,fluxes)
!!Settlement of POCD, explicit (Euler), central-differences scheme
!conc(1) = gas%POCD(1) - dt*x*(gas%POCD(2) - gas%POCD(1))/(ls%h1*gsp%ddz(1))
!do i = 2, gs%M
! conc(i) = gas%POCD(i) - dt*x*(gas%POCD(i+1) - gas%POCD(i-1))/(ls%h1*gsp%ddz2(i-1))
!enddo
!conc(gs%M+1) = gas%POCD(gs%M+1) - dt*x*(gas%POCD(gs%M+1) - gas%POCD(gs%M))/(ls%h1*gsp%ddz(gs%M))
!gas%POCD(1:gs%M+1) = conc(1:gs%M+1)
endif
! Horizontally integrated CO_2 bubble flux in case with single soil column
......
......@@ -43,6 +43,7 @@ call VERTICAL_ADVECTION(wst%v1 )
!endif
call VERTICAL_ADVECTION(gas%qwater(:,1))
call VERTICAL_ADVECTION(gas%DIC (:,1))
call VERTICAL_ADVECTION(gas%DOC (:,1))
!call VERTICAL_ADVECTION(gas%phosph)
!! Conservative filtering
......@@ -84,7 +85,7 @@ integer(kind=iintegers) :: i
!MUSCL scheme
allocate(work(1:M+1),fluxes(0:M+1))
forall(i=1:M+1) work(i) = var(i)*bathymw(i)%area_int
call MUSCL_FLUXES(work,fluxes)
call MUSCL_FLUXES(M,work,wst%w,fluxes)
do i = 1, M+1
var(i) = var(i) - dt / bathymw(i)%area_int * &
& (fluxes(i) - fluxes(i-1)) / (h1*gsp%ddz05(i-1))
......@@ -93,17 +94,58 @@ deallocate(work,fluxes)
END SUBROUTINE VERTICAL_ADVECTION
!!> Subroutine MUSCL_FLUXES computes flux following MUSCL scheme
!!! assuming equidistant grid
!SUBROUTINE MUSCL_FLUXES(var,flux)
!implicit none
!real(kind=ireals), intent(in ) :: var (1:M+1)
!real(kind=ireals), intent(out) :: flux(0:M+1)
!real(kind=ireals) :: r, varL, varR
!integer(kind=iintegers) :: i
!
!flux(0 ) = wst%w(0 )*var(1 )
!flux(M+1) = wst%w(M+1)*var(M+1)
!
!!Interpolation scheme below assumes the speed is positive
!do i = 2, M-1
! r = (var(i) - var(i-1) + small_value)/(var(i+1) - var(i) + small_value)
! varL = var(i) + 0.5*FLUX_LIMITER(r)*(var(i+1) - var(i))
! r = (var(i+1) - var(i) + small_value)/(var(i+2) - var(i+1) + small_value)
! varR = var(i+1) - 0.5*FLUX_LIMITER(r)*(var(i+2) - var(i+1))
! if (varR*varL < 0.) then
! flux(i) = wst%w(i)*0.5*(varR + varL) !sqrt(varR*varL)
! else
! flux(i) = wst%w(i)*sqrt(varR*varL)*sign(1._ireals,varR)
! endif
!enddo
!
!i = 1
!r = (var(i+1) - var(i) + small_value)/(var(i+2) - var(i+1) + small_value)
!varR = var(i+1) - 0.5*FLUX_LIMITER(r)*(var(i+2) - var(i+1))
!flux(i) = wst%w(i)*varR
!
!i = M
!r = (var(i) - var(i-1) + small_value)/(var(i+1) - var(i) + small_value)
!varL = var(i) + 0.5*FLUX_LIMITER(r)*(var(i+1) - var(i))
!flux(i) = wst%w(i)*varL
!
!END SUBROUTINE MUSCL_FLUXES
END SUBROUTINE VERTADV_UPDATE
!> Subroutine MUSCL_FLUXES computes flux following MUSCL scheme
!! assuming equidistant grid
SUBROUTINE MUSCL_FLUXES(var,flux)
SUBROUTINE MUSCL_FLUXES(M,var,w,flux)
implicit none
real(kind=ireals), intent(in ) :: var (1:M+1)
real(kind=ireals), intent(out) :: flux(0:M+1)
integer(kind=iintegers), intent(in ) :: M !> Number of layers
real(kind=ireals), intent(in ) :: var (1:M+1) !>Variable to be advected
real(kind=ireals), intent(in ) :: w (0:M+1) !>Vertical speed
real(kind=ireals), intent(out) :: flux(0:M+1) !>Variable flux due to advection
real(kind=ireals) :: r, varL, varR
integer(kind=iintegers) :: i
flux(0 ) = wst%w(0 )*var(1 )
flux(M+1) = wst%w(M+1)*var(M+1)
flux(0 ) = w(0 )*var(1 )
flux(M+1) = w(M+1)*var(M+1)
!Interpolation scheme below assumes the speed is positive
do i = 2, M-1
......@@ -112,25 +154,24 @@ do i = 2, M-1
r = (var(i+1) - var(i) + small_value)/(var(i+2) - var(i+1) + small_value)
varR = var(i+1) - 0.5*FLUX_LIMITER(r)*(var(i+2) - var(i+1))
if (varR*varL < 0.) then
flux(i) = wst%w(i)*0.5*(varR + varL) !sqrt(varR*varL)
flux(i) = w(i)*0.5*(varR + varL) !sqrt(varR*varL)
else
flux(i) = wst%w(i)*sqrt(varR*varL)*sign(1._ireals,varR)
flux(i) = w(i)*sqrt(varR*varL)*sign(1._ireals,varR)
endif
enddo
i = 1
r = (var(i+1) - var(i) + small_value)/(var(i+2) - var(i+1) + small_value)
varR = var(i+1) - 0.5*FLUX_LIMITER(r)*(var(i+2) - var(i+1))
flux(i) = wst%w(i)*varR
flux(i) = w(i)*varR
i = M
r = (var(i) - var(i-1) + small_value)/(var(i+1) - var(i) + small_value)
varL = var(i) + 0.5*FLUX_LIMITER(r)*(var(i+1) - var(i))
flux(i) = wst%w(i)*varL
flux(i) = w(i)*varL
END SUBROUTINE MUSCL_FLUXES
!> Flux limiter for TVD scheme
FUNCTION FLUX_LIMITER(r) result(fl)
implicit none
......@@ -143,7 +184,4 @@ fl = max(0.,min(1.,r))
END FUNCTION FLUX_LIMITER
END SUBROUTINE VERTADV_UPDATE
END MODULE VERTADV
# -*- coding: utf-8 -*-
class Circ(list):
def __getitem__(self, idx):
return super(Circ, self).__getitem__(idx % len(self))
import numpy as nmp
import matplotlib.pyplot as plt # for plotting
import matplotlib.ticker as mtick
......@@ -31,6 +35,9 @@ nd0 = int(input("Enter first day \n"))
ndlast = int(input("Enter last day \n"))
nh0 = int(input("Enter first hour \n"))
nhlast = int(input("Enter last hour \n"))
if nd0 == ndlast and nh0 == nhlast:
nhstep = 0
else:
nhstep = int(input("Enter time step, hours \n"))
#nh0 = 0
......@@ -67,8 +74,8 @@ time0_ = dts.date2num(time0_date)
fonts = 20 #font size
markersize = 15
linestyles = ['-','--','-.',':']
linecolors = ['b','g','r','c','m','y','k']
linestyles = Circ(['-','--','-.',':'])
linecolors = Circ(['b','g','r','c','m','y','k'])
ncolors = len(linecolors)
varname = []
......@@ -98,39 +105,39 @@ scale = []
##xmin = -8.E-2
##xmax = 8.E-2
varname.append(u'O$_2$, model')
ylabel = u'Depth, m' #$\circ C$
xlabel = varname[len(varname)-1] + u',~mg/l'
basefilename.append('Profiles')
ncol.append(15) #The column number to display, 2 for temperature
ncols.append(17) #Total number of columns
nheader.append(17) #Number of lines in the header
scale.append(1.)
#xmin = -8.E-2
#xmax = 8.E-2
#varname.append(u'CO$_2$, model')
#varname.append(u'O$_2$, model')
#ylabel = u'Depth, m' #$\circ C$
#xlabel = varname[len(varname)-1] + u',~mcmol/l'
#xlabel = varname[len(varname)-1] + u',~mg/l'
#basefilename.append('Profiles')
#ncol.append(16) #The column number to display, 2 for temperature
#ncol.append(15) #The column number to display, 2 for temperature
#ncols.append(17) #Total number of columns
#nheader.append(17) #Number of lines in the header
#scale.append(1.)
##xmin = -8.E-2
##xmax = 8.E-2
#varname.append(u'CH$_4$, model')
#varname.append(u'CO$_2$, model')
#ylabel = u'Depth, m' #$\circ C$
#xlabel = varname[len(varname)-1] + u',~mcmol/l'
#basefilename.append('Profiles')
#ncol.append(17) #The column number to display, 2 for temperature
#ncol.append(16) #The column number to display, 2 for temperature
#ncols.append(17) #Total number of columns
#nheader.append(17) #Number of lines in the header
#scale.append(1.)
##xmin = -8.E-2
##xmax = 8.E-2
varname.append(u'CH$_4$, model')
ylabel = u'Depth, m' #$\circ C$
xlabel = varname[len(varname)-1] + u',~mcmol/l'
basefilename.append('Profiles')
ncol.append(17) #The column number to display, 2 for temperature
ncols.append(17) #Total number of columns
nheader.append(17) #Number of lines in the header
scale.append(1.)
#xmin = -8.E-2
#xmax = 8.E-2
#varname.append(u'Компонента скорости по оси $x$')
#ylabel = u'Глубина, м' #$\circ C$
#xlabel = varname[len(varname)-1] + u',~м/с'
......@@ -301,9 +308,9 @@ if 'pathobs' in locals():
nprec = 3
nvars = 3
legobs = [u'T obs',u'O$_2$ obs',u'CH$_4$ obs']
nlocs = 5 #number of locations
legobssuf = [u'section~V',u'section~IV',u'section~III',u'section~II',u'section~I']
nvar = [2]
nlocs = 6 #number of locations
legobssuf = [u'section~V',u'section~IV',u'section~III',u'section~II',u'section~I',u'averaged']
nvar = [3]
nvarsd = len(nvar)
fobs = open(pathobs,'r')
#print(line)
......@@ -315,7 +322,7 @@ if 'pathobs' in locals():
line = filter(None,line)
#print line
nrecs = (len(line) - nprec)/(nvars + 1)
for i in range(nvarsd): legenditems.append(legobs[nvar[i]-1]+legobssuf[k])
for i in range(nvarsd): legenditems.append(legobs[nvar[i]-1]+', '+legobssuf[k])
varobs = nmp.zeros((nrecs,nvarsd))
z = nmp.zeros(nrecs)
for i in range(nrecs):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment