Update the Rothc model by new information from soil scientists

data/obs_data_Rostov.csv

 date,1,8,11
-1974-01-01,8.9088,8.9088,8.9088
+1975-01-01,8.9088,8.9088,8.9088
 1983-09-15,8.595600000000001,8.873999999999999,9.013200000000001
 1992-09-07,8.282399999999999,8.9088,9.1176
 1998-08-02,8.0388,8.9088,9.1872

initial_value/Evpar.txt

@@ -13,18 +13,6 @@
 0
 0
 0
-115
-184
-262
-171
-165
-125
-0
-0
-0
-0
-0
-0
 84
 167
 224

initial_value/Rostov_1.txt

@@ -14,17 +14,6 @@
 0
 0
 0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
-0
 0.0250587
 0.0501174
 0.08770545

initial_value/Rostov_2.txt

@@ -13,18 +13,6 @@
 0
 0
 0
-0
-0
-0
-0
-0.18
-0
-0
-0
-0
-0
-0
-0
 0.0253395
 0.050679
 0.08868825

initial_value/Rostov_3.txt

@@ -13,18 +13,6 @@
 0
 0
 0
-0
-0
-0
-0
-0.36
-0
-0
-0
-0
-0
-0
-0
 0.0253314
 0.0506628
 0.0886599

initial_value/Rostov_MOI.txt

-0.412552965583605
-0.424120203859102
-0.417746586157707
-0.419412464455061
-0.419073572546531
-0.334889556039574
-0.29103367664772
-0.256135976139186
-0.274615792741693
-0.344031471499806
-0.409927574051255
-0.432412848030488
 0.426525111499154
 0.419310388578998
 0.421111007032755

initial_value/Rostov_Temp.txt

--7
--2.6
-2.2
-7.8
-15.5
-20.5
-22
-21.4
-17.7
-14
-4
-0.7
 -1.1
 -2.9
 3

initial_value/Rostov_Veg.txt

@@ -10,18 +10,6 @@
 1
 1
 1
-0
-0
-0
-0
-0
-0
-0
-0
-1
-1
-1
-1
 1
 1
 1

initial_value/rainfall.txt

@@ -10,18 +10,6 @@
 15
 24
 40
-12
-8
-38
-29
-26
-7
-45
-5
-62
-15
-24
-40
 74
 4
 12

nonlin.f90

@@ -117,7 +117,7 @@ program nonlin
 
 
      Flux(0) = 0. !w*Us(0)*C0
-     Flux(1) = w*Us(1)*C0*h(1)/(2.*h(1)+h(2)) !/1000.
+     Flux(1) = -C0*w*h(1)/(2.*h(1) + h(2)) !/1000. w*Us(1)*C0*h(1)/(2.*h(1)+h(2))
 
       do i = 2, N
         Cz(i) = M10(i)*C(i-2) + M11(i)*C(i-1) + M12(i)*C(i)

source/carbon/carbon_model_to_core_arg_kit.f90

@@ -25,7 +25,7 @@ module carbon_model_to_core_arg_kit
     integer, parameter :: year_min = 1900      !< предусмотренный диапазон лет
     integer, parameter :: year_max = 2025      !<
 	
-    ! (44)*12 = 528 ! Для случая Rostov:  1974 - 2018
+    ! (44)*12 = 528 ! Для случая Rostov:  1975 - 2018
     ! (73)*12 = 876 ! Для случая DAO3:    1939 - 2012
     ! (75)*12 = 900 ! Для случая DAO4:    1937 - 2012
     ! (62)*12 = 744 ! Для случая TRGK:    1956 - 2018

source/carbon/carbon_models/carbon_model_ROTHC.f90

@@ -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_ij = (0.46*bettar/100))
+            call set_mult(n_F, 'lin', x = Rs, k_ij = BIO_prop) !  (0.46*bettar/100)
 
             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_ij = (0.54*bettar/100))
+            call set_mult(n_F, 'lin', x = Rs, k_ij = HUM_prop) !  (0.54*bettar/100)
 
 
         end subroutine carbon_model_assembly

source/carbon/carbon_models/carbon_model_ROTHC_aux.f90

@@ -5,7 +5,7 @@ module carbon_model_rothc_aux
  	use environment_core, only: Tsoil, Wsoil, lambd, veg, sw, bettar
 	use config, only: station_name, station_opt, carbon_model_type
     use const,  only: yrs, pi, nmonth
