ATOP-The Advanced Taiwan Ocean Prediction System Based on the mpiPOM . Part 1 : Model Descriptions , Analyses and Results

A data-assimilated Taiwan Ocean Prediction (ATOP) system is being developed at the National Central University, Taiwan. The model simulates sea-surface height, three-dimensional currents, temperature and salinity and turbulent mixing. The model has options for tracer and particle-tracking algorithms, as well as for wave-induced Stokes drift and wave-enhanced mixing and bottom drag. Two different forecast domains have been tested: a large-grid domain that encompasses the entire North Pacific Ocean at 0.1° × 0.1° horizontal resolution and 41 vertical sigma levels, and a smaller western North Pacific domain which at present also has the same horizontal resolution. In both domains, 25-year spin-up runs from 1988 2011 were first conducted, forced by six-hourly Cross-Calibrated Multi-Platform (CCMP) and NCEP reanalysis Global Forecast System (GSF) winds. The results are then used as initial conditions to conduct ocean analyses from January 2012 through February 2012, when updated hindcasts and real-time forecasts begin using the GFS winds. This paper describes the ATOP system and compares the forecast results against satellite altimetry data for assessing model skills. The model results are also shown to compare well with observations of (i) the Kuroshio intrusion in the northern South China Sea, and (ii) subtropical counter current. Review and comparison with other models in the literature of “(i)” are also given.


InTRODucTIOn
The western North Pacific region is home to over a quarter of the World's population.Understanding and predicting the regional weather and climate variability are clearly of interest.The oceans play an indispensable role in weather and climate variability.Through coupling of the sea-surface temperatures (SST's) and upper-ocean processes with the atmosphere, the seas moderate winter winds and make summer heat more tolerable.The seas are also spawning grounds for storms that affect hydrological cycles upon which life depends, and extreme surges (due to e.g., storms, Tsunamis) cause human sufferings of enormous proportions.Understanding and the ability to predict the ocean and its interaction with the atmosphere and land are necessary.Over long time scales: interannual, decadal and longer, improved understanding and prediction can lead to better planning and management of the earth's increasingly scarce resources.
Major boundary currents of the western North Pacific are the Mindanao Current and the Kuroshio (Fig. 1).These are the southern and northern branches, respectively, of the westward-flowing North Equatorial Current (NEC) which carries approximately 30 ~ 60 Sv (1 Sv = 10 6 m 3 s -1 ) in the upper thermocline and which near the surface bifurcates at approximately the 12 ~ 13°N latitude off the eastern coast of the Philippines (Nitani 1972;Toole et al. 1988Toole et al. , 1990;;Qu and Lukas 2003;Kim et al. 2004;Jensen 2011).The southward-flowing Mindanao Current transports approximately 13 ~ 33 Sv of the NEC water in the upper layer (600 ~ 1000 m); it forms the Mindanao eddy, and contributes to flow through the Indonesian archipelagos (the Indonesian Throughflow; e.g., Lukas et al. 1991).The remaining transport, approximately 15 ~ 35 Sv, is carried northward by the Kuroshio along the eastern coast of the Philippines (Nitani 1972).Some Kuroshio water recirculates in relatively tight anticyclonic cells east of the current in the Philippines Sea (12 ~ 20°N), and approximately 1 ~ 5 Sv intrudes into South China Sea (SCS) through the Luzon Strait (Qu et al. 2000).The main Kuroshio continues northward along the eastern coast of Taiwan and the continental slope of the East China Sea.Estimates of the transport are 15 ~ 35 Sv east of Luzon at 18°N (Sheu et al. 2010) just south of the Strait, approximately 21 ~ 22 Sv northeast of Taiwan (Johns et al. 2001), and about 25 Sv at the PN-line (the East China Sea -Okinawa section near 29°N; Kagimoto and Yamagata 1997).Localized transport increases are due to recirculating anticyclonic cells or eddies formed along the eastern side of the Kuroshio (Nitani 1972).
Eddies of both signs are abundant in the subtropical gyre of the western North Pacific east of Philippines and Taiwan (Chang andOey 2011, 2012;Chelton et al. 2011).These eddies are generally elliptically elongated roughly in the zonal direction with sizes of 200 ~ 500 km, and are believed to be the result of instability of the subtropical counter current (STCC; and its various branches) (Hasunuma andYoshida 1978;Qiu 1999).The STCC front and eddies are potentially important in modifying atmospheric winds (e.g., Kobashi et al. 2008); through their large heat contents, warm eddies can affect the intensity and spawning of typhoons (e.g., Lin et al. 2009).
Warm and salty Kuroshio waters intrude into the South China Sea and East China Sea, moderating the regional climate in winter and modifying the shelf-seas' heat, salt and biogeochemical (e.g., nutrient and carbon) balances (Liu et al. 2000).These variations, together with tidal, wind and river-driven currents, have important consequences to the circulation and mixing of the shelves and seas (Qu 2000;Guo et al. 2006;Isobe 2008).Through air-sea-wave coupling, they also alter winds, sea-surface heat and mass (evaporative and precipitation) fluxes, as well as wind waves and swells, thereby affecting continental coastal areas as well (Chang et al. 2010;Chiang et al. 2011).
In addition to research that can lead to basic understanding of processes, models of currents and water mass distributions are also useful for practical applications -such as tracking of oil-spill and pollutants, search and rescue, estimates of the flushing rates of coastal seas and bays, and identification of locations of strong upwelling (and downwelling) frontal regions, e.g., for fish-catch.The ability to model and forecast the circulation and mixing of the western North Pacific Ocean and marginal seas is thus seen to be of interest for research and applications.This paper presents our efforts at the National Central University (NCU) to develop an ocean prediction system -the Advanced Taiwan Ocean Prediction (ATOP) system.The goals are: (i) to provide regional ocean forecasts of currents, sea levels, SST's that are accessible by the public as well as the scientific community; and (ii) to conduct basic research of oceanic, air-sea coupled and other related processes.This paper focuses on goal#1, and describes our attempt to build the ATOP system.To provide a more realistic estimate of the sea state, observational data are assimilated into the model, and a preliminary validation of the forecasts against observations is conducted to provide estimates of the forecast error and uncertainty; more extensive model validations are planned in order to build confidence.The "alpha" (test) version of the ATOP system is now operational; once a day it provides -7 day analysis and +7 day forecast of the wind-and eddy-driven sea states for the North Pacific Ocean, focusing on the western North Pacific (http://mpipom.ihs.ncu.edu.tw/index.php).
Section 2 describes mpiPOM -the ocean model used for ATOP, and section 3 the forcing and observational data.Section 4 details the specifics pertaining to ATOP: initial and boundary conditions and various experiments and modeling strategies.Section 5 describes the data-assimilation method and ATOP's ±7 day analysis and forecast cycle.Section 6 describes the methods used to access the forecast skills.Section 7 presents the results of the forecast and skillassessment.Section 8 presents the simulated general circulation of the northwestern Pacific, and compares the seasonal variability of the model STCC eddies against satellite data.Finally, section 9 outlines future plans and extensions, and concludes the paper.

