From 78979f42afd8a692b316b93232808cc805213a1b Mon Sep 17 00:00:00 2001 From: sumbel <sumbel.enikeeva@gmail.com> Date: Wed, 11 Jun 2025 10:55:08 +0300 Subject: [PATCH] minor refactoring --- build.sh | 9 ++-- source/main_4d_var.f90 | 25 +++++---- .../{carbon/carbon_run.f90 => model_run.f90} | 8 +-- .../carbon_adjoint.f90 | 15 ++---- .../carbon_obs_operator.f90 | 54 +++++++------------ .../carbon_optimize.f90 | 22 +++----- .../optimization/optimization_variables.f90 | 45 ++++++++++++++++ 7 files changed, 99 insertions(+), 79 deletions(-) rename source/{carbon/carbon_run.f90 => model_run.f90} (92%) rename source/{carbon => optimization}/carbon_adjoint.f90 (78%) rename source/{carbon => optimization}/carbon_obs_operator.f90 (70%) rename source/{carbon => optimization}/carbon_optimize.f90 (79%) create mode 100644 source/optimization/optimization_variables.f90 diff --git a/build.sh b/build.sh index bd58474..ae0f53c 100755 --- a/build.sh +++ b/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 diff --git a/source/main_4d_var.f90 b/source/main_4d_var.f90 index 6805b09..512c832 100644 --- a/source/main_4d_var.f90 +++ b/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 diff --git a/source/carbon/carbon_run.f90 b/source/model_run.f90 similarity index 92% rename from source/carbon/carbon_run.f90 rename to source/model_run.f90 index 72f33ce..738e7b5 100644 --- a/source/carbon/carbon_run.f90 +++ b/source/model_run.f90 @@ -1,4 +1,4 @@ -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 diff --git a/source/carbon/carbon_adjoint.f90 b/source/optimization/carbon_adjoint.f90 similarity index 78% rename from source/carbon/carbon_adjoint.f90 rename to source/optimization/carbon_adjoint.f90 index cdc8074..b16f247 100644 --- a/source/carbon/carbon_adjoint.f90 +++ b/source/optimization/carbon_adjoint.f90 @@ -1,14 +1,10 @@ module carbon_adjoint - use grid, only : i0, i1, j0, j1, dt, ntime, date_fst - use carbon_core, only : ntile, npool + 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 diff --git a/source/carbon/carbon_obs_operator.f90 b/source/optimization/carbon_obs_operator.f90 similarity index 70% rename from source/carbon/carbon_obs_operator.f90 rename to source/optimization/carbon_obs_operator.f90 index df57d48..872fb74 100644 --- a/source/carbon/carbon_obs_operator.f90 +++ b/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,13 +72,18 @@ 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)) - cost_o = cost_o + 0.5*(model_obs_diff)**2 * R_inv + 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 end subroutine observation_operator @@ -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 diff --git a/source/carbon/carbon_optimize.f90 b/source/optimization/carbon_optimize.f90 similarity index 79% rename from source/carbon/carbon_optimize.f90 rename to source/optimization/carbon_optimize.f90 index 1da43b7..67d4aca 100644 --- a/source/carbon/carbon_optimize.f90 +++ b/source/optimization/carbon_optimize.f90 @@ -1,36 +1,31 @@ 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 diff --git a/source/optimization/optimization_variables.f90 b/source/optimization/optimization_variables.f90 new file mode 100644 index 0000000..df74779 --- /dev/null +++ b/source/optimization/optimization_variables.f90 @@ -0,0 +1,45 @@ +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 -- GitLab