表面强迫2：使用NCEP reanalysis II大气强迫来创建表面强迫

NCEP的表面强迫

NCEP Reanalysis II数据提供了许多产品，可用于驱动FVCOM中的表面施力。他们的网页是http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.html。可通过OPeNDAP服务器访问数据。MATLAB fvcom-toolbox包含许多功能，可将数据从NCEP下载并插值到给定的FVCOM非结构化网格上并输出到netCDF。

下载和内插NCEP数据

NCEP Reanalysis II数据提供以下参数：

变量名	变量简称	单位
u风分量	uwnd	多发性硬化症
v风分量	vwnd	多发性硬化症
向下的长波辐射面	dlwrs	瓦/米2
净短波辐射面	nswrs	瓦/米2
气温	air	摄氏温度
相对湿度	rhum	％
降水率	prate	多发性硬化症
海平面压力	pres	帕
潜热通量	lhtfl	瓦/米2
显热通量	shtfl	瓦/米2
潜在蒸发速率	pevpr	瓦/米2

下载数据的过程如下（也请参见下面的示例脚本）。

1. 读入SMS网格（read_sms_mesh）。
2. 提取开放边界节点（add_obc_nodes_list），add_obc_nodes_list如果不应用平均流，则将的参数中的ObcType赋值为1 ，否则为ObcType选择2。有关ObcType选项，请参见FVCOM手册中的表6.1。
3. 使用get_NCEP_forcing下载NCEP II再分析数据模型域的程度和模型的时间。
4. grid2fvcom使用三角形插值将规则网格化的NCEP数据插值到非结构化FVCOM网格上（有关更多信息，请参见MATLAB的TriScatteredInterp函数）。
5. write_FVCOM_forcing 输出必要的netCDF文件以进行表面施力，其中包括加热，风和降水/蒸发。

示例脚本

 % 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);