Skip to content
Snippets Groups Projects
Commit b955e20c authored by Sumbel Shangareeva's avatar Sumbel Shangareeva
Browse files

working version of SOCS (+dstep, dPieWise2)

parent 78979f42
No related branches found
No related tags found
No related merge requests found
......@@ -44,10 +44,12 @@ contains
fun_kit%fun => fun_mm
case('step')
fun_kit%fun => fun_step
fun_kit%dfun => dfun_step
case('PieWise1')
fun_kit%fun => fun_PieWise1
case('PieWise2')
fun_kit%fun => fun_PieWise2
fun_kit%dfun => dfun_PieWise2
case default
stop "check failed : unknown funtype at set_fun"
end select
......@@ -153,6 +155,14 @@ contains
y = pars(1) + pars(3) * theta
end function fun_step
function dfun_step(args, pars) result(dy)
implicit none
real, intent(in) :: args(narg_default)
real, intent(in) :: pars(npar_default)
real :: dy
dy = 0
end function dfun_step
!> Exponent function: y = amp*a**[k*(x-x0)] + y0
function fun_exp(args, pars) result(y)
implicit none
......@@ -186,6 +196,19 @@ contains
& theta3*(pars(5)*args(1) + pars(6))
end function fun_PieWise2
function dfun_PieWise2(args, pars) result(dy)
implicit none
real, intent(in) :: args(narg_default)
real, intent(in) :: pars(npar_default)
real :: dy, theta1, theta2, theta3
theta1 = 0.5*(1. + sign(1.,pars(2) - args(1)))
theta2 = 0.5*(1. + sign(1.,(args(1) - pars(2))*(pars(1) - args(1)) ))
theta3 = 0.5*(1. + sign(1.,args(1) - pars(1)))
dy = theta1*pars(3) + theta2*pars(4) + theta3*pars(5)
end function dfun_PieWise2
! args(1) - x
! pars(1) - y0
! pars(2) - x0
......
......@@ -54,11 +54,30 @@ module carbon_model_socs
call set_flux(fid = n_F, pid_out = n_Csoil1, pid_in = n_Csoil2, name = 'mineralization')
call set_mult(n_F, 'lin', x_ijn = C2, k = -1./Cm, y0 = 1., pool_num=n_Csoil2)
call set_mult(n_F, 'const', c = rare)
call set_mult(n_F, 'lin', x = kirk)
call set_mult(n_F, 'lin', x_ijn = C1, pool_num=n_Csoil1)
call set_mult(n_F, 'lin', x = Ft)
call set_mult(n_F, 'step', x_ij = Tsoil(:,:,1), x0 = -5., k = 1.0)
! piecewise of 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 = kirk)
call set_flux(fid = n_F, pid_out = n_Csoil1, pid_in = n_Catm, name = 'microbal_respiration')
call set_mult(n_F, 'lin', x_ijn = C1, k = (1. - rare), pool_num=n_Csoil1)
! piecewise of Ft
call set_mult(n_F, 'lin', x = Ft)
call set_mult(n_F, 'step', x_ij = Tsoil(:,:,1), x0 = -5., k = 1.0)
! piecewise of 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 = kirk)
end subroutine carbon_model_assembly
......
......@@ -71,7 +71,7 @@ contains
integer, intent(in) :: ii, jj, nn
! 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)))
end subroutine
......@@ -81,7 +81,7 @@ contains
use grid, only : date_c
open(unit=500, file='results/'//trim(carbon_model_type)//'/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown')
write(500,*) date_c%timestamp,';',C1(:,:,1),';',C2(:,:,1),';',Catm(:,:,1),';',Cveg(:,:,1)
write(500,*) date_c%timestamp,';',C1(:,:,1),';',C2(:,:,1) !,';',Catm(:,:,1),';',Cveg(:,:,1)
end subroutine
......
......@@ -53,7 +53,7 @@ program main
! 4D-Var
! ---------------------------------------------------------------------------------
call fourd_var(n_iter = 100, lr = 2.5, tol = 1.0e-4, cost_hist = cost_hist)
call fourd_var(n_iter = 20, lr = 2.5, tol = 1.0e-4, cost_hist = cost_hist)
write(*,'("Оптимальные IC: ",2F12.6)') theta_opt
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment