Description of get_NCEP

0001 function data = get_NCEP_forcing(Mobj, modelTime)
0002 % Get the required parameters from NCEP OPeNDAP data to force FVCOM
0003 % (through any of Casename_wnd.nc, Casename_sst.nc, Casename_hfx.nc
0004 % or Casename_pre_evap.nc).
0005 %
0006 % data = get_NCEP_forcing(Mobj, modelTime)
0007 %
0008 % DESCRIPTION:
0009 %   Using OPeNDAP, extract the necessary parameters to create an FVCOM
0010 %   forcing file. Requires the air_sea toolbox and the OPeNDAP toolbox (see
0011 %   below for where to get them).
0012 %
0013 % INPUT:
0014 %   Mobj - MATLAB mesh object
0015 %   modelTime - Modified Julian Date start and end times
0016 %
0017 % OUTPUT:
0018 %   data - struct of the data necessary to force FVCOM. These can be
0019 %   interpolated onto an unstructured grid in Mobj using
0020 %   grid2fvcom_U10V10.m.
0021 %
0022 % The parameters which can be obtained from the NCEP data are:
0023 %     - u wind component (uwnd)
0024 %     - v wind component (vwnd)
0025 %     - Net longwave radiation surface (nlwrs)
0026 %     - Net shortwave radiation surface (nswrs)
0027 %     - Air temperature (air)
0028 %     - Relative humidity (rhum)
0029 %     - Precipitation rate (prate)
0030 %     - Sea level pressure (slp)
0031 %     - Latent heat flux (lhtfl)
0032 %     - Surface heat flux (shtfl)
0033 %     - Potential evaporation rate (pevpr)
0034 %
0035 % In addition to these, the momentum flux is calculated from wind data.
0036 % Precipitation is converted from kg/m^2/s to m/s. Evaporation is
0037 % calculated from the mean daily latent heat net flux (lhtfl) at the
0038 % surface.
0039 %
0040 % REQUIRES:
0041 %   The air_sea toolbox:
0042 %       http://woodshole.er.usgs.gov/operations/sea-mat/air_sea-html/index.html
0043 %   The OPeNDAP toolbox:
0044 %       http://www.opendap.org/pub/contributed/source/ml-toolbox/
0045 %
0046 %
0047 % Author(s)
0048 %   Pierre Cazenave (Plymouth Marine Laboratory)
0049 %   Ricardo Torres (Plymouth Marine Laboratory)
0050 %
0051 % Revision history:
0052 %   2012-10-31 First version based on get_NCEP_L4.m.
0053 %
0054 %==========================================================================
0055 
0056 subname = 'get_NCEP_forcing';
0057 
0058 global ftbverbose;
0059 if(ftbverbose);
0060   fprintf('\n')
0061   fprintf(['begin : ' subname '\n'])
0062 end
0063 
0064 % Get the extent of the model domain (in spherical)
0065 if ~Mobj.have_lonlat
0066     error('Need spherical coordinates to extract the forcing data')
0067 else
0068     % Add a buffer of one grid cell in latitude and two in longitude to
0069     % make sure the model domain is fully covered by the extracted data.
0070     [dx, dy] = deal(2.5, 2.5); % NCEP resolution in degrees
0071     extents = [min(Mobj.lon(:))-(2*dx), max(Mobj.lon(:))+(2*dx), min(Mobj.lat(:))-dy, max(Mobj.lat(:))+dy];
0072 end
0073 
0074 if modelTime(end) - modelTime(1) > 365
0075     error('Can''t (yet) process more than a year at a time.')
0076 end
0077 
0078 yearStart = mjulian2greg(modelTime(1));
0079 yearEnd = mjulian2greg(modelTime(end));
0080 if yearEnd ~= yearStart
0081     error('Can''t (yet) process across a year boundary.')
0082 else
0083     year = yearEnd;
0084 end
0085 
0086 % Set up a struct of the NCEP remote locations in which we're interested.
0087 ncep.uwnd   = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/uwnd.10m.gauss.',num2str(year),'.nc'];
0088 ncep.vwnd   = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/vwnd.10m.gauss.',num2str(year),'.nc'];
0089 ncep.nlwrs  = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/nlwrs.sfc.gauss.',num2str(year),'.nc'];
0090 ncep.nswrs  = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/nswrs.sfc.gauss.',num2str(year),'.nc'];
0091 ncep.air    = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/air.2m.gauss.',num2str(year),'.nc'];
0092 ncep.rhum   = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface/rhum.sig995.',num2str(year),'.nc'];
0093 ncep.prate  = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/prate.sfc.gauss.',num2str(year),'.nc'];
0094 ncep.slp    = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface/slp.',num2str(year),'.nc'];
0095 ncep.lhtfl  = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/lhtfl.sfc.gauss.',num2str(year),'.nc'];
0096 ncep.shtfl  = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/shtfl.sfc.gauss.',num2str(year),'.nc'];
0097 ncep.pevpr  = ['http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/pevpr.sfc.gauss.',num2str(year),'.nc'];
0098 
0099 fields = fieldnames(ncep);
0100 
0101 for aa=1:length(fields)
0102     data.(fields{aa}).data = [];
0103     data.(fields{aa}).time = [];
0104     data.(fields{aa}).lat = [];
0105     data.(fields{aa}).lon = [];
0106     data_attributes.(fields{aa}) = [];
0107 
0108     % Get the data time and convert to Modified Julian Day.
0109     data_time = loaddap([ncep.(fields{aa}),'?time']);
0110     data_attributes.(fields{aa}) = loaddap('-A',[ncep.(fields{aa})]);
0111     timevec = datevec((data_time.time)/24+365);
0112     data.time = greg2mjulian(timevec(:,1), timevec(:,2), timevec(:,3), ...
0113         timevec(:,4), timevec(:,5), timevec(:,6));
0114     % Clip the time to the given range
0115     data_time_mask = data.time >= modelTime(1) & data.time <= modelTime(end);
0116     data_time_idx = 1:size(data.time,1);
0117     data_time_idx = data_time_idx(data_time_mask);
0118     data.time = data.time(data_time_mask);
0119 
0120     % Check the times
0121     %[yyyy,mm,dd,hh,MM,ss] = mjulian2greg(data.time(1))
0122     %[yyyy,mm,dd,hh,MM,ss] = mjulian2greg(data.time(end))
0123 
0124     % Clip the data to the model domain
0125     data_lon = loaddap([ncep.(fields{aa}),'?lon']);
0126     % If the extents are negative in longitude, we need to extract the NCEP
0127     % data in two goes, once for the end of the grid (west of Greenwich),
0128     % once for the beginning (east of Greenwich), and then stick the two
0129     % bits together.
0130     clear index_lon index_lat
0131     if extents(1) < 0 && extents(2) < 0
0132         % This is OK, we can just shunt the values by 360.
0133         extents(1) = extents(1) + 360;
0134         extents(2) = extents(2) + 360;
0135         index_lon = find(data_lon.lon > extents(1) & data_lon.lon < extents(2));
0136     elseif extents(1) < 0 && extents(2) > 0
0137         % This is the tricky one. We'll do two passes to extract the
0138         % western chunk first (extents(1)+360 to 360), then the eastern
0139         % chunk (0-extents(2))
0140         index_lon{1} = find(data_lon.lon >= extents(1) + 360);
0141         index_lon{2} = find(data_lon.lon <= extents(2));
0142     else
0143         % Dead easy, we're in the eastern hemisphere, so nothing too
0144         % strenuous here
0145         index_lon = find(data_lon.lon > extents(1) & data_lon.lon < extents(2));
0146     end
0147 
0148     % Latitude is much more straightforward
0149     data_lat = loaddap([ncep.(fields{aa}),'?lat']);
0150     index_lat = find(data_lat.lat > extents(3) & data_lat.lat < extents(4));
0151 
0152     % Get the data
0153     if iscell(index_lon)
0154         % We need to do each half and merge them
0155         eval(['data1_west.(fields{aa}) = loaddap(''', ncep.(fields{aa}),'?',...
0156             fields{aa},'[', num2str(min(data_time_idx)-1),':',...
0157             num2str(max(data_time_idx)-1), '][',...
0158             num2str(min(index_lat)-1), ':', num2str(max(index_lat)-1),...
0159             '][', num2str(min(index_lon{1})-1), ':',...
0160             num2str(length(data_lon.lon)-1), ']'');']);
0161         eval(['data1_east.(fields{aa}) = loaddap(''', ncep.(fields{aa}),'?',...
0162             fields{aa}, '[', num2str(min(data_time_idx)-1),':',...
0163             num2str(max(data_time_idx)-1), '][',...
0164             num2str(min(index_lat)-1), ':', num2str(max(index_lat)-1),...
0165             '][', '0', ':', num2str(max(index_lon{2})-1), ']'');']);
0166         % Merge the two sets of data together
0167         structfields = fieldnames(data1_west.(fields{aa}).(fields{aa}));
0168         for ii=1:length(structfields)
0169             switch structfields{ii}
0170                 case 'lon'
0171                     % Only the longitude and the actual data need sticking
0172                     % together, but each must be done along a different
0173                     % axis (lon is a vector, the data is an array).
0174                     data1.(fields{aa}).(fields{aa}).(structfields{ii}) = ...
0175                         [data1_west.(fields{aa}).(fields{aa}).(structfields{ii});data1_east.(fields{aa}).(fields{aa}).(structfields{ii})];
0176                 case fields{aa}
0177                     % This is the actual data
0178                     data1.(fields{aa}).(fields{aa}).(structfields{ii}) = ...
0179                         [data1_west.(fields{aa}).(fields{aa}).(structfields{ii}),data1_east.(fields{aa}).(fields{aa}).(structfields{ii})];
0180                 otherwise
0181                     % Assume the data are the same in both arrays. A simple
0182                     % check of the range of values in the difference
0183                     % between the two arrays should show whether they're
0184                     % the same or not. If they are, use the western values,
0185                     % otherwise, warn about the differences. It might be
0186                     % the data are relatively unimportant anyway (i.e. not
0187                     % used later on).
0188                     try
0189                         tdata = data1_west.(fields{aa}).(fields{aa}).(structfields{ii}) - data1_east.(fields{aa}).(fields{aa}).(structfields{ii});
0190                         if range(tdata(:)) == 0
0191                             % They're the same data
0192                             data1.(fields{aa}).(fields{aa}).(structfields{ii}) = ...
0193                                 data1_west.(fields{aa}).(fields{aa}).(structfields{ii});
0194                         else
0195                             warning('Unexpected data field and the west and east halves don''t match. Skipping.')
0196                         end
0197                     catch
0198                         warning('Unexpected data field and the west and east halves don''t match. Skipping.')
0199                     end
0200                     clear tdata
0201             end
0202         end
0203     else
0204         % We have a straightforward data extraction
0205         eval(['data1.(fields{aa}) = loaddap(''', ncep.(fields{aa}),'?',...
0206             fields{aa}, '[', num2str(min(data_time_idx)-1),':',...
0207             num2str(max(data_time_idx)-1), '][',...
0208             num2str(min(index_lat)-1), ':', num2str(max(index_lat)-1),...
0209             '][', num2str(min(index_lon)-1), ':',...
0210             num2str(max(index_lon)-1), ']'');']);
0211     end
0212 
0213     datatmp = squeeze(data1.(fields{aa}).(fields{aa}).(fields{aa}));
0214     datatmp = (datatmp * data_attributes.(fields{aa}).(fields{aa}).scale_factor) + data_attributes.(fields{aa}).(fields{aa}).add_offset;
0215 
0216     data.(fields{aa}).data = cat(1, data.(fields{aa}).data, datatmp);
0217     data.(fields{aa}).time = cat(1, data.(fields{aa}).time, squeeze(data1.(fields{aa}).(fields{aa}).time));
0218     data.(fields{aa}).lat = squeeze(data1.(fields{aa}).(fields{aa}).lat);
0219     data.(fields{aa}).lon = squeeze(data1.(fields{aa}).(fields{aa}).lon);
0220 end
0221 
0222 % Now we have some data, we need to create some additional parameters
0223 % required by FVCOM.
0224 
0225 % Convert precipitation from kg/m^2/s to m/s (required by FVCOM) by
0226 % dividing by freshwater density (kg/m^3).
0227 data.prate.data = data.prate.data/1000;
0228 
0229 % Evaporation can be approximated by:
0230 %
0231 %   E(m/s) = lhtfl/Llv/rho
0232 %
0233 % where:
0234 %
0235 %   lhtfl   = "Mean daily latent heat net flux at the surface"
0236 %   Llv     = Latent heat of vaporization (approx to 2.5*10^6 J kg^-1)
0237 %   rho     = 1025 kg/m^3
0238 %
0239 Llv = 2.5*10^6;
0240 rho = 1025; % using a typical value for seawater.
0241 Et = data.lhtfl.data/Llv/rho;
0242 data.P_E.data = data.prate.data-Et;
0243 
0244 % Calculate the momentum flux
0245 WW = data.uwnd.data + data.vwnd.data * 1i;
0246 data.tau.data = stresslp(abs(WW),10);
0247 [data.tx.data,data.ty.data] = wstress(data.uwnd.data,data.vwnd.data,10);
0248 data.tx.data=reshape(data.tx.data*0.1, size(data.uwnd.data)); % dyn/cm^2 to N/m^2
0249 data.ty.data=reshape(data.ty.data*0.1, size(data.uwnd.data)); % dyn/cm^2 to N/m^2
0250 
0251 % Get the fields we need for the subsequent interpolation
0252 data.lon = data.uwnd.lon;
0253 data.lon(data.lon > 180) = data.lon(data.lon > 180) - 360;
0254 data.lat = data.uwnd.lat;
0255 
0256 % Have a look at some data.
0257 % [X, Y] = meshgrid(data.lon, data.lat);
0258 % for i=1:size(data.uwnd.data, 3)
0259 %     figure(1)
0260 %     clf
0261 %     uv = sqrt(data.uwnd.data(:, :, i).^2 + data.vwnd.data(:, :, i).^2);
0262 %     pcolor(X, Y, uv)
0263 %     shading flat
0264 %     axis('equal','tight')
0265 %     pause(0.1)
0266 % end
get_NCEP_forcing

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
