Document for Stage0 (WRF model)

in Background Error Statistics CalculationFinal Report for 2004 CWB project

Yong-Run Guo

Mesoscale and Microscale Meteorology Division

National Center for Atmospherics Research

Boulder, Colorado, USA

Submitted to Central Weather Bureau, Taiwan, ROC

26 February 2005 November 2004

1. Overview Introduction

The Background Error Statistics (BES) file is one of main input files for WRF 3D-Var system. To a certain extent, the BES determined the performance of a 3D-var system. There are several methods proposed to derive the BES (Derber, Fisher). Especially for the WRF 3D-Var application, a program for BES calculation has been developed. This program includes several parts of the codes: Stage0, Stage1, Stage2, Stage3, Stage4, and Diags. The tasks of each of the parts are

Stage0: Read in the user’s model output data, either a series of the forecasts initiated at the consecutive times with a regular interval or a set of ensemble forecasts, and compute the difference fields for streamfunction, potential velocity, temperature, relative humidity, and surface pressure.

Stage1: Create the bins depending on the latitude of the Y-direction, and remove the bias from the difference fields for each of the bins.

Stage2: Calculate the regression coefficients for balanced potential velocity, temperature, and surface pressure with the streamfunction for each of bins, and produce the unbalanced potential velocity, the unbalanced temperature, and the unbalanced surface pressure.

Satge3: Compute the local (depending on the bins) and global eigenvectors/eigenvalues for the control variables: streamfunction, unbalanced potential velocity, unbalanced temperature, relative humidity, and unbalanced surface pressure.

Stage4: Compute the horizontal scale length used in recursive filter for the regional 3D-Var applications or the power spectrum for the global 3D-Var applications for each of the control variables.

Diags: Gather the results from Stage2, Stage3, and Stage4, and write out a BES file.

This program has the following features:

1. 1. The BES produced by this program is applied to the cv_options = 5 for WRF 3D-Var, which means i) the control variables are same as the cv_options = 3 but by using the coefficients projected on the vertical modes; ii) the horizontal covariance is modeled by the first order recursive filter with multiple passes; and iii) the vertical covariance is expressed by the eigenvector derived from the vertical covariance matrix, same as cv_options = 2.
2. 2. The program can be started at any of the stages when the input files have been produced by the previous stages.
3. 3. For Stage3 and Stage4 (the most expensive parts), the codes can be applied to the control variables one by one. So the computation could be completed on several machines simultaneously. It’s equivalent to the parallel runs and speed-up the BES file produced.
4. 4. Only the code for Stage0 is an interface for the user’s applications. The codes from Stage1 to Diags are transparency for users, and in general, users do not need to touch them.

Radio occultation (RO) of the terrestrial atmosphere is now possible through the use of signals transmitted by satellites of the Global Positioning System (GPS) and received by one or more satellites in low Earth orbit. The Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) mission will be launched in late 2005, and about 2,500 ~ 3,000 RO soundings will be provided on a daily basis. This is an important source of data for numerical weather prediction (NWP) in the near future. Meanwhile, NCAR/MMM (National Center for Atmospheric Research, Mesoscale and Microscale Meteorology Division) has developed an advanced 3D-Var (3-Dimensional Variational) data assimilation system over the past few years. Currently, this system has several main features: 1) assimilation of 19 types of observations, including GPS RO local refractivity, radar radial velocity and reflectivity, etc.; 2) multiple choices of control variables (CV) and the corresponding background error statistics (BES): cv_options=2 (NCAR approach) and cv_options=3 (NCEP approach), and the global BES files are provided as default; 3) choices of the minimization algorithms: Quasi-Newton and conjugate-gradient; 4) outer-loop implemented for the incremental approach accounting for the adaptive quality control (QC) and nonlinear effects of the observation operators; 5) ability of being applied to both MM5 model and the WRF model; and 6) code developed in the WRF framework and parallelized with high computing efficiency, and can be run on multiple platforms: IBM-SP, DEC, PC-Linux, Cray X1, Fujitsu VPP5000(?). [Kuo: Why do we put a question mark here?]

