ROMS Forcing Package

This document contains information about how to generate surface and bottom boundary conditions for ROMS. The forcing package can be used to generate the forcing NetCDF file. It also can be used to modify and/or append fields to an existing forcing NetCDF file. The source code can be obtained as tar file from ROMS web site.


Basically, ROMS needs surface and bottom boundary conditions for momentum and tracers-type variables (temperature, salinity, etc). Other fields are also needed when using boundary layers formulations at the bottom, ocean-atmosphere, and sea-ice interfaces. However, it is customary to provide only forcing data for momentum and tracers surface boundary conditions. The data may include the net surface flux or all necessary fields to compute such flux. The bottom fluxes are usually set or computed internally.


Currently, the forcing data can be of two types: point or grided data time series. If point data, the same value is applied to each horizontal grid point in the model during execution. This is useful in idealized simulations or when there is data from only one buoy or meteorological tower station. The input file format can be also of two types: ASCII or NetCDF. If ASCII, the User needs to provide a file for each forcing field; vector components should be included in the same file. The input grided ASCII data is assumed to be in a regular grid (spacing in each direction can be different) and this package will interpolate to the application grid. The ASCII format is described below. If NetCDF, the data is assumed already mapped to the User application grid. Usually, the objective analysis package is used to map data to the model horizontal curvilinear grid using statistical correlation functions. For more details, see OA package documentation. The input forcing data can be contained in one or several NetCDF files.


In this package, the reading and writing of data is carried out in in routine wrt_out.F. The generic routine get_grdfld.F reads and interpolates input forcing grided data from ASCII files whereas get_ncfld.F reads forcing data from NetCDF file(s). All the information for each forcing field is managed in this package by several arrays with a unique index descriptor, as shown in Table 1.


Table 1
ID Index Forcing Field Variable Time Units
1 idUsms surface U-momentum stress sustr sms_time Pa (N/m2)
2 idVsms surface V-momentum stress svstr sms_time Pa (N/m2)
3 idUair surface U-wind component Uair wind_time m/s
4 idVair surface V-wind component Vair wind_time m/s
5 idWamp wind-induced wave amplitude Awave wave_time m
6 idWdir wind-induced wave direction Dwave wave_time degree
7 idWper wind-induced wave Period Pwave wave_time s
8 idSden bottom sediment grain density Ssize - m
9 idSsiz bottom sediment grain diameter size Sdens - kg/m3
10 idTdew surface dew-point temperature Tdew tdew_time Kelvin
11 idTair surface air temperature Tair tair_time Celsius
12 idPair surface air pressure Pair pair_time mb
13 idQair surface air relative humidity Qair qair_time percentage
14 idCfra cloud fraction cloud cloud_time 0 to 1
15 idrain rain fall rate rain rain_time kg/m2/s
16 idSrad net shortwave radiation flux swrad srf_time Watts/m2
17 idLrad net longwave radiation flux lwrad lrf_time Watts/m2
18 idShea net sensible heat flux sensible sen_time Watts/m2
19 idfSSS sea surface salinity SSS sss_time PSU
20 idfSST sea surface temperature SST sst_time Celsius
21 iddQdT surface net heat flux sensitivity to SST dQdSST sst_time Watts/m2/Celsius
22 idTsur(itemp) net surface heat flux shflux shf_time Watts/m2
23 idTsur(isalt) net surface freshwater flux swflux swf_time cm/day


These indices are initialized in routine init_vars.F. This design allows the easy inclusion of new fields and update of this package in the future. This list of indices also shows the forcing fields presently managed by this package.

Input ASCII Format

The input ASCII format for each forcing field is as follows: 

      NTIME(id....)   Length of time series (number of records).
      NLON (id....)   Number of longitude points.
      NLAT (id....)   Number of latitude points.
      TIME            Time of the data (days, Julian days).
      LON             Longitude of the grided data, LON(NLON).
      LAT             Latitude  of the grided data, LAT(NLAT).
      VAL             Scaler forcing field, VAL(NLON*NLAT).
      VALX            X-component forcing field, VALX(NLON*NLAT).
      VALY            Y-component forcing field, VALY(NLON*NLAT).
Point Time-Series Data 
         read (iunit(id....),*) NTIME(id....)
         do n = 1, NTIME(id....)
           read (iunit(id....),*) TIME, VAL            scalar data
           ...
      or
           read (iunit(id....),*) TIME, VALX, VALY     vector data
           ...
         enddo
Grided Time-Series Data 
         read (iunit(id....),*) NTIME(id....),
     &                          NLON (id....), NLAT (id....)
         read (iunit(id....),*) (LON(i), i = 1, NLON(id....))
         read (iunit(id....),*) (LAT(j), j = 1, NLAT(id....))
         do n = 1, NTIME(id....)
           read (iunit(id....),*) TIME
           NP = NLON(id....) * NLAT(id....)
           read (iunut(id....),*) (VAL(i), i = 1, NP)      scalar data
      or
           read (iunit(id....),*) TIME
           NP = NLON(id....) * NLAT(id....)
           read (iunit(id....),*) (VALX(i), i = NP)
                                                           vector data
           read (iunit(id....),*) TIME
           NP = NLON(id....) * NLAT(id....)
           read (iunit(id....),*) (VALY(i), i = NP)
        enddo
Notice that all reading is free formatted and the input grid must be ascending order: 
         LON(i) < LON(i+1)  and   LAT(j) < LAT(j+1)
and data must not have flagged values. Otherwise, the interpolation will not work correctly!. Users can modify this input ASCII format easily. The one provided above is just a guideline. Forcing Fields

The names of input and output NetCDF variables are initialized in routine init_vars.F. See array Vname. The variable names for each forcing field and its associated time are shown in Table 1.

Other names for time variables are allowed; provided that the field contains the time attribute specifying its associate time variable name. For example, a ncdump of a forcing file will give you: 

  dimensions:

        sms_time = 12 ;

  variables:

        double sms_time(sms_time) ;
                sms_time:long_name = "surface momentum stress time" ;
                sms_time:units = "day" ;
                sms_time:field = "sms_time, scalar, series" ;
                sms_time:cycle_length = 360.0 ;

        float sustr(sms_time, eta_u, xi_u) ;
                sustr:long_name = "surface u-momentum stress" ;
                sustr:units = "Newton meter-2" ;
                sustr:field = "surface u-momentum stress, scalar, series" ;
                sustr:time = "sms_time" ;

        float svstr(sms_time, eta_v, xi_v) ;
                svstr:long_name = "surface v-momentum stress" ;
                svstr:units = "Newton meter-2" ;
                svstr:field = "surface v-momentum stress, scalar, series" ;
                svstr:time = "sms_time" ;

  data:

        sms_time = 15, 45, 75, 105, 135, 165, 195, 225, 255, 285, 315, 345 ;


Notice that the variable sms_time associated with surface momentum stress has the attribute cycle_length with a value of 360.0 which indicates to the model that data must be used in a 360 days perpetual cycle. This is an example for monthly climatological forcing: every month is assumed to be 30 days long and there are 360 days in a year. Each monthly data is time centered at the 15th of each month. The model will interpolate between available records to get data at each timestep.

The net surface heat flux is the sum of all its constituents: 

            Q_net = shortwave + longwave + sensible + latent

                    (downward flux is positive: heating)
                    (upward   flux is negative: cooling)


The User can provide the net surface heat flux at input or activate the ocean-atmosphere boundary layer of the model to compute both the net heat flux and wind stress from atmospheric fields. This boundary layer is based on Liu et al. (1979) air-sea exchange parameterization and the model routine bulk_flux.F was adapted from the COARE code described in Fairall et al. (1996a, 1996b). The turbulent fluxes for wind, heat and moisture are computed using Monin-Obukhov similarity theory (Liu et al., 1979).