-    use grid,   only: date_c, i0, i1, j0, j1, ml
+    use grid,   only: date_c, i0, i1, j0, j1, ml, dt
 ! ---------------------------------------------- Variables ------------------------------------------------------------
     implicit none
 
@@ -37,10 +37,18 @@ module carbon_model_rothc_aux
     real :: Ft, Fs, Fv                       !< Функции: температуры, влажности, фракции растительного покрова [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)/)
+    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-6), 9.65*(1.E-9), 2.12*(1.E-8), 6.43*(1.E-9)/)
+    !real :: ks(ntiles) = (/ 10., 0.3, 0.66, 0.02/)    
     !< Скорость дыхания единицы массы каждого пула в стандартных условиях, [1-s]
     real, parameter :: alphadr = 1.44        !< Определяет соотношение поступления опада между DPM и RPM, [dim],
     ! Для деревьев, для кустарников, для натуральной травы, для СХ (/ 0.25, 0.33, 0.67, 1.44 /)
+    
+    real, allocatable :: SA(:,:)       ! Определяет распределение выдыхаемой части пула по остальным пулам
+    real, allocatable :: CO2_prop(:,:) ! Определяет распределение выдыхаемой части пула по остальным пулам
+    real, allocatable :: BIO_prop(:,:) ! Определяет распределение выдыхаемой части пула по остальным пулам
+    real, allocatable :: HUM_prop(:,:) ! Определяет распределение выдыхаемой части пула по остальным пулам
+    
     ! ------ Climate variables --------
     ! ------ Determined within
     real, allocatable :: s0(:,:,:)           !< Оптимальная влажность почвы, [dim]
@@ -52,17 +60,21 @@ module carbon_model_rothc_aux
     real :: cs                               ! Сумма пулов
     
     ! ----- The topsoil moisture deficit -------
-   !  real, allocatable :: TSMD(:,:)       ! the topsoil moisture deficit
-   !  real, allocatable :: rainfall(:,:)   ! осадки (из внешних данных)
-   !  real, allocatable :: Evpar(:,:)      ! Evaporation
-   !  real, allocatable :: evpar_part(:,:) ! 75% from evaporation
-   !  real, allocatable :: moi_bal(:,:)    ! rainfall - evaporation
-   !  real, allocatable :: aTSMD(:,:)      ! accumulated topsoil mositure deficit
-   !  real :: b                            ! moisture coefficient    
-   !  integer :: mncX    
-   !  integer :: mnc = 528 ! Случай для Ростова
-   !  real, allocatable :: in_rainfall(:)
-   !  real, allocatable :: in_Evpar(:)
+     real, allocatable :: TSMD(:,:)       ! the topsoil moisture deficit
+     real :: deep = 30.                   ! глубина слоя
+     real, allocatable :: SMD(:,:)        ! TSMD на актуальной глубине
+     real, allocatable :: bSMD(:,:)       ! Случай TSMD для голой почвы (bare soil)
+     real, allocatable :: moi_bal(:,:)    ! rainfall minus evaporation
+     real, allocatable :: aTSMD(:,:)      ! accumulated topsoil mositure deficit 
+
+     real, allocatable :: in_rainfall(:)
+     real, allocatable :: in_Evpar(:)
+     real, allocatable :: rainfall(:,:)   ! осадки (из внешних данных)
+     real, allocatable :: Evpar(:,:)      ! Evaporation
+     real, allocatable :: evpar_part(:,:) ! 75% from evaporation  
+     integer :: mncX    
+     integer :: mnc = 516 ! Случай для Ростова
+
 contains
 
     subroutine carbon_model_init()
@@ -72,26 +84,36 @@ contains
       allocate(s0(i0:i1,j0:j1,ml))
       allocate(smin(i0:i1,j0:j1,ml))
 
+     allocate(SA(i0:i1,j0:j1))
+     allocate(CO2_prop(i0:i1,j0:j1))
+     allocate(BIO_prop(i0:i1,j0:j1))
+     allocate(HUM_prop(i0:i1,j0:j1))
+
      ! init of TPMD
