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