converting formulas to a new format

source/carbon/carbon_model_to_core_fun_kit.f90

@@ -84,7 +84,7 @@ contains
         real, intent(in) :: args(narg_default)
         real, intent(in) :: pars(npar_default)
         real :: y
-        y = args(1) * pars(3) + pars(1)
+        y = pars(3) * (args(1) - pars(2)) + pars(1)
     end function fun_lin
 
     !> Hyperbolic function: y = k/(x-x0) + y0
@@ -93,7 +93,7 @@ contains
         real, intent(in) :: args(narg_default)
         real, intent(in) :: pars(npar_default)
         real :: y
-        y = pars(3) / (args(1) - pars(2))
+        y = pars(3) / (args(1) - pars(2)) + pars(1)
     end function fun_hyp
 
     !> Michaelis-Menthen function: y = k*x/(x-x0) + y0
@@ -102,7 +102,7 @@ contains
         real, intent(in) :: args(narg_default)
         real, intent(in) :: pars(npar_default)
         real :: y
-        y = pars(3) * args(1) / (args(1) - pars(2))
+        y = pars(3) * args(1) / (args(1) - pars(2)) + pars(1)
     end function fun_mm
 
     !> Step function: y = k*θ(x-x0) + y0

source/carbon/carbon_models/carbon_model_INMCM_aux.f90

@@ -169,8 +169,6 @@ contains
         integer :: k
         real :: work
 
-            !Wsoil(ii,jj,:) = 1. !0.21447610589527216609628319344309
-            !Tsoil(ii,jj,:) = 40.
             work = 0.
             do k = ms+1, ml-1
                 if (z(k) <= hint) then
@@ -188,7 +186,6 @@ contains
         use grid, only : date_c
 
         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

@@ -34,26 +34,105 @@ module carbon_model_rothc
 
             ! 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_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(1))
+            call set_mult(n_F, 'lin', x_ijn = CDPM)
+
+
+            !call set_mult(n_F, 'const', c_n = 47.91)
+            !call set_mult(n_F, 'exp',   x_ij = Tsoil, a = exp, k = 106.06, x0 = -18.27, y0 = 1.)
+            !call set_mult(n_F, 'step',  x_ij = Tsoil, x0 = -4.99)
+
 
             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_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(2))
+            call set_mult(n_F, 'lin', x_ijn = CRPM)
 
             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_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(3))
+            call set_mult(n_F, 'lin', x_ijn = CBIO)
 
             call set_flux(fid = n_F, pid_out = n_CHUM, pid_in = n_Catm, name = 'RHUM')
-            call set_mult(n_F, 'lin', x = RHUM)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(4))
+            call set_mult(n_F, 'lin', x_ijn = CHUM)
 
             ! atmosphere fall
-            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*Rs')
-            call set_mult(n_F, 'lin', x = Rs, k_ij = BIO_prop)
-
-            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*Rs')
-            call set_mult(n_F, 'lin', x = Rs, k_ij = HUM_prop)
-
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*RDPM')
+            call set_mult(n_F, 'lin', x_ij = BIO_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(1))
+            call set_mult(n_F, 'lin', x_ijn = CDPM)
+            
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*RRPM')
+            call set_mult(n_F, 'lin', x_ij = BIO_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(2))
+            call set_mult(n_F, 'lin', x_ijn = CRPM)
+            
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*RBIO')
+            call set_mult(n_F, 'lin', x_ij = BIO_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(3))
+            call set_mult(n_F, 'lin', x_ijn = CBIO)
+
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*RHUM')
+            call set_mult(n_F, 'lin', x_ij = BIO_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(4))
+            call set_mult(n_F, 'lin', x_ijn = CHUM)
+
+
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RDPM')
+            call set_mult(n_F, 'lin', x_ij = HUM_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(1))
+            call set_mult(n_F, 'lin', x_ijn = CDPM)
+
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RRPM')
+            call set_mult(n_F, 'lin', x_ij = HUM_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(2))
+            call set_mult(n_F, 'lin', x_ijn = CRPM)
+
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RBIO')
+            call set_mult(n_F, 'lin', x_ij = HUM_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(3))
+            call set_mult(n_F, 'lin', x_ijn = CBIO)
+
+            call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RHUM')
+            call set_mult(n_F, 'lin', x_ij = HUM_prop)
+            call set_mult(n_F, 'lin', x = Ft)
+            call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
+            call set_mult(n_F, 'lin', x = Fs)
+            call set_mult(n_F, 'lin', x = ks(4))
+            call set_mult(n_F, 'lin', x_ijn = CHUM)
 
         end subroutine carbon_model_assembly
 
 end module
