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

Add 2 new basic function Piecewise1, Piecewise2

parent 57a18ad5
No related branches found
No related tags found
No related merge requests found
...@@ -19,7 +19,7 @@ module carbon_core ...@@ -19,7 +19,7 @@ module carbon_core
integer, parameter :: nflux_default = 20 !< макс число потоков integer, parameter :: nflux_default = 20 !< макс число потоков
integer, parameter :: nmult_default = 10 !< макс число множителей в выражении для потока integer, parameter :: nmult_default = 10 !< макс число множителей в выражении для потока
integer, parameter :: narg_default = 1 !< макс число аргументов функций в множителях integer, parameter :: narg_default = 1 !< макс число аргументов функций в множителях
integer, parameter :: npar_default = 5!8 !< макс число параметров функций в множителях integer, parameter :: npar_default = 8 !< макс число параметров функций в множителях
! производные типы данных ! производные типы данных
! --------------------------------------------------------------------------------- ! ---------------------------------------------------------------------------------
......
...@@ -129,7 +129,9 @@ contains ...@@ -129,7 +129,9 @@ contains
& x0, x0_n, x0_ij, & & x0, x0_n, x0_ij, &
& k, k_n, k_ij, & & k, k_n, k_ij, &
& a, a_n, a_ij, & & a, a_n, a_ij, &
& amp, amp_n, amp_ij) & amp, amp_n, amp_ij, &
& c01, c01_n, c01_ij, &
& c02, c02_n, c02_ij)
! --------------------------------------- ! ---------------------------------------
!< @brief добавить множитель в выражение для потока !< @brief добавить множитель в выражение для потока
...@@ -149,6 +151,8 @@ contains ...@@ -149,6 +151,8 @@ contains
real, optional, intent(in) :: k, k_n(:), k_ij(:,:) real, optional, intent(in) :: k, k_n(:), k_ij(:,:)
real, optional, intent(in) :: a, a_n(:), a_ij(:,:) real, optional, intent(in) :: a, a_n(:), a_ij(:,:)
real, optional, intent(in) :: amp, amp_n(:), amp_ij(:,:) real, optional, intent(in) :: amp, amp_n(:), amp_ij(:,:)
real, optional, intent(in) :: c01, c01_n(:), c01_ij(:,:)
real, optional, intent(in) :: c02, c02_n(:), c02_ij(:,:)
integer :: m integer :: m
...@@ -164,7 +168,9 @@ contains ...@@ -164,7 +168,9 @@ contains
& x0, x0_n, x0_ij, & & x0, x0_n, x0_ij, &
& k, k_n, k_ij, & & k, k_n, k_ij, &
& a, a_n, a_ij, & & a, a_n, a_ij, &
& amp, amp_n, amp_ij) & amp, amp_n, amp_ij, &
& c01, c01_n, c01_ij, &
& c02, c02_n, c02_ij)
call set_fun(flux(fid)%mult_kit(m)%fun_kit, & call set_fun(flux(fid)%mult_kit(m)%fun_kit, &
& funtype) & funtype)
......
...@@ -41,10 +41,10 @@ contains ...@@ -41,10 +41,10 @@ contains
fun_kit%fun => fun_mm fun_kit%fun => fun_mm
case('step') case('step')
fun_kit%fun => fun_step fun_kit%fun => fun_step
case('spline') case('PieWise1')
fun_kit%fun => fun_spline fun_kit%fun => fun_PieWise1
! case('inter') case('PieWise2')
! fun_kit%fun => fun_inter fun_kit%fun => fun_PieWise2
case default case default
stop "check failed : unknown funtype at set_fun" stop "check failed : unknown funtype at set_fun"
end select end select
...@@ -67,7 +67,7 @@ contains ...@@ -67,7 +67,7 @@ contains
! --------------------------------------------------------------------------------- ! ---------------------------------------------------------------------------------
! args(1) - x ! args(1) - x
! pars(1) - y0 или c ! pars(1) - y0
! pars(2) - x0 ! pars(2) - x0
! pars(3) - k ! pars(3) - k
! pars(4) - a ! pars(4) - a
...@@ -82,7 +82,7 @@ contains ...@@ -82,7 +82,7 @@ contains
real, intent(in) :: args(narg_default) real, intent(in) :: args(narg_default)
real, intent(in) :: pars(npar_default) real, intent(in) :: pars(npar_default)
real :: y real :: y
y = pars(1) y = pars(6)
end function fun_const end function fun_const
!> Linear function: y = k*(x-x0) + y0 !> Linear function: y = k*(x-x0) + y0
...@@ -131,14 +131,41 @@ contains ...@@ -131,14 +131,41 @@ contains
y = pars(5) * pars(4)**( pars(3) * (args(1) - pars(2)) ) + pars(1) y = pars(5) * pars(4)**( pars(3) * (args(1) - pars(2)) ) + pars(1)
end function fun_exp end function fun_exp
!> Spline (piecewise1) function: y = (a+k)*x + (y0+c)/2 + (a-k)/2*|x-x0| !> Piecewise1 function: y = θ1(x1-x)*(k1*x + b1) + θ2(x-x1)*(k2*x + b2)
function fun_spline(args, pars) result(y) function fun_PieWise1(args, pars) result(y)
implicit none implicit none
real, intent(in) :: args(narg_default) real, intent(in) :: args(narg_default)
real, intent(in) :: pars(npar_default) real, intent(in) :: pars(npar_default)
real :: y real :: y, theta1, theta2
y = (pars(4)+pars(3))*args(1) + (pars(1)+pars(5))/2 + (pars(4)-pars(3))/2*abs(args(1)-pars(2)) theta1 = 0.5*(1. + sign(1.,pars(2) - args(1)))
end function fun_spline theta2 = 0.5*(1. + sign(1.,args(1) - pars(2)))
y = theta1*(pars(3)*args(1) + pars(7)) + theta2*(pars(4)*args(1) + pars(8))
end function fun_PieWise1
!> Piecewise2 function: y = θ1(x1-x)*(k1*x + b1) + θ2(x-x1)(x2-x)*(k2*x + b2) + θ3(x-x2)*(k3*x + b3)
function fun_PieWise2(args, pars) result(y)
implicit none
real, intent(in) :: args(narg_default)
real, intent(in) :: pars(npar_default)
real :: y, theta1, theta2, theta3
theta1 = 0.5*(1. + sign(1.,pars(2) - args(1)))
theta2 = 0.5*(1. + sign(1.,(args(1) - pars(2))*(pars(5) - args(1)) ))
theta3 = 0.5*(1. + sign(1.,args(1) - pars(5)))
y = theta1*(pars(3)*args(1) + pars(7)) + theta2*(pars(4)*args(1) + pars(8)) + &
& theta3*(pars(6)*args(1) + pars(1))
end function fun_PieWise2
! args(1) - x
! pars(1) - y0
! pars(2) - x0
! pars(3) - k
! pars(4) - a
! pars(5) - amp
! pars(6) - c
! pars(7) - c1
! pars(8) - c2
! !> Inerpolation function: y = (k+amp)*x + b + (a-k)/2*|x-x0| + (amp-a)*|x-y0| ! !> Inerpolation function: y = (k+amp)*x + b + (a-k)/2*|x-x0| + (amp-a)*|x-y0|
! function fun_inter(args, pars) result(y) ! function fun_inter(args, pars) result(y)
......
...@@ -22,6 +22,9 @@ module carbon_model_to_core_par_kit ...@@ -22,6 +22,9 @@ module carbon_model_to_core_par_kit
real, parameter :: k_default = 1. real, parameter :: k_default = 1.
real, parameter :: a_default = exp(1.) real, parameter :: a_default = exp(1.)
real, parameter :: amp_default = 1. real, parameter :: amp_default = 1.
real, parameter :: c_default = 0.
real, parameter :: c01_default = 0.
real, parameter :: c02_default = 0.
contains contains
...@@ -36,7 +39,9 @@ contains ...@@ -36,7 +39,9 @@ contains
& x0, x0_n, x0_ij, & & x0, x0_n, x0_ij, &
& k, k_n, k_ij, & & k, k_n, k_ij, &
& a, a_n, a_ij, & & a, a_n, a_ij, &
& amp, amp_n, amp_ij) & amp, amp_n, amp_ij, &
& c01, c01_n, c01_ij, &
& c02, c02_n, c02_ij)
! --------------------------------------- ! ---------------------------------------
!< @brief установка вида функций и их параметров !< @brief установка вида функций и их параметров
...@@ -52,6 +57,8 @@ contains ...@@ -52,6 +57,8 @@ contains
real, optional, intent(in) :: k, k_n(ntile), k_ij(i0:i1,j0:j1) real, optional, intent(in) :: k, k_n(ntile), k_ij(i0:i1,j0:j1)
real, optional, intent(in) :: a, a_n(ntile), a_ij(i0:i1,j0:j1) real, optional, intent(in) :: a, a_n(ntile), a_ij(i0:i1,j0:j1)
real, optional, intent(in) :: amp, amp_n(ntile), amp_ij(i0:i1,j0:j1) real, optional, intent(in) :: amp, amp_n(ntile), amp_ij(i0:i1,j0:j1)
real, optional, intent(in) :: c01, c01_n(ntile), c01_ij(i0:i1,j0:j1)
real, optional, intent(in) :: c02, c02_n(ntile), c02_ij(i0:i1,j0:j1)
integer :: n integer :: n
...@@ -63,11 +70,14 @@ contains ...@@ -63,11 +70,14 @@ contains
par_kit(n)%par3 = miss_v par_kit(n)%par3 = miss_v
enddo enddo
! par_kit(1) - y0 или c ! par_kit(1) - y0
! par_kit(2) - x0 ! par_kit(2) - x0
! par_kit(3) - k ! par_kit(3) - k
! par_kit(4) - a ! par_kit(4) - a
! par_kit(5) - amp ! par_kit(5) - amp
! par_kit(6) - c
! par_kit(7) - c01
! par_kit(8) - c02
if (present(y0)) then if (present(y0)) then
par_kit(1)%par1 = y0 par_kit(1)%par1 = y0
...@@ -78,15 +88,6 @@ contains ...@@ -78,15 +88,6 @@ contains
elseif (present(y0_ij)) then elseif (present(y0_ij)) then
par_kit(1)%par3 = y0_ij par_kit(1)%par3 = y0_ij
par_kit(1)%num = 3 par_kit(1)%num = 3
elseif (present(c)) then
par_kit(1)%par1 = c
par_kit(1)%num = 1
elseif (present(c_n)) then
par_kit(1)%par2 = c_n
par_kit(1)%num = 2
elseif (present(c_ij)) then
par_kit(1)%par3 = c_ij
par_kit(1)%num = 3
else else
par_kit(1)%par1 = y0_default par_kit(1)%par1 = y0_default
par_kit(1)%num = 1 par_kit(1)%num = 1
...@@ -148,6 +149,49 @@ contains ...@@ -148,6 +149,49 @@ contains
par_kit(5)%num = 1 par_kit(5)%num = 1
endif endif
if (present(c)) then
par_kit(6)%par1 = c
par_kit(6)%num = 1
elseif (present(c_n)) then
par_kit(6)%par2 = c_n
par_kit(6)%num = 2
elseif (present(c_ij)) then
par_kit(6)%par3 = c_ij
par_kit(6)%num = 3
else
par_kit(6)%par1 = c_default
par_kit(6)%num = 1
endif
if (present(c01)) then
par_kit(7)%par1 = c01
par_kit(7)%num = 1
elseif (present(c01_n)) then
par_kit(7)%par2 = c01_n
par_kit(7)%num = 2
elseif (present(c01_ij)) then
par_kit(7)%par3 = c01_ij
par_kit(7)%num = 3
else
par_kit(7)%par1 = c01_default
par_kit(7)%num = 1
endif
if (present(c02)) then
par_kit(8)%par1 = c02
par_kit(8)%num = 1
elseif (present(c02_n)) then
par_kit(8)%par2 = c02_n
par_kit(8)%num = 2
elseif (present(c02_ij)) then
par_kit(8)%par3 = c02_ij
par_kit(8)%num = 3
else
par_kit(8)%par1 = c02_default
par_kit(8)%num = 1
endif
end subroutine end subroutine
......
...@@ -162,31 +162,18 @@ contains ...@@ -162,31 +162,18 @@ contains
! Fs = 0.2 ! Fs = 0.2
! end if ! end if
!a = Fs
if (aTSMD(ii,jj) > 0.444*SMD(ii,jj)) then if (aTSMD(ii,jj) > 0.444*SMD(ii,jj)) then
Fs = 1. Fs = 1.
else else
Fs = 0.2 + (1. - 0.2)*(SMD(ii,jj) - aTSMD(ii,jj))/(SMD(ii,jj) - 0.444*SMD(ii,jj)) Fs = 0.2 + (1. - 0.2)*(SMD(ii,jj) - aTSMD(ii,jj))/(SMD(ii,jj) - 0.444*SMD(ii,jj))
end if end if
!b = Fs
if (Tsoil(ii,jj,1) >= -5.) then
Ft = 47.91/(1. + exp(106.06/(Tsoil(ii,jj,1) + 18.27))) Ft = 47.91/(1. + exp(106.06/(Tsoil(ii,jj,1) + 18.27)))
else
Ft = 0
end if
! ---- Initialization ------ ! ---- Initialization ------
Cs = CDPM(ii,jj,1) + CRPM(ii,jj,1) + CBIO(ii,jj,1) + CHUM(ii,jj,1) Cs = CDPM(ii,jj,1) + CRPM(ii,jj,1) + CBIO(ii,jj,1) + CHUM(ii,jj,1)
print*, 'moi_bal:', moi_bal(ii,jj)
print*, 'pr', pr(ii,jj)
print*, 'Evpar:', Evpar(ii,jj)
print*, 'evpar_part:', evpar_part(ii,jj)
end subroutine end subroutine
......
...@@ -59,7 +59,12 @@ module carbon_model_socs ...@@ -59,7 +59,12 @@ module carbon_model_socs
call set_mult(n_F, 'step', x_ij = Tsoil(:,:,1), x0 = -5., k = 1.0) call set_mult(n_F, 'step', x_ij = Tsoil(:,:,1), x0 = -5., k = 1.0)
! piecewise of Fs ! piecewise of Fs
call set_mult(n_F, 'lin', x = Fs) !call set_mult(n_F, 'lin', x = Fs)
call set_mult(n_F, 'PieWise2', x_ij = Wsoil(:,:,1), x0_ij = smin(:,:,1), y0_ij = s0(:,:,1), k = 0., &
& a_ij = 0.8/(s0(:,:,1) - smin(:,:,1)), amp = -0.8, c01 = 0.2, &
& c02_ij = 0.2 - 0.8*smin(:,:,1)/(s0(:,:,1) - smin(:,:,1)), c_ij = 1 + 0.8*s0(:,:,1))
call set_mult(n_F, 'lin', x = 7.5/yrs) call set_mult(n_F, 'lin', x = 7.5/yrs)
call set_flux(fid = n_F, pid_out = n_Csoil1, pid_in = n_Catm, name = 'microbal_respiration') call set_flux(fid = n_F, pid_out = n_Csoil1, pid_in = n_Catm, name = 'microbal_respiration')
...@@ -68,8 +73,13 @@ module carbon_model_socs ...@@ -68,8 +73,13 @@ module carbon_model_socs
! piecewise of Ft ! piecewise of Ft
call set_mult(n_F, 'lin', x = Ft) call set_mult(n_F, 'lin', x = Ft)
call set_mult(n_F, 'step', x_ij = Tsoil(:,:,1), x0 = -5., k = 1.0) call set_mult(n_F, 'step', x_ij = Tsoil(:,:,1), x0 = -5., k = 1.0)
! piecewise of Fs ! piecewise of Fs
call set_mult(n_F, 'lin', x = Fs) !call set_mult(n_F, 'lin', x = Fs)
call set_mult(n_F, 'PieWise2', x_ij = Wsoil(:,:,1), x0_ij = smin(:,:,1), y0_ij = s0(:,:,1), k = 0., &
& a_ij = 0.8/(s0(:,:,1) - smin(:,:,1)), amp = -0.8, c01 = 0.2, &
& c02_ij = 0.2 - 0.8*smin(:,:,1)/(s0(:,:,1) - smin(:,:,1)), c_ij = 1 + 0.8*s0(:,:,1))
call set_mult(n_F, 'lin', x = 7.5/yrs) call set_mult(n_F, 'lin', x = 7.5/yrs)
......
...@@ -33,6 +33,7 @@ module carbon_model_socs_aux ...@@ -33,6 +33,7 @@ module carbon_model_socs_aux
real :: kirk = 0. !< Коэф. скорости разложения C1 real :: kirk = 0. !< Коэф. скорости разложения C1
!real :: kd !< Коэф. перехода углерода из C2 в C1, десорбация, разрушение агрегатов !real :: kd !< Коэф. перехода углерода из C2 в C1, десорбация, разрушение агрегатов
real, allocatable :: C_atm1(:,:,:) , C_atm2(:,:,:) real, allocatable :: C_atm1(:,:,:) , C_atm2(:,:,:)
real, allocatable :: thet1(:,:,:) , thet2(:,:,:)
contains contains
...@@ -43,6 +44,9 @@ contains ...@@ -43,6 +44,9 @@ contains
allocate(C_atm1(i0:i1,j0:j1,ml)) allocate(C_atm1(i0:i1,j0:j1,ml))
allocate(C_atm2(i0:i1,j0:j1,ml)) allocate(C_atm2(i0:i1,j0:j1,ml))
allocate(thet1(i0:i1,j0:j1,ml))
allocate(thet2(i0:i1,j0:j1,ml))
allocate(s0(i0:i1,j0:j1,ml)) allocate(s0(i0:i1,j0:j1,ml))
allocate(smin(i0:i1,j0:j1,ml)) allocate(smin(i0:i1,j0:j1,ml))
...@@ -75,24 +79,16 @@ contains ...@@ -75,24 +79,16 @@ contains
integer, intent(in) :: ii, jj, nn integer, intent(in) :: ii, jj, nn
if (Wsoil(ii,jj,1) >= s0(ii,jj,1)) then ! if (Wsoil(ii,jj,1) >= s0(ii,jj,1)) then
Fs = 1. - 0.8*(Wsoil(ii,jj,1) - s0(ii,jj,1)) ! 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 ! 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)) ! 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 ! else if (Wsoil(ii,jj,1) <= smin(ii,jj,1)) then
Fs = 0.2 ! Fs = 0.2
end if
! 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 ! end if
Ft = 47.91/(1. + exp(106.06/(Tsoil(ii,jj,1) + 18.27))) Ft = 47.91/(1. + exp(106.06/(Tsoil(ii,jj,1) + 18.27)))
!call set_mult(n_F, 'spline', x_ij = Wsoil - smin, a = 0.8/(s0-smin), k = 0., y0 = 0.2, amp = 0.2, x0 = smin)
end subroutine end subroutine
......
...@@ -41,7 +41,7 @@ ...@@ -41,7 +41,7 @@
! 'Rostov' - Станция ФАНЦ ! 'Rostov' - Станция ФАНЦ
! !
!> Тип подачи удобрения: !> Тип подачи удобрения:
station_opt = '3' station_opt = '1'
! '1' - контрольный случай, без дополнительных удобрений ! '1' - контрольный случай, без дополнительных удобрений
! '2' - 2 способ подачи удобрений ! '2' - 2 способ подачи удобрений
! '3' - 3 способ подачи удобрений ! '3' - 3 способ подачи удобрений
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment