minor refactoring

build.sh

@@ -65,14 +65,15 @@ gfortran $keys_compile -c source/carbon/carbon_models/carbon_model_SOCS_aux.f90
 gfortran $keys_compile -c source/carbon/carbon_models/carbon_model_SOCS.f90 -o bin/carbon_model_SOCS.o -Jbin
 gfortran $keys_compile -c source/carbon/carbon_models/carbon_model_ROTHC_aux.f90 -o bin/carbon_model_ROTHC_aux.o -Jbin
 gfortran $keys_compile -c source/carbon/carbon_models/carbon_model_ROTHC.f90 -o bin/carbon_model_ROTHC.o -Jbin
-gfortran $keys_compile -c source/carbon/carbon_obs_operator.f90 -o bin/carbon_obs_operator.o -Jbin
-gfortran $keys_compile -c source/carbon/carbon_adjoint.f90 -o bin/carbon_adjoint.o -Jbin
+gfortran $keys_compile -c source/optimization/optimization_variables.f90 -o bin/optimization_variables.o -Jbin
+gfortran $keys_compile -c source/optimization/carbon_obs_operator.f90 -o bin/carbon_obs_operator.o -Jbin
+gfortran $keys_compile -c source/optimization/carbon_adjoint.f90 -o bin/carbon_adjoint.o -Jbin
 gfortran $keys_compile -c source/carbon/carbon_solver.f90 -o bin/carbon_solver.o -Jbin
 gfortran $keys_compile -c source/carbon/carbon_postprocessing.f90 -o bin/carbon_postprocessing.o -Jbin
 gfortran $keys_compile -c source/carbon/carbon.f90 -o bin/carbon.o -Jbin
 #gfortran $keys_compile -c source/main.f90 -o bin/main.o -Jbin