Although the GPS refractivity data are not yet available from the COSMIC mission, GPS RO data are available from two existing missions (CHAMP and SAC-C) through CDAAC (COSMIC Data Analysis and Archival Center) in Boulder. The Refractivity data from CHAMP are available from 19 May 2001 to present, and the data from SAC-C are available from 2 August 2001 to 14 November 2002. In preparation for the arrival of the COSMIC data, CWB (Central Weather Bureau) established collaboration with NCAR MMM Division to start the initial experimentation on the assimilation of GPS occultation with the versatile WRF 3D-Var data assimilation system, using CWB observation database and NWP model. It is anticipated that this system will be implemented into CWB operational environment. During the first year (2004) of the CWB-NCAR collaboration, six tasks were performed. This includes establishing the WRF 3D-Var system on the CWB computers, followed by a series of preliminary experiments with CWB operational model domain, CWB observation data base, and the NCEP AVN first guess for an interesting case.

In this documentreport, we will give the detailed description of the Stage0 code, especially for WRF model’s application, and the user’s guide to run the Stage0 program.

the accomplishments for the CWB-NCAR collaborative project (through joint efforts between NCAR and CWB staff) are summarized task by task below. The preliminary results from a series of the 3D-Var experiments with WRF 3D-Var/WRF model for Typhoon Dujuan case (30 to 31 August 2003) are presented. Based on these results, the future directions of the work are also suggested.

2. Description of the Stage0 code for WRF modelReview of the project tasks and accomplishments

A, Task of the stage0 1: NCAR/MMM provided the web site to access the WRF 3D-Var code and test data

There are 5 tasks in Stage0:

The WRF 3D-Var code and test data can be downloaded from

http://www.mmm.ucar.edu/wrf/users/download/get_source2.html

From this address, users can also find many links to get the information on the NCAR 3D-Var system, such as documentation, tutorial presentation, and activities, etc.

Ms. Hui-Chuan Lin, CWB staff, has successfully downloaded, implemented, and run 3D-Var system on CWB computers: PC-Linux, Fujitsu VPP5000(?). [Kuo: Why question mark?]

B, Task 2: NCAR/MMM provided the technical consulting support for CWB staff

NCAR staff has had many discussions with Ms. Hui-Chuan Lin during her visit to NCAR from June to September 2004, and helped her run the 3D-Var system: both the observation preprocessor (3DVAR_OBSPROC) and the 3D-Var code (wrf3dvar) itself. In addition to the public release version on the web site, Ms. Hui-Chuan Lin also received the updated versions of the 3D-Var system from NCAR staff.

C, Task 3: Assimilation of GPS radio occultation data

From the COSMIC web site

http://www.cosmic.ucar.edu

the GPS refractivity data (level2/wetPrf) can be downloaded from CDAAC. For the Typhoon Dujuan case, the CHAMP data on the Julian day 242 (30 August 2003, 170 soundings) and 243 (31 August 2003, 168 soundings) has been downloaded and processed. A decoder program: wetprf_decoder, which converts the downloaded file to the ASCII LITTLE_R format file, has been delivered to CWB staff, and it works properly.

D, task 4: 3D-Var observation preprocessor (3DVAR_OBSPROC) upgraded

1) Read in the forecast fields at the specific times from a series of the model forecast files initiated at the different dates.

At a specific date, for use of the NMC approach, two sets of the forecast fields at the same valid time from two forecast files should be read in while for use of the ENS (ensemble forecasts) approach, one set of the forecast fields from each of the ensemble members.

2) Convert the original model forecast fields to WRF 3D-Var analysis fields at A-grid: u and v components of wind, temperature, relative humidity, and surface pressure.

3) Compute the difference fields.

For the NMC approach, the difference fields between two forecasts at the same valid time are computed for each of the forecast dates.

For the ENS approach, the ensemble mean over the ensemble members is computed first, the difference between the ensemble forecasts and the ensemble mean are computed.

4) Compute the difference fields of the streamfunction and potential velocity based on the u and v difference fields by solving the Poisson equation using SIN_FFT method.

