Replace the interpolation of a single velocity with the interpolation of the...

Replace the interpolation of a single velocity with the interpolation of the original u and v components and instead calculate the velocity at the end. This means the ASCII files needed for the mean flow can more easily be generated

fvcom_prepro/get_POLCOMS_meanflow.m

@@ -45,6 +45,8 @@ function Mobj = get_POLCOMS_meanflow(Mobj, files)
 %
 % Revision history
 %    2013-02-20 First version.
+%    2013-02-26 Add interpolation of the u and v vectors separately and
+%    then calculate the interpolated velocity at the end.
 %
 %==========================================================================
 
@@ -70,7 +72,9 @@ fvlat = Mobj.lat(Mobj.obc_nodes(Mobj.obc_nodes ~= 0));
 % Number of boundary nodes
 nf = sum(Mobj.nObcNodes);
 
-fvmf = nan(nf, nt); % FVCOM interpolated mean flow
+% fvmf = nan(nf, nt); % FVCOM interpolated mean flow
+fvmfu = nan(nf, nt); % FVCOM interpolated mean flow vectors
+fvmfv = nan(nf, nt); % FVCOM interpolated mean flow vectors
 
 if ftbverbose
     tic
@@ -87,7 +91,7 @@ for t = 1:nt
     % Calculate the velocity and direction (?) from the u and v vectors.
     % This approach means we don't have to interpolate u and v separately
     % (which is expensive) but just do velocity on its own.
-    pcvel3 = sqrt(pcu3.^2 + pcv3.^2);
+%     pcvel3 = sqrt(pcu3.^2 + pcv3.^2);
        
     % For the velocity, it appears we don't need to have a depth varying
     % value (instead mean flow is scaled according to the values in MFDIST,
@@ -95,17 +99,23 @@ for t = 1:nt
     % we'll calculate a depth average velocity here for the time being.
     % Transpose the array here so it's (x, y) rather than the original (y,
     % x).
-    pcvel3mean = mean(pcvel3, 3)';
+%     pcvel3mean = mean(pcvel3, 3)';
+    pcu3mean = mean(pcu3, 3)';
+    pcv3mean = mean(pcv3, 3)';
 
     % We need to create a mask to eliminate land values and apply it to the
     % depth averaged values.
     mask = squeeze(pc.depth.data(:, :, end, t))' >= 0;
     
-    tpcvel3mean = pcvel3mean;
+%     tpcvel3mean = pcvel3mean;
+    tpcu3mean = pcu3mean;
+    tpcv3mean = pcv3mean;
     tlon = lon;
     tlat = lat;
     % Clear out the land values and flatten the arrays.
-    tpcvel3mean(mask) = [];
+%     tpcvel3mean(mask) = [];
+    tpcu3mean(mask) = [];
+    tpcv3mean(mask) = [];
     tlon(mask) = [];
     tlat(mask) = [];
     
@@ -125,93 +135,33 @@ for t = 1:nt
         % parallelisation.
         plon = tlon(ixy);
         plat = tlat(ixy);
-        pvel = tpcvel3mean(ixy);
+%         pvel = tpcvel3mean(ixy);
+        pu = tpcu3mean(ixy);
+        pv = tpcv3mean(ixy);
 
         % Use a triangulation to do the horizontal interpolation.