-gfortran $keys_compile -c source/carbon/carbon_run.f90 -o bin/carbon_run.o -Jbin
-gfortran $keys_compile -c source/carbon/carbon_optimize.f90 -o bin/carbon_optimize.o -Jbin
+gfortran $keys_compile -c source/model_run.f90 -o bin/model_run.o -Jbin
+gfortran $keys_compile -c source/optimization/carbon_optimize.f90 -o bin/carbon_optimize.o -Jbin
 gfortran $keys_compile -c source/main_4d_var.f90 -o bin/main_4d_var.o -Jbin
 
 gfortran -o run.exe bin/*.o $keys_compile $keys_linking

source/main_4d_var.f90

@@ -8,18 +8,25 @@ program main
                    & dlon, dlat, dlon_nc, dlat_nc, dt_nc, date_fst, date_lst, ncell_tot, ncell_mask
     use carbon_obs_operator, only: obs_operator_init
 
-    use carbon_optimize, only: fourd_var, nparam
+    use carbon_optimize, only: fourd_var
+    use model_run, only: run_model
+    use optimization_variables, only: nparam, opt_pool_num, theta_b, theta_opt, B_inv, R_inv
 
     implicit none
 
-    integer, parameter :: opt_pool_num(nparam) = [3, 4] !< id of pools to optimize
-    real, parameter :: theta_b(nparam) = [  0., 6. ] !< a priori parameters
-    real, parameter :: B_inv(nparam,nparam) = reshape([1.0,0.0,0.0,0.1],[nparam,nparam]) !< error covariance of a priori
-
     real, allocatable :: cost_hist(:)
-    real :: theta_opt(nparam) !< optimal parameters
+
+
 
     call config_init()
+
+    nparam = 2
+    allocate(opt_pool_num(nparam), theta_b(nparam), theta_opt(nparam), B_inv(nparam, nparam))
+    opt_pool_num = [3, 4] !< id of pools to optimize
+    theta_b = [  0., 6. ] !< a priori parameters
+    B_inv = reshape([1.0,0.0,0.0,0.1],[nparam,nparam]) !< error covariance of a priori
+    R_inv = 1.0/0.1
+
     call grid_init()
     call obs_operator_init()
 
@@ -46,10 +53,8 @@ program main
 
     ! 4D-Var
     ! ---------------------------------------------------------------------------------
-    call fourd_var(theta_b, B_inv, opt_pool_num, &
-            n_iter = 100, lr = 2.5, tol = 1.0e-4, &
-            theta_opt = theta_opt, cost_hist = cost_hist)
-
+    call fourd_var(n_iter = 100, lr = 2.5, tol = 1.0e-4, cost_hist = cost_hist)
     write(*,'("Оптимальные IC: ",2F12.6)') theta_opt
 
+
 end program main

source/carbon/carbon_run.f90 → source/model_run.f90

-module carbon_run
+module model_run
 
     use environment, only : environment_calc_at_timestep, &
                           & environment_calc_at_cell, &
@@ -9,9 +9,9 @@ module carbon_run
                      & carbon_postprocessing
     use grid, only : i0, i1, j0, j1, ntime, ii, jj, tt, nn, mask, date_c, date_fst ,dt
     use carbon_adjoint, only : adjoint_run, &
-                             & adjoint_postprocessing, J_right_side, J_right_side_t
+                             & adjoint_postprocessing
     use carbon_core, only : ntile, tile_weight, nflux, flux, npool
-    use carbon_obs_operator,  only : record_direct
+    use optimization_variables, only: J_right_side, J_right_side_t, record_direct
 
     implicit none
 
@@ -74,4 +74,4 @@ module carbon_run
     end subroutine run_model
 
 
-end module carbon_run
\ No newline at end of file
+end module model_run
\ No newline at end of file

source/carbon/carbon_adjoint.f90 → source/optimization/carbon_adjoint.f90

@@ -2,13 +2,9 @@ module carbon_adjoint
     use grid,      only: i0, i1, j0, j1, dt, ntime, date_fst
     use carbon_core, only: ntile, npool
     use carbon_obs_operator, only: calc_forcing
+    use optimization_variables, only: adj, J_right_side, J_right_side_t
     implicit none
 
-
-    real, allocatable :: adj(:,:,:,:,:) !< C^*
-    real, allocatable :: J_right_side(:,:,:,:,:,:) !< dF/dC(t)
-    real, allocatable :: J_right_side_t(:,:,:,:,:) !< dF/dC
-
 contains
 
     subroutine adjoint_init()
@@ -53,9 +49,4 @@ contains
 
     end subroutine adjoint_postprocessing
 
-
-    subroutine adjoint_deallocate()
-        deallocate(J_right_side, J_right_side_t, adj)
-    end subroutine adjoint_deallocate
-
 end module carbon_adjoint
\ No newline at end of file

source/carbon/carbon_obs_operator.f90 → source/optimization/carbon_obs_operator.f90

@@ -4,13 +4,9 @@ module carbon_obs_operator
     use carbon_core, only: npool, ntile, pool
     use const, only: yrs, day2sec
     use config, only: station_name
+    use optimization_variables, only: nobs, obs_val, obs_idx, state_ts, R_inv, opt_pool_num
 
-    integer :: nobs
-    real, parameter :: R_inv = 1.0/0.1
 
-    integer, allocatable :: obs_idx(:)
-    real,    allocatable :: obs_val(:)
-    real, allocatable :: state_ts(:,:,:,:,:)
 
 contains
 
@@ -52,17 +48,21 @@ contains
         integer, intent(in) :: nn
         integer, intent(in) :: t
 
-        integer :: m
+        integer :: m, k
         real:: forcing(npool)
-        real :: residual
+        real :: residual, residual1
 
         forcing = 0.
         residual = 0.
         do m = 1, nobs
             if (obs_idx(m) == t-1) then
-                residual = (state_ts(3,ii,jj,nn,t-1) + state_ts(4,ii,jj,nn,t-1) - obs_val(m))/yrs
-                forcing(3) = residual * R_inv
-                forcing(4) = residual * R_inv
+                do k = 1, size(opt_pool_num)
+                    residual = residual + state_ts(opt_pool_num(k),ii,jj,nn,t-1)
+                end do
+                residual = (residual - obs_val(m))/yrs
+                do k = 1, size(opt_pool_num)
+                    forcing(opt_pool_num(k)) = residual * R_inv
+                end do
             end if
         end do
     end function calc_forcing
@@ -72,12 +72,17 @@ contains
         use grid, only : i0,i1,j0,j1
 
         real, intent(inout) :: cost_o
-        real :: model_obs_diff
-        integer :: k
+        real :: model_obs_diff, model_obs_diff1
+        integer :: k, m
 
         cost_o = 0.0
+
         do k = 1, nobs
-          model_obs_diff = (state_ts(3,i0,j0,1,obs_idx(k)) + state_ts(4,i0,j0,1,obs_idx(k)) - obs_val(k))
+            model_obs_diff = 0.
+            do m = 1, size(opt_pool_num)
+                model_obs_diff = model_obs_diff + state_ts(opt_pool_num(m),i0,j0,1,obs_idx(k))
+            end do
+            model_obs_diff = model_obs_diff - obs_val(k)
             cost_o = cost_o + 0.5*(model_obs_diff)**2 * R_inv
         end do
 
@@ -108,25 +113,4 @@ contains
         jdn = d + (153*(m + 12*a -3)+2)/5 + 365*b + b/4 - b/100 + b/400 - 32045
     end subroutine ymd_to_jdn
 
-    subroutine record_direct(t)
-
-        integer, intent(in) :: t
-        integer :: p, ii, jj, nn
-
-        do p = 1, npool
-            do jj = j0, j1
-                do ii = i0, i1
-                    do nn = 1, ntile
-                        state_ts(p,ii,jj,nn,t) = pool(p)%val(ii,jj,nn)
-                    end do
-                end do
-            end do
-        end do
-    end subroutine record_direct
-
-
-    subroutine carbon_obs_operator_deallocate()
-        deallocate(state_ts)
-    end subroutine carbon_obs_operator_deallocate
-
 end module carbon_obs_operator
\ No newline at end of file

source/carbon/carbon_optimize.f90 → source/optimization/carbon_optimize.f90

 module carbon_optimize
   use carbon_core,     only : pool, ntile, npool
-  use carbon_run,      only : run_model
-  use carbon_adjoint,  only : adjoint_init, adj, J_right_side, adjoint_deallocate
+  use model_run,      only : run_model
+  use carbon_adjoint,  only : adjoint_init
   use environment, only : environment_init
   use carbon, only : carbon_init, carbon_deallocate
-  use carbon_obs_operator, only: carbon_obs_operator_deallocate, state_ts, observation_operator
+  use carbon_obs_operator, only: observation_operator
+  use optimization_variables, only: nparam, adj, J_right_side, state_ts, theta_b, B_inv, &
+                                    opt_pool_num, theta_opt, optimization_variables_deallocate
 
   implicit none
   private
   public :: fourd_var
-  public :: nparam
 
-  integer, parameter :: nparam = 2
+
   real,    parameter :: eps = 1.0e-12
 
 contains
 
-subroutine fourd_var(theta_b, B_inv, opt_pool_num, &
-                               n_iter, lr, tol,                       &
-                               theta_opt, cost_hist)
+subroutine fourd_var(n_iter, lr, tol, cost_hist)
 
   use grid, only : i0,i1,j0,j1, ntime
 
-  real,    intent(in)  :: theta_b(nparam)      !< a priori parameters
-  real,    intent(in)  :: B_inv(nparam,nparam) !< error covariance of a priori
   integer, intent(in)  :: n_iter               !< number of iterations of opt algo
   real,    intent(in)  :: lr, tol              !< learning rate for Adagrad
-  real,    intent(out) :: theta_opt(nparam)    !< optimal parameters
   real,    allocatable, intent(out) :: cost_hist(:)
 
   integer :: k, it
-  integer :: opt_pool_num(nparam)
   real :: theta(nparam), theta_new(nparam)
   real :: cache(nparam), grad(nparam)
   real :: cost, cost_prev, diff
@@ -129,8 +124,7 @@ subroutine deallocation()
   call carbon_deallocate()
   call carbon_model_to_core_deallocate()
 
-  call adjoint_deallocate()
-  call carbon_obs_operator_deallocate()
+  call optimization_variables_deallocate()
 
 end subroutine deallocation
 

source/optimization/optimization_variables.f90

+module optimization_variables
+    integer :: nparam
+    real, allocatable :: adj(:,:,:,:,:) !< C^*
+    real, allocatable :: J_right_side(:,:,:,:,:,:) !< dF/dC(t)
+    real, allocatable :: J_right_side_t(:,:,:,:,:) !< dF/dC
+
+    integer :: nobs
+    real :: R_inv
+    integer, allocatable :: obs_idx(:)
+    real,    allocatable :: obs_val(:)
+
+    real, allocatable :: state_ts(:,:,:,:,:)
+    real, allocatable :: B_inv(:,:)
+    real, allocatable :: theta_b(:)
+    integer, allocatable :: opt_pool_num(:)
+
+    real, allocatable :: theta_opt(:)
+
+    contains
+
+    subroutine record_direct(t)
+
+        use carbon_core, only: pool, npool, ntile
+        use grid,      only: i0, i1, j0, j1
+
+        integer, intent(in) :: t
+        integer :: p, ii, jj, nn
+
+        do p = 1, npool
+            do jj = j0, j1
+                do ii = i0, i1
+                    do nn = 1, ntile
+                        state_ts(p,ii,jj,nn,t) = pool(p)%val(ii,jj,nn)
+                    end do
+                end do
+            end do
+        end do
+    end subroutine record_direct
+
+    subroutine optimization_variables_deallocate()
+        deallocate(J_right_side, J_right_side_t, adj, state_ts)
+    end subroutine optimization_variables_deallocate
+
+
+end module optimization_variables
\ No newline at end of file