Note that i) all the operations for the conversion from wind components to the streamfunction and potential velocity are linear, so theoretically doing the conversion first and difference second should be same as the difference first and conversion second; ii) the solutions of the streamfunction and potential velocity derived from the u and v components in a limited area are not unique but depend on the lateral boundary specification. Using the SIN FFT to solving the Poisson equation implies that the zero lateral boundaries are used.

5) Write out the difference fields for each of the dates.Generalization of the domain definition

Originally, the 3DVAR_OBSPROC was developed for MM5 application, i.e. the domains are defined by a coarse domain and its nests. The conversion of (longitude, latitude) to domain (x, y) is completed by subroutine LLXY. However, the CWB 45-km model domain is a WRF-type domain with the central longitude different from the standard longitude, no information on the coarse domain and nesting were provided. To obtain the correct conversion between (lon, lat) to (x, y), a module map_utils adapted from wrfsi was introduced. With the help of this module, the correct (x, y) coordinates, the coarse domain and nested information were obtained (Fig. 1).

2) Observation errors for GPS refractivity

WMO code FM=116 is assigned to GPS refractivity N. In order to minimize the code changes, the data, N, are placed in the field dew_point in both LITTLE_R and 3D-Var OBS format files. As the initial implementation, the observation error specification is based on Huang et al (2004)

where N_bottom = 10 and N_top = 3 N-units are N error at the bottom and top of the atmosphere. p₀ = 1000 hPa is the pressure at the bottom and p_t = 10 hPa is the pressure at the top of atmosphere. p is the pressure at the observed level. The observation errors are placed in the field dew_point%error. Figure 2 shows the vertical profile of the N errors.

This single profile of the N observation errors is not a good global representation of the RO errors (see below). Kuo et al. (2004) showed that the GPS RO refractivity observation errors can vary with latitudes and seasons. Next year, we may refine the N error specification based on the observation error estimates described in Kuo et al (2004). A data quality control procedure may also be introduced to remove or truncate bad data

3) GPS refractivity data distribution for Typhoon Dujuan case

The distributions of the CHAMP level2 wetPrf data within 6-h time window centered at 1200 UTC 30 and 31 August 2003 are plotted by 3DVAR_OBSPROC (Fig. 3).

BE, Details of the code structureTask 5: Use of the WRF 3D-Var analysis for CWB’s regional forecast system

1) Main program: wrf3dvar/main/gen_be/gen_be_stage0.F

The structure of this program is very similar to wrf3dvar/main/da_3dvar.F but to replace “call da_3dvar_interface” by “call da_gen_be_stats_interface”. The module “module_da_gen_be_stats_interface.F” is added under the directory “wrf3dvar/da_3dvar”. With this module, the subroutine “da_3dvar/src/da_stats_be/da_stats_be.F” is called, and then in “da_stats_be.F”, subroutines: “da_stats_be/da_init_be”, “DA_Gen_Be_Stats/DA_Stats_Namelist”, and “DA_Gen_Be_Stats/DA_Statistics_Step0” are called sequentially. Actually, the subroutine “DA_Statistics_Step0” is the core part of the stage0 program,program; other parts of the code are only for matching with the wrf3dvar framework, such as the netCDF_IO, mpp capabilities, etc. The flowchart is summarized as below. If users develop their own stage0 program with their own model forecasts, the procedure can be much simplified, directly starting from the code of DA_Statistics_Step0.F with the necessary modifications.

The flowchart is summarized as below.



2) Subroutine: DA_Statistics_Step0

This subroutine includes two loops: loop over the file_date and loop over the ensemble member. Based on the file_date and ensemble member index, the specific file name will be formed by calling subroutine “DA_Make_Filename”, then the needed data are read by “call med_initialdata_input_3dvar”, a standard WRF 3DVar subroutine to read the netCDF file, and transfer to u, v, t, rh, and Psfc by calling WRF 3DVar subroutine “DA_Setup_FirstGuess”, and assign these variables to the variables named by xb24 and xb12. For the BACKGROUND_OPTION = 2 (ensemble approach), only variable xb24 need to be read in and xb12 is stored the ensemble mean computed by subroutine “DA_Make_Ensemble_Mean”.