-    ! allocate(TSMD(i0:i1,j0:j1))
-    ! allocate(rainfall(i0:i1,j0:j1))
-    ! allocate(Evpar(i0:i1,j0:j1))
-    ! allocate(evpar_part(i0:i1,j0:j1))
-    ! allocate(moi_bal(i0:i1,j0:j1))
-    ! allocate(aTSMD(i0:i1,j0:j1))
-    ! allocate(in_rainfall(0:mnc))
-    ! allocate(in_Evpar(0:mnc))
-
-    !  TSMD(:,:) = -(20. + 1.3*bettar(:,:) - 0.01*bettar(:,:)**2)
-
-    !  open (unit = 1, file = 'initial_value/rainfall.txt', status='unknown')
-    !  read(1,*) in_rainfall(1:mnc)
-    !  close (1)
-    !  open (unit = 1, file = 'initial_value/Evpar.txt', status='unknown')
-    !  read(1,*) in_Evpar(1:mnc)
-    !  close (1)
-
-    !  mncX = 1
+     allocate(TSMD(i0:i1,j0:j1))
+     allocate(SMD(i0:i1,j0:j1))
+     allocate(bSMD(i0:i1,j0:j1))
+     allocate(moi_bal(i0:i1,j0:j1))
+     allocate(aTSMD(i0:i1,j0:j1))
+
+     allocate(in_rainfall(0:mnc))
+     allocate(in_Evpar(0:mnc))
+     allocate(rainfall(i0:i1,j0:j1))
+     allocate(Evpar(i0:i1,j0:j1))
+     allocate(evpar_part(i0:i1,j0:j1))
+
+      TSMD(:,:) = -(20. + 1.3*bettar(:,:) - 0.01*bettar(:,:)**2)
+      SMD(:,:) = TSMD(:,:)/23.*deep
+      bSMD(:,:) = SMD(:,:)/1.8
+
+      open (unit = 1, file = 'initial_value/rainfall.txt', status='unknown')
+      read(1,*) in_rainfall(1:mnc)
+      close (1)
+      open (unit = 1, file = 'initial_value/Evpar.txt', status='unknown')
+      read(1,*) in_Evpar(1:mnc)
+      close (1)
+
+      mncX = 1
 
       ! main init part of Rothc
       do i = 1, ml
@@ -101,6 +123,16 @@ 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.)
+      HUM_prop(:,:) = 0.54/(SA(:,:) + 1.)
+
     end subroutine
 
 
@@ -117,37 +149,21 @@ contains
       integer, intent(in) :: ii, jj
 
      ! ---- TPMD_calc_at_cell ----
-     ! rainfall(ii,jj) = in_rainfall(mncX)
-     ! Evpar(ii,jj) = in_Evpar(mncX)
-
-     ! !rainfall(ii,jj) = 48.2
-     ! !Evpar(ii,jj) = 81.8
-
-     ! evpar_part(ii,jj) = 0.75*Evpar(ii,jj)
-     ! moi_bal(ii,jj) = rainfall(ii,jj) - evpar_part(ii,jj)
-
-     ! if (moi_bal(ii,jj) < 0. .or. aTSMD(ii,jj) < 0.) then
-     !   aTSMD(ii,jj) = moi_bal(ii,jj) + aTSMD(ii,jj)
-    !    if (aTSMD(ii,jj) < TSMD(ii,jj)) aTSMD(ii,jj) = TSMD(ii,jj)
-     !   if (aTSMD(ii,jj) > 0.) aTSMD(ii,jj) = 0.
-     ! end if
-
-     ! if (aTSMD(ii,jj) > 0.444*TSMD(ii,jj)) then
-     !   b = 1.0
-     ! else
-     !   b = 0.2 + (1.0 - 0.2)*(TSMD(ii,jj) - aTSMD(ii,jj))/(TSMD(ii,jj) - 0.444*TSMD(ii,jj))
-     ! end if
-     ! mncX = mncX + 1
-! --------------
-     ! main part of rothc_calc_at_cell
+      rainfall(ii,jj) = in_rainfall(mncX)
+      Evpar(ii,jj) = in_Evpar(mncX)
 
-     ! Tsoil(ii,jj,1) = 10.2
-     ! Wsoil(ii,jj,1) = 0.355704637
+      evpar_part(ii,jj) = 0.75*Evpar(ii,jj)
+      moi_bal(ii,jj) = rainfall(ii,jj) - evpar_part(ii,jj)
+
+      if (veg(ii,jj) == 1) then
+        aTSMD(ii,jj) = max(SMD(ii,jj), min(0.,aTSMD(ii,jj) + moi_bal(ii,jj)))
+      else
+        aTSMD(ii,jj) = max(min(bSMD(ii,jj),aTSMD(ii,jj)),min(0.,aTSMD(ii,jj) + moi_bal(ii,jj)))
+      end if
 
