Get the required parameters from NAE2 data to force FVCOM (through
Casename_wnd.nc).
Mobj = get_NAE2_forcing(Mobj, inputConf)
DESCRIPTION:
Extract meteorological forcing data from Met Office NAE2 data to create
an FVCOM forcing file.
INPUT:
Mobj - MATLAB mesh object
OUTPUT:
Mobj - MATLAB mesh object containing meteorological forcing data for
FVCOM, interpolated onto an unstructured grid.
The parameters which can be obtained from the AMM/S12 model output are:
- u wind component (uwnd)
- v wind component (vwnd)
- Sea level pressure (slp)
In addition to these, the momentum flux is calculated from wind data.
KJT note: took this out from Pierre's version. Implement this?
Algorithm is based on Phil Hall's 'produce_netcdf_input_data.py' script,
adapted to handle NAE2 data by Karen Thurston
('produce_netcdf_input_dataR.py'. This Python script is, in turn, based
on Jenny Brown's 'MET_INT.f' Fortran script.
Author(s)
Karen Thurston (National Oceanography Centre Liverpool)
Revision history:
2012-12-05 First version
==========================================================================0001 function Mobj = get_NAE2_forcing(Mobj, inputConf) 0002 % Get the required parameters from NAE2 data to force FVCOM (through 0003 % Casename_wnd.nc). 0004 % 0005 % Mobj = get_NAE2_forcing(Mobj, inputConf) 0006 % 0007 % DESCRIPTION: 0008 % Extract meteorological forcing data from Met Office NAE2 data to create 0009 % an FVCOM forcing file. 0010 % 0011 % INPUT: 0012 % Mobj - MATLAB mesh object 0013 % 0014 % OUTPUT: 0015 % Mobj - MATLAB mesh object containing meteorological forcing data for 0016 % FVCOM, interpolated onto an unstructured grid. 0017 % 0018 % The parameters which can be obtained from the AMM/S12 model output are: 0019 % - u wind component (uwnd) 0020 % - v wind component (vwnd) 0021 % - Sea level pressure (slp) 0022 % 0023 % In addition to these, the momentum flux is calculated from wind data. 0024 % KJT note: took this out from Pierre's version. Implement this? 0025 % 0026 % Algorithm is based on Phil Hall's 'produce_netcdf_input_data.py' script, 0027 % adapted to handle NAE2 data by Karen Thurston 0028 % ('produce_netcdf_input_dataR.py'. This Python script is, in turn, based 0029 % on Jenny Brown's 'MET_INT.f' Fortran script. 0030 % 0031 % Author(s) 0032 % Karen Thurston (National Oceanography Centre Liverpool) 0033 % 0034 % Revision history: 0035 % 2012-12-05 First version 0036 % 0037 %========================================================================== 0038 0039 subname = 'get_NAE2_forcing'; 0040 0041 global ftbverbose; 0042 if(ftbverbose); 0043 fprintf('\n') 0044 fprintf(['begin : ' subname '\n']) 0045 end 0046 0047 %% Where are the files we'll need? 0048 if isunix % Unix? 0049 metfname = ['/bank/jane/met/',datestr(inputConf.startDate,'YYYY'),... 0050 '/',lower(datestr(inputConf.startDate,'mmmYY')),'nae10R.dat']; 0051 comprfname = '/login/jane/NAE2/metintco.cs3x.nae2.compress.2'; 0052 setupfname = '/work/jane/cs3x/prep/setupcs3xSGIL.uda'; 0053 % elevfname = ['/bank/jane/cs3x/sarray.uda.',... 0054 % datestr(inputConf.startDate,'YYYY')]; 0055 elseif ispc % Or Windows? 0056 metfname = ['\\store\bank\jane\met\',datestr(inputConf.startDate,'YYYY'),... 0057 '\',lower(datestr(inputConf.startDate,'mmmYY')),'nae10R.dat']; 0058 comprfname = '\\store\kthurs\from_Jane\metintco.cs3x.nae2.compress.2'; 0059 setupfname = '\\store\work\jane\cs3x\prep\setupcs3xSGIL.uda'; 0060 % elevfname = ['\\store\bank\jane\cs3x\sarray.uda.',... 0061 % datestr(inputConf.startDate,'YYYY')]; 0062 end 0063 0064 % Define the parameters of the cs3x grid 0065 cs3x_nx = 198; % number of grid cells in the x direction 0066 cs3x_ny = 207; % number of grid cells in the y direction 0067 cs3x_lonstart = -19-(5/6); % longitude of the bottom left corner 0068 cs3x_latstart = 40+(1/9); % latitude of the bottom left corner 0069 cs3x_loninc = 1/6; % grid resolution in the x direction (degrees) 0070 cs3x_latinc = 1/9; % grid resolution in the y direction (degrees) 0071 cs3x_lonfin = cs3x_lonstart+(cs3x_loninc*(cs3x_nx-1)); % longitude of the top right corner 0072 cs3x_latfin = cs3x_latstart+(cs3x_latinc*(cs3x_ny-1)); % latitude of the top right corner 0073 cs3x_lon = cs3x_lonstart:cs3x_loninc:cs3x_lonfin; % array of grid lon points 0074 cs3x_lat = cs3x_latstart:cs3x_latinc:cs3x_latfin; % array of grid lat points 0075 % cs3x_area = cs3x_loninc*cs3x_latinc; 0076 0077 % Sanity check. Does our FVCOM grid fit within the cs3x domain? 0078 if min(Mobj.lon) < cs3x_lonstart || max(Mobj.lon) > cs3x_lonfin || ... 0079 min(Mobj.lat) < cs3x_latstart || max(Mobj.lat) > cs3x_latfin 0080 error('Your FVCOM grid is bigger than the available cs3x grid. Choose another met forcing option or crop your FVCOM grid.') 0081 end 0082 0083 % Open the met compression info file 0084 fid = fopen(comprfname); 0085 tline = fgets(fid); 0086 count=0; 0087 while ischar(tline) 0088 count=count+1; 0089 compr_info{count}=tline; 0090 tline = fgets(fid); 0091 end 0092 fclose(fid); 0093 0094 % This section reads info from the met data compression file. The first 0095 % line is a set of constants (NCCP, NRRP, NCCW, NRRW, NMOD/NTOTI). 0096 0097 % Read the first line of the compression info file 0098 temp = sscanf(char(compr_info{1}),'%u'); 0099 NCCP = temp(1); % Number of columns in the pressure data 0100 NRRP = temp(2); % Number of rows in the pressure data 0101 NCCW = temp(3); % Number of columns in the wind data 0102 NRRW = temp(4); % Number of rows in the wind data 0103 NTOTI = temp(5); % Also referred to as 'NMOD' 0104 0105 % Initialise compr_array 0106 compr_array = nan(10,length(compr_info)-1); 0107 0108 for i=2:length(compr_info) 0109 temp = sscanf(char(compr_info{i}),'%8f'); 0110 compr_array(:,i-1) = temp; 0111 end 0112 0113 nelf = max(compr_array(1:4*NTOTI))+1; % No idea what this is for 0114 0115 % Index of specified corner of source gridbox (pressure) 0116 IBLP = compr_array(1:NTOTI)-1; % Bottom left point 0117 IBRP = compr_array(NTOTI+1:2*NTOTI)-1; % Bottom right point 0118 0119 % Index of specified corner of source gridbox (wind) 0120 IBLW = compr_array(2*NTOTI+1:3*NTOTI)-1; % Bottom left point 0121 IBRW = compr_array(3*NTOTI+1:4*NTOTI)-1; % Bottom right point 0122 0123 % Weight applied to value at specified corner of source gridbox (pressure) 0124 WTRP = compr_array(4*NTOTI+1:5*NTOTI); % Top right point 0125 WBRP = compr_array(5*NTOTI+1:6*NTOTI); % Bottom right point 0126 WTLP = compr_array(6*NTOTI+1:7*NTOTI); % Top left point 0127 WBLP = compr_array(7*NTOTI+1:8*NTOTI); % Bottom left point 0128 0129 % Weight applied to value at specified corner of source gridbox (wind) 0130 WTRW = compr_array(8*NTOTI+1:9*NTOTI); % Top right point 0131 WBRW = compr_array(9*NTOTI+1:10*NTOTI); % Bottom right point 0132 WTLW = compr_array(10*NTOTI+1:11*NTOTI); % Top left point 0133 WBLW = compr_array(11*NTOTI+1:12*NTOTI); % Bottom left point 0134 0135 % Coefficients to rotate winds onto equatorial lat/lon grid from ELF model 0136 COEFF1 = compr_array(12*NTOTI+1:13*NTOTI); 0137 COEFF2 = compr_array(13*NTOTI+1:14*NTOTI); 0138 0139 clear compr_info compr_array 0140 0141 % This section reads cs3x model control data from the setup file. 0142 fid=fopen(setupfname,'r','b'); 0143 temp=fread(fid,inf,'*int32','b'); 0144 fclose(fid); 0145 0146 NRR = temp(2); % Number of rows in full rectangle 0147 NCC = temp(3); % Number of columns in full rectangle 0148 ITOT = temp(4); % Total number of points in compact arrays 0149 IINZ = temp(5); % Number of internal Z-points 0150 0151 NTRNS = temp(8:NRR+7); % Transformation matrix for compact addressing 0152 NTRNT = temp(NRR+10:2*NRR+9); % Transformation matrix for compact addressing 0153 NSUM = temp(2*NRR+12:3*NRR+11); % Transformation matrix for compact addressing 0154 0155 clear temp 0156 0157 % The array section numbers (as in MET_INT.F). 0158 hw0 = 1; 0159 hw1 = 326; 0160 hw2 = 100; 0161 hw3 = 521; 0162 hw4 = 600; 0163 hw5 = (hw1-hw0)*(hw3-hw2); 0164 0165 % As far as I can tell, this section generates the cross-indexing between 0166 % the met data and the cs3x grid. 0167 ELFI = zeros(hw5,1); 0168 FIEJ = zeros(hw5,1); 0169 0170 count = 1; 0171 j = 1; 0172 0173 for iy = hw0:(hw1-1) 0174 for ix = hw2:(hw3-1) 0175 i = ix+iy*hw4; 0176 j = j+1; 0177 ELFI(count) = i-1; 0178 FIEJ(count) = j-1; 0179 count = count+1; 0180 end 0181 end 0182 0183 if max(ELFI)+1 > nelf 0184 nelf = max(ELFI)+1; 0185 end 0186 0187 if max(FIEJ)+1 > nelf 0188 nelf = max(FIEJ)+1; 0189 end 0190 0191 % I may need to find a different method for storing the data. I'm not sure 0192 % if I can do a whole month of data on my laptop. 0193 PASTIT = false; 0194 kline = 0; 0195 0196 % Open the met file 0197 fid = fopen(metfname); 0198 0199 % calculate the number of lines per data segment 0200 I = fgets(fid); 0201 m1 = sscanf(I,'%u'); 0202 data_seg = ceil(m1(1)/m1(2)); 0203 0204 frewind(fid); 0205 0206 % Initialise the met data array 0207 met_temp = cell(1); 0208 0209 % Get met data from the met file and write to temporary cell array 0210 while PASTIT == false 0211 kline = kline + 2; 0212 I = fgets(fid); 0213 if I==-1 0214 % break the loop, it's the EOF 0215 break 0216 end 0217 J = fgets(fid); 0218 if J==-1 0219 % break the loop, it's the EOF 0220 break 0221 end 0222 0223 % find the date in the met file 0224 datem = sscanf(J,'%u'); 0225 datem = datenum(datem(4),datem(3),datem(2),datem(1),0,0); 0226 0227 % Compare met file date with inputConf.startDate (is it bigger?) and 0228 % with inputConf.endDate (is it smaller?) If yes and yes, then we want 0229 % this. Write data to a temporary array. 0230 if (datem >= datenum(inputConf.startDate)) && ... 0231 (datem <= datenum(inputConf.endDate)) 0232 % Initialise temp cell array 0233 temp = cell(data_seg+2,1); 0234 % write I and J 0235 temp{1} = sscanf(I,'%u'); 0236 temp{2} = sscanf(J,'%u'); 0237 % Get the actual data 0238 for m = 1:data_seg 0239 I = fgetl(fid); 0240 temp{m+2} = sscanf(I,'%10f'); 0241 end 0242 met_temp{end+1}=temp; 0243 clear temp; 0244 else 0245 PASTIT = true; 0246 end 0247 end 0248 0249 fclose(fid); 0250 0251 met_temp = met_temp(2:end); 0252 0253 % Not sure what these numbers are for. 0254 JBLP = IBLP-NCCP; 0255 JBRP = IBRP-NCCP; 0256 JBLW = IBLW-NCCW; 0257 JBRW = IBRW-NCCW; 0258 0259 elf = zeros(nelf,1); 0260 0261 % Temporary arrays to hold pressure/u-wind/v-wind from extracted met data 0262 P2 = zeros(NTOTI,1); 0263 tmpa = zeros(NTOTI,1); 0264 tmpb = zeros(NTOTI,1); 0265 0266 % i1j1 = zeros(4,Mobj.nElems); 0267 % A1 = zeros(2,Mobj.nElems); 0268 % A2 = zeros(2,Mobj.nElems); 0269 % A3 = zeros(2,Mobj.nElems); 0270 % A4 = zeros(2,Mobj.nElems); 0271 0272 % FirstInt = true; 0273 kline = 1; 0274 0275 % Arrays to hold the pressure/u-wind/v-wind data when interpolated to the 0276 % FVCOM grid 0277 slp = zeros(Mobj.nVerts,size(met_temp,2)/3); 0278 uwnd = zeros(Mobj.nElems,size(met_temp,2)/3); 0279 vwnd = zeros(Mobj.nElems,size(met_temp,2)/3); 0280 time = zeros(size(met_temp,2)/3,1); 0281 0282 % Take our forcing data, translate it to the cs3x grid, then interpolate 0283 % onto the FVCOM grid. 0284 for m=0:3:size(met_temp,2)-1 0285 for ISW = 1:3 0286 kline = kline+2; 0287 ipts = met_temp{m+ISW}{1}(1); 0288 kpts = met_temp{m+ISW}{1}(2); 0289 0290 % Store the date 0291 if ISW == 1 0292 time((m/3)+1) = datenum(met_temp{m+ISW}{2}(4),... 0293 met_temp{m+ISW}{2}(3),met_temp{m+ISW}{2}(2),... 0294 met_temp{m+ISW}{2}(1),0,0); 0295 end 0296 0297 % Preallocate 'field' for speed 0298 field = NaN(8,data_seg); 0299 0300 % Get the met data we're interested in 0301 for n=1:data_seg 0302 field(:,n)=met_temp{m+ISW}{n+2}; 0303 end 0304 field = reshape(field,data_seg*8,1); 0305 0306 % if the last column of data is not full, adjust for that 0307 if size(met_temp{m+ISW}{end},1)<8 0308 too_many = 8-size(met_temp{m+ISW}{end},1); 0309 field = field(1:end-too_many); 0310 end 0311 0312 elf(ELFI)=field(FIEJ); 0313 0314 % Interpolate the data to the sea model grid 0315 if ISW == 1 0316 P2 = (WBLP.*elf(IBLP)')+(WBRP.*elf(IBRP)')+(WTLP.*elf(JBLP)')+... 0317 (WTRP.*elf(JBRP)'); 0318 elseif ISW == 2 0319 tmpa = (WBLW.*elf(IBLW)')+(WBRW.*elf(IBRW)')+(WTLW.*elf(JBLW)')+... 0320 (WTRW.*elf(JBRW)'); 0321 elseif ISW == 3 0322 tmpb = (WBLW.*elf(IBLW)')+(WBRW.*elf(IBRW)')+(WTLW.*elf(JBLW)')+... 0323 (WTRW.*elf(JBRW)'); 0324 end 0325 0326 kline = kline+data_seg; 0327 end 0328 0329 % Rotate the winds to equatorial lat/lon 0330 tmpc = (COEFF1.*tmpa)+(COEFF2.*tmpb); 0331 tmpd = (COEFF1.*tmpb)-(COEFF2.*tmpa); 0332 0333 k = 0; 0334 0335 % Convert the data to a gridded array instead of a vector 0336 for j = 0:NRR-1 0337 I1 = NTRNT(j+1); 0338 I2 = NTRNS(j+1); 0339 for i = I1+1:I2 0340 k = k+1; 0341 PRESS(i,NRR-j) = P2(k); 0342 U10E(i,NRR-j) = tmpc(k); 0343 U10N(i,NRR-j) = tmpd(k); 0344 end 0345 end 0346 0347 % Interpolate the pressure data onto the FVCOM grid 0348 slp(:,(m/3)+1) = interp2(cs3x_lon,cs3x_lat,PRESS',Mobj.lon,Mobj.lat); 0349 0350 % Calculate the coordinates of the centre of each FVCOM grid element 0351 xelement = mean(Mobj.lon(Mobj.tri(:,:)),2); 0352 yelement = mean(Mobj.lat(Mobj.tri(:,:)),2); 0353 0354 % Interpolate the wind data onto the FVCOM grid 0355 uwnd(:,(m/3)+1) = interp2(cs3x_lon,cs3x_lat,U10E',xelement,yelement); 0356 vwnd(:,(m/3)+1) = interp2(cs3x_lon,cs3x_lat,U10N',xelement,yelement); 0357 0358 0359 % if FirstInt==1 0360 % % Interpolate the pressure onto the model grid nodes 0361 % for i = 1:Mobj.nVerts 0362 % xmet = cs3x_lonstart; 0363 % ymet = cs3x_latstart; 0364 % i1 = 1; 0365 % i2 = 2; 0366 % j1 = 1; 0367 % j2 = 2; 0368 % 0369 % while Mobj.lon(i) > xmet+cs3x_loninc 0370 % i1 = i1+1; 0371 % i2 = i2+1; 0372 % xmet = xmet+cs3x_loninc; 0373 % end 0374 % 0375 % while Mobj.lat(i) > ymet+cs3x_latinc 0376 % j1 = j1+1; 0377 % j2 = j2+1; 0378 % ymet = ymet+cs3x_latinc; 0379 % end 0380 % 0381 % a1 = (xmet + cs3x_loninc - Mobj.lon(i)).*... 0382 % (ymet + cs3x_latinc - Mobj.lat(i)); 0383 % a2 = (xmet + cs3x_loninc - Mobj.lon(i)) .* (Mobj.lat(i)-ymet); 0384 % a3 = (Mobj.lon(i) - xmet).*(ymet + cs3x_latinc - Mobj.lat(i)); 0385 % a4 = (Mobj.lon(i) - xmet).*(Mobj.lat(i) - ymet); 0386 % 0387 % P(i,(m/3)+1) = (PRESS(i1,j1).*a1 + PRESS(i1,j2).*a2 + PRESS(i2,j1).*a3 ... 0388 % + PRESS(i2,j2).*a4) ./cs3x_area; 0389 % 0390 % i1j1(1,i) = i1; 0391 % i1j1(2,i) = j1; 0392 % A1(1,i) = a1; 0393 % A2(1,i) = a2; 0394 % A3(1,i) = a3; 0395 % A4(1,i) = a4; 0396 % 0397 % % tidy up 0398 % clear xmet ymet i1 i2 j1 j2 a1 a2 a3 a4 0399 % end 0400 % 0401 % % Interpolate the winds onto the model grid nodes 0402 % for i = 1:Mobj.nElems 0403 % xmet = cs3x_lonstart; 0404 % ymet = cs3x_latstart; 0405 % i1 = 0; 0406 % i2 = 1; 0407 % j1 = 0; 0408 % j2 = 1; 0409 % 0410 % % Calculate the coordinates of the centre of the element 0411 % xelement = mean(Mobj.lon(Mobj.tri(i,:))); 0412 % yelement = mean(Mobj.lat(Mobj.tri(i,:))); 0413 % 0414 % while xelement > xmet+cs3x_loninc 0415 % i1 = i1+1; 0416 % i2 = i2+1; 0417 % xmet = xmet+cs3x_loninc; 0418 % end 0419 % 0420 % while yelement > ymet+cs3x_latinc 0421 % j1 = j1+1; 0422 % j2 = j2+1; 0423 % ymet = ymet+cs3x_latinc; 0424 % end 0425 % 0426 % a1 = (xmet + cs3x_loninc - xelement).*... 0427 % (ymet + cs3x_latinc - yelement); 0428 % a2 = (xmet + cs3x_loninc - xelement) .* (yelement-ymet); 0429 % a3 = (xelement - xmet).*(ymet + cs3x_latinc - yelement); 0430 % a4 = (xelement - xmet).*(yelement - ymet); 0431 % 0432 % Uwind(i,(m/3)+1) = (U10E(i1,j1).*a1 + U10E(i1,j2).*a2 + U10E(i2,j1).*a3 ... 0433 % + U10E(i2,j2).*a4) ./cs3x_area; 0434 % Vwind(i,(m/3)+1) = (U10N(i1,j1).*a1 + U10N(i1,j2).*a2 + U10N(i2,j1).*a3 ... 0435 % + U10N(i2,j2).*a4) ./cs3x_area; 0436 % 0437 % i1j1(3,i) = i1; 0438 % i1j1(4,i) = j1; 0439 % A1(2,i) = a1; 0440 % A2(2,i) = a2; 0441 % A3(2,i) = a3; 0442 % A4(2,i) = a4; 0443 % 0444 % % tidy up 0445 % clear xmet ymet i1 i2 j1 j2 a1 a2 a3 a4 xelement y element 0446 % end 0447 % else 0448 % % If it's not the first pass, there's no need to re-calculate i1j1 0449 % % or A1, A2, A3, A4. 0450 % for i = 1:Mobj.nVerts 0451 % P(i,(m/3)+1) = (PRESS(i1j1(1,i),i1j1(2,i)).*A1(1,i) + ... 0452 % PRESS(i1j1(1,i),i1j1(2,i)+1).*A2(1,i) + ... 0453 % PRESS(i1j1(1,i)+1,i1j1(2,i)).*A3(1,i) + ... 0454 % PRESS(i1j1(1,i)+1,i1j1(2,i)+1).*A4(1,i))./cs3x_area; 0455 % end 0456 % 0457 % for i = 1:Mobj.nElems 0458 % Uwind(i,(m/3)+1) = (U10E(i1j1(3,i),i1j1(4,i)).*A1(2,i) + ... 0459 % U10E(i1j1(3,i),i1j1(4,i)+1).*A2(2,i) + ... 0460 % U10E(i1j1(3,i)+1,i1j1(4,i)).*A3(2,i) + ... 0461 % U10E(i1j1(3,i)+1,i1j1(4,i)+1).*A4(2,i))./cs3x_area; 0462 % Vwind(i,(m/3)+1) = (U10N(i1j1(3,i),i1j1(4,i)).*A1(2,i) + ... 0463 % U10N(i1j1(3,i),i1j1(4,i)+1).*A2(2,i) + ... 0464 % U10N(i1j1(3,i)+1,i1j1(4,i)).*A3(2,i) + ... 0465 % U10N(i1j1(3,i)+1,i1j1(4,i)+1).*A4(2,i))./cs3x_area; 0466 % end 0467 % end 0468 0469 % FirstInt = false; 0470 end 0471 0472 % Convert data time to Modified Julian Day time for FVCOM 0473 time = datevec(time); 0474 MJD_time = greg2mjulian(time(:,1),time(:,2),time(:,3),time(:,4),time(:,5),... 0475 time(:,6)); 0476 0477 % Add slp, uwnd, vwnd and time to Mobj 0478 Mobj.Met.slp.node = slp; 0479 Mobj.Met.uwnd.data = uwnd; 0480 Mobj.Met.vwnd.data = vwnd; 0481 Mobj.Met.time = MJD_time; 0482 0483 %% mjd0 = nml file start time (inputConf.startDate) 0484 % mjd1 = nml file finish time (inputConf.endDate) 0485 % mjd2 = variable time from met data file 0486 % mjd3 = variable time from temporary met data file 0487 % mjd4 = variable time from surge data file 0488 % mjd5 = time at start of surge data file, the reference start point of the 0489 % harmonic analysis 0490 0491 % What are the surge parameters? (Don't know if I need these in this file. 0492 % Let's keep them here for now) 0493 % surge_lonstart = -19.916667; 0494 % surge_latstart = 40.0555555; 0495 % surge_loninc = 1/6; 0496 % surge_latinc = 1/9;