-        trivel = TriScatteredInterp(plon', plat', pvel', 'natural');
-        ivelobc(i) = trivel(fx, fy);
+%         trivel = TriScatteredInterp(plon', plat', pvel', 'natural');
+        triu = TriScatteredInterp(plon', plat', pu', 'natural');
+        triv = TriScatteredInterp(plon', plat', pv', 'natural');
+%         ivelobc(i) = trivel(fx, fy);
+        iuobc(i) = triu(fx, fy);
+        ivobc(i) = triv(fx, fy);
 
         % Check if we have NaNs (mostly if the position is outside the
         % model domain).
-        if isnan(ivelobc(i))
+        if isnan(ivelobc(i)) || isnan(iuobc(i)) || isnan(ivobc(i))
             warning('FVCOM boundary node at %f, %f is outside the PML POLCOMS-ERSEM domain. Setting to the closest PML POLCOMS-ERSEM value.', fx, fy)
-            ivelobc(i) = tpcvel2(ii(1));
+%             ivelobc(i) = tpcvel2(ii(1));
+            iuobc(i) = tpcu2(ii(1));
+            ivobc(i) = tpcv2(ii(1));
         end
     end
 
-  
-%     for j = 1:nz
-%         % Now extract the relevant layer from the 3D subsets. Transpose the
-%         % data to be (x, y) rather than (y, x).
-%         pcvel2 = pcvel3(:, :, j)';
-%         pcdepth2 = squeeze(pc.depth.data(:, :, j, t))';
-% 
-%         % Create new arrays which will be flattened when masking (below).
-%         tpcvel2 = pcvel2;
-%         tpcdepth2 = pcdepth2;
-%         tlon = lon;
-%         tlat = lat;
-%         
-%         % Create and apply a mask to remove values outside the domain. This
-%         % inevitably flattens the arrays, but it shouldn't be a problem
-%         % since we'll be searching for the closest values in such a manner
-%         % as is appropriate for an unstructured grid (i.e. we're assuming
-%         % the PML POLCOMS-ERSEM data is irregularly spaced).
-%         mask = tpcdepth2 > 20000;
-%         tpcvel2(mask) = [];
-%         tpcdepth2(mask) = [];
-%         % Also apply the masks to the position arrays so we can't even find
-%         % positions outside the domain, effectively meaning if a value is
-%         % outside the domain, the nearest value to the boundary node will
-%         % be used.
-%         tlon(mask) = [];
-%         tlat(mask) = [];
-% 
-%         % Preallocate the intermediate results arrays.
-%         ivelobc = nan(nf, 1);
-%         idepthobc = nan(nf, 1);
-%         
-%         % Speed up the tightest loop with a parallelized loop.
-%         parfor i = 1:nf
-%             % Now we can do each position within the 2D layer.
-% 
-%             fx = fvlon(i);
-%             fy = fvlat(i);
-% 
-%             [~, ii] = sort(sqrt((tlon - fx).^2 + (tlat - fy).^2));
-%             % Get the n nearest nodes from PML POLCOMS-ERSEM data (more?
-%             % fewer?).
-%             ixy = ii(1:16);
-% 
-%             % Get the variables into static variables for the
-%             % parallelisation.
-%             plon = tlon(ixy);
-%             plat = tlat(ixy);
-%             pvel = tpcvel2(ixy);
-%             pdepth = tpcdepth2(ixy);
-%             
-%             % Use a triangulation to do the horizontal interpolation.
-%             trivel = TriScatteredInterp(plon', plat', pvel', 'natural');
-%             triz = TriScatteredInterp(plon', plat', pdepth', 'natural');
-%             ivelobc(i) = trivel(fx, fy);
-%             idepthobc(i) = triz(fx, fy);
-%             
-%             % Check both, though if one is NaN, they both will be.
-%             if isnan(ivelobc(i)) || isnan(idepthobc(i))
-%                 warning('FVCOM boundary node at %f, %f is outside the PML POLCOMS-ERSEM domain. Setting to the closest PML POLCOMS-ERSEM value.', fx, fy)
-%                 ivelobc(i) = tpcvel2(ii(1));
-%                 idepthobc(i) = tpcdepth2(ii(1));
-%             end
-%         end
-%         
-%         % Put the results in this intermediate array.
-%         ivelz(:, j) = ivelobc;
-%         idepthz(:, j) = idepthobc;
-%     end
-
     % The horizontally-interpolated values in the final FVCOM results
     % array.
-    fvmf(:, t) = ivelobc;
+%     fvmf(:, t) = ivelobc;
+    fvmfu(:, t) = iuobc;
+    fvmfv(:, t) = ivobc;
 
     if ftbverbose
         fprintf('done.\n')
@@ -222,7 +172,11 @@ if ftbverbose
 end
 
 % Stick the values in the mesh structure.
-Mobj.velocity = fvmf;
+Mobj.meanflow_u = fvmfu;
+Mobj.meanflow_v = fvmfv;
+% Mobj.velocity = fvmf;
+Mobj.velocity = sqrt(Mobj.meanflow_u.^2 + Mobj.meanflow_v.^2);
+
 
 % Convert the current times to Modified Julian Day (this is a bit ugly).
 pc.time.all = strtrim(regexp(pc.time.units, 'since', 'split'));