Until now, the two sets of the fields, u, v, t, rh, and Psfc, for the difference calculation are ready. Then the difference fields are created by calling “DA_Difference”, and the difference of the u, v components will be converted to streamfunction psi and potential velocity chi by calling “DA_New_Statistics_Variable”.

Finally, the difference file will be written out for the specific date and ensemble member by calling:”DA_Write_Diff”. Each of the binary file contains a header record: time (character*10), ide, jde, and kde (the dimensions of domain, integer), and 7 data records, difference of psi, chi, t, rh, Psfc, and the full fields of h and lat.

3) Convention of the input and output file names

· Input file:

The input file name is composed of three parts: i) directory_name and file_head, ii) file_date; and iii) ensemble member index. The directory_name and file_head are provided through the namelist file: namelist.stats. The file_date has 19 characters, like ccyy-mm-dd_hh:00:00, and the ensemble member index are one or two digits, i.e. currently only maximum of 99 members are allowed for each of the times So the input file name is trim(directory_name)/trim(file_head)_ccyy-mm-dd_hh:00:00.trim(index), such as “GEN_BE_data/wrfout_d01_2002-01-01_00:00:00.5”. For NMC approach, there is no “index” part, i.e. “GEN_BE_data/wrfout_d01_2002-01-01_00:00:00”.

· Output file name:

The output file names are formed by the program stage0 and recognized by the stage1 program. The file names are hardwired by program as wrf.diff_ccyy-mm-dd_hh:00:00.index, such as “wrf.diff_2002-01-01:00:00.5”, and for NMC approach, there is no “index” part, i.e. “wrf.diff_2002-01-01:00:00”.



4) Directories special for program stage0

There are two directories special for program stage0: da_3dvar/src/DA_Gen_Be_Stats and da_3dvar/src/da_stats_be. As mentioned before, the main program gen_be_stage0.F and interface module da_gen_be_stats_interface” are resided in main/gen_be and da_3dvar, respectively.

3. User guide for program Stage0

1) Compile and running

The compiling procedure is similar to rgw wrf3dvar but need two compiling steps:

configure

compile 3dvar

compile gen_be

The executables for gen_be_stage0, 1, 2,…, are under the directory main/gen_be.

2) Shell script for running Stage0

The shell scripts for running Stage0 program is /run/gen_be/da_wrf_stage0.csh

In this script, there are 4 parts need to be modified by users for their applications

i) setup the directories for WRF 3DVar source code, input data, and job running;

ii) setup the domain dimension and grid size;

iii) Setup the information for job running, including the start and end times, which type of the BES produced, the prefix of the filename, file interval, and first forecast time in hours.

iv) Namelist.stats. In general use, users do not need to touch this part since most of namelist variables can obtain the values from the above part 3). Only in case of debugging, you may need to set a non-zero value to PRINT_DETAIL, and .true. to TEST_TRANSFORMS.

Two more namelist variables: from_mss and mss_directory, are used to acquire the input data files from the NCAR Massive Storage System (MSS) directly. This may be useful for use of the ensemble forecast approach with the large number of ensemble for many dates. If users do not have NCAR MSS available, ignore these variables.

3) Namelist.stats

There are 3 records in the namelist file:

&FILERECD

DIRECTORY ; the directory for input data file (character (len=120))

FILE_HEAD ; the prefix of the input file names (character (len=120)). For example, the file_head is ‘wrfout_d01’ for input file: ‘wrfout_01
_2002-01-01_00:00:00’.

BGN_DATE ; the starting time for BES calculation (character (len=24)): ccyy-mm-dd_hh:00:00

END_DATE ; the ending time for BES calculation.

TEST_TRANSFORMS ; .FALSE. not do the transform test, .TRUE. doing the test of the transform between (u,v) and (psi,chi).

PRINT_DETAIL ; =0 no details printed, >0 print more things for debugging

From_mss ;.FALSE. not use NCAR MSS, .TRUE. acquire the input data from NCAR MSS during the execution.

