Met Office Unified Model surface forcing

MATLAB toolbox steps

The MATLAB Tools:fvcom-toolbox provides a suite of functions to process Met Office Unified Model outputs to drive surface forcing in FVCOM. The necessary steps are outlined below.

1. Load the model grid (read_sms_mesh).
2. Convert to lat/long (if necessary).
3. Use get_MetUM_pp to download the Met Office PP files from the BADC FTP server. This requires a username and password from here.
4. Convert the downloaded PP files to netCDF with pp2nc_subset to cover the model domain.
5. Extract the variables of interest from the newly created netCDF files with read_MetUM_forcing.
6. Export the variables to a WRF format netCDF file with write_WRF_forcing which can be used directly in FVCOM (no separate interpolation necessary).

N.B. If your model is running in cartesian coordinates, then you need to compile FVCOM with the PROJ libraries (FLAG_6 = -DPROJ, PROJLIBS = ... and PROJINCS = ... uncommented in make.inc).

Sample MATLAB script

    %%%------------------------------------------------------------------------
    %%%                          INPUT CONFIGURATION
    %%%
    %%% Define the model input parameters here e.g. forcing source etc.
    %%%
    %%%------------------------------------------------------------------------
    
    
    conf.base = '/data/medusa/pica/models/FVCOM/irish_sea/run/';
    
    % 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 = 2011;
    conf.startDate = [conf.modelYear,01,01,00,00,00];
    conf.endDate = [conf.modelYear,02,01,00,00,00];
    
    % Surface forcing
    conf.doForcing = 'MetUM';
    conf.convsh = '/users/modellers/pica/Software/bin/convsh';
    conf.tcl = '/users/modellers/pica/MATLAB/fvcom-toolbox/utilities/subset.tcl';
    
    %%%------------------------------------------------------------------------
    %%%                      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.
    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);
    
    % Convert UTM to lat/lon.
    if Mobj.have_lonlat == 0
        utmZones = cellfun(@(x) repmat(x, length(Mobj.x), 1), conf.utmZone, 'uni', false);
        [Mobj.lat, Mobj.lon] = utm2deg(Mobj.x, Mobj.y, utmZones{1});
        Mobj.have_lonlat = true;
    end
    
    % Model domain +/- 3 degrees. Interpolation resolution = 0.11 degrees
    % ~= Met Office Unified Model resolution.
    xextents = [min(Mobj.lon) - 3, max(Mobj.lon) + 3, 0.11];
    yextents = [min(Mobj.lat) - 3, max(Mobj.lat) + 3, 0.11];
    
    if strcmpi(conf.doForcing, 'MetUM')
        % Get and convert PP to netCDF.
        met_pp_files = get_MetUM_pp([conf.startDateMJD, conf.endDateMJD], {'user', 'password'});
        met_files = pp2nc_subset(met_pp_files, conf.convsh, conf.tcl, xextents, yextents);
    
        % Extract the variables of interest from the new netCDF files.
        varlist = {'longitude', 'latitude', 't_1', 'sh', 'x-wind', ...
            'y-wind', 'lh', 'rh', 'precip', 'solar', 'longwave', ...
            'p_4', 'temp_3'}
        MetUM = read_MetUM_forcing(met_files, varlist);
    
    
        %%%--------------------------------------------------------------------
        % After that, the net heat flux has to be created as the sum of the
        % latent, sensible, shortwave and longwave heat fluxes. The 
        % evaporation can also be estimated from the latent heat (see the 
        % example create_files.m file in the fvcom-toolbox). The field
        % names in the MetUM struct also need to be changed to NCEP versions.
        % I've used the following little bit of code to do that:
        %%%--------------------------------------------------------------------
    
    
        replacements.longitude = 'lon';
        replacements.latitude = 'lat';
        replacements.time = 'time';
        replacements.p_4 = 'pres';
        replacements.longwave = 'nlwrf';
        replacements.solar = 'nswrf';
        replacements.nshf = 'nshf';
        replacements.rh = 'rhum';
        replacements.lh = 'lhtfl';
        replacements.sh = 'shtfl';
        replacements.xwind = 'u10';
        replacements.ywind = 'v10';
        replacements.precip = 'P_E';
        replacements.Et = 'evap';
    
        fields = fieldnames(MetUM);
        for ff = 1:length(fields)
            if isfield(replacements, fields{ff}) && ~isfield(MetUM, replacements.(fields{ff}))
                MetUM.(replacements.(fields{ff})) = MetUM.(fields{ff});
                MetUM = rmfield(MetUM, fields{ff});
            end
        end
    end
    
    conf.outbase = fullfile(conf.base,'input/configs/',conf.casename);
    
    % Make the output directory if it doesn't exist.
    if exist(conf.outbase, 'dir') ~= 7
        mkdir(conf.outbase)
    end
    
    % Output to a WRF format forcing file.
    if strcmpi(conf.doForcing, 'MetUM')
        write_WRF_forcing(MetUM, ...
            fullfile(conf.outbase, [conf.casename, '_wnd.nc']));
    end