\ No newline at end of file
-

source/carbon/carbon_models/carbon_model_ROTHC_aux.f90

@@ -34,10 +34,12 @@ module carbon_model_rothc_aux
     real :: RBIO = 0.0                                 !< Дыхание пула BIO
     real :: RHUM = 0.0                                 !< Дыхание пула HUM
     ! ------- Functions ------------
-    real :: Ft, Fs, Fv                       !< Функции: температуры, влажности, фракции растительного покрова [dim]
-    real :: fdpm                             !< Функция определяющая разлагаемую часть опада, [dim]
+    real :: Ft = 0.                                    !< Функция температуры [dim]
+    real :: Fs = 0.                                    !< Функция влажности [dim]
+    real :: Fv = 0.                                    !< Функция фракции растительного покрова [dim]           
+    real :: fdpm = 0.                                  !< Функция определяющая разлагаемую часть опада, [dim]
     ! ------- Coefficients ------------
-    real :: ks(ntiles) = (/ 3.22*(1.E-7), 9.65*(1.E-9), 2.12*(1.E-8), 6.43*(1.E-10)/)   
+    real :: 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 /)
@@ -55,9 +57,9 @@ module carbon_model_rothc_aux
     !real :: bettar(i0:i1,j0:j1)             !< Fraction of soil respiration, [%] from 0 to 1
     !real :: sw                              !< Влажность увядания
     ! ------ Initialization ------
-    real :: cs                               !< Сумма пулов
+    real :: cs = 0.                          !< Сумма пулов
     ! ----- The topsoil moisture deficit -------
-     real :: deep = 30.                      ! глубина слоя (@todo сопоставить с переменной Ml)
+    real :: deep = 30.                      ! глубина слоя, см (@todo сопоставить с переменной Ml)
     real, allocatable :: TSMD      (:,:)    !< the topsoil moisture deficit
     real, allocatable :: SMD       (:,:)    !< TSMD на актуальной глубине
     real, allocatable :: bSMD      (:,:)    !< Случай TSMD для голой почвы (bare soil)
@@ -87,6 +89,21 @@ contains
       allocate(aTSMD(i0:i1,j0:j1))
       allocate(evpar_part(i0:i1,j0:j1))
 
+      TSMD = 0.
+      SMD = 0.
+      bSMD = 0.
+      moi_bal = 0.
+      aTSMD = 0.
+      evpar_part = 0.
+     
+      SA = 0.
+      CO2_prop = 0.
+      BIO_prop = 0.
+      HUM_prop = 0.
+      
+      s0 = 0.
+      smin = 0.
+
       TSMD(:,:) = -(20. + 1.3*bettar(:,:) - 0.01*bettar(:,:)**2)
       SMD(:,:) = TSMD(:,:)/23.*deep
       bSMD(:,:) = SMD(:,:)/1.8
@@ -99,11 +116,6 @@ contains
 
       fdpm = alphadr/(1 + alphadr)
 
-      SA = 0.
-      CO2_prop = 0.
-      BIO_prop = 0.
-      HUM_prop = 0.
-
       SA(:,:) = 1.67*(1.85 + 1.6*exp(-0.0786*bettar(:,:)))
       CO2_prop(:,:) = SA(:,:)/(SA(:,:) + 1.)
       BIO_prop(:,:) = 0.46/(SA(:,:) + 1.)
@@ -146,11 +158,11 @@ contains
       end if
 
 
-      if (aTSMD(ii,jj) > 0.444*SMD(ii,jj)) then
-        Fs = 1.
-      else
-        Fs = 0.2 + (1. - 0.2)*(SMD(ii,jj) - aTSMD(ii,jj))/(SMD(ii,jj) - 0.444*SMD(ii,jj))
-      end if
+      !if (aTSMD(ii,jj) > 0.444*SMD(ii,jj)) then
+      !  Fs = 1.
+      !else
+      !  Fs = 0.2 + (1. - 0.2)*(SMD(ii,jj) - aTSMD(ii,jj))/(SMD(ii,jj) - 0.444*SMD(ii,jj))
+      !end if
 
 
       if (Tsoil(ii,jj,1) >= -5.) then
@@ -159,14 +171,19 @@ contains
         Ft = 0
       end if
 
-      Fv = 0.6 + 0.4*(1. - veg(ii,jj))
+      !Fv = 0.6 + 0.4*(1. - veg(ii,jj))
+
+      !RDPM = CDPM(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(1)*dt))/dt
+      !RRPM = CRPM(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(2)*dt))/dt
+      !RBIO = CBIO(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(3)*dt))/dt
+      !RHUM = CHUM(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(4)*dt))/dt
 
