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

Big upgrade: add humidity of soil, carbon intake, initialization

parent 024fe212
No related branches found
No related tags found
No related merge requests found
Showing
with 507 additions and 21 deletions
0.071019834989088
0.0610112990645545
0.0503732443249746
0.0514795620943083
0.0502044758900219
0.0469449347498228
0.0598391677136819
0.0696076051984264
0.127847393851185
0.105430074931238
0.0941460890763334
0.0776258919910274
\ No newline at end of file
0.344667980428586
0.306241937762506
0.259597392525599
0.25839805237134
0.243648486597415
0.217886398136598
0.26104169976804
0.276927816333597
0.436763320410713
0.419444403285626
0.410701820529836
0.360286495192741
\ No newline at end of file
0.0805681635225855
0.0698134109936384
0.058022
0.0589065284851548
0.0569542638854971
0.052634083099385
0.0659615308894217
0.0747096379133999
0.130782587714246
0.113157789308251
0.103778518318311
0.0870331947736925
\ No newline at end of file
0.365476061634786
0.331176598946869
0.285299233380516
0.279437195258345
0.258357881675492
0.225451483277588
0.261766010614518
0.266375070095089
0.396457983332979
0.403460850576339
0.411841392420236
0.372795757058519
\ No newline at end of file
0.0865488194689593
0.0756517716655732
0.0633010280575017
0.0638327676465571
0.0611820362720759
0.0558871154447951
0.0689036005034001
0.0761453855979163
0.127903742203471
0.115332422173861
0.108407331828413
0.0924122530008228
\ No newline at end of file
0.371136442326107
0.339726564688768
0.295172984611121
0.286651408019961
0.262359248983505
0.226168085403905
0.258692827220753
0.258331429526321
0.375200464410152
0.391277675685208
0.407006295132092
0.373980695956709
\ No newline at end of file
0.321186494247689
0.28688637556055
0.244230185318244
0.242066391162298
0.227049240668277
0.201692013182339
0.239541965481311
0.251118516336518
0.389326309286388
0.38035274907833
0.376875883829632
0.33350823713317
\ No newline at end of file
0.310158932772697
0.277398764886249
0.236404835732225
0.23406102084036
0.21925377135633
0.194447560067398
0.230447623914992
0.240894860910728
0.371965494966236
0.364867648642256
0.362567585036274
0.321529157003837
\ No newline at end of file
0.262430650395193
0.235124435913182
0.200665607546636
0.198391169899117
0.185514275872262
0.16416376561568
0.194006268323069
0.202033953983955
0.310301906867695
0.306007581304697
0.305233714250582
0.271453327317091
\ No newline at end of file
0.0837991551907178
0.0722950458377144
0.0598823381244564
0.0610004598881861
0.0592382771238908
0.0550718259943483
0.0696083929045549
0.079893477053301
0.143175566206182
0.121009410512475
0.109516194981339
0.0910641294778255
\ No newline at end of file
0.0463741084720905
0.04072393782247
0.0342
0.0343616759190522
0.032782219376686
0.0297628292233768
0.0363881436104635
0.0397226787272637
0.0654288425843253
0.0601653365712916
0.0572501513735094
0.0492143865776693
\ No newline at end of file
0.0956100457997356
0.0819003956975829
0.0674730963333333
0.0691051751152043
0.0675874879170618
0.0634508722097931
0.0813530919562843
0.0955333898857609
0.178640303103547
0.144698160860139
0.127994351101248
0.104919318327788
\ No newline at end of file
-0.9
3.3
1.5
12.8
15.2
23.5
24.2
26.5
17.2
11.5
5.0
0.1
0
0
0
1
1
1
1
1
1
1
1
0
\ No newline at end of file
......@@ -6,6 +6,11 @@ import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
#nml = f90nml.read(self.nml_path)
#start_lon = nml['drivpar_basic']['lam_min_']
'''
Script for plotting.
Prerequisites:
......@@ -17,17 +22,35 @@ Usage:
set model_data_paths, obs_data_path, pools in plot.py
'''
#if enviromental = station
model_data_paths = [
'results/Fast1.txt',
#'results/Fast2.txt',
#'results/Fast3.txt'
'results/VLDMR_1.txt', #Rostov_obs_1.1
'results/VLDMR_2.txt',
'results/VLDMR_3.txt'
]
obs_data_path = 'data/obs_data_fast.csv'
obs_data_path = 'data/obs_data_VLDMR.csv' #obs_data_fast_1.1
#obs_data_path = f'data/obs_data_{station}.csv'
#else
#model_data_paths = [
#'results/common.txt'
#]
#if carbon_model_type = INMCM
pools = {'Csoil' : 'Почва',
'Csoilb': 'Почва типа b'
}
#elif carbon_model_type = ROTHC
#pools = {'CDPM': 'Разлагаемый растительный материал',
# 'CRPM': 'Устойчивый растительный материал',
# 'CBIO': 'Пул микробной биомассы',
# 'CHUM': 'Долгоживущий гумусовый пул'}
#elif carbon_model_type = SOCS
#pools = {'csoil1': 'Огранический углерод в почве',
# 'csoil2': 'Минерализованный углерод в почве',}
#else
#pools = {'C1': 'Pool 1',
# 'C2': 'Pool 2',}
pools = {'CDPM': 'Разлагаемый растительный материал',
'CRPM': 'Устойчивый растительный материал',
'CBIO': 'Пул микробной биомассы',
'CHUM': 'Долгоживущий гумусовый пул'}
def get_experiment_name(path):
......@@ -120,7 +143,7 @@ def pools_sums_plots(file_names: List[str], obs_data: pd.DataFrame, out_path: st
for i, file_name in enumerate(file_names):
df, experiment_name, opt, _ = read_data(file_name)
pools_sum = df[df.columns[1:]].sum(axis=1)
ax.plot(df['date'], pools_sum, label=f'Результаты расчетов, {experiment_name}, {int(opt)}', linewidth=2,
ax.plot(df['date'], pools_sum, label=f'Результаты расчетов {experiment_name} {int(opt)}', linewidth=2,
color=colors[color_index])
date_min = min(date_min, df['date'].min())
date_max = max(date_max, df['date'].max())
......@@ -131,7 +154,7 @@ def pools_sums_plots(file_names: List[str], obs_data: pd.DataFrame, out_path: st
ax.scatter(obs_data['date'],
obs_data[col],
color=colors[color_index],
label=f'Данные наблюдений, {i + 1}',
label=f'Данные наблюдений {i + 1}',
s=50)
color_index += 1
if i + 1 >= len(file_names):
......
......@@ -13,14 +13,19 @@ module carbon_model_to_core_arg_kit
private
public :: set_args
public :: get_args
public :: year_min, year_max, nmonth
! параметры
! ---------------------------------------------------------------------------------
integer, parameter :: nmonth = 12
integer, parameter :: year_min = 2022 !< предусмотренный диапазон лет
integer, parameter :: year_max = 2023
integer, parameter :: year_max = 2523 !<
! (43+1)*12 = 528 ! Для случая Rostov: 1974 - 2017
! (72+1)*12 = 876 ! Для случая DAO3: 1939 - 2011
! (74+1)*12 = 900 ! Для случая DAO4: 1937 - 2011
! (61+1)*12 = 744 ! Для случая TRGK: 1956 - 2017
! (49+1)*12 = 600 ! Для случая VLDMR: 1968 - 2017
contains
......
module carbon_model_rothc
use carbon_model_to_core, only : set_tiles, set_pool, set_flux, set_mult
use carbon_model_rothc_aux
implicit none
private
public :: carbon_model_assembly
contains
subroutine carbon_model_assembly()
integer :: n_Catm, n_Cveg, n_CDPM, n_CRPM, n_CBIO, n_CHUM !< id пулов
integer :: n_F !< id потока
call set_tiles(ntiles)
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.03589581800064, alias = CDPM)
call set_pool(pid = n_CRPM, name = 'CRPM', initial_value = 0.5878742278272 , alias = CRPM)
call set_pool(pid = n_CBIO, name = 'CBIO', initial_value = 0.025696504154112 , alias = CBIO)
call set_pool(pid = n_CHUM, name = 'CHUM', initial_value = 1.05686011826944 , 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) !x_year = lambdac x = lambd
call set_mult(n_F, 'lin', x = 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) !x_year lambd
call set_mult(n_F, 'lin', x = fdpm, k = -1., y0 = 1.)
! to atmosphere
call set_flux(fid = n_F, pid_out = n_CDPM, pid_in = n_Catm, name = 'RDPM')
call set_mult(n_F, 'lin', x = RDPM)
call set_flux(fid = n_F, pid_out = n_CRPM, pid_in = n_Catm, name = 'RRPM')
call set_mult(n_F, 'lin', x = RRPM)
call set_flux(fid = n_F, pid_out = n_CBIO, pid_in = n_Catm, name = 'RBIO')
call set_mult(n_F, 'lin', x = RBIO)
call set_flux(fid = n_F, pid_out = n_CHUM, pid_in = n_Catm, name = 'RHUM')
call set_mult(n_F, 'lin', x = RHUM)
! 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_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))
end subroutine carbon_model_assembly
end module
module carbon_model_rothc_aux
use const, only : yrs, day2sec, pi
use carbon_model_to_core_arg_kit, only : year_min, year_max, nmonth
use grid, only: date_c, date_fst, date_lst, dt
implicit none
! Интерфейс
! -------------------------------------------------------------------------------------------------------------------
! ------- Station of observation ------
character(len=10) :: station = 'VLDMR' !< Станция наблюдения за климатом ! Нужно указать название
character(len=2) :: opt = '1' !< Имя варианта подачи удобрения ! Нужно указать номер
! ------- Serve value -------
integer, parameter :: mnc = (year_max - year_min + 1)*nmonth !Как часть параметров для задания климатических данных
!< Колличество месяцев в расчете
integer :: k !< Номер месяца с начала работы программы !Как часть параметров для шага по времени
! -------------------------------------------------------------------------------------------------------------------
! Main part of rothc
! ------- Pools -------
real, dimension(:,:,:), pointer :: Catm !< атмосфера
real, dimension(:,:,:), pointer :: Cveg !< растительность
real, dimension(:,:,:), pointer :: CDPM !< разлагаемый растительный материал / decomposable plant material
real, dimension(:,:,:), pointer :: CRPM !< устойчивый растительный материал / resistant plant material
real, dimension(:,:,:), pointer :: CBIO !< пул микробной биомассы / microbial biomass pool
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 :: atmosphere1 !< Между атмосферой и почвой BIO
real, dimension(:,:,:), pointer :: atmosphere2 !< Между атмосферой и почвой HUM
real, dimension(:,:,:), pointer :: respiration1 !< Между почвой DPM и атмосферой ! 1 - DPM
real, dimension(:,:,:), pointer :: respiration2 !< Между почвой RPM и атмосферой ! 2 - RPM
real, dimension(:,:,:), pointer :: respiration3 !< Между почвой BIO и атмосферой ! 3 - BIO
real, dimension(:,:,:), pointer :: respiration4 !< Между почвой HUM и атмосферой ! 4 - HUM
! ------- Serve value -------
integer, parameter :: ntiles = 4 !< В данном случае, номер пула
! ------- Initial value ------
real :: Rs = 0.0 !< Функции дыхания
real :: RDPM = 0.0 !< Дыхание пула DPM
real :: RRPM = 0.0 !< Дыхание пула RPM
real :: RBIO = 0.0 !< Дыхание пула BIO
real :: RHUM = 0.0 !< Дыхание пула HUM
! ------- Functions ------------
real :: a1,a2,a3,a4 !< Функции коэффициентов для пулов
real :: Ft, Fs, Fv, Fa !< Функции: температуры, влажности, фракции растительного покрова, общая функция, [dim]
real :: fdpm !< Функция определяющая разлагаемую часть опада, [dim]
! ------- Coefficients ------------
real, parameter :: ks(ntiles) = (/ 3.22*(1.E-7), 9.65*(1.E-9), 2.12*(1.E-8), 6.43*(1.E-10)/)
!< Скорость дыхания единицы массы каждого пула в стандартных условиях, [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]
! ------ Determined externally
real :: in_temp(0:mnc) !< Поступление извне данных cредней температуры воздуха в месяц, [Celsius]
real :: Temp !< Средняя температура воздуха в интервал времени dt, [Celsius]
real :: Tsoil !< Средняя температура почвы в интервал времени dt, [Celsius]
real :: in_veg(0:mnc) !< Поступление извне данных средней концентрации растений в месяц, [%]
real :: veg !< Средняя концентрация растений в интервал времени dt, [%]
real :: in_wsoil(0:mnc) !< Поступление извне данных влажности почвы, [dim]
real :: Wsoil !< Влажность почвы в интервал времени dt, [dim]
real :: sw !< Влажность увядания DAO: 0.120 Rostov: 0.146
! ------ Litterfall --------
real :: in_lambd(0:mnc) !< Поступление извне данных по поступлению углерода в почву, [kg/m**3 / year]
real :: lambd !< Поступлени еуглерода в почву в интервал времени dt (lambdac)
contains
subroutine carbon_model_init()
! ---- Part of rothc carbon_model_init() ----
s0 = 0.5*(1 + sw)
smin = 1.7*sw
fdpm = alphadr/(1 + alphadr)
! ---- Part of enviromental carbon_model_init() ----
open (unit = 1, file = 'initial_value/'//trim(station)//'_Temp.txt', status='unknown')
read(1,*) in_temp(1:mnc)
close (1)
open (unit = 1, file = 'initial_value/'//trim(station)//'_Veg.txt', status='unknown')
read(1,*) in_veg(1:mnc)
close (1)
open (unit = 1, file = 'initial_value/'//trim(station)//'_MOI.txt', status='unknown')
read(1,*) in_Wsoil(1:mnc)
close (1)
open (unit = 1, file = 'initial_value/'//trim(station)//'_'//trim(opt)//'.txt', status='unknown')
read(1,*) in_lambd(1:mnc)
close (1)
select case (station)
case('Rostov')
sw = 0.146
bettar = 0.41
case('DAO3', 'DAO4')
sw = 0.120
bettar = 0.25
case('VLDMR')
sw = 0.120
bettar = 0.08
case('TRGK')
sw = 0.120
bettar = 0.07
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 ! Для шага по времени день:
k = 0
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 ! Для шага по времени месяц:
k = 1
end if
end subroutine
subroutine carbon_model_calc_at_timestep() !Пусто
end subroutine
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) ! Количество дней в месяце
!k ! Номер месяца с начала работы программы
! -------------
if (dt == 86400) then ! Для шага по времени день:
Temp = in_temp(k)*sin(pi*(j-1)/(N(i)-1)) + in_temp(k+1)*(1-sin(pi*(j-1)/(N(i)-1)))
Wsoil = in_wsoil(k)*sin(pi*(j-1)/(N(i)-1)) + in_wsoil(k+1)*(1-sin(pi*(j-1)/(N(i)-1)))
veg = in_veg(k)
lambd = in_lambd(k)/N(i)
if (j == (N(i)/2)) k = k+1
else ! Для шага по времени месяц:
Temp = in_temp(k)
Wsoil = in_wsoil(k)
veg = in_veg(k)
lambd = in_lambd(k)
k = k+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
Fs = 0.2
end if
if (Tsoil >= -18.27) then
Ft = 47.91/(1 + exp(106.06/(Tsoil + 18.27)))
else
Ft = 0
end if
Fv = 0.6 + 0.4*(1 - veg)
Fa = Ft*Fv*Fs
a1 = ks(1)*Fa
a2 = ks(2)*Fa
a3 = ks(3)*Fa
a4 = ks(4)*Fa
RDPM = a1*CDPM(ii,jj,1)
RRPM = a2*CRPM(ii,jj,1)
RBIO = a3*CBIO(ii,jj,1)
RHUM = a4*CHUM(ii,jj,1)
Rs = RDPM + RRPM + RBIO + RHUM
!print*, 'date_fst, ', date_fst%y
!print*, 'date_lst, ', date_lst%y
!print*, 'timestamp_fst, ', timestamp_fst
!print*, 'important thing 4, ', Temp(k)
end subroutine
subroutine carbon_model_calc_at_tile(ii,jj,nn) !Пусто
integer, intent(in) :: ii, jj, nn
end subroutine
subroutine carbon_model_postprocessing() !Вывод переменных
use grid, only : date_c
open(unit=500, file='results/'//trim(station)//'_'//trim(opt)//'.txt', status='unknown')
write(500,*) date_c%timestamp,';',CDPM(:,:,1),';',CRPM(:,:,1),';',CBIO(:,:,1),';',CHUM(:,:,1)
end subroutine
end module
......@@ -130,7 +130,7 @@ module config
! -----------------------------------------------------------------------------
select case(carbon_model_type)
case('inmcm', 'socs', 'other')
case('inmcm', 'socs', 'other', 'rothc')
case default
stop "check failed : unknown carbon_model_type"
end select
......
......@@ -321,19 +321,19 @@ contains
j1 = j1_nc
end select
allocate(lon(i0:i1))
do i = i0, i1
allocate(lon(nlon))
do i = 1, nlon
lon(i) = lon_ref + (i-1)*dlon
enddo
allocate(lat(j0:j1))
do j = j0, j1
allocate(lat(nlat))
do j = 1, nlat
lat(j) = lat_ref + (j-1)*dlat
enddo
allocate(area(i0:i1,j0:j1))
do j = j0, j1
do i = i0, i1
allocate(area(nlon,nlat))
do j = 1, nlat
do i = 1, nlon
area(i,j) = r_earth*abs(dlat)*deg2rad * r_earth*cos(lat(j)*deg2rad)*abs(dlon)*deg2rad
enddo
enddo
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment