Surface forcing from NCEP

NCEP Reanalysis II data provide a number of products which can be used to drive surface forcing in FVCOM. Their webpage is http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.html. Data is accessible through an OPeNDAP server. The MATLAB fvcom-toolbox includes a number of functions to download and interpolate data from NCEP onto a given FVCOM unstructured grid and output to netCDF.

The NCEP Reanalysis II data provide the following parameters:

Variable nameVariable shortnameUnits
u wind component uwnd m/s
v wind component vwnd m/s
Downward longwave radiation surface dlwrs W/m2
Net shortwave radiation surface nswrs W/m2
Air temperature air celsius
Relative humidity rhum %
Precipitation rate prate m/s
Sea level pressure pres Pa
Latent heat flux lhtfl W/m2
Sensible heat flux shtfl W/m2
Potential evaporation rate pevpr W/m2

The process for downloading the data is as follows (see example script below too).

  1. Read in an SMS mesh (read_sms_mesh).
  2. Extract open boundary nodes (add_obc_nodes_list), assigning a value of 1 to the ObcType in add_obc_nodes_list's arguments if no mean flow will be applied, otherwise select 2 for the ObcType. See Table 6.1 in the FVCOM manual for the ObcType options.
  3. Use get_NCEP_forcing to download the NCEP Reanalysis II data for the extent of the model domain and for the duration of the model.
  4. grid2fvcom interpolates the regularly gridded NCEP data onto the unstructured FVCOM grid with a triangular interpolation (see MATLAB's TriScatteredInterp function for more information).
  5. write_FVCOM_forcing outputs the necessary netCDF file for the surface forcing, which includes heating, wind and precipitation/evaporation.

Example script

snippet.matlab
    % Script to load an FVCOM unstructured grid and interpolate
    % NCEP Renalysis 2 forcing data onto it.
 
 
    %%%------------------------------------------------------------------------
    %%%                          INPUT CONFIGURATION
    %%%
    %%% Define the model input parameters here e.g. forcing source etc.
    %%%
    %%%------------------------------------------------------------------------
 
 
    conf.base = '/data/medusa/pica/models/FVCOM/irish_sea/run/';
 
    % Which version of FVCOM are we using (for the forcing file formats)?
    conf.FVCOM_version = '3.1.6';
 
    % Case name for the model inputs and outputs
    conf.casename = 'irish_sea_v10';
 
    conf.coordType = 'cartesian'; % 'cartesian' or 'spherical'
    % Input grid UTM Zone (if applicable)
    conf.utmZone = {'30 U'}; % syntax for utm2deg
 
    % Sponge layer parameters
    conf.spongeRadius = 20000; % in metres, or -1 for variable width.
    conf.spongeCoeff = 0.00001;
 
    % Model time ([Y,M,D,h,m,s])
    conf.modelYear = 2000;
    conf.startDate = [conf.modelYear,01,01,00,00,00];
    conf.endDate = [conf.modelYear + 1,01,01,00,00,00];
 
    % NCEP surface forcing
    conf.doForcing = 'NCEP';
 
    %%%------------------------------------------------------------------------
    %%%                      END OF INPUT CONFIGURATION
    %%%------------------------------------------------------------------------
 
    % Convert times to Modified Julian Date.
    conf.startDateMJD = greg2mjulian(conf.startDate(1), ...
        conf.startDate(2), ...
        conf.startDate(3), ...
        conf.startDate(4), ...
        conf.startDate(5), ...
        conf.startDate(6));
    conf.endDateMJD = greg2mjulian(conf.endDate(1), ...
        conf.endDate(2), ...
        conf.endDate(3), ...
        conf.endDate(4), ...
        conf.endDate(5), ...
        conf.endDate(6));
    conf.inputTimeTS = conf.startDateMJD:conf.dtTS:conf.endDateMJD;
 
    % Read the input mesh and bathymetry. Also creates the data necessary for
    % the Coriolis correction in FVCOM.
    Mobj = read_sms_mesh(...
        '2dm',fullfile(conf.base,'raw_data/',[conf.casename,'.2dm']),...
        'bath',fullfile(conf.base,'raw_data/',[conf.casename,'.dat']),...
        'coordinate',conf.coordType,'addCoriolis',true);
 
    if strcmpi(conf.doForcing, 'NCEP')
        % Use the OPeNDAP NCEP script (get_NCEP_forcing.m) to get the following
        % parameters:
        %     - u wind component (uwnd) [m/s]
        %     - v wind component (vwnd) [m/s]
        %     - Downward longwave radiation surface (dlwrs) [W/m^2]
        %     - Net shortwave radiation surface (nswrs) [W/m^2]
        %     - Air temperature (air) [celsius]
        %     - Relative humidity (rhum) [%]
        %     - Precipitation rate (prate) [m/s]
        %     - Sea level pressure (pres) [Pa]
        %     - Latent heat flux (lhtfl) [W/m^2]
        %     - Sensible heat flux (shtfl) [W/m^2]
        %     - Potential evaporation rate (pevpr) [W/m^2]
        %     - Momentum flux (tx, ty) [Ns/m^2/s?]
        %     - Precipitation-evaporation (P_E) [m/s]
        %     - Evaporation (Et) [m/s]
        %
        % The script converts the NCEP data from the OPeNDAP sever from
        % longitudes in the 0 to 360 range to the -180 to 180 range. It also
        % subsets for the right region (defined by Mobj.lon and Mobj.lat).
        %
        forcing = get_NCEP_forcing(Mobj, [conf.startDateMJD, conf.endDateMJD]);
 
        forcing.domain_cols = length(forcing.lon);
        forcing.domain_rows = length(forcing.lat);
        forcing.domain_cols_alt = length(forcing.rhum.lon);
        forcing.domain_rows_alt = length(forcing.rhum.lat);
 
        % Convert the small subdomain into cartesian coordinates. We need to do
        % this twice because some of the NCEP data are on different grids (e.g.
        % sea level pressure, relative humidity etc.).
        tmpZone = regexpi(conf.utmZone,'\ ','split');
        [tmpLon, tmpLat] = meshgrid(forcing.lon, forcing.lat);
        [tmpLon2, tmpLat2] = meshgrid(forcing.rhum.lon, forcing.rhum.lat);
        [forcing.x, forcing.y] = wgs2utm(tmpLat(:), tmpLon(:), str2double(char(tmpZone{1}(1))), 'N');
        [forcing.xalt, forcing.yalt] = wgs2utm(tmpLat2(:), tmpLon2(:), str2double(char(tmpZone{1}(1))), 'N');
        clear tmpLon tmpLat tmpLon2 tmpLat2
        % Create arrays of the x and y positions.
        forcing.x = reshape(forcing.x, forcing.domain_rows, forcing.domain_cols);
        forcing.y = reshape(forcing.y, forcing.domain_rows, forcing.domain_cols);
        forcing.xalt = reshape(forcing.xalt, forcing.domain_rows_alt, forcing.domain_cols_alt);
        forcing.yalt = reshape(forcing.yalt, forcing.domain_rows_alt, forcing.domain_cols_alt);
 
        [forcing.lon, forcing.lat] = meshgrid(forcing.lon, forcing.lat);
 
        forcing = rmfield(forcing, {'domain_rows', 'domain_cols', ...
            'domain_rows_alt', 'domain_cols_alt'});
 
        tic
 
        % Interpolate the NCEP data to the model grid and output to netCDF.
        interpfields = {'uwnd', 'vwnd', 'pres', 'nswrs', 'prate', 'Et', 'time', 'lon', 'lat', 'x', 'y'};
        forcing_interp = grid2fvcom(Mobj, interpfields, forcing);
        if ftbverbose
            fprintf('Elapsed interpolation time: %.2f minutes\n', toc/60)
        end
    end
 
 
    %% Write out all the required files.
 
    conf.outbase = fullfile(conf.base,'input',conf.casename);
 
    % Make the output directory if it doesn't exist.
    if exist(conf.outbase, 'dir')~=7
        mkdir(conf.outbase)
    end
 
    forcingBase = fullfile(conf.outbase,conf.casename);
    write_FVCOM_forcing(Mobj, ...
        forcingBase, ...
        forcing_interp, ...
        [conf.doForcing, ' atmospheric forcing data'], ...
        conf.FVCOM_version);