Mss_directory ; the input data directory in NCAR MSS (character (len=120)).

&TIMERECD

T_FORECAST1 ; the earlier forecast time for BES calculation

T_FORECAST2 ; the later forecast time for BES calculation, i.e. T_FORECAST2 = T_FORECAST1 + FILE_INTERVAL. So when the T_FORECAST1 and FILE_INTERVAL are specified, you don’t need to specify this variable to avoid inconsistent.

FILE_INTERVAL ; Time interval in hour of the input data files.

&ANALTYPE

BACKGOUND_OPTION ; BES type: 1 = NMC, 2 = ENS

MEMBERS ; If BACKGROUND_OPTION = 1, MEMBERS = 1,

If BACKGROUND_OPTION = 2, set the number of ensemble members for each of the times to MEMBERS.

4) Shell script for running Stage1 to Diags

The script for running all other stages is run/gen_be/gen_be_sample.csh

This shell script can be used to run all stages or the select stage.

5) Print and plot the background error statistics

· ./gen_be_diags_read.exe can be used to print the BES in ASCII format to fort.171, 172,…..

· The shell script da_3dvar/utl/plot_eigen_in_be.csh can be used to plot the global the first five eigenvalues, eigenvectors, and the first five local eigenvalues, and the scale lengths. You just need to copy that shell script to your working directory and edit it for your application. The BES plots are easy to be produced for your check.

· Also there are some NCL plotting routines for plotting the correlation coefficients, eigenvalues and eigenvectors.

2) At the request of CWB, NCAR developed the code to output the WRF 3D-Var analysis increments file, including the variables of u, v, t, p, q, and ph, for converting to CWB model initial condition.

3) To reduce the conversion errors, NCAR proposed a procedure through adjustment of water vapor to ensure that the surface pressure P_sfc (total air mass) will not be changed before and after the conversion (this is not yet fully implemented).

NFS surface pressure and water vapor

Water vapor pressure Pv in a column is defined as

where r is the total air density, q is the specific humidity. Note in CWB NFS model, the vertical coordinate s is defined by the total pressure p and p_sfc,

WRF surface pressure and water vapor

In WRF the vertical coordinate h is defined by the dry hydrostatic pressure. The water vapor pressure in a column is defined as

where r_d is the dry air density, r is the mixing ratio; p_d and pd_sfc are the dry pressure and dry surface pressure.

Procedure of the conversion

a) To obtain m for WRF from NFS data

b) To obtain the dry air pressure p_d at the model h levels



c) To obtain the NFS dry pressure p_d(nfs) at the NFS s levels.

The water vapor pressure Pv_k at the level s_k can be calculated as:



d) To obtained the mixing ratio r at the WRF h levels based on p_dk, p_dk(nfs) and r_k(nfs)

Use the linear or the log-p interpolation in the vertical.

e) To obtain the WRF water vapor pressure Pv_wrf in a column based on the mixing ratio r and m, dh, and the WRF surface pressure:

At this point, we may say that the NFS to WRF conversion was completed. But due to interpolation truncation errors, the total water vapor pressure Pv from NFS may be different from Pv_wrf, thus the surface pressure P_sfc from NFS may be slightly different from P_sfc(WRF).

Water vapor adjustment

Because Pv_wrf only depends on r, m, and h, and since m and h cannot be modified, the only way to make Pv_wrf=Pv is to adjust the mixing ratio r at the WRF model levels.






Because Pv is a linear function of the mixing ratio, with the new mixing ratio r_adj, the new Pv_wrf will be equal to Pv, and P_sfc(WRF) = P_sfc.

The advantages of the procedures of “NFS->WRF->xb” and “xa->WRF->NFS” are that the analysis from WRF 3D-Var analysis can be used to initialize both the CWB NFS model and the WRF model. However, it will introduce some truncation errors from the conversion between the NFS and WRF. In the future, we should consider the procedures of “NFS->xb” and “xa->NFS” directly.

F, Task 6: Analysis of the WRF 3D-Var experimental results

To ensure the WRF 3D-Var system working properly, the following tasks were conducted.

