Skip to content
Snippets Groups Projects
Commit fab2c539 authored by Georgiy Faikin's avatar Georgiy Faikin
Browse files

Add new environmenmt data type: environment_data_station; ROTHC approved

parent 5d54d7c7
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 233 additions and 251 deletions
data/obs_data_DAO4.csv 0 → 100644
+ 9
− 0
View file @ fab2c539
date,Control,NPK1.5+FYM,NPK3,NPK3+FYM
10-05-1937,3.0997,3.2387,3.2248,3.0441
10-05-1949,2.65026,3.07653,2.5854,2.75
10-05-1952,2.53906,2.919,2.5298,2.43713
10-05-1958,2.53906,3.06726,2.63173,2.5298
10-05-1966,2.10353,2.8726,2.2796,2.34446
10-05-1972,2.09426,2.8356,2.2796,2.21473
10-05-1977,2.0294,2.7336,2.24253,2.15913
10-05-2009,2.74108,2.910382,2.90232,2.9399426
nonlin.f90
+ 106
− 49
View file @ fab2c539
program nonlin
implicit none
integer, parameter :: N = 11 ! количество элементов сетки
! ----------------------------------------------- Variables -----------------------------------------------------------
! общие переменные
integer :: i,j !< count
integer, parameter :: N = 5 ! количество элементов сетки 1
integer, parameter :: M = 5 ! количество элементов сетки 2
logical :: odn = .false. !false
! Величины для шага по сетке
integer, parameter :: a = 12 ! коэффициент сгущения сетки
real :: x(N) ! Значение аргумента, аналог x(i-1/2) в статье
real :: E(N) ! интерполяция
real :: h(N) ! Шаг сетки
real :: Cor(2:N-1) ! Оригинальное значение поступающей функции
real :: z(N) ! Уровень сетки 1 (глубина)
real :: x(N) ! Значение аргумента, промежуточное значение уровня сетки, аналог x(i-1/2) в статье
real :: h(N) ! Шаг сетки 1
real :: t(M) ! Уровень сетки 2 (время)
real :: dt = 1. ! Шаг сетки 2
! Внешняя функция, и её усреднение
real :: Cor(N)!Cor(2:N-1) ! Оригинальное значение поступающей функции
real :: c0 = 3. ! начальное значение функции
real :: w = 0.05 ! const функции
real :: Cr1(N) ! Среднее значение функции (функция разбивается на небольшие интервалы)
real :: InCor ! Интеграл от поступающей из вне функции
real :: M10(2:N-1), M10p1(2:N-1), M10p2(2:N-1)
real :: M11(2:N-1), M11p1(2:N-1), M11p2(2:N-1)
real :: M12(2:N-1), M12p1(2:N-1), M12p2(2:N-1)
real :: Ct(M,N) ! Функция от времени
! Коэффициенты для финального вычисления, вычисления реконструкции слева и справа
real :: M10(2:N-1), M10p1(2:N-1), M10p2(2:N-1) ! Коэффициент С10 и его составляющие
real :: M11(2:N-1), M11p1(2:N-1), M11p2(2:N-1) ! Коэффициент С11 и его составляющие
real :: M12(2:N-1), M12p1(2:N-1), M12p2(2:N-1) ! Коэффициент С12 и его составляющие
real :: M00p1(N-1), M00p2(N-1), M00(N-1)
real :: M01p1(N-1), M01p2(N-1), M01(N-1)
real :: M02p1(N-1), M02p2(N-1), M02(N-1)
real :: M00p1(N-2), M00p2(N-2), M00(N-2) ! Коэффициент С00 и его составляющие
real :: M01p1(N-2), M01p2(N-2), M01(N-2) ! Коэффициент С01 и его составляющие
real :: M02p1(N-2), M02p2(N-2), M02(N-2) ! Коэффициент С02 и его составляющие
real :: CL(2:N-1), CR(1:N-2) ! Реконструкция слева и справа
! Финальное вычисление, реконструкция слева и справа
real :: CzL(2:N-1), CzR(1:N-2) ! Реконструкция слева и справа
real :: InCor
integer :: i !< count
! ----------------------------------------------- Main programm -------------------------------------------------------
! ---- initial value ----
do i = 2, N-1
Cor(i) = 0.003 + 0.001*i
t = 1.
Ct = 0.
! ---- Вычисление аргумента сетки, шага по сетке ----
if (odn .eqv. .true.) then ! ---- Случай однородной сетки ----
z(1) = 1.
h(:) = 1.
do i = 2, N
z(i) = z(i-1) + h(i)
end do
! ---- вычисление аргумента сетки, шага по сетке ----
E = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1./)
else ! ---- Неоднородная сетка ----
E = (/0.1,0.2,0.3,0.4,0.5/)!,0.6,0.7,0.8,0.9,1.,1.1/)
do i = 1, N
!E(i) = (i - 1.)/(N - 1.)*1.0
x(i) = ((1 + a)**E(i) - 1.)/a
z(i) = ((1 + a)**E(i) - 1.)/a
end do
h(1) = z(1)
do i = 2, N
h(i) = z(i) - z(i-1)
end do
h(1) = x(1)
end if
x(1) = z(1)/2.
do i = 2, N
h(i) = x(i) - x(i-1)
x(i) = (z(i) + z(i-1))/2
Cor(i) = C0*(x(i) - w)
end do
! ---- Вычисление среднего значения поступающей на вход функции Cor(x,t) в ячейке с номером i, ----
!! На выход Cr1(N-1)
do i = 2, N-1
InCor = InCor + Cor(i)*0.5*(x(i+1) - x(i-1)) !Функция Cor(i) интегрируется в пределах x(i-1) и x(i)
end do
h(1) = 0.01
do i = 1, N
Cr1(i) = 1./h(i)*InCor
end do
......@@ -69,29 +105,50 @@ program nonlin
M00(i) = M00p1(i)/M00p2(i)
M01p1(i) = h(i)*(h(i+1) + h(i+2))/((h(i) + h(i+1))*h(i+2))
M01p2(i) = h(i)*(h(i+1))**2/(h(i+2)*(h(i) + h(i+1) + h(i+2))*(h(i+1) + h(i+2)))
M01(i) = M01p1(i) + M01p2(i)
M01p2(i) = h(i)*(h(i+1)**2)/(h(i+2)*(h(i) + h(i+1) + h(i+2))*(h(i+1) + h(i+2)))
M01(i) = M01p1(i) - M01p2(i)
M02p1(i) = h(i)*h(i+1)
M02p1(i) = - h(i)*h(i+1)
M02p2(i) = (h(i) + h(i+1) + h(i+2))*(h(i+1) + h(i+2))
M02(i) = M02p1(i)/M02p2(i)
end do
! ---- Реконструкция ----
do i = 2, N-1
CL(i) = M10(i)*Cr1(i-1) + M11(i)*Cr1(i) + M12(i)*Cr1(i+1) ! Слева
CzL(i) = M10(i)*Cr1(i-1) + M11(i)*Cr1(i) + M12(i)*Cr1(i+1) ! Слева
end do
do i = 1, N-2
CR(i) = M00(i)*Cr1(i) + M01(i)*Cr1(i+1) + M02(i)*Cr1(i+2) ! Справа
CzR(i) = M00(i)*Cr1(i) + M01(i)*Cr1(i+1) + M02(i)*Cr1(i+2) ! Справа
end do
print*, 'CL', CL
print*, 'CR', CR
print*, 'Cr1:', Cr1
! ---- учет диффиренцирования по времени ----
do j = 2, M
Ct(j,2) = (CzL(2))/h(i)*dt + Ct(j-1,2)
do i = 3, N-1
Ct(j,i) = (CzL(i) - CzL(i-1))/h(i)*dt + Ct(j-1,i)
end do
end do
!print*, 'M10:', M10
!print*, 'M11:', M11
!print*, 'M12:', M12
print*, 'x:', x
print*, 'h:', h
print*, 'InCor:', InCor
!print*, 'M01:', M01
!print*, 'M02:', M02
print*,'CzL', CzL
print*,'Ct(M,2)', Ct(M,2)
print*,'Ct(M,3)', Ct(M,3)
print*,'Ct(M,4)', Ct(M,4)
print*,'Ct(M,5)', Ct(M,5)
!print*,'Ct(M,1)', Ct(M,4)
!print*, 'InCor:', InCor
!print*, 'Cr1:', Cr1
end program nonlin
......
plot.py
+ 14
− 6
View file @ fab2c539
......@@ -21,15 +21,23 @@ Usage:
nml = f90nml.read('ui1_config.nml')
carbon_model_type = nml['config_namelist']['carbon_model_type']
station = nml['config_namelist']['station_name']
station = 'Rostov' # Rostov DAO3 DAO4 VLDMR TRGK
model_data_paths = [
f'results/{carbon_model_type}/{station}_1.txt',
f'results/{carbon_model_type}/{station}_2.txt',
f'results/{carbon_model_type}/{station}_3.txt'
]
if station == 'DAO3' or 'DAO4':
model_data_paths = [
f'results/{station}_1.txt',
f'results/{station}_2.txt',
f'results/{station}_3.txt'
#f'results/{station}_4.txt'
f'results/{carbon_model_type}/{station}_1.txt',
f'results/{carbon_model_type}/{station}_2.txt',
f'results/{carbon_model_type}/{station}_3.txt',
f'results/{carbon_model_type}/{station}_4.txt'
]
#if enviromental == station
obs_data_path = f'data/obs_data_{station}.csv'
#else
......
source/carbon/carbon_model_to_core_arg_kit.f90
+ 2
− 2
View file @ fab2c539
......@@ -22,8 +22,8 @@ module carbon_model_to_core_arg_kit
integer, parameter :: nmonth = 12
integer, parameter :: year_min = 1974 !< предусмотренный диапазон лет
integer, parameter :: year_max = 2018 !<
integer, parameter :: year_min = 1937 !< предусмотренный диапазон лет
integer, parameter :: year_max = 2937 !<
! (44)*12 = 528 ! Для случая Rostov: 1974 - 2018
! (73)*12 = 876 ! Для случая DAO3: 1939 - 2012
......
source/carbon/carbon_models/carbon_model_INMCM.f90
+ 4
− 4
View file @ fab2c539
......@@ -34,8 +34,8 @@ contains
call set_pool(pid = n_Catm, name = 'catm' )
call set_pool(pid = n_Cveg, name = 'cveg' , alias = Cveg )
call set_pool(pid = n_Csoil, name = 'csoil' ,initial_value = 8.730624, alias = Csoil ) !8.730624 2.78124 2.806524 2.983512 2.941372
call set_pool(pid = n_Csoilb, name = 'csoilb' ,initial_value = 0.178176, alias = Csoilb ) !0.178176 0.05676 0.057276 0.060888 0.060028
call set_pool(pid = n_Csoil, name = 'csoil' ,initial_value = 2.941372, alias = Csoil ) !8.730624 2.78124 2.806524 2.983512 2.941372
call set_pool(pid = n_Csoilb, name = 'csoilb' ,initial_value = 0.060028, alias = Csoilb ) !0.178176 0.05676 0.057276 0.060888 0.060028
call set_flux(fid = n_F, pid_out = n_Csoil, pid_in = n_Catm, name = 'fmicr', alias = Fmicr)
call set_mult(n_F, 'lin', x_ijn = Csoil)
......@@ -44,7 +44,7 @@ contains
call set_mult(n_F, 'step', x_ij = Tsoil(:,:,ms+1), x0 = tmin_soil, k = 0.5, y0 = 0.5)
call set_mult(n_F, 'mm', x_n = rsw, k = 1., x0 = -0.20)
call set_mult(n_F, 'hyp', x_n = rsw, k = 0.23, x0 = -0.23)
call set_mult(n_F, 'const', c = 1.5) !1.5
call set_mult(n_F, 'const', c = 1.5)
call set_flux(fid = n_F, pid_out = n_Csoilb, pid_in = n_Catm, name = 'fmicrb', alias = Fmicrb)
call set_mult(n_F, 'lin', x_ijn = Csoilb)
......@@ -53,7 +53,7 @@ contains
call set_mult(n_F, 'step', x_ij = Tsoil(:,:,ms+1), x0 = tmin_soil, k = 0.5, y0 = 0.5)
call set_mult(n_F, 'mm', x_n = rsw, k = 1., x0 = -0.20)
call set_mult(n_F, 'hyp', x_n = rsw, k = 0.23, x0 = -0.23)
call set_mult(n_F, 'const', c = 1.5) !1.5
call set_mult(n_F, 'const', c = 1.5)
call set_mult(n_F, 'const', c = cv81b)
call set_flux(fid = n_F, pid_out = n_Cveg, pid_in = n_Csoil, name = 'dv68', alias = Flit)
......
source/carbon/carbon_models/carbon_model_INMCM_aux.f90
+ 9
− 8
View file @ fab2c539
......@@ -4,8 +4,8 @@ module carbon_model_inmcm_aux
! ----------------------------------------------- Use pack ------------------------------------------------------------
use environment_model_inmcm, only: nv2
use environment_core, only: Tsoil, Temp, Wsoil, lambd
use grid, only : date_c, date_fst, date_lst, dt
use config, only : station_name, station_opt
use grid, only: date_c
use config, only: station_name, station_opt, carbon_model_type
use const, only: pi, yrs, nmonth
! ---------------------------------------------- Variables ------------------------------------------------------------
implicit none
......@@ -169,8 +169,8 @@ contains
integer :: k
real :: work
!Wsoil(ii,jj,:) = 1.
!Tsoil(ii,jj,:) = 45.
!Wsoil(ii,jj,:) = 1. !0.21447610589527216609628319344309
!Tsoil(ii,jj,:) = 40.
work = 0.
do k = ms+1, ml-1
if (z(k) <= hint) then
......@@ -187,7 +187,8 @@ contains
use grid, only : date_c
open(unit=500, file='results/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
open(unit=500, file='results/'//trim(carbon_model_type)//'/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
!open(unit=500, file='results/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
write(500,*) date_c%timestamp,';',Csoil(:,:,12),';',Csoilb(:,:,12)
end subroutine
......
source/carbon/carbon_models/carbon_model_ROTHC.f90
+ 8
− 8
View file @ fab2c539
......@@ -20,17 +20,17 @@ module carbon_model_rothc
call set_pool(pid = n_Catm, name = 'catm')
call set_pool(pid = n_Cveg, name = 'cveg')
call set_pool(pid = n_CDPM, name = 'CDPM', initial_value = 0.043915506788352, alias = CDPM)
call set_pool(pid = n_CRPM, name = 'CRPM', initial_value = 0.57587725621248 , alias = CRPM)
call set_pool(pid = n_CBIO, name = 'CBIO', initial_value = 0.20176158117888 , alias = CBIO)
call set_pool(pid = n_CHUM, name = 'CHUM', initial_value = 8.11102631424 , alias = CHUM)
call set_pool(pid = n_CDPM, name = 'CDPM', initial_value = 0.012180579701505, alias = CDPM)
call set_pool(pid = n_CRPM, name = 'CRPM', initial_value = 0.3894608389488 , alias = CRPM)
call set_pool(pid = n_CBIO, name = 'CBIO', initial_value = 0.06696989650323 , alias = CBIO)
call set_pool(pid = n_CHUM, name = 'CHUM', initial_value = 2.5754859888393 , alias = CHUM)
! litterfall
call set_flux(fid = n_F, pid_out = n_Cveg, pid_in = n_CDPM, name = 'fdpm*lambdac')
call set_mult(n_F, 'lin', x = lambd, k = fdpm)
call set_mult(n_F, 'lin', x_ij = lambd, k = fdpm)
call set_flux(fid = n_F, pid_out = n_Cveg, pid_in = n_CRPM, name = '(1-fdpm)*lambdac')
call set_mult(n_F, 'lin', x = lambd, k = (1 - fdpm))
call set_mult(n_F, 'lin', x_ij = lambd, k = (1 - fdpm))
! to atmosphere
call set_flux(fid = n_F, pid_out = n_CDPM, pid_in = n_Catm, name = 'RDPM')
......@@ -47,10 +47,10 @@ module carbon_model_rothc
! atmosphere fall
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = '(0.46*bettar)*Rs')
call set_mult(n_F, 'lin', x = Rs, k = (0.46*bettar))
call set_mult(n_F, 'lin', x = Rs, k_ij = (0.46*bettar))
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = '(0.54*bettar)*Rs')
call set_mult(n_F, 'lin', x = Rs, k = (0.54*bettar))
call set_mult(n_F, 'lin', x = Rs, k_ij = (0.54*bettar))
end subroutine carbon_model_assembly
......
source/carbon/carbon_models/carbon_model_ROTHC_aux.f90
+ 48
− 141
View file @ fab2c539
module carbon_model_rothc_aux
use const, only : yrs, day2sec, pi
use carbon_model_to_core_arg_kit, only : nmonth
use grid, only: date_c, date_fst, date_lst, dt
use config, only : station_name, station_opt
!use environment_core, only: !veg !Wsoil
use environment_core, only: Tsoil, Wsoil, veg, lambd, bettar, sw
use const, only: yrs, day2sec, pi, nmonth
use grid, only: date_c, i0, i1, j0, j1, ml
use config, only: station_name, station_opt, carbon_model_type
implicit none
! Интерфейс
! -------------------------------------------------------------------------------------------------------------------
integer :: mnc !< Колличество месяцев в расчете
integer :: mncX = 0 !< Номер месяца с начала работы программы !Как часть параметров для шага по времени
integer :: mnclot_fst(5) = (/1937,1939,1956,1968,1974/) !Как часть параметров для задания климатических данных
integer :: mnclot_lst(5) = (/2012,2012,2018,2018,2018/) !Как часть параметров для задания климатических данных
!DAO4 DAO3 VLDR TRGK ROST
integer :: station_n ! 1 2 3 4 5 !Как часть параметров для задания климатических данных
! -------------------------------------------------------------------------------------------------------------------
! Main part of RoTHC
! ------- Pools -------
real, dimension(:,:,:), pointer :: Catm !< атмосфера
real, dimension(:,:,:), pointer :: Cveg !< растительность
......@@ -28,7 +17,7 @@ implicit none
real, dimension(:,:,:), pointer :: CHUM !< долгоживущий гумусовый пул / long lived humified pool
! ------- Flows n_Catm, n_Cveg, n_CDPM, n_CRPM, n_CBIO, n_CHUM -------
real, dimension(:,:,:), pointer :: litterfall1 !< Между растениями и почвой DPM
real, dimension(:,:,:), pointer :: litterfall2 !< Между растением и почвой RPM
real, dimension(:,:,:), pointer :: litterfall2 !< Между растениеми и почвой RPM
real, dimension(:,:,:), pointer :: atmosphere1 !< Между атмосферой и почвой BIO
real, dimension(:,:,:), pointer :: atmosphere2 !< Между атмосферой и почвой HUM
real, dimension(:,:,:), pointer :: respiration1 !< Между почвой DPM и атмосферой ! 1 - DPM
......@@ -51,95 +40,30 @@ implicit none
!< Скорость дыхания единицы массы каждого пула в стандартных условиях, [1-s]
real, parameter :: alphadr = 1.44 !< Определяет соотношение поступления опада между DPM и RPM, [dim],
! Для деревьев, для кустарников, для натуральной травы, для СХ (/ 0.25, 0.33, 0.67, 1.44 /)
real :: bettar !< Fraction of soil respiration, [%] from 0 to 1
! ------ Climate variables --------
! ------ Determined within
real :: s0 !< Оптимальная влажность почвы, [dim]
real :: smin !< Переменные для функций, [dim]
real, allocatable :: s0(:,:,:) !< Оптимальная влажность почвы, [dim]
real, allocatable :: smin(:,:,:) !< Переменные для функций, [dim]
! ------ Determined externally
real, allocatable :: in_temp(:) !< Поступление извне данных cредней температуры воздуха в месяц, [Celsius]
real :: Temp !< Средняя температура воздуха в интервал времени dt, [Celsius]
real :: Tsoil !< Средняя температура почвы в интервал времени dt, [Celsius]
real, allocatable :: in_veg(:) !< Поступление извне данных средней концентрации растений в месяц, [%]
real :: veg !< Средняя концентрация растений в интервал времени dt, [%]
real, allocatable :: in_wsoil(:) !< Поступление извне данных влажности почвы, [dim]
real :: Wsoil !< Влажность почвы в интервал времени dt, [dim]
real :: sw !< Влажность увядания DAO: 0.120 Rostov: 0.146
! ------ Litterfall --------
real, allocatable :: in_lambd(:) !< Поступление извне данных по поступлению углерода в почву, [kg/m**3 / year]
real :: lambd !< Поступление углерода в почву в интервал времени dt (lambdac)
!real :: bettar !< Fraction of soil respiration, [%] from 0 to 1
!real :: sw !< Влажность увядания
! ------ Initialization ------
real :: cs ! Сумма пулов
contains
subroutine carbon_model_init()
! ---- Part of rothc carbon_model_init() ----
s0 = 0.5*(1 + sw)
smin = 1.7*sw
integer :: i !< count
fdpm = alphadr/(1 + alphadr)
allocate(s0(i0:i1,j0:j1,ml))
allocate(smin(i0:i1,j0:j1,ml))
! ---- Part of enviromental carbon_model_init() ----
mnc = (date_lst%y - date_fst%y)*nmonth + (date_lst%m - date_fst%m)
allocate(in_temp(0:mnc))
allocate(in_veg(0:mnc))
allocate(in_wsoil(0:mnc))
allocate(in_lambd(0:mnc))
open (unit = 1, file = 'initial_value/'//trim(station_name)//'_Temp.txt', status='unknown')
read(1,*) in_temp(1:mnc)
close (1)
open (unit = 1, file = 'initial_value/'//trim(station_name)//'_Veg.txt', status='unknown')
read(1,*) in_veg(1:mnc)
close (1)
open (unit = 1, file = 'initial_value/'//trim(station_name)//'_MOI.txt', status='unknown')
read(1,*) in_Wsoil(1:mnc)
close (1)
open (unit = 1, file = 'initial_value/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
read(1,*) in_lambd(1:mnc)
close (1)
select case (station_name)
case('Rostov')
sw = 0.146
bettar = 0.41
station_n = 5
case('DAO3', 'DAO4')
sw = 0.120
bettar = 0.25
if (station_name == 'DAO3') then
station_n = 2
else
station_n = 1
end if
case('VLDMR')
sw = 0.120
bettar = 0.08
station_n = 3
case('TRGK')
sw = 0.120
bettar = 0.07
station_n = 4
case default
stop "check failed : unknown station name"
end select
in_lambd = in_lambd/dt
! ---- Part of timesteps init ----
! @todo учесть шаг по времени 1 час
if (dt == 86400) then ! Для шага по времени день:
mncX = (date_fst%y - mnclot_fst(station_n))*nmonth + date_lst%m - date_fst%m
in_temp(0) = in_temp(1)
in_veg(0) = in_veg(1)
in_wsoil(0) = in_wsoil(1)
in_lambd(0) = in_lambd(1)
else ! Для шага по времени месяц:
mncX = (date_fst%y - mnclot_fst(station_n))*nmonth + date_lst%m
end if
do i = 1, ml
s0(:,:,i) = 0.5*(1. + sw(:,:))
smin(:,:,i) = 1.7*sw(:,:)
end do
fdpm = alphadr/(1 + alphadr)
end subroutine
......@@ -154,48 +78,23 @@ contains
subroutine carbon_model_calc_at_cell(ii,jj)
! ---- Part of environmental_calc_at_cell(ii,jj) ----
integer, intent(in) :: ii, jj
integer i, N(12),j !< count
! ------------- Переменные для распределения значений по всему временному интервалу
i = date_c%m ! Номер месяца
j = date_c%d ! Номер дня месяца
N(i) = date_c%days(i) ! Количество дней в месяце
!mncX ! Номер месяца с начала работы программы
! -------------
if (dt == 86400) then ! Для шага по времени день:
Temp = in_temp(mncX)*sin(pi*(j-1)/(N(i)-1)) + in_temp(mncX+1)*(1-sin(pi*(j-1)/(N(i)-1)))
Wsoil = in_wsoil(mncX)*sin(pi*(j-1)/(N(i)-1)) + in_wsoil(mncX+1)*(1-sin(pi*(j-1)/(N(i)-1)))
veg = in_veg(mncX)
lambd = in_lambd(mncX)/N(i)
if (j == (N(i)/2)) mncX = mncX + 1
else ! Для шага по времени месяц:
Temp = in_temp(mncX)
Wsoil = in_wsoil(mncX)
veg = in_veg(mncX)
lambd = in_lambd(mncX)
mncX = mncX + 1
end if
Tsoil = Temp ! температура почвы равна температуре воздуха
! ---- Part of rothc_calc_at_cell(ii,jj) ----
if (Wsoil >= s0) then
Fs = 1 - 0.8*(Wsoil - s0)
else if (smin < Wsoil .and. Wsoil < s0) then
Fs = 0.2 + 0.8*(Wsoil - smin)/(s0 - smin)
else if (Wsoil <= smin) then
if (Wsoil(ii,jj,1) >= s0(ii,jj,1)) then
Fs = 1 - 0.8*(Wsoil(ii,jj,1) - s0(ii,jj,1))
else if (smin(ii,jj,1) < Wsoil(ii,jj,1) .and. Wsoil(ii,jj,1) < s0(ii,jj,1)) then
Fs = 0.2 + 0.8*(Wsoil(ii,jj,1) - smin(ii,jj,1))/(s0(ii,jj,1) - smin(ii,jj,1))
else if (Wsoil(ii,jj,1) <= smin(ii,jj,1)) then
Fs = 0.2
end if
if (Tsoil >= -18.27) then
Ft = 47.91/(1 + exp(106.06/(Tsoil + 18.27)))
if (Tsoil(ii,jj,1) >= -18.27) then
Ft = 47.91/(1 + exp(106.06/(Tsoil(ii,jj,1) + 18.27)))
else
Ft = 0
end if
Fv = 0.6 + 0.4*(1 - veg)
Fv = 0.6 + 0.4*(1 - veg(ii,jj))
RDPM = Ft*Fv*Fs*ks(1)*CDPM(ii,jj,1)
RRPM = Ft*Fv*Fs*ks(2)*CRPM(ii,jj,1)
......@@ -204,6 +103,9 @@ contains
Rs = RDPM + RRPM + RBIO + RHUM
! ---- Initialization ------
Cs = CDPM(ii,jj,1) + CRPM(ii,jj,1) + CBIO(ii,jj,1) + CHUM(ii,jj,1)
end subroutine
......@@ -222,9 +124,14 @@ contains
use grid, only : date_c
open(unit=500, file='results/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
open(unit=500, file='results/'//trim(carbon_model_type)//'/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
write(500,*) date_c%timestamp,';',CDPM(:,:,1),';',CRPM(:,:,1),';',CBIO(:,:,1),';',CHUM(:,:,1)
! ---- Initialization ------
!print*, 'finish thing 1, ', CDPM(:,:,1)/cs*100,'%'
!print*, 'finish thing 2, ', CRPM(:,:,1)/cs*100,'%'
!print*, 'finish thing 3, ', CBIO(:,:,1)/cs*100,'%'
!print*, 'finish thing 4, ', CHUM(:,:,1)/cs*100,'%'
end subroutine
......
source/environment/environment_data_station.f90
+ 10
− 10
View file @ fab2c539
......@@ -37,8 +37,8 @@ module environment_data_station
real, allocatable :: in_veg(:) !< Поступление извне данных средней концентрации растений в месяц, [%]
real, allocatable :: in_wsoil(:) !< Поступление извне данных влажности почвы, [dim]
real, allocatable :: in_lambd(:) !< Поступление извне данных по поступлению углерода в почву, [kg/m**3 / year]
! ---- For INMCM ----
! station = DAO4 DAO3 TRGK VLDR ROST
! ---- For INMCM ----
real :: rhodry_in(station_max) = (/1.26, 1.26, 1.23, 1.30, 2.62/)
! 1.0 1.0 1.28
! ---- For Rothc ----
......@@ -74,16 +74,16 @@ contains
read(station_opt,*) opt_n
select case (station_name)
case('Rostov')
station_n = 5
case('DAO3')
station_n = 2
case('DAO4')
station_n = 1
case('VLDMR')
station_n = 4
case('DAO3')
station_n = 2
case('TRGK')
station_n = 3
case('VLDMR')
station_n = 4
case('Rostov')
station_n = 5
end select
open (unit = 1, file = 'initial_value/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
......@@ -165,7 +165,7 @@ contains
mncX = mncX + 1
end if
!if (mncX == mnc) mncX = 0
if (mncX == mnc) mncX = 0
do i = 1, ml
Tsoil(:,:,i) = Temp(:,:) ! температура почвы равна температуре воздуха
......
ui1_config.nml
+ 23
− 23
View file @ fab2c539
......@@ -12,7 +12,7 @@
! -----------------------------------------------------------------
!> Углеродная модель:
carbon_model_type = 'inmcm'
carbon_model_type = 'rothc'
!
! 'inmcm' - модель inmcm [1,2]
! 'socs' - модель SOCS [3]
......@@ -32,21 +32,22 @@
! 'lsm_online' - онлайн-каплинг с моделью деятельного слоя (@todo пока не поддерживается)
! дополнительно для environment_data_type = 'station':
station_name = 'Rostov'
!> Имя станции:
station_name = 'DAO4'
! 'Rostov' - Станция ФАНЦ
! 'DAO3' - Станция в Долгопрудном 1
! 'DAO4' - Станция в Долгопрудном 2
! 'VLDMR' - Данные Владимир
! 'TRGK' - Данные Торжок
station_opt = '1'
!
!> Тип подачи удобрения:
station_opt = '4'
! '1' - контрольный случай, без дополнительных удобрений
! '2' - 2 способ подачи удобрений
! '3' - 3 способ подачи удобрений
! '4' - 4 способ подачи удобрений
! дополнительно для environment_data_type = 'lsm_offline':
!> Путь к файлу:
lsm_datafile = 'data/inmcm_ru-fyo_fluxnet_1998-2015_0.5h.nc', UTC = +3
!lsm_datafile = 'data/inmcm_etr_era5_2016-2020_1h_ref.nc'
......@@ -100,12 +101,11 @@
!> Дата и время на момент старта:
datetime_init_mode = 'manual' !auto
datetime_init_mode = 'manual'
!
! 'auto' - получить из входного файла (для lsm_offline)
! 'manual' - задать вручную [yyyy-mm-dd hh:mm:ss] 1937 1939 1956 1968 1974
datetime_init = '1974-01-01 00:00:00'
! (для lsm_offline) 'auto' - получить из входного файла
datetime_init = '1937-01-01 00:00:00'
!> Продолжительность расчета:
ntime_mode = 'datetime_last'
......@@ -113,7 +113,7 @@
! 'ntime' - указанное число шагов
ntime = 1000
! 'datetime_last' - до достижения указанной даты 2012 2018 2975
datetime_last = '2018-01-01 00:00:00'
datetime_last = '2012-01-01 00:00:00'
! (для lsm_offline) 'auto' - до конца входного файла
! продвинутые настройки
......@@ -129,7 +129,7 @@
testing_log_mode = 'read' ! none, write, read
nv_singlecolumn = 12 ! Для модели INMCM нужно поставить 12, для rothc и SOCS 1
nv_singlecolumn = 1 ! Для модели INMCM нужно поставить 12, для rothc и SOCS 1
environment_model_type = 'inmcm' ! inmcm
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment