data has evolved. The script
was provided as pseudo-code for the user to modify and adapt. I don't have the time to look at it now. It depends on whether you are extracting instantaneous or accumulated fields, where the time interval is essential to determining the scales. I suggest you read the
documentation and fix the scales according to your needs. I believe the
changed the shortname of some of the data in the output NetCDF files. Users should invest some time ensuring that the data used to force ROMS is accurate, and always verify processing scripts provided by others, including whether I wrote such a script. I am constantly verifying all the available scripts in my applications. It is paramount to do so to analyze ROMS solutions.
Code: Select all
% Set file names. Input ERA5-Interim input data file names. The data
% was extracted for the LAKE ERIE (85W, 40N) to (75W, 45N) hourly. The
% dataset was extracted in several monthly files:
%
% era5_201001.nc
% era5_201002.nc
% ...
% era5_2012011.nc
% era5_2012012.nc
%
% ERA-5 Interim variable names: @ denotes accumulated quantity that must
% be divided by the time interval into the
% cycle 1, 3, 6, 9 or 12 hours
%
% alnip 0.0 - 1.0 neat IR albedo for direct radiation
% asn 0.0 - 1.0 snow albedo
% d2m K 2m dewpoint temperature
% @ e m evaporation
% @ ewss N m-2 s east-west surface stress
% lblt K lake bottom temperature
% lict K lake ice surface temperature
% mer kg m-2 s-1 mean evaporation rate
% metss N m-2 mean eastward turbulent surface stress
% mlspr kg m-2 s-1 mean large-scale precipitation rate
% mntss N m-2 mean northward turbulent surface stress
% mror kg m-2 s-1 mean runoff rate
% msr kg m-2 s-1 mean snowfall rate
% msdwlwrf W m-2 mean surface downward longwave radiation flux
% msdwswrf W m-2 mean surface downward shortwave radiation flux
% msl Pa mean sea level pressure
% mslhf W m-2 mean surface latent heat flux
% msnlwrf W m-2 mean surface net longwave radiation flux
% msnswrf W m-2 mean surface net shortwave radiation flux
% msshf W m-2 mean surface sensible heat flux
% mtnlwrf W m-2 mean top net longwave radiaton flux
% mtnswrf W m-2 mean top net shortwave radiation flux
% mtpr kg m-2 s-1 mean total precipitation rate
% @ nsss N m-2 s north-south surface stress
% @ par W m-2 s photosynthetically active radiation at surface
% @ ro m runoff
% sf m water equiv snowfall
% @ slhf W m-2 s surface latent heat flux
% @ sshf W m-2 s surface sensible heat flux
% @ ssr W m-2 s surface net solar radiation (shortwave)
% @ str W m-2 s surface net thermal radiation (longwave)
% @ strd W m-2 s surface thermal radiation downwards
% t2m K 2 metre temperature
% @ tcc nondimensional total cloud cover [0:1]
% @ tp m total precipitation
% v10 m s-1 10 metre U wind component
% u10n m s-1 neutral wind at 10m u-component
% vl0 m s-1 10 metre V wind component
% v10n m s-1 neutral wind at 10m v-component
%
% This dataset is written in compact way (short numbers). We need to
% convert to floating-point data and scale to ROMS units:
%
% deltat 1*3600 1-hour step
%
% Uwind (m s-1) v10
% Vwind (m s-1) v10
% sustr (N m-2) ewss / deltat
% svstr (N m-2) nsss / deltat
% shflux (W m-2) (ssr+str+sshf+slhf) / deltat
% swrad (W m-2) ssr / deltat
% lwrad_down (W m-2) strd / deltat
% latent (W m-2) slhf / deltat
% sensible (W m-2) sshf / deltat
% rain (kg m-2 s-1) tp * Rho_w / deltat
% evaporation (kg m-2 s-1) e * Rho_w / deltat
% swflux (cm day-1) (e - tp) * 100 / (3/24); 0.125 day step
% cloud (nondimesional) tcc
% Pair (mb) msl / 100; (1 mb = 100 Pa)
% Tair (Celsius) t2m - 273.15; (1 C = 273.15 K)
% Qair (percentage) 100 * (E/Es)
%
% where
%
% Rho_w = 1000 kg m-3 (water density)
%
% E = 6.11 * 10.0 ^ (7.5 * d2m / (237.7 + d2m)) vapor pressure (mb)
% d2m in Celsius
%
% Es = 6.11 * 10.0 ^ (7.5 * t2m / (237.7 + t2m)) saturation vapor
% pressure (mb)
% t2m in Celsius
%
%-------------------------------------------------------------------------
% Constants.
Rho_w = 1000; % kg m-3
FluxScale = 1.0;
StressScale = 1.0;
% Path to downloaded ERA data files.
Dir = fullfile(getenv('HOME'), ...
'ROMS/Projects/LAKE_ERIE/Data/FRC/ERA5');
% Base date for ROMS forcing files: "hours since 2000-01-01 00:00:00".
mybasedate = datenum(2000,1,1,0,0,0);
% Set ROMS and ECMWF-ERA field NetCDF variable names for each forcing
% variable, input file name prefix, output file name. Notice that output
% forcing variables follows ROMS metadata design.
%
% If the ROMS option BULK_FLUXES is NOT activatived, we just need to
% process elements 1:5 in the structure vector (F). Otherwise, we need
% process almost all the fields below.
clear S F
F( 1).Vname = {'sustr', 'metss'};
F( 1).Tname = {'sms_time', 'time'};
F( 1).input = 'era5_';
F( 1).output = 'LakeErie_sustr_era5.nc';
F( 1).scale = StressScale;
F( 2).Vname = {'svstr', 'mntss'};
F( 2).Tname = {'sms_time', 'time'};
F( 2).input = 'era5_';
F( 2).output = 'LakeErie_svstr_era5.nc';
F( 2).scale = StressScale;
F( 3).Vname = {'shflux', 'null'};
F( 3).Tname = {'shf_time', 'valid_time'};
F( 3).input = 'era5_';
F( 3).output = 'LakeErie_shflux_era5.nc';
F( 3).scale = FluxScale;
F( 4).Vname = {'swflux', 'null'};
F( 4).Tname = {'swf_time', 'time'};
F( 4).input = 'era5_';
F( 4).output = 'LakeErie_swflux_era5.nc';
F( 4).scale = Rho_w;
F( 5).Vname = {'swrad', 'msnswrf'};
F( 5).Tname = {'srf_time', 'valid_time'};
F( 5).input = 'era5_201001_201212.nc';
F( 5).output = 'LakeErie_swrad_era5.nc';
F( 5).scale = FluxScale;
F( 6).Vname = {'Uwind', 'u10'};
F( 6).Tname = {'wind_time', 'time'};
F( 6).input = 'era5_';
F( 6).output = 'LakeErie_uwind_era5.nc';
F( 6).scale = 1.0;
F( 7).Vname = {'Vwind', 'v10'};
F( 7).Tname = {'wind_time', 'time'};
F( 7).input = 'era5_';
F( 7).output = 'LakeErie_vwind_era5.nc';
F( 7).scale = 1.0;
F( 8).Vname = {'lwrad', 'msnlwrf'};
F( 8).Tname = {'lrf_time', 'valid_time'};
F( 8).input = 'era5_201001_201212.nc';
F( 8).output = 'LakeErie_lwrad_era5.nc';
F( 8).scale = FluxScale;
F( 9).Vname = {'lwrad_down', 'msdwlwrf'};
F( 9).Tname = {'lrf_time', 'time'};
F( 9).input = 'era5_';
F( 9).output = 'LakeErie_lwrad_era5.nc';
F( 9).scale = FluxScale;
F(10).Vname = {'latent', 'mslhf'};
F(10).Tname = {'lhf_time', 'time'};
F(10).input = 'era5_';
F(10).output = 'LakeErie_latent_era5.nc';
F(10).scale = FluxScale;
F(11).Vname = {'sensible', 'msshf'};
F(11).Tname = {'sen_time', 'time'};
F(11).input = 'era5_';
F(11).output = 'LakeErie_sensible_era5.nc';
F(11).scale = FluxScale;
F(12).Vname = {'cloud', 'tcc'};
F(12).Tname = {'cloud_time', 'time'};
F(12).input = 'era5_';
F(12).output = 'LakeErie_cloud_era5.nc';
F(12).scale = 1.0;
F(13).Vname = {'rain', 'mtpr'};
F(13).Tname = {'rain_time', 'time'};
F(13).input = 'era5_';
F(13).output = 'LakeErie_rain_era5.nc';
F(13).scale = Rho_w;
F(14).Vname = {'Pair', 'msl'};
F(14).Tname = {'pair_time', 'time'};
F(14).input = 'era5_';
F(14).output = 'LakeErie_Pair_era5.nc';
F(14).scale = 0.01;
F(15).Vname = {'Tair', 't2m'};
F(15).Tname = {'tair_time', 'time'};
F(15).input = 'era5_';
F(15).output = 'LakeErie_Tair_era5.nc';
F(15).scale = 1.0;
F(16).Vname = {'Qair', 't2m'}; % Use temperature (t2m) and
F(16).Tname = {'qair_time', 'time'}; % dewpoint temperature (d2m)
F(16).input = 'era5_'; % to compute relative humidity
F(16).output = 'LakeErie_Qair_era5.nc';
F(16).scale = 1.0;
F(17).Vname = {'PAR', 'par'};
F(17).Tname = {'srf_time', 'time'};
F(17).input = 'era5_';
F(17).output = 'LakeErie_PAR_era5.nc';
F(17).scale = 1.0;
% Set field element in structure F to process.
%doFields = [5 8];
doFields = [3];
% Various parameters.
spherical = true; % Spherical switch
LonMin = -85.0; % Lake Erie left-bottom longitude
LonMax = -75.0; % Lake Erie right-top longitude
LatMin = 40.0; % Lake Erie left-bottom latitude
LatMax = 45.0; % Lake Erie right-top latitude
FlipLat = true; % Flip data in ERA NetCDF files
% (see below information)
StrDay = datenum('01-Jan-2010'); % starting day to process
EndDay = datenum('31-Dec-2012'); % ending day to process
nctype = 'nc_float'; % Input data is in single precision
Unlimited = true; % time dimension is umlimited in
% output files
%--------------------------------------------------------------------------
% Create surface forcing NetCDF files: build creation parameters
% structure, S. We want to create a single file for each forcing
% field. However, the wind components "Uwind" and "Vwind" are
% saved in the same NetCDF file.
%--------------------------------------------------------------------------
Tindex = [];
ReplaceValue = NaN;
PreserveType = true;
mode = netcdf.getConstant('CLOBBER'); % overwrite!!!
mode = bitor(mode,netcdf.getConstant('64BIT_OFFSET'));
% The strategy here is to build manually the NetCDF metadata structure to
% facilitate creating several NetCDF file in a generic and compact way.
% This structure is similar to that returned by "nc_inq" or native Matlab
% function "ncinfo".
%
% Notice that we call the "roms_metadata" function to create the fields
% in the structure. Then, we call "check_metadata" for fill unassigned
% values and to check for consistency.
OutFiles = {};
for n = doFields,
ncname = char(F(n).output);
Vname = char(F(n).Vname{1});
Tname = char(F(n).Tname{1});
create_file = true;
% if (n == 1)
% create_file = true;
% else
% create_file = ~any(strcmp(OutFiles, ncname));
% end
OutFiles = [OutFiles, ncname];
% InpFile = fullfile(Dir, strcat(char(F(n).input),'201001.nc'));
InpFile = fullfile(Dir, char(F(5).input));
if (create_file),
disp(blanks(1));
disp(['** Creating ROMS NetCDF forcing file: ', ncname,' **']);
disp(blanks(1))
end
S.Filename = ncname;
lon = nc_read(InpFile, 'longitude', ...
Tindex, ReplaceValue, PreserveType);
lat = nc_read(InpFile, 'latitude', ...
Tindex, ReplaceValue, PreserveType);
% We want the longitude in the range from -180 to 180 instead of 0 to 360.
% Notice that the latitude is flip so the origin is at southwest corner
% of the extracted grid (LonMin,LatMin)
ROMS_lon = lon;
ind = find(ROMS_lon > 180);
if (~isempty(ind))
ROMS_lon(ind) = ROMS_lon(ind) - 360;
end
if (FlipLat),
ROMS_lat = flipud(lat);
else
ROMS_lat = lat;
end
Im = length(ROMS_lon);
Jm = length(ROMS_lat);
S.Attributes(1).Name = 'type';
S.Attributes(1).Value = 'FORCING file';
S.Attributes(2).Name = 'title';
S.Attributes(2).Value = ['ECMWF ERA5-Interim Dataset, ', ...
'Lake Erie'];
S.Attributes(3).Name = 'ERA5_variable';
S.Attributes(3).Value = char(F(n).Vname{2});
S.Attributes(4).Name = 'history';
S.Attributes(4).Value = ['Forcing file created with ', ...
which(mfilename) ' on ' date_stamp];
S.Dimensions(1).Name = 'lon';
S.Dimensions(1).Length = Im;
S.Dimensions(1).Unlimited = false;
S.Dimensions(2).Name = 'lat';
S.Dimensions(2).Length = Jm;
S.Dimensions(2).Unlimited = false;
S.Dimensions(3).Name = Tname;
S.Dimensions(3).Length = nc_constant('nc_unlimited');
S.Dimensions(3).Unlimited = true;
S.Variables(1) = roms_metadata('spherical');
S.Variables(2) = roms_metadata('lon');
S.Variables(3) = roms_metadata('lat');
S.Variables(4) = roms_metadata(Tname, [], [], Unlimited);
S.Variables(5) = roms_metadata(Vname, spherical, nctype, Unlimited);
% Edit the time variable 'units" attribute for the correct reference
% time and add calendar attribute.
natts = length(S.Variables(4).Attributes);
iatt = strcmp({S.Variables(4).Attributes.Name}, 'units');
S.Variables(4).Attributes(iatt).Value = ['hours since ' datestr(mybasedate,31)];
S.Variables(4).Attributes(natts+1).Name = 'calendar';
S.Variables(4).Attributes(natts+1).Value = 'gregorian';
% Change the dimensions of "shflux", "swflux", "sustr" or "svstr" to
% (lon,lat) dimensions to allow interpolation from a coarser grid in
% ROMS. As default, 'roms_metadata' uses ROMS native coordinates.
%
% This illustrates also how to manipulate the default metadata values set
% in 'roms_metadata'. It has to be done before calling 'check_metadata'.
if (strcmp(Vname, 'shflux') || ...
strcmp(Vname, 'swflux') || ...
strcmp(Vname, 'sustr') || ...
strcmp(Vname, 'svstr')),
S.Variables(5).Dimensions(1).Name = 'lon';
S.Variables(5).Dimensions(2).Name = 'lat';
end
% Check ROMS metadata structure. Fill unassigned fields.
S = check_metadata(S);
% Create forcing NetCDF files. Write our coordinates. Notice that
% "nc_append" is used here since we wand both wind components "Uwind"
% and "Vwind" to be in the same output NetCDF file for CF conventions.
if (create_file),
ncid = nc_create(S.Filename, mode, S); % create a new NetCDF file
lon = repmat(ROMS_lon, [1 Jm]);
lat = repmat(ROMS_lat', [Im 1]);
status = nc_write(S.Filename, 'spherical', int32(spherical));
status = nc_write(S.Filename, 'lon', lon);
status = nc_write(S.Filename, 'lat', lat);
else
nc_append(S.Filename, S) % append to existing NetCDF file
end
end
%--------------------------------------------------------------------------
% Convert forcing data to floating-point and scale to ROMS units.
%--------------------------------------------------------------------------
% This particular data set has a time coordinate in hours starting
% on 1-Jan-1900 00:00:00.
%epoch = datenum('01-Jan-1900');
epoch = datenum('01-Jan-1970');
% Set reference time for this application using the specified base date.
ref_time = (mybasedate - epoch)*24; % hours since 01-Jan-2000
StrYear = str2num(datestr(StrDay, 'yyyy'));
EndYear = str2num(datestr(EndDay, 'yyyy'));
MyRec = zeros([1 length(F)]);
disp(blanks(1));
for n = doFields
year = StrYear;
% while (year <= EndYear)
% for month = 1:12
% suffix = strcat(num2str(year), num2str(month,'%2.2i'), '.nc');
% InpFile = fullfile(Dir, strcat(char(F(n).input),suffix));
InpFile = fullfile(Dir, char(F(5).input));
OutFile = char(F(n).output);
Troms = char(F(n).Tname{1});
Tecmwf = char(F(n).Tname{2});
Vroms = char(F(n).Vname{1});
Vecmwf = char(F(n).Vname{2});
scale = abs(F(n).scale);
% Determine time records to process.
time = nc_read(InpFile, Tecmwf); % hours
StrRec = 1;
EndRec = length(time);
% Read and write forcing fields. The forcing field are scaled to ROMS
% units.
for Rec=StrRec:EndRec
MyRec(n) = MyRec(n) + 1;
mydate = datestr(epoch+time(Rec)/86400);
disp(['** Processing: ', Vroms, ' for ',mydate,' **']);
% frc_time = time(Rec) - ref_time; % hours since 01-Jan-2000
frc_time = time(Rec)/3600 - ref_time; % hours since 01-Jan-2000
switch Vroms
case 'Tair'
field = nc_read(InpFile, Vecmwf, Rec);
field = field - 273.15; % Kelvin to Celsius
case 'Qair'
tsur = nc_read(InpFile, 't2m', Rec); % 2m temperature
tdew = nc_read(InpFile, 'd2m', Rec); % 2m dewpoint
tsur = tsur - 273.15;
tdew = tdew - 273.15;
% E = 6.11 .* 10.0 .^ (7.5 .* tdew ./ (237.7 + tdew));
% Es = 6.11 .* 10.0 .^ (7.5 .* tsur ./ (237.7 + tsur));
E = 6.112 * exp( 17.62*tdew ./ (243.12+tdew));
Es = 6.112 * exp( 17.62*tsur ./ (243.12+tsur));
field = 100.0 .* (E ./ Es);
case 'swflux'
evap = nc_read(InpFile, 'mer' , Rec); % evaporation
prec = nc_read(InpFile, 'mtpr', Rec); % precipitation
evap = -evap; % evap < 0 where latent < 0
field = (evap - prec) .* scale;
case 'shflux'
fname = 'LakeErie_sensible_era5.nc';
sensbl = nc_read(fname, 'sensible' , Rec); % sensible
fname = 'LakeErie_latent_era5.nc';
latent = nc_read(fname, 'latent', Rec); % latent
fname = 'LakeErie_lwrad_era5.nc';
nlwrad = nc_read(fname, 'lwrad', Rec); % net longwave
fname = 'LakeErie_swrad_era5.nc';
nswrad = nc_read(fname, 'swrad', Rec); % net shortwave
field = (sensbl+latent+nlwrad+nswrad) .* scale;
otherwise
ReplaceValue = 0.0;
field = nc_read(InpFile, Vecmwf, Rec, ReplaceValue);
ind = find(field < 1.0E-13);
if (~isempty(ind))
field(ind)=0.0;
end
field = field.*scale;
end
fieldfinal = field;
% If the scale F(n).scale is set to negative, the input ECMWF data is a
% cummulative integral in forecast cycle from hour zero.
% For steps at 1, 6, 9 and 12 hours we must separate last 3 hours of
% integration from previous accumulation.
% At 3 hour step don't change anything
if (F(n).scale < 0)
step = rem(frc_time,0.5)*24;
if step == 3
fieldfinal = field;
else
fieldfinal = field - field_previous; % At other steps subtract
end % the previous accumulation
frc_time = frc_time - 1.5/24; % Center forcing time on the
% accumulation interval
field_previous = field; % Save this accumulation
% to on the next step
end
% If appropriate, flip the data so the origin (1,1) corresponds to
% the lower-left corner (LonMin,LatMin) of the extracted region.
% ERA data is written into NetCDF files with the origin at (1,end).
% That is, the origin is the left-upper corner of the extracted grid
% (LonMin,LatMax). Users need to check and plot the extracted data
% from the ECMWF server.
if (FlipLat)
fieldfinal = fliplr(fieldfinal);
end
% Write out data.
status = nc_write(OutFile, Troms, frc_time, MyRec(n));
status = nc_write(OutFile, Vroms, fieldfinal, MyRec(n));
end
% end % input monthly data loop
% Advance year counter. Recall all input files are split in annual
% files.
year = year + 1;
% end
end