-     ! if (Tsoil(ii,jj,1) <= 0.) then ! Если влажность зануляется при отрицательной температуре
-     !   Wsoil(ii,jj,1) = 0.
-     !   s0(ii,jj,1) = 0.
-     ! end if
+      mncX = mncX + 1
+! --------------
+     ! main part of rothc_calc_at_cell
 
       if (Wsoil(ii,jj,1) >= s0(ii,jj,1)) then
         Fs = 1. - 0.8*(Wsoil(ii,jj,1) - s0(ii,jj,1))
@@ -157,20 +173,37 @@ contains
         Fs = 0.2
       end if
 
-     ! Fs = b
 
-      if (Tsoil(ii,jj,1) >= -18.27) then
+      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
         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(ii,jj))
 
-      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)
+      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)*Ft*Fv*Fs*ks(1)
+      !RRPM = CRPM(ii,jj,1)*Ft*Fv*Fs*ks(2)
+      !RBIO = CBIO(ii,jj,1)*Ft*Fv*Fs*ks(3)
+      !RHUM = CHUM(ii,jj,1)*Ft*Fv*Fs*ks(4)
+      
+      print*, 'RDPM:', RDPM
+      print*, 'RRPM:', RRPM
+      print*, 'RBIO:', RBIO
+      print*, 'RHUM:', RHUM
 
       Rs = RDPM + RRPM + RBIO + RHUM
 

source/carbon/carbon_models/carbon_model_SOCS.f90

@@ -48,7 +48,7 @@ module carbon_model_socs
 
             call set_flux(fid = n_F, pid_out = n_Csoil2, pid_in = n_Csoil1, name = 'destabilization')
             call set_mult(n_F, 'lin', x_ijn = C2)
-            call set_mult(n_F, 'lin', x = kd)            
+            call set_mult(n_F, 'lin', x_ij = kd)            
 
             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.)

source/environment/environment_data_station.f90

@@ -24,7 +24,7 @@ module environment_data_station
     integer :: opt_n                       !< Номер способа подачи удобрения
 	integer :: station_n !                 DAO4 DAO3 TRGK VLDR ROST   !< Номер станции наблюдения за климатом
                          !                  1    2    3    4    5
-    integer :: mnclot_fst(station_max) = (/1937,1939,1956,1968,1974/) !< Даты начала сбора данных по климату, в рамках станции 
+    integer :: mnclot_fst(station_max) = (/1937,1939,1956,1968,1975/) !< Даты начала сбора данных по климату, в рамках станции 
     integer :: mnclot_lst(station_max) = (/2012,2012,2018,2018,2018/) !< Даты конца сбора данных по климату, в рамках станции 
 ! ---------------------------------------------------------------------------------------------------------------------
     ! ------ Determined externally  ------
@@ -38,7 +38,7 @@ module environment_data_station
 	                                                                                   ! 2.62
     ! ---- For Rothc ----
     real(kind = 16) :: sw_in    (station_max) = (/0.120,    0.120,    0.120,    0.120,   0.146/) 
-    real(kind = 16) :: bettar_in(station_max) = (/25.,      25.,      7.,       8.,      41./)
+    real(kind = 16) :: bettar_in(station_max) = (/25.,      25.,      7.,       8.,      44./)
     ! ---- For SOCs ----
     real :: kd_in    (station_max,opt_max)	
 	!          opt =    1      2      3      4

ui1_config.nml

@@ -12,7 +12,7 @@
 ! -----------------------------------------------------------------
 
     !> Углеродная модель:
-    carbon_model_type = 'socs'
+    carbon_model_type = 'rothc'
     !
     ! 'inmcm' - модель inmcm [1,2]
     ! 'socs'  - модель SOCS [3]
@@ -41,7 +41,7 @@
             ! 'TRGK'   - Данные Торжок
             !    
             !> Тип подачи удобрения:
-            station_opt = '3'
+            station_opt = '1'
             ! '1' - контрольный случай, без дополнительных удобрений
             ! '2' - 2 способ подачи удобрений
             ! '3' - 3 способ подачи удобрений
@@ -104,8 +104,8 @@
     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'
+    ! 'manual' - задать вручную [yyyy-mm-dd hh:mm:ss] 1937 1939 1956 1968 1975  
+        datetime_init = '1975-01-01 00:00:00'
 
     !> Продолжительность расчета:
     ntime_mode = 'datetime_last'