Closed HCJung-jbnu closed 4 months ago
ESM4-PDAF 예시
! *** All other filters ***
! *** SEIK, LSEIK, ETKF, LETKF, ESTKF, LESTKF, NETF, LNETF ***
! *** and generation of synthetic observations (GENOBS) ***
filter_param_i(1) = dim_state_p_ocn ! State dimension
filter_param_i(2) = dim_ens_ocn ! Size of ensemble <- **앙상블 갯수**
filter_param_r(1) = forget_ocn ! Forgetting factor
! Optional parameters; you need to re-set the number of parameters if you use them
filter_param_i(3) = dim_lag_ocn ! Size of lag in smoother
!!!------
filter_param_i(4) = incremental_ocn ! Whether to perform incremental analysis
filter_param_i(5) = type_forget_ocn ! Type of forgetting factor
filter_param_i(6) = type_trans_ocn ! Type of ensemble transformation
filter_param_i(7) = type_sqrt_ocn ! Type of transform square-root (SEIK-sub4/ESTKF)
! filter_param_i(8) = observe_ens_ocn ! Apply H to ensemble or ensmeble mean for residual (ESTKF/ETKF/SEIK)
CALL PDAF_init(filtertype_ocn, subtype_ocn, step_null_ocn, &
filter_param_i, 7, &
filter_param_r, 2, &
COMM_model_ocn, COMM_filter_ocn, COMM_couple_ocn, &
task_id_ocn, n_modeltasks_ocn, filterpe_ocn, init_seik_pdaf_ocn, &
screen_ocn, status_pdaf)
....
....
....
! ******************************'***
! *** Prepare ensemble forecasts ***
! ******************************'***
CALL PDAF_get_state(steps, timenow, doexit, next_observation_pdaf_ocn, &
distribute_state_pdaf_ocn, prepoststep_ens_pdaf_ocn, status_pdaf, 0, 0)
**dim_state_p_ocn :: 각 PE의 모델 dimension 갯수(i * j * k) 해당 값은 0가 되면 안됨 최소 값=4**
**dim_ens_ocn :: 앙상블 갯수**
**task_id_ocn :: 앙상블 rank = task number**
**n_modeltasks_ocn :: 앙상블 갯수**
**filterpe_ocn :: 첫 번째 앙상블=task 인지 여부, filter task 확인 true/false**
**init_seik_pdaf_ocn :: pdaf의 eofcovar 앙상블 생성 활용시 필요, 그렇지 않으면 사용하지 않음**
사용 안할 시 해당 코드는 껍데기만 만들고 안의 내용은 삭제
SUBROUTINE init_seik_pdaf_ocn(filtertype, dim_p, dim_ens, state_p, Uinv, & ens_p, flag)
**마지막에 PDAF_get_state 호출**
CALL PDAF_get_state(steps, timenow, doexit, next_observation_pdaf_atm, & distribute_state_pdaf_atm, prepoststep_ens_pdaf_atm, status_pdaf,0, 0)
CALL PDAF_get_state(steps, time, doexit, U_next_observation, U_distribute_state, &
U_prepoststep, outflag, assim_var_flag, varnum) !!hcjuing 20210712
....
assim_var_flag = 1 !!Temp
CALL assimilate_pdaf_ocn(nstep_ocn) !!<<--- main subroutine of data assimilation 1
....
assim_var_flag = 2 !!SALT
CALL assimilate_pdaf_ocn(nstep_ocn) !!<<--- main subroutine of data assimilation 2
....
assim_var_flag = 3 !!BIO
CALL assimilate_pdaf_ocn(nstep_ocn) !!<<--- main subroutine of data assimilation 3
....
....
IF ( assim_var_flag .EQ. 1 .OR. assim_var_flag .EQ. 2 ) THEN !!temp/salt
CALL PDAF_assimilate_lestkf(collect_state_pdaf_ocn, distribute_state_pdaf_ocn, init_dim_obs_f_pdaf_ocn, &
obs_op_f_pdaf_ocn, init_obs_f_pdaf_ocn, init_obs_l_pdaf_ocn, &
prepoststep_ens_pdaf_ocn, &
prodRinvA_l_pdaf_ocn, init_n_domains_pdaf_ocn, init_dim_l_pdaf_ocn, &
init_dim_obs_l_pdaf_ocn, g2l_state_pdaf_ocn, l2g_state_pdaf_ocn, &
g2l_obs_pdaf_ocn, init_obsvar_pdaf_ocn, init_obsvar_l_pdaf_ocn, &
next_observation_pdaf_ocn, status_pdaf, assim_var_flag, num_3dvar_ocn, dim_state_p_ocn)
! next_observation_pdaf_ocn, status_pdaf, assim_var_flag, 2) !!chl test
ELSE IF ( assim_var_flag .EQ. 3 ) THEN !chl -> biovar
CALL PDAF_assimilate_lestkf(collect_state_pdaf_ocn, distribute_state_pdaf_ocn, init_dim_obs_f_pdaf_ocn, &
obs_op_f_pdaf_ocn, init_obs_f_pdaf_ocn, init_obs_l_pdaf_ocn, &
prepoststep_ens_pdaf_ocn, &
prodRinvA_l_pdaf_ocn, init_n_domains_pdaf_ocn, init_dim_l_pdaf_ocn, &
init_dim_obs_l_pdaf_ocn, g2l_state_pdaf_ocn, l2g_state_pdaf_ocn, &
g2l_obs_pdaf_ocn, init_obsvar_pdaf_ocn, init_obsvar_l_pdaf_ocn, &
next_observation_pdaf_ocn, status_pdaf, assim_var_flag, num_3dvar_ocn, dim_state_p_bio_ocn)
! next_observation_pdaf_ocn, status_pdaf, assim_var_flag, 2) !!chl test
END IF
....
Before the analysis step is called the following is executed:
increment가 적용될 model field를 가져옴
SUBROUTINE collect_state(dim_p, state_p)
INTEGER, INTENT(in) :: dim_p ! State dimension for PE-local model sub-domain
REAL, INTENT(inout) :: state_p(dim_p) ! State vector for PE-local model sub-domain
When the ensemble integration of the forecast is completed, the analysis step is executed. Before the loop over all local analysis domains, the following routines are executed:
rmse 정보 log 에 보여줌 screen_ocn 이 1보다 커야 로그에 보여줌
SUBROUTINE prepoststep(step, dim_p, dim_ens, dim_ens_p, dim_obs_p, &
state_p, Uinv, ens_p, flag)
INTEGER, INTENT(in) :: step ! Current time step
! (When the routine is called before the analysis -step is provided.)
INTEGER, INTENT(in) :: dim_p ! PE-local state dimension
INTEGER, INTENT(in) :: dim_ens ! Size of state ensemble
INTEGER, INTENT(in) :: dim_ens_p ! PE-local size of ensemble
INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector
REAL, INTENT(inout) :: state_p(dim_p) ! PE-local forecast/analysis state
! The array 'state_p' is not generally not initialized in the case of SEIK/EnKF/ETKF/ESTKF.
! It can be used freely in this routine.
REAL, INTENT(inout) :: Uinv(dim_ens-1, dim_ens-1) ! Inverse of matrix U
REAL, INTENT(inout) :: ens_p(dim_p, dim_ens) ! PE-local state ensemble
INTEGER, INTENT(in) :: flag ! PDAF status flag
analysis domain(각 pe의 모델 격자 개수) 의 개수 전달, column 전체를 자료동화 하는 경우 모델의 수평 격자 개수가 해당됨.
SUBROUTINE init_n_domains(step, n_domains_p)
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(out) :: n_domains_p ! Number of analysis domains for local model sub-domain
AM4의 경우 NCEP2를 모델 격자로 수평내삽하여 모든 pe에 내삽된 global obs와 lat/lon 정보를 공유함.
SUBROUTINE init_dim_obs_f(step, dim_obs_f)
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(out) :: dim_obs_f ! Dimension of full observation vector
홈페이지 공식 설명에는 dim_ens 만큼 call이 된다고 하지만, 확인을 해 보면 filterpe 즉 1번 task(=ens)에서 모든 앙상블 데이터를 요구함, 그러므로 1번 task에 모든 global model field를 전달해 줘야 함.
SUBROUTINE obs_op_f(step, dim_p, dim_obs_f, state_p, m_state_f)
INTEGER, INTENT(in) :: step ! Currrent time step
INTEGER, INTENT(in) :: dim_p ! PE-local dimension of state
INTEGER, INTENT(in) :: dim_obs_f ! Dimension of the full observed state
REAL, INTENT(in) :: state_p(dim_p) ! PE-local model state
REAL, INTENT(out) :: m_state_f(dim_obs_f) ! Full observed state
AM4의 경우 timestep 날짜에 스캔된, 수평/연직 내삽된 ncep2데이터를 전체 tile 자료를 전달함.
SUBROUTINE init_obs_f(step, dim_obs_f, observation_f)
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector
REAL, INTENT(out) :: observation_f(dim_obs_f) ! Full observation vector
MOM5, AM4 모두 U_init_dim_obs_f 에서 계산하여 전달
SUBROUTINE init_obsvar(step, dim_obs_f, obs_f, meanvar_f)
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: dim_obs_f ! Full dimension of observation vector
REAL, INTENT(in) :: obs_f(dim_obs_f) ! Full observation vector
REAL, INTENT(out) :: meanvar_f ! Mean observation error variance
AM4의 경우 U_init_dim_obs_f가 아닌 해당 서브루틴에서 kdtree 활용해서 반경 내 관측자료 스캔한 정보 저장.
SUBROUTINE init_dim_l(step, domain_p, dim_l)
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: domain_p ! Current local analysis domain
INTEGER, INTENT(out) :: dim_l ! Local state dimension
각 domain 별로 scan 된 관측자료 갯수 전달
SUBROUTINE init_dim_obs_l(domain_p, step, dim_obs_f, dim_obs_l)
INTEGER, INTENT(in) :: domain_p ! Current local analysis domain
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: dim_obs_f ! Full dimension of observation vector
INTEGER, INTENT(out) :: dim_obs_l ! Local dimension of observation vector
increment가 적용 될 전체 모델 값에서 해당 domain에 해당되는 값 전달
SUBROUTINE g2l_state(step, domain_p, dim_p, state_p, dim_l, state_l)
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: domain_p ! Current local analysis domain
INTEGER, INTENT(in) :: dim_p ! State dimension for model sub-domain
INTEGER, INTENT(in) :: dim_l ! Local state dimension
REAL, INTENT(in) :: state_p(dim_p) ! State vector for model sub-domain
REAL, INTENT(out) :: state_l(dim_l) ! State vector on local analysis domain
operator에서 관측 위치로 내삽 된 모델 값에서 각 domain에 scan된 위치의 값을 전달
SUBROUTINE g2l_obs(domain_p, step, dim_obs_f, dim_obs_l, mstate_f, mstate_l)
INTEGER, INTENT(in) :: domain_p ! Current local analysis domain
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector for model sub-domain
INTEGER, INTENT(in) :: dim_obs_l ! Dimension of observation vector for local analysis domain
REAL, INTENT(in) :: mstate_f(dim_obs_f) ! Full observation vector for model sub-domain
REAL, INTENT(out) :: mstate_l(dim_obs_l) ! Observation vector for local analysis domain
각 domain에 scan된 위치의 관측 자료 값을 전달
SUBROUTINE init_obs_l(domain_p, step, dim_obs_l, observation_l)
INTEGER, INTENT(in) :: domain_p ! Current local analysis domain
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: dim_obs_l ! Local dimension of observation vector
REAL, INTENT(out) :: observation_l(dim_obs_l) ! Local observation vector
SUBROUTINE prodRinvA_l(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l)
INTEGER, INTENT(in) :: domain_p ! Current local analysis domain
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: dim_obs_l ! Dimension of local observation vector
INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix
REAL, INTENT(in) :: obs_l(dim_obs_l) ! Local vector of observations
REAL, INTENT(inout) :: A_l(dim_obs_l, rank) ! Input matrix from analysis routine
REAL, INTENT(out) :: C_l(dim_obs_l, rank) ! Output matrix
increment가 적용된 domain별 local 모델 값을 전체 field에 전달
SUBROUTINE l2g_state(step, domain_p, dim_l, state_l, dim_p, state_p)
INTEGER, INTENT(in) :: step ! Current time step
INTEGER, INTENT(in) :: domain_p ! Current local analysis domain
INTEGER, INTENT(in) :: dim_p ! State dimension for model sub-domain
INTEGER, INTENT(in) :: dim_l ! Local state dimension
REAL, INTENT(in) :: state_p(dim_p) ! State vector for model sub-domain
REAL, INTENT(out) :: state_l(dim_l) ! State vector on local analysis domain
자료동화 통계값 계산
SUBROUTINE prepoststep(step, dim_p, dim_ens, dim_ens_p, dim_obs_p, &
state_p, Uinv, ens_p, flag)
INTEGER, INTENT(in) :: step ! Current time step
! (When the routine is called before the analysis -step is provided.)
INTEGER, INTENT(in) :: dim_p ! PE-local state dimension
INTEGER, INTENT(in) :: dim_ens ! Size of state ensemble
INTEGER, INTENT(in) :: dim_ens_p ! PE-local size of ensemble
INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector
REAL, INTENT(inout) :: state_p(dim_p) ! PE-local forecast/analysis state
! The array 'state_p' is not generally not initialized in the case of SEIK/EnKF/ETKF/ESTKF.
! It can be used freely in this routine.
REAL, INTENT(inout) :: Uinv(dim_ens-1, dim_ens-1) ! Inverse of matrix U
REAL, INTENT(inout) :: ens_p(dim_p, dim_ens) ! PE-local state ensemble
INTEGER, INTENT(in) :: flag ! PDAF status flag
SUBROUTINE next_observation_pdaf_ocn(stepnow, nsteps, doexit, time)
! !ARGUMENTS:
INTEGER, INTENT(in) :: stepnow ! Number of the current time step
INTEGER, INTENT(out) :: nsteps ! Number of time steps until next obs
INTEGER, INTENT(out) :: doexit ! Whether to exit forecasting (1 for exit)
REAL, INTENT(out) :: time ! Current model (physical) time
COMM_model, COMM_coupler의 rank는 0부터 시작하지만 COMM_couple의 rank는 1부터 시작함n_model_task=1 일때,
n_model_task=4 일때,