If the ocean-atmosphere boundary layer is activated, the User needs to provide the following atmospheric fields:



Here, surface means near surface sea level and it is usually at 10 m. If some of the above fields are not available, the User can set them internally in the model using the analytical options. For instance, if the daily heating cycle is not well resolved, the User can employ the provided global equation to compute incoming shortwave radiation from the albedo formulation, see subroutine ana_srfluxin ROMS. The net longwave radiation is always computed internally via the Berliand (1952) formula.

Some knowledge in ocean-atmosphere interactions is required before activating this boundary layer in the model. It is not recommended to use it as a black box!


Sometimes, very long climatological solutions present heating or cooling trends. It is common practice to compute surface boundary conditions for temperature and salinity involving relaxation to its observed values. If a correction to the heat flux is activated, the User needs to provide:



Then, the model will correct the net surface heat flux as: 
       Q_mod(:,:) = Q_net(:,:) + dQdSST(:,:) * [T_mod(:,:,N) - SST(:,:)]

If a relaxation to surface salinity is activated, the User needs to provide:


and the model will compute the net surface freshwater flux as: 
       F_mod(:,:) = - Tnudg(isalt) * Hz(:,:,N) * [S_mod(:,:,N) - SSS(:,:)]
were Tnudg is a nudging time scale and Hz(:,:,N) is upper level thickness. For a description of these surface corrected fluxes see, for example, Haney (1971) and Killworth et al. (2000).

If the K-profile parameterization (KPP) for vertical mixing or the vertical heat distribution (CPP: SOLAR_SOURCE) is activated in the model, the User needs to provide:


to compute surface radiative buoyancy flux in the KPP (Large et al., 1994), or to penetrate heat in the water column using an exponential decay function in Jerlov water type.

Finally, if the wave/current/sediment boundary layer formulation is activated to compute bottom stress, the User needs to provide or set the following fields:


Tidal Forcing

ROMS (version 1.7.0 and up) can be forced at the open boundaries with tides. The forcing data for each tidal angular component consists of:




The tidal data can be appended to an existing forcing NetCDF file. See Table 2 for variable names and units.

Table 2
Forcing Field Variable Record Units
tidal angular period tide_period tide_period hour
tidal elevation amplitude tide_Eamp tide_period m
tidal elevation phase angle tide_Ephase tide_period degree
tidal current phase angle tide_Cphase tide_period degree
tidal current inclination angle tide_Cangle tide_period degree
maximum tidal current tide_Cmax tide_period m/s
minimum tidal current tide_Cmin tide_period m/s


A CDL representation of the forcing NetCDF file with tide data using ncdump will give you: 

  dimensions:
        xi_rho = 111 ;
        eta_rho = 241 ;
        tide_period = 7 ;

  variables:

        double tide_period(tide_period) ;
                tide_period:long_name = "tide angular period" ;
                tide_period:units = "hours" ;
                tide_period:field = "tide_period, scalar" ;

        double tide_Ephase(tide_period, eta_rho, xi_rho) ;
                tide_Ephase:long_name = "tidal elevation phase angle" ;
                tide_Ephase:units = "degrees";
                tide_Ephase:field = "tide_Ephase, scalar" ;

        double tide_Eamp(tide_period, eta_rho, xi_rho) ;
                tide_Eamp:long_name = "tidal elevation amplitude" ;
                tide_Eamp:units = "meter" ;
                tide_Eamp:field = "tide_Eamp, scalar" ;

        double tide_Cphase(tide_period, eta_rho, xi_rho) ;
                tide_Cphase:long_name = "tidal current phase angle" ;
                tide_Cphase:units = "degrees" ;
                tide_Cphase:field = "tide_Cphase, scalar" ;

        double tide_Cangle(tide_period, eta_rho, xi_rho) ;
                tide_Cangle:long_name = "tidal current inclination angle" ;
                tide_Cangle:units = "degrees between semi-major axis and East" ;
                tide_Cangle:field = "tide_Cangle, scalar" ;

        double tide_Cmin(tide_period, eta_rho, xi_rho) ;
                tide_Cmin:long_name = "minimum tidal current" ;
                tide_Cmin:units = "meter second-1" ;
                tide_Cmin:field = "tide_Cmin, scalar" ;

        double tide_Cmax(tide_period, eta_rho, xi_rho) ;
                tide_Cmax:long_name = "maximum tidal current" ;
                tide_Cmax:units = "meter second-1" ;
                tide_Cmax:field = "tide_Cmax, scalar" ;


  data:

        tide_period = 25.82, 23.93, 12.66, 12.42, 12.00, 6.21, 4.14 ;