1) The innovations of the GPS refractivity soundings were plotted. We found that in most cases, the innovations are less than the observation errors assigned here. As we discussed before, the observational errors may be too large, as defined in Huang et al. 2004.

2) The Background Error Statistics (BES) tuning was conducted. Since we do not have BES file from the CWB model, the global NCEP-derived BES has to be used for the 3D-Var experiments. The subjective tuning is performed through the single obs tests for GPS refractivity.

3) A series of the 3D-Var/WRF forecast experiments were carried out. Preliminary results are reported here. For Typhoon Dujuan’s track forecast, the impact of the GPS Ref. data is moderate, but 3D-Var analysis usually gave superior results over the control run.

Observation operator (Local refractivity observation operator)

The refractivity N is defined as:

where p is the pressure in hPa, T in the temperature in K, and q is the specific humidity in kg/kg (Zou et al 1995). This is a one-step formula for refractivity calculation because p, T, and q are the analysis variables in wrf3dvar.

The innovations (O-B) of GPS refractivity with NCEP analysis as the first guess and the observation errors

The (O-B) and observation errors for 7 GPS soundings within the 6-h time window centered at 1200 UTC 30 August 2003 are shown in Fig. 4. As pointed out earlier, the observation errors based on Huang et al. (2004) may be too large. The large values of (O-B) located below 8-km and the large negative values near the surface indicate potential data problems. This suggests that certain QC (Quality Control) should be applied prior to the assimilation.

Background error statistics tuning

Several steps can be taken to tune the background error statistics.

a) Tuning the BES subjectively (when collected enough (O-B) data, tuning the parameters is possible based on the Hollingsworth and Lonnberg (1986) method);

b) Increasing the variances of BES will reduce the magnitudes of the increments;

c) Reducing the scale-length will reduce the influence range of the observations.

The tuning factors for cv_options=3 and cv_options=2 are listed in Table 1.

Table 1 Tuning factors for cv_options=3 and cv_options=2

CV_OPTION = 3				CV_OPTION = 2
Control	Variance Default / Tuned	H-length Default / Tuned	V-length Default / Tuned	Control	Variance Default / Tuned	H-length Default / Tuned
y	0.25 / 0.25	0.75 / 1.0	1.5 / 1.5	y	1.0 / 1.35	1.0 / 0.2
cu	0.25 / 0.25	0.75 / 1.0	1.5 / 1.5	c	1.0 / 1.35	1.0 / 0.2
Tu	0.25 / 0.35	0.75 / 1.0	1.5 / 1.5	pu	1.0 / 1.35	1.0 / 0.2
qp	0.25 / 0.10	0.75 / 1.75	1.5 / 1.5	q	1.0 / 1.00	1.0 / 1.0
Psfcu	0.25 / 0.35	0.75 / 1.0	1.5 / 1.5

For cv_options=3, the background variance factors for the mass variables, T_u and Psfc_u are increased (background weightings decreased) but for moisture, q_p, is decreased (background weighting increased). This means that increments for T_u and Psfc_u will be increased and for q_p decreased. The horizontal scale-lengths are increased for all of the control variables, especially for moisture (1.75), because the default value (0.75) is too small. As a result, the influence range for the mass variables (T and P) is only about 1,500-km (Fig. 5d). With the tuned value (1.0), the influence range became more reasonable, about 2,300-km (Fig. 5b). This helps to spread the GPS data impact more widely in a 45-km resolution CWB domain (Fig. 6a).

For cv_options=2, the default background variance factors are too small, and the scale-length are too large for wind and pressure (Fig. 5c). With the tuned values, more reasonable increments are obtained (Fig. 5a).

Single OBS tests for GPS refractivity

The cross-sections of the temperature and pressure increments for single GPS refractivity OBS tests for Typhoon Dujuan case in CWB domain are shown in Fig. 5. After tuning, the influence range of an OBS is more reasonable by subjective inspection as mentioned above. A single GPS refractivity OBS (data=1.0 N-unit with the error=1.0 N-unit) is located at the middle level in the center of the domain (see explanation in Fig. 5).

