From 11335c2679fe9a608d32124c0e53156ebf5a6605 Mon Sep 17 00:00:00 2001 From: George Faykin <pofnutsy@gmail.com> Date: Sat, 30 Nov 2024 21:34:49 +0300 Subject: [PATCH] Add new environmenmt data type: environment_data_station; SOCs approved --- build.sh | 1 + initial_value/VLDMR_4.txt | 0 nonlin.f90 | 160 +++++++++++++++ .../carbon/carbon_model_to_core_arg_kit.f90 | 2 +- .../carbon_models/carbon_model_INMCM_aux.f90 | 9 +- .../carbon_models/carbon_model_SOCS.f90 | 8 +- .../carbon_models/carbon_model_SOCS_aux.f90 | 132 ++---------- source/config.f90 | 73 +++++-- source/const.f90 | 1 + source/environment/environment.f90 | 17 +- source/environment/environment_core.f90 | 22 +- .../environment_data_generator.f90 | 1 + .../environment/environment_data_station.f90 | 191 ++++++++++++++++++ .../environment/environment_model_INMCM.f90 | 18 +- ui1_config.nml | 27 ++- 15 files changed, 503 insertions(+), 159 deletions(-) create mode 100644 initial_value/VLDMR_4.txt create mode 100644 nonlin.f90 create mode 100644 source/environment/environment_data_station.f90 diff --git a/build.sh b/build.sh index fa41a45..b02f493 100755 --- a/build.sh +++ b/build.sh @@ -48,6 +48,7 @@ gfortran $keys_compile -c source/datetime.f90 -o bin/datetime.o -Jbin gfortran $keys_compile -c source/grid.f90 -o bin/grid.o -Jbin gfortran $keys_compile -c source/environment/environment_core.f90 -o bin/environment_core.o -Jbin gfortran $keys_compile -c source/environment/environment_data_generator.f90 -o bin/environment_data_generator.o -Jbin +gfortran $keys_compile -c source/environment/environment_data_station.f90 -o bin/environment_data_station.o -Jbin gfortran $keys_compile -c source/environment/environment_data_lsm_offline.f90 -o bin/environment_data_lsm_offline.o -Jbin gfortran $keys_compile -c source/environment/environment_model_INMCM.f90 -o bin/environment_model_INMCM.o -Jbin gfortran $keys_compile -c source/environment/environment.f90 -o bin/environment.o -Jbin diff --git a/initial_value/VLDMR_4.txt b/initial_value/VLDMR_4.txt new file mode 100644 index 0000000..e69de29 diff --git a/nonlin.f90 b/nonlin.f90 new file mode 100644 index 0000000..6da8093 --- /dev/null +++ b/nonlin.f90 @@ -0,0 +1,160 @@ +program nonlin + + implicit none + integer, parameter :: N = 11 ! количество элементов сетки + integer, parameter :: a = 12 ! коэффициент сгущения сетки + real :: x(N) ! Значение аргумента, аналог x(i-1/2) в статье + real :: E(N) ! интерполяция + real :: h(N) ! Шаг сетки + + real :: Cor(2:N-1) ! Оригинальное значение поступающей функции + real :: Cr1(N) ! Среднее значение функции (функция разбивается на небольшие интервалы) + + real :: M10(2:N-1), M10p1(2:N-1), M10p2(2:N-1) + real :: M11(2:N-1), M11p1(2:N-1), M11p2(2:N-1) + real :: M12(2:N-1), M12p1(2:N-1), M12p2(2:N-1) + + real :: M00p1(N-1), M00p2(N-1), M00(N-1) + real :: M01p1(N-1), M01p2(N-1), M01(N-1) + real :: M02p1(N-1), M02p2(N-1), M02(N-1) + + real :: CL(2:N-1), CR(1:N-2) ! Реконструкция слева и справа + + real :: InCor + integer :: i !< count + ! ---- initial value ---- + do i = 2, N-1 + Cor(i) = 0.003 + 0.001*i + end do + ! ---- вычисление аргумента сетки, шага по сетке ---- + E = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1./) + do i = 1, N + !E(i) = (i - 1.)/(N - 1.)*1.0 + x(i) = ((1 + a)**E(i) - 1.)/a + end do + + h(1) = x(1) + do i = 2, N + h(i) = x(i) - x(i-1) + end do + ! ---- Вычисление среднего значения поступающей на вход функции Cor(x,t) в ячейке с номером i, ---- + !! На выход Cr1(N-1) + + do i = 2, N-1 + InCor = InCor + Cor(i)*0.5*(x(i+1)-x(i-1)) !Функция Cor(i) интегрируется в пределах x(i-1) и x(i) + end do + h(1) = 0.01 + do i = 1, N + Cr1(i) = 1./h(i)*InCor + end do + + ! ---- Вычисление коэффициентов при значениях среднего значения функции C(x,t) + do i = 2, N-1 + M10p1(i) = -h(i)*h(i+1) + M10p2(i) = (h(i) + h(i-1))*(h(i) + h(i-1) + h(i+1)) + M10(i) = M10p1(i)/M10p2(i) + + M11p1(i) = h(i+1)*(3*h(i)**2 + 3*h(i)*h(i-1) + 2*h(i)*h(i+1) + h(i-1)**2 + h(i-1)*h(i+1)) + M11p2(i) = (h(i)**2 + h(i)*h(i-1) + 2*h(i)*h(i+1) + h(i-1)*h(i+1) + h(i+1)**2)*(h(i) + h(i-1)) + M11(i) = M11p1(i)/M11p2(i) + + M12p1(i) = h(i)*(h(i) + h(i-1)) + M12p2(i) = h(i)**2 + h(i)*h(i-1) + 2*h(i)*h(i+1) + h(i-1)*h(i+1) + h(i+1)**2 + M12(i) = M12p1(i)/M12p2(i) + end do + + do i = 1, N-2 + M00p1(i) = h(i+1)*(h(i+1) + h(i+2)) + M00p2(i) = (h(i) + h(i+1) + h(i+2))*(h(i) + h(i+1)) + M00(i) = M00p1(i)/M00p2(i) + + M01p1(i) = h(i)*(h(i+1) + h(i+2))/((h(i) + h(i+1))*h(i+2)) + M01p2(i) = h(i)*(h(i+1))**2/(h(i+2)*(h(i) + h(i+1) + h(i+2))*(h(i+1) + h(i+2))) + M01(i) = M01p1(i) + M01p2(i) + + M02p1(i) = h(i)*h(i+1) + M02p2(i) = (h(i) + h(i+1) + h(i+2))*(h(i+1) + h(i+2)) + M02(i) = M02p1(i)/M02p2(i) + end do + + ! ---- Реконструкция ---- + do i = 2, N-1 + CL(i) = M10(i)*Cr1(i-1) + M11(i)*Cr1(i) + M12(i)*Cr1(i+1) ! Слева + end do + + do i = 1, N-2 + CR(i) = M00(i)*Cr1(i) + M01(i)*Cr1(i+1) + M02(i)*Cr1(i+2) ! Справа + end do + + print*, 'CL', CL + print*, 'CR', CR + print*, 'Cr1:', Cr1 + print*, 'x:', x + print*, 'h:', h + print*, 'InCor:', InCor + +end program nonlin + + + + + + + +!MODULE MUSCL_MOD + +!часть с use + ! implicit none !некоторые переменные + + +!contains + + ! !> Function action_1 - implements MINMOD limiter after Roe (1986) + ! FUNCTION action_1(a,b) + ! implicit none + ! real, intent(in) :: a, b + ! real :: action_1 + ! action_1 = 0.5*(sign(1.,a) + sign(1.,b))*min(abs(a),abs(b)) + ! END FUNCTION action_1 + + + + !SUBROUTINE Non_lineare_grid(N,array) + + !integer, intent(in) :: N !< Number of grid cells + !real, intent(in ) :: array(1:N) !< Grid variable to be interpolated + + + +! do i = 1, N-1 +! fun(i) = call(f, x(i)) !!!! +! call integral(x(i-1), x(i),n,intgr) +! srdzn(i) = intgr(i)/h(i) +! end do +! +! intgr = h(i)*(sum(fun) - (fun(0) + fun(n))/2 ) +! +! +! function integral(x0,xn,n,f, intgr) +! +! implicit none +! integer, intent(in) :: n +! real, intent(in) :: x0, xn, f +! real :: hop, x(0:n), intgr +! +! hop = (xn - x0)/n +! x(0) = x0 +! x(n) = xn +! +! do i = 1, N-1 +! x(i) = x(i - 1) + hop +! fun(i) = call(f, x(i)) !!!! +! end do +! +! intgr = hop*(sum(fun) - (fun(0) + fun(n))/2 ) +! +! end function integral + + !END SUBROUTINE Non_lineare_grid + +!END MODULE MUSCL_MOD \ No newline at end of file diff --git a/source/carbon/carbon_model_to_core_arg_kit.f90 b/source/carbon/carbon_model_to_core_arg_kit.f90 index dec00b0..d098686 100644 --- a/source/carbon/carbon_model_to_core_arg_kit.f90 +++ b/source/carbon/carbon_model_to_core_arg_kit.f90 @@ -22,7 +22,7 @@ module carbon_model_to_core_arg_kit integer, parameter :: nmonth = 12 - integer, parameter :: year_min = 1974 !< предусмотренный диапазон лет + integer, parameter :: year_min = 1956 !< предусмотренный диапазон лет integer, parameter :: year_max = 2018 !< ! (44)*12 = 528 ! Для случая Rostov: 1974 - 2018 diff --git a/source/carbon/carbon_models/carbon_model_INMCM_aux.f90 b/source/carbon/carbon_models/carbon_model_INMCM_aux.f90 index 024b756..a605f22 100644 --- a/source/carbon/carbon_models/carbon_model_INMCM_aux.f90 +++ b/source/carbon/carbon_models/carbon_model_INMCM_aux.f90 @@ -136,8 +136,6 @@ contains mncX = (date_fst%y - mnclot_fst(station_n))*nmonth + date_lst%m end if - print*, 'mncX', mncX - ! ---- Part of inmcm carbon_model_init() ---- ! ---- Variables from INMCM ---- @@ -221,18 +219,19 @@ contains do jj = j0, j1 do n = 1, nv2 ers_weight(ii,jj,n) = Csoil(ii,jj,n) * soer(n) / sc8 + print*, area(ii,jj), sc8, ers_weight(ii,jj,n) end do end do end do end if - + end subroutine subroutine carbon_model_calc_at_cell(ii,jj) ! ---- Part of environmental_calc_at_cell(ii,jj) ---- - use grid, only: date_c + use grid, only: date_c, area integer, intent(in) :: ii, jj integer i, N(12), j !< count @@ -262,8 +261,8 @@ contains !print*, 'Flit: ', Flit(:,:,12) !print*, 'Fers', Fers(:,:,12) !write(333,*) date_c%timestamp, Fmicr(:,:,12), Fmicrb(:,:,12), Flit(:,:,12), Fers(:,:,12) - !write(334,*) date_c%timestamp, rsw(12), Temp(ii,jj), Csoil(:,:,12) + end subroutine diff --git a/source/carbon/carbon_models/carbon_model_SOCS.f90 b/source/carbon/carbon_models/carbon_model_SOCS.f90 index 63d3f55..e7740b6 100644 --- a/source/carbon/carbon_models/carbon_model_SOCS.f90 +++ b/source/carbon/carbon_models/carbon_model_SOCS.f90 @@ -25,14 +25,14 @@ module carbon_model_socs call set_pool(pid = n_Catm, name = 'catm') call set_pool(pid = n_Cveg, name = 'cveg') - call set_pool(pid = n_Csoil1, name = 'csoil1',initial_value = 0.178176, alias = C1) - call set_pool(pid = n_Csoil2, name = 'csoil2',initial_value = 8.730624, alias = C2) + call set_pool(pid = n_Csoil1, name = 'csoil1',initial_value = 0.060028, alias = C1) !0.056761 + call set_pool(pid = n_Csoil2, name = 'csoil2',initial_value = 2.941372, alias = C2) !2.78124 call set_flux(fid = n_F, pid_out = n_Cveg, pid_in = n_Csoil1, name = 'litterfall') - call set_mult(n_F, 'lin', x = lambd) + call set_mult(n_F, 'lin', x_ij = lambd) 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, k = kd) + call set_mult(n_F, 'lin', x_ijn = C2, k_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.) diff --git a/source/carbon/carbon_models/carbon_model_SOCS_aux.f90 b/source/carbon/carbon_models/carbon_model_SOCS_aux.f90 index 2b9ba4c..df7d777 100644 --- a/source/carbon/carbon_models/carbon_model_SOCS_aux.f90 +++ b/source/carbon/carbon_models/carbon_model_SOCS_aux.f90 @@ -1,26 +1,13 @@ module carbon_model_socs_aux - use carbon_model_to_core_arg_kit, only : year_min, year_max, nmonth + use environment_core, only: kd, lambd use const, only : yrs - use grid, only: date_c, date_fst, date_lst, dt + use grid, only: date_c + use config, only : station_name, station_opt + ! ---------------------------------------------- Environmental factors ------------------------------------------------ implicit none - ! Интерфейс -! ------------------------------------------------------------------------------------------------------------------- - ! ------- Station of observation ------ - character(len=10) :: station = 'Rostov' !< Станция наблюдения за климатом ! Нужно указать название - character(len=2) :: opt = '1' !< Имя варианта подачи удобрения ! Нужно указать номер - ! ------- Serve value ------- - integer :: mnc !< Колличество месяцев в расчете - integer :: mncX = 0 !< Номер месяца с начала работы программы !Как часть параметров для шага по времени - integer :: mnclot_fst(5) = (/1937,1939,1956,1968,1974/) !Как часть параметров для задания климатических данных - integer :: mnclot_lst(5) = (/2012,2012,2018,2018,2018/) !Как часть параметров для задания климатических данных - !DAO4 DAO3 VLDR TRGK ROST - integer :: station_n ! 1 2 3 4 5 !Как часть параметров для задания климатических данных -! ------------------------------------------------------------------------------------------------------------------- - ! Main part of SOCS - ! ------- Pools ------- real, dimension(:,:,:), pointer :: Catm !< атмосфера real, dimension(:,:,:), pointer :: Cveg !< растительность @@ -37,123 +24,42 @@ implicit none real :: rare = 0.45 !< Доля углерода переходящего в пул защищенного C1 real :: Cm = 12. !< Max кол-во органического углерода которое может быть защищено в почве real :: kirk = 4./yrs !< Коэф. скорости разложения C1 - real :: kd !< Коэф. перехода углерода из C2 в C1, десорбация, разрушение агрегатов - ! ------ Litterfall -------- - real, allocatable :: in_lambd(:) !< Поступление извне данных по поступлению углерода в почву, [kg/m**3 / year] - real :: lambd !< Поступление углерода в почву в интервал времени dt (lambdac) + !real :: kd !< Коэф. перехода углерода из C2 в C1, десорбация, разрушение агрегатов + contains subroutine carbon_model_init() -! ---- Part of enviromental carbon_model_init() ---- - mnc = (date_lst%y - date_fst%y)*nmonth + (date_lst%m - date_fst%m) - allocate(in_lambd(0:mnc)) - - open (unit = 1, file = 'initial_value/'//trim(station)//'_'//trim(opt)//'.txt', status='unknown') - read(1,*) in_lambd(1:mnc) - close (1) - - select case (station) - case('Rostov') - station_n = 5 - select case (opt) - case ('1') - kd = 0.012/yrs - case ('2') - kd = 0.007/yrs - case ('3') - kd = 0.007/yrs - case default - stop "check failed : unknown opt name" - end select - case('DAO3', 'DAO4') - if (station == 'DAO3') then - station_n = 2 - else - station_n = 1 - end if - select case (opt) - case ('1') - kd = 0.03/yrs - case ('2') - kd = 0.03/yrs - case ('3') - kd = 0.035/yrs - case ('4') - kd = 0.03/yrs - case default - stop "check failed : unknown opt name" - end select - case('VLDMR') - stop "check failed : can't work with this station name" - case('TRGK') - stop "check failed : can't work with this station name" - case default - stop "check failed : unknown station name" - end select - - in_lambd = in_lambd/dt - -! ---- Part of timesteps init ---- -! @todo учесть шаг по времени 1 час - if (dt == 86400) then ! Для шага по времени день: - mncX = (date_fst%y - mnclot_fst(station_n))*nmonth + date_lst%m - date_fst%m - in_lambd(0) = in_lambd(1) - else ! Для шага по времени месяц: - mncX = (date_fst%y - mnclot_fst(station_n))*nmonth + date_lst%m - end if - + end subroutine subroutine carbon_model_calc_at_timestep() - - + + end subroutine - - + + subroutine carbon_model_calc_at_cell(ii,jj) - + integer, intent(in) :: ii, jj -! ---- Part of environmental_calc_at_cell(ii,jj) ---- - integer i, N(12),j !< count - ! ------------- Переменные для распределения значений по всему временному интервалу - i = date_c%m ! Номер месяца - j = date_c%d ! Номер дня месяца - N(i) = date_c%days(i) ! Количество дней в месяце - !mncX ! Номер месяца с начала работы программы - ! ------------- - - if (dt == 86400) then ! Для шага по времени день: - lambd = in_lambd(mncX)/N(i) - if (j == (N(i)/2)) mncX = mncX + 1 - else ! Для шага по времени месяц: - lambd = in_lambd(mncX) - mncX = mncX + 1 - end if - - !print*, 'date_fst%y, ', date_fst%y - !print*, 'mnclot_fst(station_n), ', mnclot_fst(station_n) - !print*, 'date_lst%m, ', date_lst%m - !print*, 'station_n, ', station_n end subroutine - - + + subroutine carbon_model_calc_at_tile(ii,jj,nn) - + integer, intent(in) :: ii, jj, nn - + end subroutine - - - + + subroutine carbon_model_postprocessing() use grid, only : date_c - open(unit=500, file='results/'//trim(station)//'_'//trim(opt)//'.txt', status='unknown') + open(unit=500, file='results/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown') write(500,*) date_c%timestamp,';',C1(:,:,1),';',C2(:,:,1) end subroutine diff --git a/source/config.f90 b/source/config.f90 index a1a7776..f4d90da 100644 --- a/source/config.f90 +++ b/source/config.f90 @@ -6,11 +6,17 @@ module config public private :: config_check - + ! конфигурация character(len_default) :: carbon_model_type + + ! входные данные character(len_default) :: environment_data_type character(len_default) :: lsm_datafile character(len_default) :: lsm_dataformat + + character(len_default) :: station_name + character(len_default) :: station_opt +! пространственная сетка character(len_default) :: spatial_grid_mode real :: spatial_grid_res(2) character(len_default) :: spatial_sample_mode @@ -18,6 +24,7 @@ module config integer :: ich real :: point(2) real :: polygon(4) +! продолжительность расчета character(len_default) :: dt_mode integer :: dt character(len_default) :: datetime_init_mode @@ -26,6 +33,7 @@ module config integer :: ntime character(len_default) :: datetime_last integer :: UTC = 0 +! продвинутые настройки logical :: if_datafile_read_at_first logical :: if_out_yearly logical :: if_standard_print @@ -33,12 +41,16 @@ module config character(len_default) :: testing_log_mode integer :: nv_singlecolumn character(len_default) :: environment_model_type + + namelist /config_namelist/ & & carbon_model_type, & & environment_data_type, & & lsm_datafile, & & lsm_dataformat, & + &station_name, & !Добавил + &station_opt, & !Добавил & spatial_grid_mode, & & spatial_grid_res, & & spatial_sample_mode, & @@ -98,26 +110,26 @@ module config select case (lsm_dataformat) case('inmcm') - nc_namespace%lon = 'lon' - nc_namespace%lat = 'lat' - nc_namespace%lev = 'lev' - nc_namespace%time = 'time' - nc_namespace%us = 'UPBL' - nc_namespace%vs = 'VPBL' - nc_namespace%ta = 'TPBL' - nc_namespace%qa = 'QPBL' - nc_namespace%ps = 'PGR' - nc_namespace%p = 'PRECIP' - nc_namespace%rsw = 'TABL2_SDWSW' - nc_namespace%rlw = 'TABL2_SDWLW' + nc_namespace%lon = 'lon' + nc_namespace%lat = 'lat' + nc_namespace%lev = 'lev' + nc_namespace%time = 'time' + nc_namespace%us = 'UPBL' + nc_namespace%vs = 'VPBL' + nc_namespace%ta = 'TPBL' + nc_namespace%qa = 'QPBL' + nc_namespace%ps = 'PGR' + nc_namespace%p = 'PRECIP' + nc_namespace%rsw = 'TABL2_SDWSW' + nc_namespace%rlw = 'TABL2_SDWLW' nc_namespace%tsrf1 = 'TSRFOL1' - nc_namespace%tgr = 'TGROLD' - nc_namespace%s = 'SNOLD' - nc_namespace%ga = 'CTxAUS' + nc_namespace%tgr = 'TGROLD' + nc_namespace%s = 'SNOLD' + nc_namespace%ga = 'CTxAUS' nc_namespace%tsoil = 'TO' nc_namespace%wsoil = 'WLO' nc_namespace%isoil = 'WIO' - nc_namespace%mask = 'mask' + nc_namespace%mask = 'mask' ! см также перевод единиц в environment_data_lsm_offline.f90 end select @@ -136,12 +148,26 @@ module config end select select case(environment_data_type) - case('none', 'generator', 'lsm_offline') + case('none', 'generator', 'lsm_offline', 'station') case('lsm_online') stop "check failed : environment_data_type == 'lsm_online' is not supported yet" case default stop "check failed : unknown environment_data_type" end select + + select case(station_name) + case('Rostov', 'DAO3','DAO4', 'VLDMR', 'TRGK') + case default + stop "check failed : unknown station name" + end select + + select case(station_opt) + case('1', '2','3', '4') + case default + stop "check failed : unknown station_opt" + end select + + select case(lsm_dataformat) case('inmcm') @@ -188,8 +214,17 @@ module config end select ! тест на сочетания настроек - ! ----------------------------------------------------------------------------- + ! ----------------------------------------------------------------------------- + + if (carbon_model_type == 'socs' .and. environment_data_type == 'station' & + & .and. (station_name == 'VLDMR' .or. station_name == 'TRGK' )) then + stop "check failed : this carbon_model_type can't work with this station name" + endif + if ((station_name == 'Rostov' .or. station_name == 'VLDMR' .or. station_name == 'TRGK' ) .and. station_opt == '4') then + stop "check failed : this station name can't work with this station_opt" + endif + if (carbon_model_type == 'inmcm' .and. environment_data_type == 'none') then stop "check failed : carbon_model_type == 'inmcm' requires environment_data_type /= 'none'" endif diff --git a/source/const.f90 b/source/const.f90 index 2bc9f78..6285f4a 100644 --- a/source/const.f90 +++ b/source/const.f90 @@ -23,6 +23,7 @@ module const real, parameter :: r_earth = 6371000. !< радиус Земли, м real, parameter :: rho_w = 1. !< плотность воды [г/см3] real, parameter :: deg2rad = pi/180. + integer, parameter :: nmonth = 12 !< количество месяцев в году ! технические константы: integer, parameter :: miss_v = -999 !< значение по умолчанию diff --git a/source/environment/environment.f90 b/source/environment/environment.f90 index cf1a506..aa75d10 100644 --- a/source/environment/environment.f90 +++ b/source/environment/environment.f90 @@ -16,6 +16,13 @@ module environment & environment_model_calc_at_timestep_inmcm => environment_model_calc_at_timestep, & & environment_model_calc_at_cell_inmcm => environment_model_calc_at_cell, & & environment_model_calc_at_tile_inmcm => environment_model_calc_at_tile + use environment_data_station, only : environment_data_init_station => environment_data_init, & + & environment_data_calc_at_timestep_station => environment_data_calc_at_timestep, & + & environment_data_calc_at_cell_station => environment_data_calc_at_cell, & + & environment_data_calc_at_tile_station => environment_data_calc_at_tile + + + use config, only : environment_data_type, & & environment_model_type @@ -25,7 +32,7 @@ module environment public :: environment_init public :: environment_calc_at_timestep public :: environment_calc_at_cell - public :: environment_calc_at_tile + public :: environment_calc_at_tile interface subroutine interface_init() @@ -50,7 +57,6 @@ module environment procedure(interface_calc_at_cell), pointer :: environment_model_calc_at_cell procedure(interface_calc_at_tile), pointer :: environment_model_calc_at_tile - contains @@ -72,7 +78,14 @@ contains environment_data_calc_at_timestep => environment_data_calc_at_timestep_lsm_offline environment_data_calc_at_cell => environment_data_calc_at_cell_lsm_offline environment_data_calc_at_tile => environment_data_calc_at_tile_lsm_offline + + case('station') + environment_data_init => environment_data_init_station + environment_data_calc_at_timestep => environment_data_calc_at_timestep_station + environment_data_calc_at_cell => environment_data_calc_at_cell_station + environment_data_calc_at_tile => environment_data_calc_at_tile_station end select + select case(environment_model_type) case('inmcm') diff --git a/source/environment/environment_core.f90 b/source/environment/environment_core.f90 index 9e1a786..9eabd6a 100644 --- a/source/environment/environment_core.f90 +++ b/source/environment/environment_core.f90 @@ -8,6 +8,7 @@ module environment_core private public :: Temp, e, Rswd, Rlwd, p, pr, Wind, Tsrf, Tgr, snow, ra, Tsoil, Wsoil, Isoil + public :: sw, bettar, veg, kd, lambd public :: environment_core_init ! переменные @@ -27,7 +28,12 @@ module environment_core real, allocatable, dimension(:,:,:), target :: Tsoil !< Temerature at soil on depth levels, [Celsius] real, allocatable, dimension(:,:,:), target :: Wsoil !< Mass water content at soil on depth levels, [g/g] real, allocatable, dimension(:,:,:), target :: Isoil !< Mass ice content at soil on depth levels, [g/g] - + + real, allocatable, dimension(:,:), target :: sw !< Влажность увядания + real, allocatable, dimension(:,:), target :: bettar !< Fraction of soil respiration, [%] from 0 to 1 + real, allocatable, dimension(:,:), target :: veg !< Средняя концентрация растений в интервал времени dt, [%] + real, allocatable, dimension(:,:), target :: kd !< Коэф. перехода углерода из C2 в C1, десорбация, разрушение агрегатов + real, allocatable, dimension(:,:), target :: lambd !< Поступление углерода в почву в интервал времени dt (lambdac) contains @@ -54,7 +60,19 @@ contains allocate(Tsoil(i0:i1,j0:j1,ml)) allocate(Wsoil(i0:i1,j0:j1,ml)) allocate(Isoil(i0:i1,j0:j1,ml)) - + + allocate(sw(i0:i1,j0:j1)) + allocate(bettar(i0:i1,j0:j1)) + allocate(veg(i0:i1,j0:j1)) + allocate(kd(i0:i1,j0:j1)) + allocate(lambd(i0:i1,j0:j1)) + + sw = miss_v + bettar = miss_v + veg = miss_v + kd = miss_v + lambd = miss_v + Temp = miss_v e = miss_v Rswd = miss_v diff --git a/source/environment/environment_data_generator.f90 b/source/environment/environment_data_generator.f90 index a83a546..42c16d3 100644 --- a/source/environment/environment_data_generator.f90 +++ b/source/environment/environment_data_generator.f90 @@ -46,6 +46,7 @@ contains pr(i,j) = 0. !< precipitation ra(i,j) = 20. !< aerodynamical resistance snow(i,j) = 0. + ! veg(i,j) = 0.7 ! --- Tsrf(i,j) = 20. + 8.*cos(w_day*(hour-12.5)) !< temperature of surface Tgr(i,j) = Tsrf(i,j) diff --git a/source/environment/environment_data_station.f90 b/source/environment/environment_data_station.f90 new file mode 100644 index 0000000..2ad061f --- /dev/null +++ b/source/environment/environment_data_station.f90 @@ -0,0 +1,191 @@ + +module environment_data_station + + ! интерфейс + ! --------------------------------------------------------------------------------- + + use environment_core, only : Temp, e, Rlwd, Rswd, p, pr, Wind, ra, Tsrf, Tsoil, Wsoil, Isoil, snow, Tgr, & + & sw,bettar,kd,veg,lambd + use grid, only: date_c, date_fst, date_lst, dt, i0, i1, j0, j1 + use config, only : station_name, station_opt + use const, only : yrs, nmonth, pi + + implicit none + + private + public :: environment_data_init + public :: environment_data_calc_at_timestep + public :: environment_data_calc_at_cell + public :: environment_data_calc_at_tile + +! ------------------------------------------------------------------------------------------------------------------- + ! ------- Serve value ------- + integer :: mnc !< Количество месяцев в расчете + integer :: mncX = 0 !< Номер месяца с начала работы программы !Как часть параметров для шага по времени + + integer, parameter :: station_max = 5 !< Максимальное число станций + integer :: opt_n !< Номер способа подачи удобрения + integer, parameter :: opt_max = 4 !< Максимальное количество способов подачи удобрения + ! DAO4 DAO3 TRGK VLDR ROST + integer :: station_n ! 1 2 3 4 5 !< Номер станции наблюдения за климатом + integer :: mnclot_fst(station_max) = (/1937,1939,1956,1968,1974/) !< Даты начала сбора данных по климату в рамках станции + integer :: mnclot_lst(station_max) = (/2012,2012,2018,2018,2018/) !< Даты конца сбора данных по климату в рамках станции +! ------------------------------------------------------------------------------------------------------------------- +! real :: Temp !< Средняя температура воздуха в интервал времени dt, [Celsius] +! real :: Tsoil !< Средняя температура почвы в интервал времени dt, [Celsius] +! real :: Wsoil !< Влажность почвы в интервал времени dt, [dim] + ! ------ Determined externally + real, allocatable :: in_temp(:) !< Поступление извне данных cредней температуры воздуха в месяц, [Celsius] + real, allocatable :: in_veg(:) !< Поступление извне данных средней концентрации растений в месяц, [%] + real, allocatable :: in_wsoil(:) !< Поступление извне данных влажности почвы, [dim] + real, allocatable :: in_lambd(:) !< Поступление извне данных по поступлению углерода в почву, [kg/m**3 / year] + + ! station = DAO4 DAO3 TRGK VLDR ROST + real :: sw_in(station_max) = (/0.120, 0.120, 0.120, 0.120, 0.146/) + real :: bettar_in(station_max) = (/0.25, 0.25, 0.07, 0.08, 0.41/) + real :: kd_in(station_max,opt_max) + ! opt = 1 2 3 4 + data kd_in(01,:) /0.03, 0.03, 0.035, 0.03/ + data kd_in(02,:) /0.03, 0.03, 0.035, 0.03/ + data kd_in(03,:) /0., 0., 0., 0./ + data kd_in(04,:) /0., 0., 0., 0./ + data kd_in(05,:) /0.012, 0.007, 0.007, 0./ + +contains + + ! внешние процедуры + ! --------------------------------------------------------------------------------- + + subroutine environment_data_init() + +! ---- Part of enviromental carbon_model_init() ---- + mnc = (date_lst%y - date_fst%y)*nmonth + (date_lst%m - date_fst%m) + allocate(in_lambd(0:mnc)) + allocate(in_temp(0:mnc)) + allocate(in_veg(0:mnc)) + allocate(in_wsoil(0:mnc)) + + print*, mnc, station_name + + kd_in = kd_in/yrs + + read(station_opt,*) opt_n + + select case (station_name) + case('Rostov') + station_n = 5 + case('DAO3') + station_n = 2 + case('DAO4') + station_n = 1 + case('VLDMR') + station_n = 4 + case('TRGK') + station_n = 3 + end select + + open (unit = 1, file = 'initial_value/'//trim(station_name)//'_'//trim(station_opt)//'.txt', status='unknown') + read(1,*) in_lambd(1:mnc) + close (1) + + open (unit = 1, file = 'initial_value/'//trim(station_name)//'_Temp.txt', status='unknown') + read(1,*) in_temp(1:mnc) + close (1) + + open (unit = 1, file = 'initial_value/'//trim(station_name)//'_MOI.txt', status='unknown') + read(1,*) in_wsoil(1:mnc) + close (1) + + open (unit = 1, file = 'initial_value/'//trim(station_name)//'_Veg.txt', status='unknown') + read(1,*) in_veg(1:mnc) + close (1) + + in_lambd = in_lambd/dt + + sw(:,:) = sw_in(station_n) + bettar(:,:) = bettar_in(station_n) + kd(:,:) = kd_in(station_n, opt_n) + +! ---- Part of timesteps init ---- +! @todo учесть шаг по времени 1 час + if (dt == 86400) then ! Для шага по времени день: + mncX = (date_fst%y - mnclot_fst(station_n))*nmonth + date_lst%m - date_fst%m + in_temp(0) = in_temp(1) + in_veg(0) = in_veg(1) + in_wsoil(0) = in_wsoil(1) + in_lambd(0) = in_lambd(1) + else ! Для шага по времени месяц: + mncX = (date_fst%y - mnclot_fst(station_n))*nmonth + date_lst%m + end if + +! Набор переменных которые должны быть, но не данны через качественные внешние данные + p(:,:) = 980. !< pressure + Rlwd(:,:) = 350. !< longwave radiation + pr(:,:) = 0. !< precipitation + ra(:,:) = 20. !< aerodynamical resistance + snow(:,:) = 0. !< snow cover + Isoil(:,:,:) = 0. !< part of ice in the soil + + end subroutine + + + subroutine environment_data_calc_at_timestep() + ! --------------------------------------- + use const, only : w_day, miss_v + use grid, only : z, ms, ml + + integer nom,nod, N(12) + integer :: i, j !< count + real :: hour + ! ------------- Переменные для распределения значений по всему временному интервалу + nom = date_c%m ! Номер месяца + nod = date_c%d ! Номер дня месяца + N(nom) = date_c%days(nom) ! Количество дней в месяце + !mncX ! Номер месяца с начала работы программы + ! ------------- + if (dt == 86400) then ! Для шага по времени день: + Temp(:,:) = in_temp(mncX)*sin(pi*(nod-1)/(N(nom)-1)) + in_temp(mncX+1)*(1-sin(pi*(nod-1)/(N(nom)-1))) + Wsoil(:,:,:) = in_wsoil(mncX)*sin(pi*(nod-1)/(N(nom)-1)) + in_wsoil(mncX+1)*(1-sin(pi*(nod-1)/(N(nom)-1))) + veg(:,:) = in_veg(mncX) + lambd(:,:) = in_lambd(mncX)/N(nom) + if (nod == (N(nom)/2)) mncX = mncX + 1 + else ! Для шага по времени месяц: + Temp(:,:) = in_temp(mncX) + Wsoil(:,:,:) = in_wsoil(mncX) + veg(:,:) = in_veg(mncX) + lambd(:,:) = in_lambd(mncX) + mncX = mncX + 1 + end if + + do i = 1, ml + Tsoil(:,:,i) = Temp(:,:) ! температура почвы равна температуре воздуха + end do + + hour = date_c%h +! Набор переменных которые должны быть, но не данны через качественные внешние данные + do i = i0, i1 + do j = j0, j1 + Wind(i,j) = 2 + 1.*cos(w_day*(hour-13)) !< wind + e(i,j) = 15. + 2.*cos(2*w_day*(hour-7)) !< humidity + Rswd(i,j) = max(0.,800.*cos(w_day*(hour-12))) !< shortwave radiation + ! --- + Tsrf(i,j) = 20. + 8.*cos(w_day*(hour-12.5)) !< temperature of surface + Tgr(i,j) = Tsrf(i,j) + end do + end do + + end subroutine + + + subroutine environment_data_calc_at_cell(ii,jj) + + integer, intent(in) :: ii, jj + + end subroutine + + subroutine environment_data_calc_at_tile(ii,jj,nn) + ! --------------------------------------- + integer, intent(in) :: ii, jj, nn + end subroutine + +end module diff --git a/source/environment/environment_model_INMCM.f90 b/source/environment/environment_model_INMCM.f90 index 097f747..ffa0e25 100644 --- a/source/environment/environment_model_INMCM.f90 +++ b/source/environment/environment_model_INMCM.f90 @@ -198,16 +198,20 @@ contains deallocate(work1i) close(1) - call special_soilpar_read(path_inmcm_data//'soil_par_average.nc', 'theta_s', sPORu, fill_v = 0.4) - call special_soilpar_read(path_inmcm_data//'CH_par_average.nc', 'psi_s', sPSIMAXu, fill_v = -17.) - call special_soilpar_read(path_inmcm_data//'CH_par_average.nc', 'lambda', sBHu, fill_v = 0.15) + ! call special_soilpar_read(path_inmcm_data//'soil_par_average.nc', 'theta_s', sPORu, fill_v = 0.4) + ! call special_soilpar_read(path_inmcm_data//'CH_par_average.nc', 'psi_s', sPSIMAXu, fill_v = -17.) + ! call special_soilpar_read(path_inmcm_data//'CH_par_average.nc', 'lambda', sBHu, fill_v = 0.15) + sPORu = 0.4 + sPSIMAXu = -17. + sBHu = 0.15 sBHu(:,:) = 1./sBHu(:,:) allocate(MCD15A2H_LAI(720,360,nv2,12)) - call nc_errhand( nf90_open(path_inmcm_data//'mcd15a2h_to_inmcm.nc', nf90_nowrite, ncid) ) - call nc_errhand( nf90_inq_varid(ncid, 'lai', varid) ) - call nc_errhand( nf90_get_var(ncid, varid, MCD15A2H_LAI(:,:,:,:), (/1,1,1,1/), (/720,360,nv2,12/)) ) - call nc_errhand( nf90_close(ncid) ) + ! call nc_errhand( nf90_open(path_inmcm_data//'mcd15a2h_to_inmcm.nc', nf90_nowrite, ncid) ) + ! call nc_errhand( nf90_inq_varid(ncid, 'lai', varid) ) + ! call nc_errhand( nf90_get_var(ncid, varid, MCD15A2H_LAI(:,:,:,:), (/1,1,1,1/), (/720,360,nv2,12/)) ) + ! call nc_errhand( nf90_close(ncid) ) + MCD15A2H_LAI = 2. if (nv_singlecolumn /= miss_v) then vegg_nv1(:,:,:) = 0. diff --git a/ui1_config.nml b/ui1_config.nml index ddc752d..fd094e2 100644 --- a/ui1_config.nml +++ b/ui1_config.nml @@ -14,7 +14,7 @@ !> Углеродная модель: carbon_model_type = 'socs' ! - ! 'inmcm' - модель INMCM [1,2] + ! 'inmcm' - модель inmcm [1,2] ! 'socs' - модель SOCS [3] ! 'rothc' - модель ROTHC ! 'other' - шаблон @@ -23,13 +23,28 @@ ! ----------------------------------------------------------------- !> Входные данные об окружающей среде: - environment_data_type = 'generator' + environment_data_type = 'station' ! ! 'none' - не требуются ! 'generator' - генерация суточных синусоид + ! 'station' - выбор готовых данных, со станций наблюдения ! 'lsm_offline' - netcdf-файл с выходными данными модели деятельного слоя ! 'lsm_online' - онлайн-каплинг с моделью деятельного слоя (@todo пока не поддерживается) + ! дополнительно для environment_data_type = 'station': + station_name = 'DAO4' + ! 'Rostov' - Станция ФАНЦ + ! 'DAO3' - Станция в Долгопрудном 1 + ! 'DAO4' - Станция в Долгопрудном 2 + ! 'VLDMR' - Данные Владимир + ! 'TRGK' - Данные Торжок + + station_opt = '4' + ! '1' - контрольный случай, без дополнительных удобрений + ! '2' - 2 способ подачи удобрений + ! '3' - 3 способ подачи удобрений + ! '4' - 4 способ подачи удобрений + ! дополнительно для environment_data_type = 'lsm_offline': !> Путь к файлу: @@ -64,8 +79,8 @@ ich = 46967 ! Федоровское, 05x05 ! 'point' - ячейка по координатам [lon, lat] ! point(:) = 32.75, 56.25 ! Федоровское - point(:) = 39.888735, 47.364103 ! ФАНЦ (Ростов) - ! point(:) = 37.535197, 55.942416 ! ДАОС (Долгопрудный) + ! point(:) = 39.888735, 47.364103 ! ФАНЦ (Ростов) + point(:) = 37.535197, 55.942416 ! ДАОС (Долгопрудный) ! 'polygon' - область по координатам [lon_west, lon_east, lat_south, lat_north] polygon(:) = 26.25, 69.75, 50.75, 69.75 ! лесная зона ЕТР ! (для lsm_offline) 'all' - вся область nectdf-файла @@ -88,7 +103,7 @@ datetime_init_mode = 'manual' !auto ! ! 'manual' - задать вручную [yyyy-mm-dd hh:mm:ss] 1937 1939 1956 1968 1974 - datetime_init = '1974-01-01 00:00:00' + datetime_init = '1937-01-01 00:00:00' ! (для lsm_offline) 'auto' - получить из входного файла @@ -98,7 +113,7 @@ ! 'ntime' - указанное число шагов ntime = 1000 ! 'datetime_last' - до достижения указанной даты 2012 2018 2975 - datetime_last = '2018-01-01 00:00:00' + datetime_last = '2012-01-01 00:00:00' ! (для lsm_offline) 'auto' - до конца входного файла ! продвинутые настройки -- GitLab