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

1.Morphometry added to salinity transport in water 2.Switch for salinity transport in soil 3.Misc

parent 7d7b863d
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 125 additions and 98 deletions
docs/users_guide.pdf
+ 0
− 0
View file @ 18f6eedf
No preview for this file type
docs/users_guide.tex
+ 2
− 0
View file @ 18f6eedf
......@@ -352,6 +352,8 @@ give significant speedup \\
\hline
\emph{sedim} & Gravitational sedimentation of tracer in water: 0 - off, 1 - on \\
\hline
\emph{salsoil} & Switch for salinity transport in soil: 0 - off (zero flux at the bottom), 1 - on \\
\hline
\emph{dyn\_pgrad} & Dynamic pressure gradient: 0 - off, 1 - on \\
\hline
\emph{zero\_model} & Switch for zero-dimensional model: 0 -off, 1 - on \\
......
source/model/driving_params_mod.f90
+ 4
− 1
View file @ 18f6eedf
......@@ -24,7 +24,7 @@ type, public :: intpar
character(len=30) :: name
end type
integer(kind=iintegers), parameter :: numder = 35 ! Number derived-type model integer controls handled in unified manner
integer(kind=iintegers), parameter :: numder = 36 ! Number derived-type model integer controls handled in unified manner
integer(kind=iintegers), parameter :: numder_real = 13 ! The same, but real controls
real(kind=ireals), parameter :: missing_value = -999.
......@@ -110,6 +110,7 @@ type(intpar), target :: skin
type(intpar), target :: massflux
type(intpar), target :: ifrad
type(intpar), target :: sedim
type(intpar), target :: salsoil
type(intpar), target :: nstep_keps ! The number of timesteps of k-epsilon parmeterization per on model timestep
type(intpar), target :: nscreen
type(intpar), target :: SoilType
......@@ -235,6 +236,7 @@ if (firstcall) then
dertypepar(33)%p => nsoilcols_
dertypepar(34)%p => cuette
dertypepar(35)%p => N_tribin
dertypepar(36)%p => salsoil
dt_out%name = 'dt_out' ;
......@@ -268,6 +270,7 @@ if (firstcall) then
massflux%name = 'massflux' ;
ifrad%name = 'ifrad' ;
sedim%name = 'sedim' ;
salsoil%name = 'salsoil' ;
nstep_keps%name = 'nstep_keps' ;
nscreen%name = 'nscreen' ;
SoilType%name = 'soiltype' ;
......
source/model/init_var_mod.f90
+ 19
− 19
View file @ 18f6eedf
......@@ -8,7 +8,7 @@
& init_T, skin, zero_model, &
& h10, l10, hs10, ls10, tempair, Ts0, Tb0, Tm, Tbb0, &
& Sals0, Salb0, us0, vs0, h_ML0, &
& rosoil, por, depth_area, ip, &
& rosoil, rosdry, por, depth_area, ip, &
& flag_snow, flag_snow_init, itop, nstep, &
& h1, l1, hs1, ls1, &
......@@ -27,23 +27,23 @@
! Subroutine INIT_VAR initializes the prognostic variables of the model
use PHYS_FUNC!, only : &
!& MELTPNT, &
!& UNFRWAT, &
!& WL_MAX, &
!& WI_MAX, &
!& HENRY_CONST, &
!& MELTINGPOINT
use PHYS_FUNC, only : &
& MELTPNT, &
& UNFRWAT, &
& WL_MAX, &
& WI_MAX, &
& HENRY_CONST, &
& MELTINGPOINT
use PHYS_CONSTANTS!, only : &
!& row0, g, Kelvin0, R_univ
use PHYS_CONSTANTS, only : &
& row0, g, Kelvin0, R_univ
use METH_OXYG_CONSTANTS!, only : &
!& ch4_atm0, o2_atm0, co2_atm0, &
!& rel_conc_ebul_crit, &
!& Henry_const0_ch4, &
!& Henry_temp_dep_ch4, &
!& Henry_temp_ref
use METH_OXYG_CONSTANTS, only : &
& ch4_atm0, o2_atm0, co2_atm0, &
& rel_conc_ebul_crit, &
& Henry_const0_ch4, &
& Henry_temp_dep_ch4, &
& Henry_temp_ref
use DRIVING_PARAMS, only : &
& tricemethhydr, &
......@@ -64,7 +64,7 @@
real(kind=ireals), intent(in) :: Sals0, Salb0
real(kind=ireals), intent(in) :: us0, vs0
real(kind=ireals), intent(in) :: h_ML0
real(kind=ireals), intent(in) :: rosoil(1:ns), por(1:ns)
real(kind=ireals), intent(in) :: rosoil(1:ns), rosdry(1:ns), por(1:ns)
real(kind=ireals), intent(in) :: depth_area(1:ndatamax,1:2) ! Data for lake bathymetry
type(initprof), intent(in) :: ip
real(kind=ireals), intent(in) :: zsoil(1:ns)
......@@ -151,7 +151,7 @@
endif
Sals1(1:ns,1:nsoilcols) = Salb0 ! assuming, that salinity in ground is the same
Sals1(1:ns,1:nsoilcols) = Salb0*row0/( rosdry(1)*(1 - por(1)) ) ! assuming, that salinity in ground is the same
! as one in near bottom layer of water
allocate(pressoil(1:ns))
......
source/model/lake.f90
+ 2
− 1
View file @ 18f6eedf
......@@ -592,7 +592,7 @@ if (init(ix,iy) == 0) then
& init_T, skin%par, zero_model%par, &
& h10, l10, hs10, ls10, tempair, Ts0, Tb0, Tm, Tbb0, &
& Sals0, Salb0, us0, vs0, h_ML0, &
& rosoil, por, depth_area, ip, &
& rosoil, rosdry, por, depth_area, ip, &
& flag_snow, flag_snow_init, itop, nstep, &
& h1, l1, hs1, ls1, &
......@@ -712,6 +712,7 @@ if1: IF (layer_case == 1) THEN
!---------------------------- CASE 1: LIQUID WATER LAYER AT THE TOP ("Summer")-----------------------
!Check the bottom layer water temperature
x = Meltpnt(Sal1(M+1),preswat(M+1),nmeltpoint%par)
if (WATER_FREEZE_MELT(Tw1(M+1), 0.5*ddz(M)*h1, x, +1) .and. &
......
source/model/salinity_mod.f90
+ 32
− 32
View file @ 18f6eedf
......@@ -7,30 +7,30 @@ use LAKE_DATATYPES, only : ireals, iintegers
contains
SUBROUTINE S_DIFF(dt,dhilow)
use ARRAYS!, only : &
!& deepice, water, &
!& water_salinity_indic, soil_salinity_indic, &
!& nstep
use ARRAYS_WATERSTATE!, only : lamsal, Sal1, Sal2
use ARRAYS_BATHYM!, only : h1, l1, dhw0, dhw
use ARRAYS_GRID!, only : ddz, dzs, nsoilcols
use ARRAYS_SOIL!, only : wsoil, Sals1, Sals2,rosdry,por
use ARRAYS, only : &
& deepice, water, &
& water_salinity_indic, soil_salinity_indic, &
& nstep
use ARRAYS_WATERSTATE, only : lamsal, Sal1, Sal2
use ARRAYS_BATHYM, only : h1, l1, dhw0, dhw, area_int, area_half
use ARRAYS_GRID, only : ddz, nsoilcols
use ARRAYS_SOIL, only : wsoil, Sals1, Sals2,rosdry,por,dzs
use DRIVING_PARAMS, only : sedim, soilswitch, M, ns
use DRIVING_PARAMS, only : sedim, soilswitch, M, ns, salsoil
use ATMOS!, only: &
!& Sflux0
use ATMOS, only: &
& Sflux0
use PHYS_FUNC!, only: &
!& w_sedim
use PHYS_FUNC, only: &
& w_sedim
use PHYS_CONSTANTS!, only : &
!& roi_d_row0, &
!& salice0, &
!& row0
use PHYS_CONSTANTS, only : &
& roi_d_row0, &
& salice0, &
& row0
use T_SOLVER_MOD!, only : &
!& DIFF_COEF
use T_SOLVER_MOD, only : &
& DIFF_COEF
implicit none
......@@ -77,21 +77,21 @@ endif
if (water == 1) then
call DIFF_COEF(a,b,c,d,2,M,2,M,water_salinity_indic,dt)
x = 0.5*( - lamsal(1)/(h1*ddz(1)) + dhw0/(2.d0*dt) )
x = 0.5*( - area_half(1)*lamsal(1)/(h1*ddz(1)*area_int(1)) + dhw0/(2.d0*dt) )
c(1) = x - ddz(1)*h1/(2.d0*dt) ! - dhw0/(2*dt) is wrong sign!
b(1) = x
d(1) = - ddz(1)*h1*Sal1(1)/(2.d0*dt) - Sflux0 + x*(Sal1(2) - Sal1(1))
if (deepice == 1) then
! Case water,deepice and soil; upper ice and snow are allowed
! Salinity diffusion in water
x = - lamsal(M)/(ddz(M)*h1) + (dhw-dhw0)/(2.d0*dt)
x = - area_half(M)*lamsal(M)/(ddz(M)*h1*area_int(M+1)) + (dhw-dhw0)/(2.d0*dt)
c(M+1) = x - ddz(M)*h1/(2.d0*dt)
a(M+1) = x
d(M+1) = - Sal1(M+1)*ddz(M)*h1/(2.d0*dt) + Sflux1
call PROGONKA (a,b,c,d,Sal,1,M+1)
Sal (1:M+1) = max(Sal(1:M+1),0._ireals)
Sal2(1:M+1) = Sal(1:M+1)
if (soilswitch%par == 1) then
if (soilswitch%par == 1 .and. salsoil%par == 1) then
! Salinity diffusion in soil
call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
c(1) = - 1.d0 - dt*wsoil(1)/dzs(1)
......@@ -109,13 +109,13 @@ if (water==1) then
endif
else
! Case water, soil; upper ice and snow are allowed
if (soilswitch%par == 1) then
if (soilswitch%par == 1 .and. salsoil%par == 1) then
!------WATER-SOIL INTERFACE-------------------------
x = 0.5*(dhw-dhw0)/dt - lamsal(M)/(ddz(M)*h1)
x = 0.5*(dhw-dhw0)/dt - area_half(M)*lamsal(M)/(ddz(M)*h1*area_int(M+1))
y = rosdry(1)*(1 - por(1))/row0
a(M+1) = x
b(M+1) = 0.5*wsoil(1)*y
c(M+1) = - ( 0.5*(dzs(1)*y + ddz(M)*h1)/dt + 0.5*wsoil(1)*y - x ) ! sign of (dhw-dhw0) corrected
c(M+1) = - ( 0.5*(dzs(1)*y + ddz(M)*h1)/dt + 0.5*wsoil(1)*y - x )
d(M+1) = - Sal1(M+1)*0.5*(dzs(1)*y + ddz(M)*h1)/dt
call DIFF_COEF(a,b,c,d,2,ns-1,M+2,M+ns-1,soil_salinity_indic,dt)
c(M+ns) = - 1.d0 + wsoil(ns-1)*dt/(dzs(ns-1))
......@@ -135,7 +135,7 @@ if (water==1) then
Sals2(1:ns,nsoilcols) = Sal(M+1:M+ns)
Sal2(1:M+1) = Sal(1:M+1)
else
x = - lamsal(M)/(ddz(M)*h1) + (dhw-dhw0)/(2.d0*dt)
x = - area_half(M)*lamsal(M)/(ddz(M)*h1*area_int(M+1)) + (dhw-dhw0)/(2.d0*dt)
c(M+1) = x - ddz(M)*h1/(2.d0*dt)
a(M+1) = x
d(M+1) = - Sal1(M+1)*ddz(M)*h1/(2.d0*dt) + Sflux1
......@@ -147,7 +147,7 @@ if (water==1) then
if (sedim%par == 1) call SAL_SEDIM(ddz,h1,dt,Sal2)
else
! Case of the soil and the ice above
if (soilswitch%par == 1) then
if (soilswitch%par == 1 .and. salsoil%par == 1) then
call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
c(1) = - 1.d0 - dt*wsoil(1)/dzs(1)
b(1) = dt*wsoil(1)/dzs(1)
......
source/model/soil_mod.f90
+ 18
− 16
View file @ 18f6eedf
......@@ -200,6 +200,8 @@
!end do
end if
rosdry(:) = 1.2E+3
if (firstcall) firstcall=.false.
RETURN
END SUBROUTINE COMSOILFORLAKE
......@@ -689,21 +691,21 @@ END SUBROUTINE SOILFORLAKE
SUBROUTINE SOIL_COND_HEAT_COEF(is)
use ARRAYS_SOIL!, only : rosdry,csoil,rosoil,lamsoil,
!lammoist,filtr,wsoil,por,bh,wl1,wi1,psimax,flwmax,dlmax
use ARRAYS_GRID!, only : dzs, nsoilcols
use ARRAYS_SOIL, only : rosdry,csoil,rosoil,lamsoil, &
lammoist,filtr,wsoil,por,bh,wl1,wi1,psimax,flwmax,dlmax,dzs
use ARRAYS_GRID, only : nsoilcols
use DRIVING_PARAMS, only : ns, tricemethhydr
use PHYS_CONSTANTS!, only : &
!& row0, &
!& cw, &
!& ci, &
!& roi, &
!& row0_d_roi
use PHYS_FUNC!, only : &
!& WL_MAX, &
!& SOIL_COND_JOHANSEN
use METH_OXYG_CONSTANTS!, only : &
!& cpmethhydr
use PHYS_CONSTANTS, only : &
& row0, &
& cw, &
& ci, &
& roi, &
& row0_d_roi
use PHYS_FUNC, only : &
& WL_MAX, &
& SOIL_COND_JOHANSEN
use METH_OXYG_CONSTANTS, only : &
& cpmethhydr
implicit none
......@@ -733,7 +735,7 @@ if (firstcall) then
endif
endif
rosdry(:) = 1200. !2400.
!rosdry(:) = 1200. !2400.
do i = 1, ns
wlmax_i = POR(i)*ROW0/(rosdry(i)*(1-por(i)))
......
source/model/t_solver_mod.f90
+ 4
− 4
View file @ 18f6eedf
......@@ -755,10 +755,10 @@ SELECT CASE (substr)
case (water_salinity_indic)
do i = m0, m1
j = i - m0 + n0
a(i) = 0.5d0*( - lamsal(j-1)/(ddz(j-1)*h1) + dhw*dzeta_int(j)*dt_inv05 - &
& dhw0*dt_inv05)
b(i) = 0.5d0*( - lamsal(j)/(ddz(j)*h1) - dhw*dzeta_int(j)*dt_inv05 + &
& dhw0*dt_inv05)
a(i) = 0.5d0*( - area_half(j-1)*lamsal(j-1)/(ddz(j-1)*h1*area_int(j)) + &
& dhw*dzeta_int(j)*dt_inv05 - dhw0*dt_inv05)
b(i) = 0.5d0*( - area_half(j) *lamsal(j) /(ddz(j) *h1*area_int(j)) - &
& dhw*dzeta_int(j)*dt_inv05 + dhw0*dt_inv05)
c(i) = a(i) + b(i) - dt_inv*ddz05(j-1)*h1
d(i) = - Sal1(j)*dt_inv*ddz05(j-1)*h1 + &
& a(i)*Sal1(j-1) + b(i)*Sal1(j+1) - &
......
tools/plot_contour.py
+ 44
− 25
View file @ 18f6eedf
......@@ -7,18 +7,19 @@ from plot_timser import *
nplot = input("Choose a variable to plot \n \
1 - temperature \n \
2 - CH_4 \n \
3 - O_2 \n \
4 - CO_2 \n \
5 - CH_4 bubble flux \n \
6 - Richardson number \n \
7 - Heat conductance \n \
8 - temperature measured \n \
9 - oxygen measured \n \
10 - methane measured \n \
11 - carbon dioxide measured \n ")
2 - salinity \n \
3 - CH_4 \n \
4 - O_2 \n \
5 - CO_2 \n \
6 - CH_4 bubble flux \n \
7 - Richardson number \n \
8 - Heat conductance \n \
9 - temperature measured \n \
10 - oxygen measured \n \
11 - methane measured \n \
12 - carbon dioxide measured \n ")
path = '../results/BolshoiVilui2015_base/time_series/' # path to files
path = '../results/BolshoiVilui2015_ext05/time_series/' # path to files
pathobs = os.path.expanduser('~/Files/main/observdata/Kuivajarvi/')
exp = 'base'
......@@ -44,7 +45,25 @@ if int(nplot) == 1: # Temperature
levmax = 20.
nlevplot = 15 # Number of countour levels
pathsave = path
elif int(nplot) == 2: # Methane
if int(nplot) == 2: # Salinity
basename = 'sal_water 1 1'
fileout = 'sal_wat_cont'
nheader = 6
labelx = 'Time, months'
labely = 'Depth, m'
timescale = 365/12 # 365/12*24 - number of hours in a month
ctime = 5
cstart = 7
time0 = 7.48 #Initial time of data
#t0disp = 7.50 # Initial time for display
#t1disp = 7.95 # Final time for display
#nvardisp = nlevs_water
title = "Salinity, $kg/kg$ "
#levmin = 0.
#levmax = 20.
nlevplot = 15 # Number of countour levels
pathsave = path
elif int(nplot) == 3: # Methane
basename = 'methane_water 1 1'
fileout = 'meth_wat_cont'
nheader = 6 #Number of lines in the header
......@@ -63,7 +82,7 @@ elif int(nplot) == 2: # Methane
levmax = 600.
nlevplot = 20 # Number of countour levels
pathsave = path
elif int(nplot) == 3: # Oxygen
elif int(nplot) == 4: # Oxygen
basename = 'oxygen_water 1 1'
fileout = 'oxyg_wat_cont'
nheader = 6 #Number of lines in the header
......@@ -81,7 +100,7 @@ elif int(nplot) == 3: # Oxygen
levmax = 10.
nlevplot = 20 # Number of countour levels
pathsave = path
elif int(nplot) == 4: # Carbon dioxide
elif int(nplot) == 5: # Carbon dioxide
basename = 'co2_water 1 1'
fileout = 'co2_wat_cont'
nheader = 6 #Number of lines in the header
......@@ -99,7 +118,7 @@ elif int(nplot) == 4: # Carbon dioxide
levmax = 20.
nlevplot = 20 # Number of countour levels
pathsave = path
elif int(nplot) == 5: # Methane bubble flux
elif int(nplot) == 6: # Methane bubble flux
basename = 'fbblflx_ch4 1 1'
fileout = 'methbubbleflx_cont'
nheader = 6 #Number of lines in the header
......@@ -114,7 +133,7 @@ elif int(nplot) == 5: # Methane bubble flux
#nvardisp = nlevs_water
title = "Methane bubble flux, $mol/(m^2*s)$"
pathsave = path
elif int(nplot) == 6: # Richardson number
elif int(nplot) == 7: # Richardson number
basename = 'Ri 1 1'
fileout = 'Ri_cont'
nheader = 6 #Number of lines in the header
......@@ -124,13 +143,13 @@ elif int(nplot) == 6: # Richardson number
ctime = 5
cstart = 7
time0 = 4 #Initial time
t0disp = 4 # Initial time for display
t1disp = 10 # Final time for display
#t0disp = 4 # Initial time for display
#t1disp = 10 # Final time for display
levmax = 1.
#nvardisp = nlevs_water
title = "Richardson number, $n/d$"
pathsave = path
elif int(nplot) == 7: # Heat conductance
elif int(nplot) == 8: # Heat conductance
basename = 'lamw 1 1'
fileout = 'lamw_cont'
nheader = 6 #Number of lines in the header
......@@ -140,12 +159,12 @@ elif int(nplot) == 7: # Heat conductance
ctime = 5
cstart = 7
time0 = 4 #Initial time
t0disp = 4 # Initial time for display
t1disp = 10 # Final time for display
#t0disp = 4 # Initial time for display
#t1disp = 10 # Final time for display
#nvardisp = nlevs_water
title = "Heat conductance, $m^2/s$"
pathsave = path
elif int(nplot) == 8: # Temperature measured
elif int(nplot) == 9: # Temperature measured
basename = 'obs'
fileout = 'tempobs_cont'
nheader = 0 #Number of lines in the header
......@@ -164,7 +183,7 @@ elif int(nplot) == 8: # Temperature measured
nlevplot = 15 # Number of countour levels
pathsave = path
path = pathobs
elif int(nplot) == 9: # Oxygen measured
elif int(nplot) == 10: # Oxygen measured
basename = 'o2man2013cont'
fileout = 'o2obsman_cont'
nheader = 1 #Number of lines in the header
......@@ -183,7 +202,7 @@ elif int(nplot) == 9: # Oxygen measured
nlevplot = 20 # Number of countour levels
pathsave = path
path = pathobs
elif int(nplot) == 10: # Methane measured
elif int(nplot) == 11: # Methane measured
basename = 'ch4man2013cont'
fileout = 'ch4obsman_cont'
nheader = 1 #Number of lines in the header
......@@ -202,7 +221,7 @@ elif int(nplot) == 10: # Methane measured
nlevplot = 20 # Number of countour levels
pathsave = path
path = pathobs
elif int(nplot) == 11: # Carbon dioxide measured
elif int(nplot) == 12: # Carbon dioxide measured
basename = 'co2man2013cont'
fileout = 'co2obsman_cont'
nheader = 1 #Number of lines in the header
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment