Milija Zupanski CIRA Colorado State University (E-mail: ZupanskiM@CIRA.colostate)

1. ARMY Research Lab and CIRA/CSU Collaborationon Assessing the value of ensemble based data assimilation/forecasting system in Army microscale meteorological applications Milija Zupanski CIRA Colorado State University (E-mail: ZupanskiM@CIRA.colostate.edu)

2. ARMY Research Lab and CIRA/CSU Collaboration on Assessing the value of ensemble based data assimilation/forecasting system in Army microscale meteorological applications White paper overview Proposed work: i- Structure and correction of model error (bias) at the microscale - understand the errors and try to correct them ii- Impact of mesoscale model forecast on initialization of microscale model - transfer of uncertainty information between scales iii-Optimal use of computer resources, including HPC, with new ensemble data assimilation methodology - HPC performance of ensemble data assimilation

3. Extended Kalman Filter Direct solution: Equivalent to solving a quadratic minimization problem: Subject to Pf- forecast error covariance Pa- analysis error covariance R - observation error covariance K - nonlinear observation operator K - linearized observation operator (Jacobian)

4. Maximum Likelihood Ensemble Filter (MLEF) • Maximize posterior probability to obtain the conditional mode • Equivalent to minimizing general non-quadratic cost function: No assumption regarding the (non)linearity of observation operator K ! • Optimal state + uncertainty • Embedded in optimization algorithm (2-3 minimization iterations!) • Analysis perturbations only (no observation perturbations) • Math library subroutines for matrix operations (BLAS, LAPACK) • Mathematically consistent feed-back between data assimilation and ensemble forecasting

5. Maximum Likelihood Ensemble Filter (MLEF)Forecast error covariance x – initial conditions r – model error (including bias) pi – i-th column of the square-root analysis erorr covariance matrix bi – i-th column of the square-root forecast erorr covariance matrix

6. Maximum Likelihood Ensemble Filter (MLEF)Change of variable (Hessian preconditioing) Inversion via an eigenvalue decomposition of A:

7. Maximum Likelihood Ensemble Filter (MLEF)Gradient calculation Use: Need only forecast difference at observation location !

8. Maximum Likelihood Ensemble Filter (MLEF)Analysis error covariance • By-product of minimization algorithm L-BFGS: C-G:

9. Ensemble Data Assimilation (EnsDA) Algorithm Background: a- NCEP operational algorithm structure - driven by UNIX scripts b- Weather Research and Forecasting (WRF) 3DVAR programming - definition of control, first-guess, as ‘type’ variables c- 4DVAR experience (NCEP Eta-model, CIRA/CSU RAMS model) - flow of data, user-friendly approach

10. Ensemble Data Assimilation (EnsDA) Algorithm prep_ensda.sh WARM start: Copy files from previously completed cycle COLD start: Run randomly-perturbed ensemble forecasts to initialize fcst err cov cycle_ensda.sh icycle < N_cycles_max fcsterr_cov.sh - Prepare first-guess (background) vector - Prepare forecast error covariance (from ensembles) prep_obs.sh Given ‘OBSTYPE’ and ‘delobs’, select and copy available obs files assimilation.sh Iterative minimization of cost function, save current cycle output post.sh (post-processing)

11. Ensemble Data Assimilation (EnsDA) Algorithm assimilation.sh script assimilation.sh iter < ioutmax forward.sh: Transformation from model space to observation space Hessian preconditioning (only for iter=1) Gradient calculation (ensembles) Cost function calculation (diagnostic) Minimization (ensemble subspace) Step-length (line-search) Control variable update (transformation from ensemble subspace to model (physical) space) - Analysis error covariance calculation - Save current cycle output files

12. Ensemble Data Assimilation (EnsDA) Algorithm Organization chart ../ensda/include: Include files (control variable list, first-guess list) ../ensda/makefiles: Makefile_ensda, Makefile_’MODEL’ ../ensda/runscripts: Run EnsDA program (runensda, envir.sh) ../ensda/scripts: Scripts used to run the EnsDA program ../ensda/src: Source programs ../’MODEL’/src: Model source programs (general) ../’MODEL’/src_ensda: Model source programs (changes due to EnsDA) ../exec/ensda: Executable directory ../’Obsdir’: Observations directory ../ensda/work: Working (temporary) directory ../hold.’print’/current: General EnsDA output ../hold.’print’/cycle1, cycle2, . . . : Cycle specific EnsDA output

13. Ensemble Data Assimilation (EnsDA) Algorithm Options – envir.sh MODEL: WRF,UW-NMS,PU-NTU OBS: OBSTYPE, delobs,OBSFLAG STAT: Mode (max likelihood), Mean (ensemble mean) DA_TYPE: Filter, Smoother COVAR: Localized, Non-localized forecast error covariance N_CYCLES: Number of data assimilation cycles CYCLE_INTERVAL: Time interval between cycles ENS_SIZE: Number of ensemble members START: Cold (from scratch), Warm (from previous cycle) IOUTMAX: Number of minimization iterations MINIM_ALG: Minimization algorithm (C-G, L-BFGS) MPI_RUN: Parallel run (YES, NO) nCPU_ensfcst: No. of CPUs for ensemble fcst nCPU_da: No. of CPUs for data assimilation

14. Ensemble Data Assimilation (EnsDA) Algorithm Control variables Optional control variable components: - initial conditions - model bias - empirical parameters ../include/Cntrl_vrbl_list.h !--------------------------------------- max_num_of_cntrl_vrbls = 5 if(.not.allocated (cvar_list)) then allocate (cvar_list(1:max_num_of_cntrl_vrbls)) end if !--------------------------------------- ncv= 1 cvar_list(ncv)%ndim = 3 cvar_list(ncv)%start_index(1) = 1 cvar_list(ncv)%start_index(2) = 1 cvar_list(ncv)%start_index(3) = 1 cvar_list(ncv)%start_index(4) = 1 cvar_list(ncv)%end_index(1) = NNXP(1) cvar_list(ncv)%end_index(2) = NNYP(1) cvar_list(ncv)%end_index(3) = NNZP(1) cvar_list(ncv)%end_index(4) = 1 cvar_list(ncv)%name = 'T ' cvar_list(ncv)%stddev = 1.5 cvar_list(ncv)%description = 'Temperature ‘ !--------------------------------------- ../namelists/Cntrl_vrbl_’MODEL’.h 6 p T F F 0 t T F T 1 q T F T 3 u T F F 0 v T F F 0 param1 F T F 0 param2 F T F 0 !-- first line number is the number of control variables defined --! !-- vrbl name, ic flag, param flag, bias flag, number of biases ----! !--- recall that vrbl name has 9 characters !!!!! !--- All inputs have to be separated by at least one blank space !!!

15. Ensemble Data Assimilation (EnsDA) Algorithm Observations # Precip obs export iobs=02 #---------- specified by user --------- export OBSTYPE=precip export OBSFLAG=NOT let "delobs=1*60" ; export delobs #----------- end user ------------------ -------------------------- OBSTYPE – algorithm recognizable obs type delobs – time interval between two successive observation files OBSFLAG – Specifies if the OBSTYPE obs should be assimilated (YES, NOT) interface_’MODEL’_’OBSTYPE’.sh (model produced first-guess file => obs operators compatible first-guess file) runobs_’OBSTYPE’.sh (forward-only observation operator: first-guess => innovation vector R-1/2[y-K(x)] )

16. Ensemble Data Assimilation (EnsDA) Algorithm Model convert_vrbls_from_1d.sh run_’MODEL’.sh convert_vrbls_to_1d.sh - Data assimilation is defined using 1-d control variable - Model is defined using its own 1-d, 2-d, or 3-d, or 4-d variable - Need ‘convert’ as an interface between the 1-d and model specific variables ../’MODEL’/src: (‘MODEL’_fcst.F) read namelists: ‘MODEL_DIMENSION’,’LENGTH’,’MODEL_CNSTS’ call: ‘read/write_model_ic’,’read/write_model_bias’ ,’read/write_model_param’ ../’MODEL’/src_ensda: ‘MODEL’_get_init_fcst.F –cold start only, get model_dimension file ‘read/write_model_ic.F’,’read/write_model_bias.F’ ,’read/write_model_param.F’

17. Ensemble Data Assimilation (EnsDA) Algorithm Including model bias DOn=1,N_time !! Main time-loop . . . if(num_bias.gt.0) then nerr=N_time/max(1,num_bias) nbias=0 if((n.ne.N_time).and.(mod(n,nerr).eq.0.or.n.eq.1)) then nbias=nbias+1 call read_model_bias( . . . ,u_bias,v_bias,t_bias, . . . ,nbias,num_bias) endif endif . . . if(num_bias.gt.0) then u_phi(:,:,:)=aa*u_phi(:,:,:)+(1.-aa)*u_bias(:,:,:) v_phi(:,:,:)=aa*v_phi(:,:,:)+(1.-aa)*v_bias(:,:,:) t_phi(:,:,:)=aa*t_phi(:,:,:)+(1.-aa)*t_bias(:,:,:) . u(:,:,:)=u(:,:,:)+u_phi(:,:,:) v(:,:,:)=v(:,:,:)+v_phi(:,:,:) t(:,:,:)=t(:,:,:)+t_phi(:,:,:) endif . . . END DO !! End main time loop

18. Ensemble Data Assimilation (EnsDA) Algorithm MPI !-------------------------------------------------------------------- #ifdef MPI_USE include "mpif.h" #endif . . . . . . . . . !==============start calculation=================== #ifdef MPI_USE CALL MPI_INIT(IERR) IF (IERR .NE. MPI_SUCCESS) STOP 'FAILED TO INIT MPI' CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPROC, IERR) CALL MPI_COMM_RANK(MPI_COMM_WORLD, MPIRANK, IERR) CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #else write(*,*) "This is a NO MPI run" MPIRANK=0 NPROC=1 #endif !================================================ !-- read ensemble size (NENS) …… !== define local dimensions and indexes allocate(jdisp(0:NPROC-1)) allocate(jlen(0:NPROC-1)) do irank=0,NPROC-1 CALL para_range(1,NENS,NPROC,irank,jsta,jend) jlen(irank)=jend-jsta+1 jdisp(irank)=jsta-1 end do CALL para_range(1,NENS,NPROC,MPIRANK,jsta,jend) #ifdef MPI_USE CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #endif N_LOC=jlen(MPIRANK)