The above example shows data for seven angular periods: O1, K1, N2, M2, S2, M4, and M6. The tidal data is stored by angular period in any order since the net tidal forcing is a summation of all tidal components. Notice that the tidal fields are provided for the entire model grid and are specified at RHO-points.

In ROMS, the tidal forcing is supper-imposed on the sea surface elevation (ssh) and vertically integrated momentum (ubarclm, vbarclm) climatology arrays which are used in the Flather open boundary conditions. The tidal contributions to the climatology arrays are added in subroutine clm_tides. For better results, both the tidal elevation and currents need to be specified at the open boundary. This can be achieved by activating the CPP options SSH_TIDES and UV_TIDES, and the climatology CPP options ZCLIMATOLOGY and M2CLIMATOLOGY. The tidal forcing requires special open boundary conditions. The Flather CPP options need to be activated for 2D momentum and the Chapman conditions for free-surface.

The tidal forcing data in Table 2 can be obtained from tidal models. If this is the case, a conversion is needed from tidal current amplitude and phase lag fields and current ellipse fields before it can be used in ROMS. This processing is done using several Matlab scripts written by Zhigang Xu (Bedford Institute of Oceanography, Canada):



He has written a nice manual dicussing the theory and programs for the rotary decomposition of tidal ellipse parameters. A tar file containing the above and other Matlab scripts is provided. The script wrt_tides.m can be used to append tidal data to an existing forcing NetCDF file. The scripts kagon_tides.m and njb_tides.m show how to proccess tidal data from ADCIRC tidal data in the US East Coast. The ADCIRC data and ploting scripts was kindly provided by Rich Signell (USGS). The triangular interpolation package was written by Charles Denham (USGS). This package is used to interporlate ADCIRC data to the application grid.

References
  1. Fairall, C.W., E.F. Bradley, D.P. Rogers, J.B. Edson and G.S. Young, 1996a: Bulk parameterization of air-sea fluxes for tropical ocean global atmosphere Coupled-Ocean Atmosphere Response Experiment, J. Geophys. Res., 101, 3747-3764.


  2. Fairall, C.W., E.F. Bradley, J.S. Godfrey, G.A. Wick, J.B. Edson, and G.S. Young, 1996b: Cool-skin and warm-layer effects on sea surface temperature, J. Geophys. Res., 101, 1295-1308.


  3. Haney, R.L., 1971: Surface thermal boundary conditions for ocean circulation models, J. Phys. Oceanogr., 1, 241-248.


  4. Killworth, P.D., D.A. Smeed, and A.J. George Nurser, 2000: The effects on ocean models of relaxation toward observations at the surface, J. Phys. Oceanogr., 30, 160-174.


  5. Large, W.G., J.C. McWilliams, and S.C. Doney, 1994: A Review and model with a nonlocal boundary layer parameterization, Reviews of Geophysics, 32,363-403.


  6. Liu, W.T., K.B. Katsaros, and J.A. Businger, 1979: Bulk parameterization of the air-sea exchange of heat and water vapor including the molecular constraints at the interface, J. Atmos. Sci., 36, 1722-1735.




IMCS, Ocean Modeling Group, by Webmaster.