The MPIPOM
ATOP uses the mpi-version (Message Passing Interface) of the Princeton Ocean Model (Blumberg and Mellor 1987); Dr. Toni Jordi implemented MPI's into the model (http://www.imedea.uib-csic.es/users/toni/sbpom/index.php).The mpiPOM is an efficient parallel code and retains most of the original POM physics.It has some new features and incorporates some of the Princeton Regional Ocean Forecast System (PROFS) -http://www.aos.princeton.edu/WWWPUBLIC/PROFS/publications.html) modules, outlined below and in subsequent sections.Second-order Leapfrog time-differencing and centered-space differencing are defaults.The centered-space differencing may be replaced by the Smolarkiewicz (1984) iterative upstream scheme to reduce or eliminate 2Δ-grid oscillations; 4 iterations (FOR-TRAN variable nitera = 4) and Smolarkiewicz weight > 0.95 (1 > sw > 0.95) work well.Terrain-following sigma-coordinate is used in the vertical and there is a fourth-order scheme option to reduce the internal pressure-gradient errors (Berntsen and Oey 2010).The Mellor and Yamada's (1982) turbulence closure scheme is modified to include turbulence energy due to breaking waves near the surface (Craig and Banner 1994;Mellor and Blumberg 2004).The Smagorinsky's (1963) shear and grid-dependent horizontal viscosity is used with a nondimensional coefficient = 0.1; the corresponding horizontal diffusivity is (generally) made 5 ~ 10 times smaller.The Yin and Oey's (2007) wind-drag formula with high wind-speed limit is used.Various open boundary condition options are available; those work well are Flather (1976) radiation and flow-relaxation schemes (see review in Oey 2010; http://www.aos.princeton.edu/WWWPUBLIC/PROFS/PUBLICATION/OFES20101104.pdf).Periodic and double-periodic boundary conditions are also available for idealized simulations.Data-assimilation subroutines are also included, to assimilate satellite altimetry (SSH anomaly) data, surface drifters, ADCP's and temperature and salinity data.These use simple (and fast) statistical and/or optimum interpolation (OI) methods (see section 4).A more sophisticated Localized Ensemble Transformation Kalman Filter (LETKF; Hunt et al. 2007;Miyoshi et al. 2010;Miyazawa et al. 2012) scheme is presently being tested.

Additional mpiPOM Modules
Wave-Enhanced Bottom Drag: A Grant and Madsen (1979) type bottom-drag scheme is used in the bottom boundary layer to empirically account for increased bottom roughness due to surface gravity waves; this is especially effective over shallow shelves and under strong wind conditions.The wave-induced velocity near the bottom is estimated as [assume linear wave; Nielson (1992) .u u where u τ = bottom friction velocity.A knowledge of k p is required.At present, the Mellor's (2008) wave model has yet to be incorporated into mpiPOM; the k p is therefore fixed at 2 200 r m thus favoring high wind conditions when the wavelengths ~ 200 m (Li et al. 2002).
Passive Tracer: This is solved in the same way as temperature or salinity, except that zero-flux condition is applied at all boundaries; modifications to allow non-zero fluxes across the air-sea interface are straightforward.The Smolarkiewicz (1984) iterative upstream scheme is used for the advection terms.Multi-component tracers may be tracked, and they are specified either as point sources (on the RHS of the tracer equation) or directly at the release locations.The tracers may be two-dimensional at the surface only or three-dimensional.
Particle Tracking: Awaji et al.'s (1980) scheme is used (Chang et al. 2011).The position of a marked particle is tracked Lagrangianly from its position x o at time t o to a new position x = x o + Δx at time t: The first integral on the right is seen to be the Eulerian part, while the {·} within the second integral is the correction to the velocity field as the particle moves across ., x u x t dt

#
Similarly as tracers, multiple marked particles can be tracked, distinguished by their initial positions.At present, the scheme is implemented to track horizontally (e.g., at the sea surface) only.
Stokes Drift: Inclusion of wave-induced drift (Stokes 1847) can be important for tracking particles or tracers especially near the ocean's surface.A simple formulation due to Kenyon (1969) is used.The author derived the Stokes drift U s as a double integral (in wave number space) of the two-dimensional wave spectrum.By using the Pierson and Moskowitz's (1964) empirical (one-dimensional) spectrum for fully-developed seas, he simplified the formula which we rewrite below in vector form for the Stokes drift velocity: where u a is the wind velocity.The Stokes drift is added to the modeled velocity when used as an option in mpiPOM for tracking tracers and/or particles.

Model Flowchart
After initializations in which model parameters and switches, the input grid, topography, initial (potential) temperature (T) and salinity (S) fields (from climatology or restart files) are specified, the code goes through the internal time-stepping loop (iint = 1, iend; Fig. 2) in which various time-dependent forcing data such as monthly climatological T/S, SST (e.g., from satellite MCSST), surface fluxes (e.g., winds), rivers, are first prepared (interpolation).The (original) POM algorithm then begins with the subroutine advance.fwhere various interactive feedback variables between external (2d depth-averaged) and internal (3d dependent) modes are calculated and the code loops through the external time-stepping (ient = 1, isplit) scheme 'isplit' times for every one internal time step 'iint'.The 3d-dependent fields such as the velocity, T, S, and various turbulence variables are then time-stepped forward.Data assimilations are then performed depending on the corresponding switches (on or off, and how often, i.e., assimilation time intervals).At the end of the model integration (iint = iend), NetCDF outputs are produced.Figure 3 gives the directory structure of a typical mpiPOM run.

Topography
Model topography is from ETOPO2 which has a 2-minute resolution (~3 km; http://www.ngdc.noaa.gov/mgg/fliers/01mgg04.html).For model resolution that is coarser than this, care is taken to retain unresolved islands and island groups if their summed areal coverage exceeds a (user-specified) percentage (say 50%) of the model's coarser cell's area.A topographic smoother is used so that at each cell the depth change is no more than a certain percentage of the total water depth at that cell, i.e., ΔH/H < SL min (= 0.2 is used here).A grid-generation FORTRAN code is used to directly read, smooth and interpolate the global ETOPO2 data onto the user's specified domain and grid sizes.See Fig. 1 for the interpolated topography in the western North Pacific domain.In the vertical, finer, logarithmically distributed sigma cells are used near the free surface and ocean's bottom, so that the first grid cell is approximately 2 ~ 3 m below (above) the surface (bottom) in 1000 m water depth, and the corresponding grid spacing in the middle of the water column is approximately 30 ~ 38 m.

WOA T/S Data
The grid-generation FORTRAN code reads and interpolates the World Ocean Atlas (WOA) data from NODC (http://www.nodc.noaa.gov/OC5/WOA05/pr_woa05.html) onto the model grid.The T and S are first horizontally interpolated onto the model's (orthogonal curvilinear) grid cells at each WOA standard level - (1, 10, 20, 30, 50, 75, ..., 150, 200, ..., 300, 400, ..., 1500, 1750, 2000, 2500, ..., 5500) m, then they are vertically interpolated onto the model's sigma cells.To avoid problems near land-(seafloor-) water edges, land (seafloor) values are first filled with values from the nearest water points before the vertical interpolation.The deepest WOA values are used for model's cells below the deepest WOA level (= -5500 m).The WOA sea-surface temperature and salinity (i.e., at z = 0 m) are also similarly interpolated (horizontally only) onto the model's grid.The resulting twelve monthly fields are then used in the mpi-POM simulation for boundary conditions and for the purposes of relaxing the T/S fields to the observed climatology especially in the deep portions (below z = -2000 m) of the modeled domain.

Wind and Surface Flux Data
Three global (near-global) wind and surface-flux data sources have been used for the various spin-up runs as well as for the forecast runs.The first set is the Cross-Calibrated Multi-Platform (CCMP; http://podaac.jpl.nasa.gov/datasetlist?search=ccmp;Atlas et al. 2009) product from July/02/1987 to December/31/2009.This is a 6-hourly gridded (1/4° × 1/4°) dataset that combines ERA-40 re-analysis with satellite surface winds from Seawinds on QuikSCAT, Seawinds on ADEOS-II, AMSR-E, TRMM TMI and SSM/I, as well as wind from ships and buoys.The second set is the 6-hourly NCEP 0.5° × 0.5° GFS (Global Forecast System; http://nomads.ncdc.noaa.gov/data.php)product from 2010 through present.This GFS dataset consists of both analysis and forecast (+7 days) data and is everyday automatically downloaded by the ATOP system.The third dataset is the global atmospheric reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF).This is a daily product at 1.5° × 1.5° resolution from 1979 through approximately the present minus 3 months; details are given in Dee et al. (2011) (http://www.ecmwf.int/products/data/archive/descriptions/ei/index.html).For spin-up simulations, the three datasets have been used to update the model through December 2011.In general, the CCMP and ECMWF give similar large-scale fields though dur-ing the time that it is available (July/02/1987 to December/ 31/2009), the higher-resolution CCMP products give more energetic currents on the shelves and in the marginal seas.The GFS product is used for simulations from January 2012 through present (hindcast and nowcast) and for forecasts through present +7 days.The dataset is interpolated horizontally onto the model grid and then also temporally during the model time-stepping.

Satellite Altimeter and McSST Data for Sea-Surfaceheight Anomaly (SShA)
Gridded sea-surface-height anomaly (SSHA) data is from AVISO (http://www.aviso.oceanobs.com/duacs/).This dataset has a temporal resolution of 7 days and spatial resolution of 1/3° × 1/3° (Le Traon et al. 1998), and is available from October/1992 through present.This dataset is used for data-assimilation (see section 4).Finally, SST data is from AVHRR MCSST (AVHRR-Advanced Very High Resolution Radiometer, Multi-Channel Sea Surface Temperature; http://gcmd.nasa.gov/records/GCMD_NAVOCEANO_MCSST.html), which has a resolution of 0.1°.The MCSST data may be used also for data assimilation using a similar statistical correlation method as that used for assimilating the SSHA in section 4. Its effects on mesoscale eddy features are relatively minor however (Yin and Oey 2007), and we use it to relax the model's surface temperatures using an e-folding time constant of 1 day -1 .

Other Data
Other data such as tides, rivers, surface heat and E-P (evaporation minus precipitation) fluxes can be specified in mpiPOM.These variables will be systematically added into the present version of ATOP in the near future.Farther along, we plan to add ice and sediment components as well into mpiPOM; in which case ice concentration data as well as coastal sediment input data will also be needed.

ATOP InITIAl & BOunDARy cOnDITIOnS AnD OTheR SPecIFIcS
Two different domains have been tested in the present work: a large-grid domain that encompasses the North Pacific Ocean, 98°E -73°W and 16°S -70°N, at 0.1° × 0.1° horizontal resolution and 41 vertical sigma levels (the "Pac10" domain), and a smaller western North Pacific China Seas (WPacCS) domain 99°E -131°W and 12°S -42°N which at present has the same 0.1° × 0.1° horizontal resolution but with 31 vertical sigma levels.In the near future, the resolution of this WPacCS model will be increased: 5/3 km resolution and 61 vertical sigma levels.The North Pacific domain is shown in Fig. 4 superimposed on an example of ATOP-simulated SST on July/06/1988.Shown are also other nested domains which we plan to implement in the future.
Prior to the ATOP hindcast and forecast runs, we have conducted over 40 non-assimilated (i.e., free-running) experiments on the "Pac10" domain.Each of these experiments uses either the ECMWF (1979-2011), the CCMP (1987-2009), or the CCMP + GFS (1987 -2009 + 2010 -2011) forcing.Some experiments use steady-state mean wind forcing, and additional idealized experiments specify zonal mean wind only.Results of these free-running experiments will be reported separately.The present paper focuses on ATOP hindcasts and forecasts beginning in January/2012.We use one of the CCMP + GFS experiments to initialize the forecast in January/2012.The free-running experiments are used, however, to obtain eddy statistics such as the correlations between modeled SSHA and subsurface temperature and salinity fluctuations (see below).These correlations are then used in the data assimilation scheme to be described in section 5.
Boundary conditions are zero normal flux (of any kind) across solid boundaries.No-slip conditions (for velocity) are specified along the land-sea boundaries. 1 At the ocean's bottom, stresses are computed by matching the near-bottom velocity to the law-of-the-wall logarithmic profile.Along the open boundaries WOA climatological T and S are specified within 3°-wide (for Pac10; 1.5°-wide for WPacCS) flow-relaxation zones (Oey and Chen 1992a, b).Transports from SODA (Giese and Ray 2011) are specified together with the Flather-radiation scheme.To prevent temperature and salinity drift in deep layers in long-term integration, the T and S for z < -1000 m are restored to monthly climatological values with a time scale of 720 days.
At the sea-surface, the SST is relaxed to monthly WOA values or MCSST when the latter is available using an e-folding time constant of 1/30 day -1 for the free-running spin-up runs and 1 day -1 for the hindcast and forecast runs (when MCSST is then used).The SSS (sea-surface salinity) is similarly relaxed to climatological WOA values.Future plan will combine these relaxations with specifications of surface fluxes (e.g., from NCEP reanalysis and GFS forecast) as described in Oey and Chen (1992a).To calculate wind stresses, we use a bulk formula with a high wind-speed limited drag coefficient that curve-fits data for low-to-moderate winds (Large and Pond 1981) and data for high wind speeds (Powell et al. 2003 where ua is the wind speed. 2 According to this formula, C d is constant at low winds, linearly increases for moderate winds, reaches a broad maximum for hurricane-force winds, ua ≈ 30 ~ 50 m s -1 , and then decreases slightly for extreme winds.It is necessary to use a C d formula that accounts for high winds in order to simulate the ocean response due to typhoons.Donelan et al. (2004) suggest that the C d -leveling 1 Because POM uses C-grid the zero tangential velocity is actually specified 1 2 -grid point inside the wall.The no-slip may be considered as being 1 2 -slip.
2 This same formula was used in Oey et al. (2006), but the coefficient for u 2 a was erroneously rounded off to 0.0002.at high wind may be caused by flow separation from steep waves.Moon et al. (2004) found that C d decreases for younger waves that predominate in hurricane-forced wave fields.Bye and Jenkins (2006) attribute the broad C d -maximum to the effect of spray, which flattens the sea surface by transferring energy to longer wavelengths.The free-running experiments are initialized from diagnostic runs in which the WOA T and S fields are held fixed and the corresponding velocity (which in the ocean's interior is in near-geostrophic and thermal-wind balances); zero surface stresses and bottom logarithmic layer are used for these diagnostic runs.The hindcast and forecast runs begin in January/2012 and are then initialized using the result from one of the free-running experiments.The hindcast is relatively insensitive to the initial condition, since with data assimilation it quickly adjusts (in about 15 ~ 20 days) to the satellite observations.The ATOP's real-time hindcast and forecast cycle, to be described in the next section, then begins on January/18/2012.

DATA ASSIMIlATIOn MeThOD AnD ATOP'S AnAlySIS-FORecAST cycle
Satellite data are assimilated into the model following the methodology given in Mellor and Ezer (1991) and Ezer and Mellor (1994).In this method, SSHA is projected into the subsurface temperature field using pre-computed correlation factors derived from the long-term (20 ~ 30 years) free-running experiments described in the previous sections (ECMWF or CCMP + GFS).These long-term integrations ensure that the model has yielded a statistical equilibrium eddy field, as judged for example from the leveling-off of the plot of basin eddy kinetic energy (EKE) with time.The last 10-year outputs are then used to compute the correlation factors.Thus the resulting temperature anomaly (δT) is (<·> is time-averaging, and T is the potential temperature): , , ,  , ,  , ,  T x y z t F x y z x y t where the correlation factor is (dh = model SSHA) and the corresponding correlation coefficient is (5.2b) Ezer and Mellor (1994) assimilate along-track o dh data assuming a linear-saturation error growth model for the firstguess error.Our experience has been that if AVISO o dh maps are assimilated the following simplified formula (see Wang et al. 2003) suffices: where T is the model (first-guess) temperature, T a denotes the analysis temperature, R A is the ratio of the assimilated time step Δt A to the de-correlation time scale Δt E of the model eddy field, and T O is the 'observed' temperature inferred from Eq. (5.1), Instead of using the model mean for < T > in Eq. ( 5.4), our past experience has been that setting < T > = T C , the observed temperature climatology, helps to control long-term (~ 10 years) drift in the model.Formula (5.3) assumes that the AVISO map errors are small compared to the model errors, and that Δt A << Δt E .We set Δt A = 1 day (Yin and Oey 2007).The Δt E is estimated from the above-mentioned free-running model runs and is ≈ 20 days in western North Pacific regions comparable to also the 20 day-value used by Ezer and Mellor's (1994) for the Gulf Stream, and the 30 day-value used by Yin and Oey (2007) for the Loop Current and eddies in the Gulf of Mexico.The assimilation (5.3) is such that T a ≈ T O in regions where the correlation is high (C T 2 ≈ 1), but T a ≈ T where the correlation is low.Yin and Oey (2007) also use a similar assimilation of SST using Eq. ( 5.3) with C T and F T replaced by the corresponding functions that use d(SST) in place of dh in Eq. (5.2).However, in this work, a simpler method is used in which the SST is relaxed to MCSST as described previously.The assimilation is turned off for isobaths < 1000 m and in the relaxation zones next to the open boundaries.Over the shelves, then, flow dynamics are purely wind-driven (and tide-and riverdriven if these are also included).
Figure 5 shows plots of C T at three z-levels in the North Pacific.The C T is large (> 0.7) in the Kuroshio and northern South China Sea regions, as well as in some parts of the North Pacific.The relatively large values near the southern boundary are an artifact of the relaxation zone used there, but they do not affect the results since there is no data assimilation in the relaxation zones.The values are larger than 0.5 in the Philippines Sea. Figure 6 shows typical profiles of area-averaged C T in the northern South China Sea and around Taiwan.The C T is larger near the surface, and then decays with depth.Its maximum is typically just below the surface, away from the wind-driven Ekman layers (c.f.Mellor and Ezer 1991).
Figure 7 shows ATOP's real-time hindcast and forecast cycle.Each day at specified times, the ATOP system automatically downloads satellite SSHA and MCSST data, as well as the GFS surface flux data.A hindcast/nowcast    analysis is then conducted from -7 day to present to update the initial condition (the "restart" file) for forecast.A +7 day ocean forecast is then conducted using the +7 day GFS atmospheric forecast wind field.Since the eddy field changes relatively little in +7 days, we plan to extend the forecast beyond the +7 day period using the climatological wind field in the near future.

MeThODS OF AnAlySIS: FOR ASSeSSIng FORecAST SkIllS
To assess model skill, we compare modeled and observed (i.e., AVISO) SSH (h) fields using various metrics.First we evaluate how well eddy fields are simulated at a given time.We split h into two components: a large-scale (i.e., "smoothed") <h> and a small-scale h': We define the smoothing operator < F > for the variable 'F' to be a simple spatial averaging over P° × Q° (longitude × latitude) rectangles: and for simplicity we have taken P° = Q°.The 'i' is a running index from 1 through N which is the total number of points in each rectangle (square).For F = h, the \ h then is the deviation from the mean SSH within each square: it represents the deformation of the sea surface due to eddies with scales smaller than (approximately) P° (Fig. 8).

RMS
The root mean squared error RMS is then: The second form shows that the squared RMS error of the total h (i.e., RMS 2 h ) is the sum of the squared differences of the mean and the eddy fields.Equation ( 6.3) suggests that we define the RMS for the eddy field as: (6.4)

correlation
We define correlation in each P° × Q° rectangles; by the definition of correlation: The smallness of RMSh and RMS \ h , and high C o m h h (≈ +1) suggest good hindcast and forecast.

FORecAST ReSulTS
An example of the forecast results (SSH and surface currents) for May/28/2012 are shown in Fig. 9, and the corresponding RMS error RMS \ h and correlation C o m h h are shown in Fig. 10.There are visual correspondences between the two fields -for example in South China Sea and off the eastern coast of Taiwan and Philippines.The model appears to underestimate the strengths of eddies and their spatial scales.These underestimations are expected as the model resolution is still relatively coarse.Some additional sensitivity experiments with model parameters may also be needed.Nonetheless, the correlations between forecast and AVISO ranges from 0.4 ~ 0.8 in deep-water regions where also the RMS errors are generally less than 0.1 m (≈ twice the standard error of the altimetry measurements).Winddriven variability, which the model simulates, dominates over the shelves, and the agreements between model and AVISO are poor over these shallower waters.This is because AVISO has large errors there and is not expected to account for the wind-driven currents.
The results of these ATOP forecasts and error analyses are updated daily on our website: http://mpipom.ihs.ncu.edu.tw/index.php.Forecasts on both the Pac10 North Pacific domain, as well as on the WPacCS (western North Pacific China Seas) domain are posted, and the reader is welcome to browse and download the corresponding plots.Forecast products which may be useful during the typhoon season are SST and the depth of 26°C isotherm (D26); the latter provides a measure of the upper-ocean heat content that is potentially available to fuel tropical cyclones (see e.g., Oey et al. 2006Oey et al. , 2007, where other references are also found).Figure 11 shows ATOP SST and D26 on August/03/2012.In the region east of the Philippines, between 10 ~ 20°N and west of the dateline, the D26 exceeds 100 m, and SST > 30°C.East of Taiwan, the D26 exceeds 50 m, and large D26-values extend along the Kuroshio over the East China Sea continental slope and outer shelf, and also into the entire SCS.These values suggest that, in 2012 typhoon season, the D26 may be anomalously deeper than climatology in most parts of the western North Pacific (see, e.g., Lin et al. 2008, their Fig. 2).

SIMulATeD geneRAl cIRculATIOn, kuRO-ShIO InTRuSIOn InTO ScS, AnD The STcc
To close this paper, we present some results obtained from the unassimilated, free-running mpiPOM.We begin by describing the model's results in the equatorial region (Fig. 4); the model response here develops very quickly, within 1 year after the model's initialization.Figure 4 shows the modeled SST in June of 1988 during the developing stage of the strong 88 -89 La Niña.Over the eastern equatorial Pacific, surface waters as cold as 19°C can be seen to emanate from the central South American (Peruvian)  coast to central Pacific as far west as the dateline.In the model, these processes are forced by the anomalously strong upwelling-favorable wind along the Peruvian coast and the equatorial trade wind (Philander et al. 1989).The simulated SST patterns and numerical values agree well with observation (http://www.pmel.noaa.gov/tao/elnino/lanina-story.html#).Tropical instability waves are also seen on both sides of the cold SST tongue, with wavelengths of approximately 1000 km (Legeckis 1977;Philander et al. 1986;Willett et al. 2006).
The unassimilated run was initialized from a 22-year (1988 -2009) forced by the CCMP wind, and then repeated for another 11 years.The last 8-year (1991 -1998) is used below.Figure 12a shows model mean (8-year: 1991 -1998) near-surface currents and temperature from which various features of the western North Pacific and equatorial Pacific  are seen (see the Introduction; also Myers and Weaver 1996;Qiu 1999;and Chang and Oey 2012 and references cited therein).These include: (i) westward-flowing NEC that bifurcates near the 12 ~ 13°N latitude off the eastern coast of the Philippines; (ii) the Mindanao (anticlockwise) and Halmahera (clockwise) "twin" eddies off the Philippines-Indonesian archipelagos; (iii) the eastward North Equatorial Counter Current (NECC) near 5°N latitude; (iv) the New Guinea Coastal Current (NGCC) that flows west-northwestward along the New Guinean coast; (v) Pacific-Indonesian Throughflow that flows into the Makassar Strait; and (vi) the Kuroshio that flows northward along the eastern coasts of Philippines and Taiwan, along the East China Sea continental slope, and finally separating off the eastern coast of Japan.
We now focus on two circulation features of the western North Pacific (i) northern South China Sea, and (ii) the STCC.We compare the model results against observations and models in the literature, as well as against satellite observations.

northern South china Sea circulation
Part of the Kuroshio intrudes into the northern SCS through the Luzon Strait at approximately the 20 ~ 21°N, which is also approximately the latitude at which nearsurface drifters tend to intrude into SCS (Centurioni et al. 2004).An excellent description is given by Hsueh and Zhong (2004).A tongue-like patch of warm water is seen (Fig. 12a) to spread westward over the upper continental slope of the northern SCS, similar to the β-plume that is seen in the model of Sheremet (2001).The intruded flow reaches westward to 118°E (east of the Dongsha Island at 20.7°N, 116.7°E), before turning northward and then northeastward following the shelfbreak.This description agrees well with Liang et al's (2003) 10-year composite of Shipboard ADCP measurements (their Figs. 4 and 12).The modeled anticyclonic eddy has a maximum speed of ≈ 0.3 m s -1 on its northern (i.e., shoreward) side, which compares well with the Liang et al.'s observed speeds (see their Fig. 4).The shelfbreak current extends west to the Hainan Island (19°N, 110°E), and Hsueh and Zhong (2004) named the portion west of Dongsha the South China Sea Warm Current (SCSWC).The authors also pointed out the existence of the westward-flowing splinter current (the South China Sea Branch of Kuroshio or SCSBK) which lies further offshore over the upper continental slope.In Fig. 12a, this SCSBK is shoreward of the northern arm (westward-flowing along 19°N) of the annual-mean SCS cyclonic gyre (Qu 2000).The SCSBK appears to feed the SCSWC, presumably by onshore flow forced by the along-slope pressure gradient resulting from the tongue-like patch of warm water mentioned above (Hsueh and Zhong 2004).In the model, the SCS gyre appears to be separate from the Kuroshio intrusion.
Consistent with the Hsueh and Zhong's (2004) account based on in situ measurements, summarized above, many models in the literature also reproduce the anticyclonic turning of the Kuroshio intruded flow, of varying intrusion distances and eddy strengths.Table 1 summarizes the results in terms of the Kuroshio's intrusion turning longitude (KuITL) and the eddy's maximum speed U eddy which usually occurs along the shelfbreak.Most models show the Kuroshio intruded flow with KuITL ≈ 118°E.Two exceptions are Gan et al's (2006) POM andJia andChassignet's (2011) dataassimilated HYCOM run which basically show no intruded flow: the modeled Kuroshio shows only a small-amplitude clockwise turn in the Luzon Strait.Their results appear to be consistent with Maximenko and Niiler's (2005) 1/2°-dataset but differ from Rio's (2009) 1/4°-dataset as well as Liang et al's (2003Liang et al's ( , 2008) ) Sb-ADCP composites.Other studies using satellite observations show intruded flow but do not show the anticyclonic eddy (e.g., Nan et al. 2011).However, satellite observations themselves are prone to large errors over the shallow shelves (which Nan et al. correctly excluded, for H < 200 m), thus making it difficult to infer the geostrophic current associated with the shelfbreak portion of the eddy.The U eddy is also fairly consistent amongst the models, approximately 0.3 m s -1 as observed by Liang et al. (2003Liang et al. ( , 2008)).The exceptions are Metzger and Hurlburt (2001) and Chern and Wang (2003) whose values are weaker ≈ 0.1 m s -1 , and Hsin et al. (2012) whose value is strong ≈ 0.9 m s -1 .

The STcc
Observations have shown that, east of Taiwan and the Luzon Strait, at approximately 19 ~ 25°N latitudes, a weak (a few cm s -1 ) and shallow (near the surface approximately 100 m) eastward current penetrates into the open Pacific for thousands of kilometers, from 130 ~ 180°E.The STCC was first predicted by Yoshida and Kidokoro (1967a, b) and subsequently confirmed by numerous observations (e.g., White et al. 1978).The STCC exists in a region where one would normally expect sluggish westward flow of the Sverdrup's anticyclonic gyre.More recent observations show the existence of multiple (3) subtropical fronts, between 19 ~ 21°N and also 21 ~ 24°N, as well as along 26°N further east of 170°E (Aoki et al. 2002;Kobashi et al. 2006).STCC and eddies play an important role in the Kuroshio transport and intrusion into the SCS; it may also play an important role in the western North Pacific climate (Chang and Oey 2012, and references quoted therein).For a circulation model of the western North Pacific to be useful, it is necessary that STCC exists in the model.
The existence of multiple-front STCC is seen in the mpiPOM simulation (Fig. 12b).Two fronts are seen west of 170°E, one at 19 ~ 21°N and another one at 21 ~ 24°N, and a third front is seen at 26°N east of 170°E, in good agreements with the observations.The model fronts are stronger in summer than in winter, as demonstrated in Figs.13a  and b for March and September respectively.Both show a meandering STCC, but the amplitude is larger in summer (Fig. 13b) with stronger speeds.The meanders are signatures of STCC eddies (Qiu 1999) which also tend to be stronger in summer.The model also indicates that the southern front is more distinct in summer; its intensity especially west of approximately 140°E appears to be related to the movement of the NEC.The northern front is also stronger in summer, and it tends to coalesce with the eastern front east of 170°E.Figures 14a and b show the sectional plots of temperature and zonal velocity from model and observation at 137°E.In observation (Fig. 14b), the STCC appears smooth with only one front, probably because only 4 surveys per year were taken, and features are smoothed out after averaging from 1993 -2008.In the model, two fronts are discernible in this 8-year mean plot.The modeled STCC's location, width and vertical extents are nonetheless consistent with the observation.Figures 14b and c show the same sectional plots for model's monthly means in March and September, respectively (corresponding to Figs. 13a and b).These show that in winter (Fig. 14c), cooler waters occupy the near-surface 100 m and the 24 ~ 26°C isotherms thicken in the region between 17 ~ 23°N.The NEC shifts northward during the positive phase of the so called seasonal PTO that we previously described (Chang and Oey 2012; see their Fig. 2 and discussions), when the windstress curl distribution is such that surface convergence tends to strengthen in the latitudes of 17 ~ 23°N.Only one front then exists.In summer (Fig. 14d), the near-surface waters become warm, NEC shifts southward, isotherms shallow, and multiple fronts form (Fig. 14d).In addition to NEC, these descriptions need to be further analyzed in a future study in connection also with the southward progression of the mode waters (Kubokawa and Inui 1999;Aoki et al. 2002), as well as with the seasonal (and interannual) eddy variability of the STCC.The point here is that realistic subtropical front(s) exist in the ATOP simulation, and this is significant given the im-Table 1. Kuroshio's intrusion (clockwise) turning-longitude (KuITL; °E) measured at 21°N in the northern South China Sea and the eddy's maximum speed (U eddy ; m s -1 ), from literature and from this study (last entry).Speeds are read from the indicated figures of the quoted references, and are very approximate.The 1 st 5 entries (shaded) are observations; remaining entries are models.Except for 1/12° HYCOM (Jia and Chassignet 2011) which assimilates altimetry SSHA (and other observations), all other models are without assimilation.MDT = Mean Dynamic Topography.Most values are long-term (years) or composite mean, but some are shorter-term mean (e.g., Liang et al. 2008).
portance of the STCC and eddies to the circulation of the western North Pacific and its marginal seas (Chang and Oey 2012).We plan in the near future to conduct more extensive sensitivity experiments to further understand and quantify their dynamics.

cOncluSIOnS AnD FuTuRe PlAnS
A Taiwan ocean prediction system (ATOP) has been developed and since February 2012 has produced real-time forecasts of the wind and eddy-driven ocean circulation fields in the North Pacific Ocean and the seas around Taiwan.Forecast products are updated daily, and are available from http://mpipom.ihs.ncu.edu.tw/index.php.This paper describes details of the ATOP system and preliminary validations against satellite observations.To further examine model's fidelity, we also compare the unassimilated model results with those available in the literature and also with AVISO data.The free-running model captures many wellknown circulation features of the western North Pacific.In particular, the northern SCS circulation that results from the Kuroshio intrusion through the Luzon Strait compares well  with available long-term observations, and is generally consistent with many existing model results.Finally, the freerunning model can simulate the STCC and eddies.
In the near future, we plan to incorporate tides and rivers into the model, and to incorporate other observational data for assimilations, e.g., ARGOs and drifters.More comprehensive skill assessments are also being planned, as well as implementations of more sophisticated data-assimilation schemes.

Fig. 1 .
Fig. 1.Western North Pacific Ocean's topography (colors, in meters) and schematic sketches of major currents and regions of eddies.

Fig. 4 .
Fig. 4. ATOP North Pacific domain at 0.1° × 0.1° (~10 km) horizontal resolution, and other planned nested domains at (i) 0.05° × 0.05° (~5 km) that includes all the western North Pacific marginal seas; (ii) 0.017° × 0.017° (~5/3 km) that includes the China Seas and straits and seas in the Indonesian and Philippines archipelagos (the "WPacCS" domain); and (iii) highest resolution 555 m for the waters surrounding Taiwan.These domains are shown superimposed on the ATOP sea-surface temperature (SST) on June/06/1988 during the developing stage of the strong 88 -89 La Niña.

Fig. 6 .
Fig. 6.Vertical profiles of the correlation coefficients C T T < > < >< > / T 2 2 1 2 d dh d dh = ^h averaged over the three regions in the northern South

Fig. 8 .
Fig. 8. Large scale field (large blue transparent circle) < > h has scales ~P (or Q), while eddy field (red and blue filled circles) \ h has scales < P, = P/5 shown here.At any point, < > \ h h h = + and \ h is therefore the deviation of SSH from a large-scale field < > h that gently varies in space.

Fig. 9 .
Fig. 9.An example of ATOP forecast SSH and surface (u, v) on May/28/2012 (a), and comparison against AVISO SSH and geostrophic (u, v) on the same date (b).For the latter, the geostrophic (u, v) vectors near the equator are omitted.On both plots, purple vectors are winds which for (a) are GFS +7 day forecast for May 28 while for (b) they are updated GFS wind analysis for the same day -they are slightly different.
Fig. 10.RMS error ( RMS \ h , color, blank areas are for values less than 0.1 (twice the standard error of altimetry measurements), and they are not plotted) between ATOP forecast SSHA and AVISO SSHA, and correlation (C o m h h , contours) between them on May/28/2012.

Fig. 12 .
Fig. 12. Modeled mean T (°C, color) and (u, v) vectors at z = -25 m for (a) southwestern North Pacific including straits in the Indonesian and Philippines archipelagos; and (b) the STCC region [rectangle in (a)] shown for vectors with eastward zonal component only; note change in vector scale from (a).

Fig. 13 .
Fig. 13.Modeled monthly mean T (°C, color) and (u, v) vectors at z = -25 m for the STCC region (rectangle in Fig. 12a) shown for vectors with eastward zonal component only, for (a) March, and (b) September.

Fig. 14 .
Fig. 14.Mean temperature (contours; °C) and zonal velocity (color; cm s -1 ) for (a) model and (b) observation from Japan Meteorological Agency repeat hydrographic surveys of 1993 -2008; 4 surveys per year at the nominal 1° latitude spatial resolution, from Qiu and Chen (2010).The geostrophic velocity in (b) is referenced to 1000 dbar and dashed lines denote the zero velocity contours.Bottom panels show model monthly mean plots for (c) March and (d) September.