These figures for the single OBS tests are easily plotted with the analysis increment output file from WRF 3D-Var by using a utility program: da_3dvar/utl/convert_to_mm5v3.f90 and MM5/GRAPH.

Background error for GPS refractivity

Because there are many transforms in 3D-Var system and the GPS refractivity is not one of the direct model variables, it is not straightforward to know what the background errors of the GPS refractivity are in the analysis system. However, based on the following derivation, the background error of the GPS refractivity can be calculated from the (O-B), (O-A), and the observation error s_o. This has already been implemented in the latest WRF 3D-Var code.

For single OBS

B is the background value, O is the observation value, and the

s_b² is the background variance and s_o² is the observation variance.

Then, the analysis A should be






When the (O-B), (O-A), and s_o² of GPS refractivity are known, the background error for GPS refractivity can be computed.

Table 2, The background error s_b for (O-B)=1.0 and s_o=1.0

BES	CV_OPTION=3		CV_OPTION=2
	(O-A)	sb	(O-A)	sb
Default	0.171	2.204	0.494	1.012
Tuned	0.324	1.444	0.487	1.027

From the single GPS refractivity tests, we obtained the background error, s_b (Table 2). It is clear that the bigger s_b is, the better fit to the observation, with smaller (O-A). With the tuned factors for cv_options=3, the (O-A) = 0.324 and s_b = 1.444 are more realistic in comparing with the case of (O-B) = 1.0 and s_o = 1.0.

Even when more sophisticate approaches (NMC-method with the ensemble forecasts, the advanced tuning techniques, etc.) are used, subjective inspection for the analysis increment responses to the observation may still be a necessary step prior to the 3D-Var experiments.

Experiment design

CNTRL : No assimilation, initiated with NCEP AVN analysis

3DREF : Assimilation of GPS refractivity only

3DCWB : Assimilation of CWB observations (SOUND, AIREP, SYNOP. SHIPS, PILOT, SATOB, SATEM, METAR).

3DCWBREF : Assimilation of CWB observations + GPS refractivity.

In the WRF 3D-Var experiments, the first guess fields are obtained from NCEP AVN, and the background error statistics are interpolated from the NCEP global statistics (192 latitudes with resolution of 0.945^o in horizontal and 42 h levels in vertical). The control variables, y, c_u, T_u, q_p and Psfc_u, are selected (CV_OPTION = 3). The initial times are 1200 UTC 30 and 1200 UTC 31 August 2003.

The analysis increments from 3DREF and 3DCWBREF are shown in Fig. 6. From Figure 6, the increments from WRF 3D-Var are reasonable as a result of the assimilation of GPS refractivity and all other types of data.

Experiment results

Table 3 the averaged track forecast errors for Typhoon Dujuan

Experiment	2003083012Z			2003083112Z
	24-h(km)	48-h(km)	72-h(km)	24-h(km)	48-h(km)	72-h(km)
CNTRL	69.2	67.1	221.1	110.2	115.7	_
3DREF	69.8	85.8	214.2	115.1	120..7	_
3DCWB	58.5	89.3	259.1	103.6	104.9	_
3DCWBREF	52.5	79.8	260.0	105.6	108.1	_

Here we just present some preliminary results of the Typhoon track forecast. Figure 7 shows the track forecasts from the different experiments for Typhoon Dujuan initiated at 1200 UTC 30 and 31 August 2003. From the figures, there are only small differences between the experiments. When the averaged track forecast errors are calculated based on the 3-hourly model output for the periods of 3~24h, 27~48h, and 51~72h, we can see that the 3D-Var experiments, in general, gave better results than the control experiment. The results are listed in Table3.

Bearing in mind that no bogus vortex was used, and the resolution is 45-km for all these experiments, and there are many weakness pointed out earlier (GPS refractivity observation error specification, global background error statistics, etc.) in the current sets of 3D-Var experiments, we cannot expect significant positive impacts from the few GPS refractivity profiles. However, it is encouraging to find that there is moderate impact of assimilating a relatively small number of GPS data. Also, the use of 3D-Var gives positive impact. There is much room for improvements for the assimilation of the GPS data with WRF 3D-Var system. As the number of GPS data is increased (with the availability of the COSMIC data), data quality is improved and the 3D-Var system is fine tuned, we believe that the assimilation of the GPS data will have a positive impact on numerical weather prediction and typhoon forecast in the vicinity of Taiwan.