-      RDPM = CDPM(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(1)*dt))/dt
-      RRPM = CRPM(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(2)*dt))/dt
-      RBIO = CBIO(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(3)*dt))/dt
-      RHUM = CHUM(ii,jj,1)*(1. - exp(-Ft*Fv*Fs*ks(4)*dt))/dt
+      !RDPM = Ft*Fv*Fs*ks(1)*CDPM(ii,jj,1)
+      !RRPM = Ft*Fv*Fs*ks(2)*CRPM(ii,jj,1)
+      !RBIO = Ft*Fv*Fs*ks(3)*CBIO(ii,jj,1)
+      !RHUM = Ft*Fv*Fs*ks(4)*CHUM(ii,jj,1)
 
-      Rs = RDPM + RRPM + RBIO + RHUM
+      !Rs = RDPM + RRPM + RBIO + RHUM
 
       ! ---- Initialization ------
       Cs = CDPM(ii,jj,1) + CRPM(ii,jj,1) + CBIO(ii,jj,1) + CHUM(ii,jj,1)

source/carbon/carbon_models/carbon_model_SOCS_aux.f90

@@ -32,6 +32,7 @@ module carbon_model_socs_aux
     real :: Cm = 12.         !< Max кол-во органического углерода которое может быть защищено в почве
     real :: kirk = 0.        !< Коэф. скорости разложения C1
 	!real :: kd               !< Коэф. перехода углерода из C2 в C1, десорбация, разрушение агрегатов
+    real, allocatable :: C_atm1(:,:,:) , C_atm2(:,:,:) 
 
 contains
 
@@ -39,6 +40,9 @@ contains
 
       integer :: i !< count
 
+      allocate(C_atm1(i0:i1,j0:j1,ml))
+      allocate(C_atm2(i0:i1,j0:j1,ml))
+
       allocate(s0(i0:i1,j0:j1,ml))
       allocate(smin(i0:i1,j0:j1,ml))
 
@@ -95,7 +99,10 @@ contains
       !kd = 0.007/yrs !0.012 007
       !print*, kd      
       !print*, 0.035/yrs
+      write(530,*) date_c%timestamp,';', 7.5/yrs, Ft*Fs*7.5/yrs
 
+     C_atm1(:,:,1) = C_atm1(:,:,1) + C1(:,:,1)*(1.-rare)*4./yrs
+     C_atm2(:,:,1) = C_atm2(:,:,1) + C1(:,:,1)*(1.-rare)*Ft*Fs*7.5/yrs
     end subroutine
 
 
@@ -106,6 +113,7 @@ contains
 
       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) 
+      print*, C_atm1(:,:,1), C_atm2(:,:,1)
 
 	end subroutine
 

ui1_config.nml

@@ -64,10 +64,10 @@
     !> Пространственная сетка:
     spatial_grid_mode = 'manual'
     !
+    ! 'auto' - использовать сетку из входного netcdf-файла (для lsm_offline) 
     ! 'none'   - не требуется
     ! 'manual' - задать вручную [dlon, dlat] (@todo в режиме lsm_offline переинтерполяция не предусмотрена)
         spatial_grid_res(:) = 0.5, 0.5
-    ! (для lsm_offline) 'auto' - использовать сетку из входного netcdf-файла
 
         ! дополнительно для spatial_grid_mode = 'manual' или 'auto':
 
@@ -92,13 +92,12 @@
     !> Шаг по времени:
     dt_mode = 'manual'
     !
+    ! 'auto' - использовать dt из входного netcdf-файла (для lsm_offline) 
     ! 'manual' - задать вручную [c] (@todo режиме lsm_offline переинтерполяция не предусмотрена)
         dt = 2629800 
         ! '3600'    - 1 час
         ! '86400'   - 1 день
         ! '2629800' - 1 месяц
-    !  (для lsm_offline) 'auto' - использовать dt из входного netcdf-файла
-
 
     !> Дата и время на момент старта:
     datetime_init_mode = 'manual' 
@@ -110,11 +109,11 @@
     !> Продолжительность расчета:
     ntime_mode = 'datetime_last' 
     !
+    ! 'auto' - до конца входного файла (для lsm_offline)
     ! 'ntime' - указанное число шагов
         ntime = 1000
     ! 'datetime_last' - до достижения указанной даты 2012 2018 2975
         datetime_last = '2018-01-01 00:00:00'
-    ! (для lsm_offline) 'auto' - до конца входного файла
 
 ! продвинутые настройки
 ! -----------------------------------------------------------------