3, Future work

The basic WRF 3D-Var system has been implemented into the CWB environment, and the preliminary results of assimilating the GPS refractivity are encouraging. In order to assess and improve the performance of the WRF 3D-Var system at CWB, the following tasks should be performed in the near future:

Perform semi-operational run with 3D-Var system in CWB

To develop the direct data converters between NFS model and WRF 3D-Var.

To improve the OBS preprocessor incorporating the sophisticate observation error specification and quality control procedure.

To improve the efficiency of 3D-Var system in CWB operational environment (Fujitsu VPP 5000 ?)

To develop the post processing (display and verification) software

2) Continuation of the Typhoon Dujuan case study

To do more experiments with both higher resolution (15-km, 5-km) CWB model and WRF model by using CWB archived data.

To introduce the Typhoon bogus technique in 3D-Var system

To improve the method of determining the forecast Typhoon track, and more diagnostic program.

Perform better estimate of the background error and observation error statistics by using more data

To establish the archive system for the semi-operational run (observations, background, innovations, etc.).

To estimate the s_b and s_o for CWB model and observation data.

To derive the CWB-specified BES for different seasons, different resolutions, etc.

Assimilation more types of observations, such as satellite radiances, etc.

To develop the BUFR format observation preprocessor for satellite radiance data.

To develop the radiance assimilation code in 3D-Var system

To assess the impact of radiance assimilation by comparing with the assimilation of the retrieved satellite data.

The selection of these tasks for 2005 will be determined jointly by CWB and NCAR based on available resources and project priorities.

References

Huang, C.-Y, Y.-H. Kuo. S.-H. Chen, and F. Vandenberghe, 2004: Improvements in Typhoon Forecasting with Assimilated GPS Occultation Refractivity. Weather and Forecasting, Submitted.

Kuo, Y.-H., T.-K. Wee, S.Sokolovskiy, C. Rocken, W. Schreiner, D. Hunt, and R. A. Anthes, 2004: Inversion and Error Estimation of GPS Radio Occultation data. Special issue of the Journal of the Meteorological Society of Japan, Vol. 82, No. 1B, 507-531.

Zou, X., Y.-H. Kuo, and Y.-R. Guo, 1995: Assimilation of Atmospheric Radio Refractivity Using a Nonhydrostatic Adjoint Model. Mon. Wea. Rev. 123, 2229-2249.



Figure 1, CWB 45-km resolution domain



Figure 2, the vertical profile of the GPS refractivity observation error.



2003083009Z to 2003083015Z 2003083109Z to 2003083115Z

Figure 3, GPS refractivity data from CHAMP level2 wetPrf for Typhoon Dujuan case



Figure 4, Innovations (O-B) of GPS refractivity with NCEP AVN analysis as the background, and the observation errors.

Figure 5, The T/P increments cross-section (West-East about 4,500 km) for Single GPS Ref OBS tests for Typhoon Dujuan case in CWB domain (a) cv_options=2 with the tuning factors: 1.35, 1.35, 1.35, 1.00 for VAR and 0.2, 0.2, 0.2, and 1.0 for LEN; (b) cv_options=3 with the tuning factors: as1=(0.25, 1.0, 1.5), as2=(0.25, 1.0, 1.5), as3=(0.35, 1.0, 1.5), as4=(0.10, 1.75, 1.5), and as5=(0.35, 1.0); (c) cv_options=2 default; (d) cv_options=3 default. The OBS value and error of the Ref. is 1.0 located at x=112, y= 65, z=15. The CWB domain is 222x128x30 with 45-km grid distance.



Figure 6, wind and pressure increments at the lowest model level for Exp. 3DREF and 3DCWBREF.



Figure 7, Typhoon Dujuan track forecast for different experiments initiated at 1200 UTC 30 and 31 August 2003.