/model/model_handles2.m
MATLAB | 1544 lines | 1135 code | 150 blank | 259 comment | 79 complexity | 32b376cbf647d8712f9e4053adb5389c MD5 | raw file
Possible License(s): GPL-2.0
- %--------------------------------------------------------------------------
- % This file is part of SHEL SHallow-water numerical modEL
- %
- % SHEL is free software: you can redistribute it and/or modify
- % it under the terms of the GNU General Public License as published by
- % the Free Software Foundation, either version 3 of the License, or
- % (at your option) any later version.
- %
- % SHEL is distributed in the hope that it will be useful,
- % but WITHOUT ANY WARRANTY; without even the implied warranty of
- % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- % GNU General Public License for more dr.etails.
- %
- % You should have received a copy of the GNU General Public License
- % along with SHEL. If not, see <http://www.gnu.org/licenses/>.
- %--------------------------------------------------------------------------
-
- function s = model_handles2(r)
- %function s = model_handles2
- %
- %Call it this way:
- % s = model_handles;
- % y = s.model_compute_input(); % --> the parenthesis ARE REQUIRED!
- % y = s.model_compute_output( handles );
-
- %globals are replaced by a big structure
- r.time = nan;
- r.timetr = nan;
- r.dt = nan;
- r.dttr = nan;
- r.ltr = nan;
- r.frame = nan;
- r.duration = nan;
- r.outputdt = nan;
- r.L = nan;
- r.outputL = nan;
-
- r.frame = nan;
- r.printG = nan;
- r.film = nan;
- r.mov = nan;
- r.myfullfile = nan;
- r.mydir = nan;
- r.optprint = nan;
-
- r.dx = nan;
- r.dy = nan;
- r.M = nan;
- r.N = nan;
- r.x0 = nan;
- r.y0 = nan;
- r.x = nan;
- r.y = nan;
- r.x_u = nan;
- r.y_u = nan;
- r.x_v = nan;
- r.y_v = nan;
- r.x_w = nan;
- r.y_w = nan;
-
- r.loadbathymetry = nan;
- r.bathymetryfile = nan;
- r.varsinfile = nan;
-
- %Test-cases
- r.tc_bump = nan;
- r.tc_geostrophic = nan;
- r.tc_isla = nan;
- r.tc_taylor = nan;
-
- %boundary condition
- r.u_closed = nan;
- r.v_closed = nan;
- %r.bc_closed;
- r.noslip = nan;
-
- %r.step bathym
- r.step = nan;
-
- %bump
- r.bump_x0 = nan;
- r.bump_y0 = nan;
- r.bump_sx = nan;
- r.bump_sy = nan;
- r.bump_d0 = nan;
-
- %island
- r.isla_x0 = nan;
- r.isla_y0 = nan;
- r.isla_R = nan;
-
- %r.tracer
- r.tracer = nan;
- r.TrXo = nan;
- r.TrYo = nan;
- r.TrR = nan;
-
- %Grid and bathymetry related
- r.land = nan;
- r.d0 = nan;
- r.d0_r = nan;
-
- %T-cells
- r.Tr = nan;
- r.eta0 = nan;
- r.eta_old = nan;
- r.eta = nan;
- r.eta_new = nan;
- r.H_old = nan;
- r.H = nan;
- r.H_new = nan;
- r.d = nan;
- r.mask = nan;
- r.masknan = nan;
- r.u_t = nan;
- r.v_t = nan;
-
- %Diagnostic quantities
- r.ke = nan;
- r.pe = nan;
- r.eke = nan;
- r.gradux = nan;
- r.graduy = nan;
- r.gradvx = nan;
- r.gradvy = nan;
- r.curl_t = nan;
- r.potvorticity_t = nan;
- r.strechrate_t = nan;
- r.shearrate_t = nan;
- r.sqstrechrate_t = nan;
- r.sqshearrate_t = nan;
- r.enstrophy_t = nan;
- r.sqstrain_t = nan;
- r.divergence_t = nan;
- r.sqdivergence_t = nan;
- r.okuboweiss_t = nan;
-
- %r.u-cells
- r.u_a = nan;
- r.u_old = nan;
- r.u = nan;
- r.u_new = nan;
- r.mask_u = nan;
-
- %r.v-cells
- r.v_a = nan;
- r.v_old = nan;
- r.v = nan;
- r.v_new = nan;
- r.mask_v = nan;
-
- %Z-cells
- r.curl_w = nan;
- r.potvorticity_w = nan;
- r.shearrate_w = nan;
- r.sqshearrate_w = nan;
- r.enstrophy_w = nan;
- r.okuboweiss_w = nan;
-
- %r.time-dependent r.properties (domain integrated)
- r.volume = nan;
- r.iKe = nan;
- r.iPe = nan;
- r.vtime = nan;
- r.iVort = nan;
- r.iEnst = nan;
- r.iSqStrech = nan;
- r.iSqShear = nan;
- r.iSqStrain = nan;
- r.iOWeiss = nan;
- r.iMomentumU = nan;
- r.iMomentumV = nan;
- r.iEke = nan;
-
- r.K = nan;
- r.gama = nan;
-
- s.model_compute_input = @ComputeModelInput;
- s.model_compute_output = @ComputeModelOutput;
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Here starts the model Inputs part
- %%Initial conditions
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- function ComputeModelInput(r)
- %function ComputeModelInputs
-
- %GR : clears figure of images
- r.frame = 1;
- r.time = 0.;
- r.timetr = r.time + 2 * r.dt;
- ComputeTime(r);
- makecoordinates(r);
- fillfields(r);
-
- function ComputeTime(r)
- %function updateTime
-
- r.L = getL(r.duration,r.dt);
- r.outputL = getL(r.outputdt, r.dt);
-
- function L_L = getL (duration_L, dt_L)
- %function L_L = getL (duration_L, dt_L)
-
- L_L = floor(duration_L / dt_L);
-
- function makecoordinates(r)
- %function makecoordinates
-
- if r.loadbathymetry
- r.varsinfile = who('-file',r.bathymetryfile);
- if ~isempty(r.varsinfile)
- load(r.bathymetryfile, 'r.d');
- siz = size(r.d);
- r.M = siz(1);
- r.N = siz(2);
- end
- end
-
- r.x0 = r.dx/2;
- r.y0 = r.dy/2;
-
- %T-cell
- r.x = (1:r.M)' * (1:r.N);
- r.y = (1:r.M)' * (1:r.N);
- for j = 1:r.N
- r.x(:,j) = r.x0 + (r.x(:,j)/j - 1) * r.dx;
- end
- for i = 1:r.M
- r.y(i,:) = r.y0 + (r.y(i,:)/i - 1) * r.dy;
- end
-
- %r.u-cell
- r.x_u = (1:r.M+1)' * (1:r.N);
- r.y_u = (1:r.M+1)' * (1:r.N);
- for j = 1:r.N
- r.x_u(:,j) = r.x0 + (r.x_u(:,j)/j - 1) * r.dx - r.dx/2;
- end
- for i = 1:r.M
- r.y_u(i,:) = r.y0 + (r.y_u(i,:)/i - 1) * r.dy;
- end
-
- %r.v-cell
- r.x_v = (1:r.M)' * (1:r.N+1);
- r.y_v = (1:r.M)' * (1:r.N+1);
- for j = 1:r.N
- r.x_v(:,j) = r.x0 + (r.x_v(:,j)/j - 1) * r.dx;
- end
- for i = 1:r.M
- r.y_v(i,:) = r.y0 + (r.y_v(i,:)/i - 1) * r.dy - r.dy/2;
- end
-
- %W-cell
- r.x_w = (1:r.M)' * (1:r.N);
- r.y_w = (1:r.M)' * (1:r.N);
- for j = 1:r.N
- r.x_w(:,j) = r.x0 + (r.x_w(:,j)/j - 1) * r.dx - r.dx/2;
- end
- for i = 1:r.M
- r.y_w(i,:) = r.y0 + (r.y_w(i,:)/i - 1) * r.dy - r.dy/2;
- end
-
- function fillfields(r)
- %function fillfields
-
- % Arakawa C staggered grid (Arakawa 1966)
- %
- % --j y
- % | --- r.v --- r.v ------- r.v r.u -- T -- r.u
- % i | | | | | |
- % x r.u T r.u T r.u T | r.v |
- % | | | | | |
- % --- r.v --- r.v ------- r.v r.u -- T -- r.u
- % T-cell r.u-cell r.v-cell
- %
- % T(r.eta, r.H, r.Tr, r.d, f, r.mask, r.u(r.u, r.mask_u r.v(r.v, r.mask_v
- % x, y, wind and bottom r.x_u, r.y_u) r.x_v, r.y_v)
- % stress)
- % --j y
- % | T--- r.u ---T
- % i | |
- % x r.v W r.v
- % | |
- % T--- r.u ---T
- % W-cell
- %
- % W(zr.eta - curl of r.v, vorticity)
-
-
- %Depth field (bathymetry)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- r.land = -99;
-
- %%Constant depth or read from bathymetry file
- if r.tc_taylor %Submarine mount equals 30% of total mean depth
- r.d = r.d0 - cylindermap(r,zeros(size(r.d)), 0.3 * r.d0, r.bump_sx, r.bump_sy, ...
- r.M * r.dx * .5 - r.bump_x0, ...
- r.N * r.dy * .5 - r.bump_y0);
- elseif r.loadbathymetry && sum(strcmp(r.varsinfile,'r.d'))
- load(r.bathymetryfile, 'r.d');
- else
- r.d = r.d0 * ones(r.M,r.N);
- end
-
- %%changing depth
- if r.step
- r.d = maker.step(r.d,r.d0_r.step);
- end
-
- %%Closed boundaries: boundaries are filled with r.land
- if r.v_closed
- for j = 1:r.N-1:r.N;
- for i = 1:r.M
- r.d(i,j) = r.land;
- end
- end
- end
-
- if r.u_closed
- for i = 1:r.M-1:r.M;
- for j = 1:r.N
- r.d(i,j) = r.land;
- end
- end
- end
-
- %%Inner islands (random)
- %%%arguments: depth, io, jo, r
- if r.tc_isla
- r.d = makeisland(r,r.d, r.M * r.dx * .5 - r.isla_x0, r.N * r.dy * .5 - r.isla_y0, r.isla_R);
- end
-
- %WATERMASK matrices%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %T-cell
- %%all-water
- r.mask = ones(r.M,r.N);
-
- %%eventual islands and closed domain
- for j=1:r.N
- for i=1:r.M
- if r.d(i,j) < 0.
- r.mask(i,j) = 0;
- end
- end
- end
-
- %r.u-Cell and r.v-cell
- r.mask_u = ones(r.M+1,r.N);
- r.mask_v = ones(r.M,r.N+1);
- for j = 1:r.N
- for i = 1:r.M
- if r.mask( i, j) == 0
- r.mask_u( i, j) = r.mask(i,j);
- r.mask_u(i+1, j) = r.mask(i,j);
- r.mask_v( i, j) = r.mask(i,j);
- r.mask_v( i,j+1) = r.mask(i,j);
- end
- end
- end
-
- if r.noslip
- %r.u-cells
- for j = 2:r.N-1
- for i = 2:r.M
- r.mask_u(i,j) = r.mask(i ,j+1) * r.mask(i ,j-1) * r.mask(i-1,j+1) * r.mask(i-1,j-1);
- end
- end
- %r.v-Cells
- for j = 2:r.N
- for i = 2:r.M-1
- r.mask_v(i,j) = r.mask(i+1, j) * r.mask(i+1,j-1) * r.mask(i-1, j) * r.mask(i-1,j-1);
- end
- end
- end
-
- %Water level fields%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if (r.tc_bump || r.tc_geostrophic)
- r.eta_old = bumpmap(r,r.bump_d0, r.bump_sx, r.bump_sy, r.M * r.dx * .5 - r.bump_x0, ...
- r.N * r.dy * .5 - r.bump_y0);
- else
- if r.loadbathymetry && sum(strcmp(r.varsinfile,'r.eta'))
- load(r.bathymetryfile, 'r.eta');
- r.eta_old = r.eta;
- else
- r.eta_old = r.eta0 * ones(r.M,r.N);
- end
- end
- r.eta = r.eta_old;
- r.eta_new = r.eta;
-
- %Bottom-to-surface depth (always positive below the water surface)
- %r.H = r.eta + r.d.
- r.H_old = r.eta_old + r.d;
- r.H =r.H_old;
- r.H_new = r.H_old;
-
- if r.tracer
- %Unitary concentration...
- r.Tr = zeros(r.M,r.N);
- r.Tr = cylindermap(r,r.Tr, 1., r.TrR, r.TrR, r.M * r.dx * .5 - r.TrXo, ...
- r.N * r.dy * .5 - r.TrYo);
- end
-
- %r.u velocity fields%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if r.tc_geostrophic
- r.u_old = ugeo(r.eta);
- elseif r.tc_taylor % 3 cm/s r.u-flow
- r.u_old = zeros(size(r.mask_u));
- %%uncomment one of the following lines
- r.u_old = 0.03 * ones(size(r.mask_u));
- %r.u_old(1:r.M:r.M+1,:) = 0.03 * ones(size(r.mask_u(1:r.M:r.M+1,:)));
- else
- if r.loadbathymetry && sum(strcmp(r.varsinfile,'r.u'))
- load(r.bathymetryfile, 'r.u');
- r.u_old = r.u;
- else
- r.u_old = zeros(r.M+1,r.N);
- end
- end
- r.u = r.u_old;
- r.u_new = r.u;
- r.u_a = r.u;
-
- %r.v velocity fields%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if r.tc_geostrophic
- r.v_old = vgeo(r.eta);
- else
- if r.loadbathymetry && sum(strcmp(r.varsinfile,'r.v'))
- load(r.bathymetryfile, 'r.v');
- r.v_old = r.v;
- else
- r.v_old = zeros(r.M,r.N+1);
- end
- end
- r.v = r.v_old;
- r.v_new = r.v;
- r.v_a = r.v;
-
- %call this function when you need
- %to manually initialize the level
- %in the Taylor column flow.
- if r.tc_taylor
- makeGeostLevelFromU(r)
- end
-
- %Z field%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- r.curl_w = zeros(r.M+1,r.N+1);
- r.potvorticity_w = zeros(r.M+1,r.N+1);
- r.enstrophy_w = zeros(r.M+1,r.N+1);
- r.shearrate_w = zeros(r.M+1,r.N+1);
- r.sqshearrate_w = zeros(r.M+1,r.N+1);
- r.okuboweiss_w = zeros(r.M+1,r.N+1);
-
- %Visualization matrices%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- r.masknan = ones(r.M,r.N);
- for j=1:r.N
- for i=1:r.M
- if r.mask(i,j) == 0
- r.masknan(i,j) = nan;
- else
- r.masknan(i,j) = 1;
- end
- end
- end
-
- r.u_t = zeros(r.M,r.N);
- r.v_t = zeros(r.M,r.N);
-
- r.ke = zeros(r.M,r.N);
- r.pe = zeros(r.M,r.N);
-
- r.curl_t = zeros(r.M,r.N);
- r.potvorticity_t = zeros(r.M,r.N);
- r.strechrate_t = zeros(r.M,r.N);
- r.shearrate_t = zeros(r.M,r.N);
- r.enstrophy_t = zeros(r.M,r.N);
- r.sqstrain_t = zeros(r.M,r.N);
- r.okuboweiss_t = zeros(r.M,r.N);
- r.sqshearrate_t = zeros(r.M,r.N);
- r.sqstrechrate_t = zeros(r.M,r.N);
- r.divergence_t = zeros(r.M,r.N);
- r.sqdivergence_t = zeros(r.M,r.N);
-
- %Eddy kinetic energy
- r.gradux=zeros(r.M,r.N);
- r.graduy=zeros(r.M,r.N);
- r.gradvx=zeros(r.M,r.N);
- r.gradvy=zeros(r.M,r.N);
- r.eke=zeros(r.M,r.N);
-
- %One-dimensional diagnostic variables
- %Energy initialization
- r.vtime = nan(1,L);
- r.iKe = nan(1,L);
- r.iPe = nan(1,L);
- r.iEke = nan(1,L);
- r.volume = nan(1,L);
- r.iVort = nan(1,L);
- r.iEnst = nan(1,L);
- r.iSqShear = nan(1,L);
- r.iSqStrech = nan(1,L);
- r.iSqStrain = nan(1,L);
- r.iOWeiss = nan(1,L);
- r.iMomentumU = nan(1,L);
- r.iMomentumV = nan(1,L);
-
- ComputeDiagnostics(1);
-
- %cylinder field
- function map_L = cylindermap(r,map_L, l, sx, sy, x_L, y_L)
-
- ind = find ( (x - x_L).^2 / sx^2 + (y - y_L).^2 / sy^2 < 1. );
- map_L(ind) = l;
-
- %bump test-case field
- function map_L = bumpmap(r,l, sx, sy, x_L, y_L)
-
- map_L = zeros(r.M,r.N);
-
- for i = 1:r.M
- for j = 1:r.N
- map_L(i,j) = l * exp( - ( ( x(i,j) - x_L )^2 / sx^2 + ( y(i,j) - y_L )^2 / sy^2 ) );
- end
- end
-
- function depth = makeisland(r,depth, x_L, y_L, r)
- %function depth = makeisland(depth, x_L, y_L, r)
-
- for i = 1:r.M
- for j = 1:r.N
- if r^2 > ( x(i,j) - x_L )^2 + ( y(i,j) - y_L )^2
- depth(i,j) = r.land;
- end
- end
- end
-
- function ugeo_L = ugeo(r, r.eta_L)
- %function ugeo_L = ugeo(r.eta_L)
-
- j_L = floor(r.N/2)+1;
- sig_L = 1.;
- ug_L = zeros(r.M,r.N);
-
- for i=1:r.M
- for j=1:r.N
- ug_L(i,j) = g / f * ( j - j_L ) / sig_L^2 * r.eta_L(i,j);
- end
- end
- ugeo_L = InterpolTU(ug_L);
-
- function vgeo_L = vgeo(r,r.eta_L)
- %function vgeo_L = vgeo(r.eta_L)
-
- i_L = floor(r.M/2)+1;
- sig_L = 1.;
- vg_L = zeros(r.M,r.N);
-
- for i=1:r.M
- for j=1:r.N
- vg_L(i,j) = - g / f * ( i - i_L ) / sig_L^2 * r.eta_L(i,j);
- end
- end
- vgeo_L = InterpolTV(vg_L);
-
- function bath = makestep(bath,depth)
- %function bath = makestep(bath,depth)
- %Creates a r.stepped bathymetry
- sz = size(bath);
- M = sz(1);
- N = sz(2);
- j0 = N/5;
-
- for i = 1:M
- for j = 1:N
- if j < 0.9 * i + j0
- bath(i,j) = depth;
- end
- end
- end
-
- function makeGeostLevelFromU(r)
-
- for j=2:r.N
- for i=1:r.M
- r.eta(i,j) = r.eta(i,j-1) * r.mask(i,j-1) - r.dy * f/g * (...
- r.mask_u(i,j) * r.u(i,j) + r.mask_u(i+1,j) * r.u(i+1,j) ...
- + r.mask_u(i,j-1) * r.u(i,j-1) + r.mask_u(i+1,j-1) * r.u(i+1,j-1) ...
- ) / ( ...
- r.mask_u(i,j) + r.mask_u(i+1,j) + r.mask_u(i,j-1) + r.mask_u(i+1,j-1) ...
- );
- r.eta(i,j) = r.mask(i,j) * r.eta(i,j);
- end
- end
-
- r.eta_old = r.eta;
- r.eta_new = r.eta;
-
- r.H_old = r.eta_old + r.d;
- r.H =r.H_old;
- r.H_new = r.H_old;
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Here starts the model Outputs part
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- function ComputeModelOutput( r, handles )
- %function ComputeModelOutputs( handles ) [old ComputeModel]
-
- %r.time parameters and output parameters
-
- s = GUI_common_handles;
-
- %Movie initialization
- clear mex;
- tic;
- ol = 0;
- movie = false;
-
- %Define the type of output
- switch optprint
- case {11, 12, 13, 14, 15}
- movie = true;
- end
-
- % Write output of initial condition to disk.
- if printG
- if movie
- %Open a new avi file.
- filename = sprintf('%s.avi', [mydir,'/',myfullfile]);
- mov = avifile(filename,'fps', 10,'quality',100,'compression','None');
- end
- s.GUI_printit(frame, handles, optprint);
- frame = frame + 1;
- end
-
- % Main loop in r.time
- r.ltr = 0;
- for l = 1 : L
-
- %Takes care of the numerical scheme
- %for momentum and for the r.tracer
- %This line can easily be replaced
- %with a switch, choosing one kernel
- %among many (mex file kernel, other
- %kernels ...).
- kernel_SHEL(r,l);
-
- % Take care of plotting the outputs and write them to disk.
- if ol == outputL
-
- ol=0;
-
- if film
- s.GUI_plotmodel(handles);
- pause(0.1); %s
- end
- if printG
- s.GUI_printit(frame, handles, optprint);
- frame = frame + 1;
- end
-
- %Display r.time
- strin = sprintf('r.time = %0.3g s', r.time);
- set(handles.timetext,'String',strin);
- pause(0.01); %s pause (required if I want to update the clock).
-
- end
- ol=ol+1;
-
- %End of main loop in r.time
- end
-
- % Take care of writing the final output on disk
- if printG
- if movie
- mov = close(mov);
- end
- end
-
- if ~film
- s.GUI_plotmodel(handles);
- end
-
- % Display simulation statistics and performance
- cput = toc;
- disp(sprintf('Elapsed r.time: %0.5f s\nTime per iteration: %0.5f s\r.N', cput, cput/L));
- if r.tracer
- disp(sprintf('Made %0.3d r.tracer iterations on a total of %0.3d iterations', r.ltr, l));
- end
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Here starts the Kernel part
- %%the numerical scheme
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- function kernel_SHEL(r,l)
- %function kernel_SHEL (old ComputeModel)
-
- %Explicit leapfrog scheme
- %combined with Asselin-Roberts filter
- %(the filter supresses the computational mode
- %and prevents decoupling between odd/even timer.steps)
- %
- % Arakawa C staggered grid (Arakawa 1966)
- %
- % --j y
- % | --- r.v --- r.v ------- r.v r.u -- T -- r.u
- % i | | | | | |
- % x r.u T r.u T r.u T | r.v |
- % | | | | | |
- % --- r.v --- r.v ------- r.v r.u -- T -- r.u
- % T-cell r.u-cell r.v-cell
- % r.M x r.N r.M+1 x r.N r.M x r.N+1
- % T(r.eta, r.H, r.Tr, r.d, f, r.mask, r.u(r.u, r.mask_u r.v(r.v, r.mask_v
- % x, y, wind and bottom r.x_u, r.y_u) r.x_v, r.y_v)
- % stress, curl(r.v,r.u,0))
- %
- % --j y
- % | --- r.u ---
- % i | |
- % x r.v W r.v
- % | |
- % --- r.u ---
- % T-cell
- % r.M+1 x r.N+1
- % W( curl(r.u,r.v,0) )
- % (deformation rate = curl(-r.u,r.v,0)
- %%%
-
- %Update the r.time
- r.time = r.time + r.dt;
-
- %Compute the r.time scheme of the momentum+continuity model
- ComputeLeapfrog;
-
- %Compute all the tracers
- if r.tracer
- if r.timetr <= r.time
- [r.Tr, r.dttr] = ComputeTracer_FT(r.Tr, r.K);
- r.timetr = r.time + r.dttr;
- r.ltr = r.ltr + 1;
- end
- end
-
- %Numerics to compute the diagnostic variables
- ComputeDiagnostics(r, l);
-
- function ComputeLeapfrog(r)
- %Shallow-water equations solver
- %
-
- %T-cells
- r.dtmask;
-
- %r.u-cells
- r.dtmask_u;
-
- %r.v-cells
- r.dtmask_v;
-
- %% Solve the spatial schemes in the interior of the domain
- % For the water elevation
- RHSr.eta = ComputeContinuity(r.H, r.eta, r.u, r.v, r.mask);
- r.eta_new = r.eta_old + 2 * r.dt * RHSr.eta .* r.mask;
- r.H_new = r.eta_new + r.d;
- ind = find(r.H_new < 0);
- r.eta_new(ind) = -r.d;
- dtmask = (r.eta_new - r.eta_old) * .5 ./ RHSr.eta .* r.mask;
- %Must build the dtmask_u and dtmask_v in order to preserve momentum
- %too.
- dtmask_u = buildmask_U(dtmask,r.mask_u);
- dtmask_v = buildmask_V(dtmask',r.mask_v')';
-
- RHSr.eta = ComputeContinuity(r.H, r.eta, r.u, r.v, dtmask);
-
- % For the r.u and r.v components of velocity
- % % The old 10%-faster way (but a lot more error-prone)
- % [RHSu, RHSv] = ComputeUV2D_old;
- % Switch r.v <-> r.u; r.v_old <-> r.u_old; r.dy <-> r.dx and -1. <-> 1. (coriolis)
- % Transpose the matrices.
- RHSu = ComputeSpaceU_CS( r.H, r.H_old, r.eta, ...
- r.u, r.u_old, r.v, r.v_old, ...
- dtmask_u, r.dx, r.dy, 1.);
- RHSv = ComputeSpaceU_CS( r.H', r.H_old', r.eta', ...
- r.v', r.v_old', r.u', r.u_old', ...
- dtmask_v', r.dy, r.dx, -1.)';
-
- %% LEAPFROG r.time-SCHEME
- % Water level r.time scheme
- r.eta_new = r.eta_old + 2 * RHSr.eta .* r.mask;
- r.H_new = r.eta_new + r.d;
- %r.u and r.v r.time schemes inside the domain.
- %Leapfrog means that 2xdt is passed as teh r.time r.step.
- r.v_new = ComputeTimeU_FT(r.v_old', r.v_new', RHSv', r.H_old', r.H_new', r.mask_v', 2.)';
- r.u_new = ComputeTimeU_FT(r.u_old, r.u_new, RHSu, r.H_old, r.H_new, r.mask_u, 2.);
-
- %% r.time scheme at the northern/southern and eastern/southern OB (2 x r.d
- %% <- LeapFrog)
-
- %Southern and northern boundaries (along r.v -> transpose inputs)
- %(transposed inputs -> transposed outputs).
- [r.eta_new, r.H_new, r.v_new, r.u_new] = ComputeOB_U( ...
- r.eta_new', r.eta_old', r.H_new', r.H_old', r.d',...
- r.v_new', r.mask_v', ...
- r.u_new', r.u_old', r.mask_u', ...
- 2*r.dt, r.dy ...
- );
- %Eastern and western boundaries (along r.u)
- %The previous outputs were transposed, so we transpose
- %them back as inputs.
- [r.eta_new, r.H_new, r.u_new, r.v_new] = ComputeOB_U( ...
- r.eta_new', r.eta_old, r.H_new', r.H_old, r.d,...
- r.u_new', r.mask_u, ...
- r.v_new', r.v_old, r.mask_v, ...
- 2*r.dt, r.dx ...
- );
-
- %% Asselin-Robert filter
- r.eta = r.eta + r.gama * ( r.eta_old - 2 * r.eta + r.eta_new );
- r.u = r.u + r.gama * ( r.u_old - 2 * r.u + r.u_new );
- r.v = r.v + r.gama * ( r.v_old - 2 * r.v + r.v_new );
-
- %% Update matrices
- r.eta_old = r.eta;
- r.eta = r.eta_new;
-
- r.H_old = r.H;
- r.H = r.H_new;
-
- r.u_old = r.u;
- r.u = r.u_new;
-
- r.v_old = r.v;
- r.v = r.v_new;
-
- function RHSr.eta = ComputeContinuity(r.H, r.eta, r.u, r.v, r.mask)
- %function ComputeContinuity
- %Shallow-water equations solver
-
- % r.H;
- % r.eta;
- % r.u;
- % r.v;
- % r.mask;
- % rHSr.eta;
- r.dx;
- r.dy;
-
- [r.M,r.N] = size(r.eta);
- RHSr.eta = zeros(r.M,r.N);
-
- %Compute along r.u then along r.v
- RHSr.eta(2:r.M-1,2:r.N-1) = - ( ...
- ComputeContinuityU(r.H,r.u,r.mask,r.dx) + ComputeContinuityU(r.H',r.v',r.mask',r.dy)' ...
- );
-
- function RHS = ComputeContinuityU(r.H,r.u,r.mask,r.dx)
- [r.M,r.N] = size(r.H);
- RHS = ( ...
- r.mask(3:r.M,2:r.N-1) .* .5 ...
- .* ( r.H(3:r.M,2:r.N-1) + r.H(2:r.M-1,2:r.N-1) ) .* r.u(3:r.M,2:r.N-1) ... %i+1, i, i+1
- - r.mask(1:r.M-2,2:r.N-1) .* .5 ...
- .* ( r.H(2:r.M-1,2:r.N-1) + r.H(1:r.M-2,2:r.N-1) ) .* r.u(2:r.M-1,2:r.N-1) ... %i, i-1, i
- )/r.dx;
-
- function RHSu = ComputeSpaceU_CS( ...
- r.H, r.H_old, r.eta, ...
- r.u, r.u_old, r.v, r.v_old, ...
- r.mask_u, r.dx, r.dy, ...
- signcoriolis ...
- )
- %% function newu = ComputeSpaceU_CS
- %Shallow-water equations solver for the r.u component
- %of the flow Centred in Space (CS).
- %
- % Centred differences (CD) or Centred in Space (CS).
- %
- %NOTE: Any viscous and frictional term
- %must be evaluated in tindex(i-1,r.M)e l-1.
- %
- % switch r.v <-> r.u; j <-> i and r.dy <-> r.dx relative to ComputeU2D
- % switch r.v_old <-> r.u_old
- % take special care for the coriolis force term
-
- %% r.variables
- r.g;
- r.f;
- r.Nu;
- r.wind;
- r.bottom;
- r.pressure;
- r.coriolis;
- r.wall;
-
- r.uwind;
- r.vwind;
- rho0;
- rho_air;
- r.lb;
- r.karman;
-
- wall = false;
-
- %% r.u spatial scheme
- [r.M, r.N] = size(r.eta);
- RHSu = zeros(size(r.u));
- r.H_u = zeros(size(r.u));
- r.H_old_u = zeros(size(r.u));
- v_old_u = zeros(size(r.u));
- v_u = zeros(size(r.u));
- r.H_u(2:r.M,2:r.N-1) = twoaverage_u(r.H,r.M,r.N);
- r.H_old_u(2:r.M,2:r.N-1) = twoaverage_u(r.H_old,r.M,r.N);
- v_old_u(2:r.M,2:r.N-1) = fouraverage_u(r.v_old(:,2:r.N),r.M, r.N-2);
- v_u(2:r.M,2:r.N-1) = fouraverage_u(r.v(:,2:r.N),r.M, r.N-2);
-
- RHSu(2:r.M,2:r.N-1) = r.mask_u(2:r.M,2:r.N-1) .* ( ...
- ...%advection along r.u ...
- - 0.2500 / r.dx * ( ... %i-1, i
- r.mask_u(3:r.M+1,2:r.N-1) .* ( r.u(3:r.M+1,2:r.N-1) + r.u(2:r.M,2:r.N-1) ).^2 .* r.H(2:r.M,2:r.N-1) ... %i+1, i, i+1, i, i
- - r.mask_u(1:r.M-1,2:r.N-1) .* ( r.u(2:r.M,2:r.N-1) + r.u(1:r.M-1,2:r.N-1) ).^2 .* r.H(1:r.M-1,2:r.N-1) ... %i, i-1, i, i-1, i-1
- ) ...
- ...%advection along r.v
- - 0.0625 / r.dy * ( ... %i-1, i
- ( ...
- r.mask_u(2:r.M,3:r.N) .* ( r.H(2:r.M,3:r.N) + r.H(2:r.M,2:r.N-1) + r.H(1:r.M-1,3:r.N) + r.H(1:r.M-1,2:r.N-1) ) ... % j+1 ; j ; i-1,j+1 ; i-1
- .* ( r.v(2:r.M,3:r.N) + r.v(1:r.M-1,3:r.N) ) .* ( r.u(2:r.M,3:r.N) + r.u(2:r.M,2:r.N-1) ) ...% j+1 ; i-1,j+1 ; j+1 ; i
- - r.mask_u(2:r.M,1:r.N-2) .* ( r.H(2:r.M,2:r.N-1) + r.H(2:r.M,1:r.N-2) + r.H(1:r.M-1,2:r.N-1) + r.H(1:r.M-1,1:r.N-2) ) ... % i ; j-1 ; i-1 ; i-1,j-1
- .* ( r.v(2:r.M,2:r.N-1) + r.v(1:r.M-1,2:r.N-1) ) .* ( r.u(2:r.M,2:r.N-1) + r.u(2:r.M,1:r.N-2) ) ...% i ; i-1,j-1 ; i ; j-1
- ) ...
- ) ...
- ...%Coriolis Force (r.u -> +; r.v -> -)
- + coriolis * ( ...
- signcoriolis * f * v_u(2:r.M,2:r.N-1) .* r.H_u(2:r.M,2:r.N-1) ... % i ; i-1 .
- )...
- ...%Pressure gradient Force
- + pressure * ( ...
- - g * r.H_u(2:r.M,2:r.N-1) .* ( r.eta(2:r.M,2:r.N-1) - r.eta(1:r.M-1,2:r.N-1) ) / r.dx ...
- ... % i ; i-1 ; i ; i-1
- )...
- ...%Viscosity along r.u
- + ( ... % i ; i-1
- + Nu * r.mask_u(3:r.M+1,2:r.N-1) .* r.H_old(2:r.M,2:r.N-1) .* ( r.u_old(3:r.M+1,2:r.N-1) - r.u_old(2:r.M,2:r.N-1) ) / r.dx ... % i ; i+1 ; i ; i
- - Nu * r.mask_u(1:r.M-1,2:r.N-1) .* r.H_old(1:r.M-1,2:r.N-1) .* ( r.u_old(2:r.M,2:r.N-1) - r.u_old(1:r.M-1,2:r.N-1) ) / r.dx ... % i-1 ; i ; i-1; i-1
- ) / r.dx ... % i ; i-1
- ...%Viscosity along r.v
- + ( ... % i ; i-1 ; i ; i-1 .
- + Nu / r.dy * r.mask_u(2:r.M,3:r.N) .* r.H_old(2:r.M,2:r.N-1) .* ( r.u_old(2:r.M,3:r.N) - r.u_old(2:r.M,2:r.N-1) ) ... % (j;j+1;i-1,j+1;i-1);j+1 ; j
- - Nu / r.dy * r.mask_u(2:r.M,1:r.N-2) .* r.H_old(2:r.M,1:r.N-2) .* ( r.u_old(2:r.M,2:r.N-1) - r.u_old(2:r.M,1:r.N-2) ) ... % (j-1;j;i-1;i-1,j-1);j ; j-1
- ) / r.dy ...
- );
-
- %Add bottom stress. Look out, so that the bottom stress doesn't
- %reverts the direction of momentum ... This wouldn't be very realistic.
- if bottom || wall
- stress = zeros(size(RHSu(2:r.M,2:r.N-1)));
- if bottom
- stress = bottomstress( r.u_old(2:r.M,2:r.N-1), v_old_u(2:r.M,2:r.N-1), ...
- r.H_old_u(2:r.M,2:r.N-1), ...
- lb, karman );
- end
- if wall
- stress = stress + bottomstress(r.u_old(2:r.M,2:r.N-1), 0, r.dy, lb, karman) ...
- .* ( ...
- (1 - r.mask(1:r.M-1,3:r.N)) .* r.H_old(1:r.M-1,3:r.N) + ...
- (1 - r.mask(2:r.M,3:r.N)) .* r.H_old(2:r.M,3:r.N) + ...
- (1 - r.mask(1:r.M-1,1:r.N-2)) .* r.H_old(1:r.M-1,1:r.N-2) + ...
- (1 - r.mask(2:r.M,1:r.N-2)) .* r.H_old(2:r.M,1:r.N-2) ...
- ) * 0.5 / r.dy;
- end
- ind = find( (RHSu(2:r.M,2:r.N-1) - r.mask_u(2:r.M,2:r.N-1) .* stress) ./ RHSu(2:r.M,2:r.N-1) > 0 );
- RHSu(ind) = RHSu(ind) - r.mask_u(ind) .* stress(ind);
- ind = find( (RHSu(2:r.M,2:r.N-1) - r.mask_u(2:r.M,2:r.N-1) .* stress) ./ RHSu(2:r.M,2:r.N-1) <= 0 );
- RHSu(ind) = 0.;
- end
-
- %Add wind stress.
- if wind
- RHSu(2:r.M,2:r.N-1) = RHSu(2:r.M,2:r.N-1) ...
- + r.mask_u(2:r.M,2:r.N-1) .* windstress(...
- uwind * ones(size(r.u(2:r.M,2:r.N-1))), ...
- vwind * ones(size(r.u(2:r.M,2:r.N-1))), ...
- rho0, rho_air ...
- );
- end
-
- function r.u_new = ComputeTimeU_FT(r.u_old, r.u_new, RHSu, r.H_old, r.H_new, r.mask_u, r.dt)
- %function r.u_new = ComputeTimeU_FT(r.u_old, r.u_new, RHSu, r.H_old, r.H_new, r.mask_u, r.dt)
- %
- %Numerical solver with the Forward in r.time (FT) Euler r.time r.stepping method.
- %
- [r.M,r.N] = size(r.H_old);
-
- r.u_new(2:r.M,2:r.N-1) = r.mask_u(2:r.M,2:r.N-1) .* ( ...
- r.u_old(2:r.M,2:r.N-1) .* ( r.H_old(2:r.M,2:r.N-1) + r.H_old(1:r.M-1,2:r.N-1) ) ...
- + 2 * r.dt * RHSu(2:r.M,2:r.N-1) ...
- ) ...
- ./ ( r.H_new(2:r.M,2:r.N-1) + r.H_new(1:r.M-1,2:r.N-1) );
-
- function [r.eta_new, r.H_new, r.u_new, r.v_new] = ComputeOB_U( ...
- r.eta_new, r.eta_old, r.H_new, r.H_old, r.d,...
- r.u_new, r.mask_u, ...
- r.v_new, r.v_old, r.mask_v, ...
- r.dt, r.dx ...
- )
- %function computeOB
- %Solve the north, south, east and west boundary condition
- %for the r.eta, r.u and r.v variables.
-
- r.g;
-
- %BC conditions
- r.flather;
- r.neumann;
- radiatr.etan;
- radiatelevel;
-
- flather = 0;
- neumann = 0;
- radiatr.etan = 0;
- radiatelevel = 1;
-
- [r.M,r.N] = size(r.d);
-
- %% Eastern (1,:) and Western (r.M,:) boundaries
- if radiatelevel
- %Waterlevel: GW radiation
- r.eta_new(1:r.M-1:r.M,:) = r.mask_u(1:r.M:r.M+1,:) .* ( ...
- r.eta_old(1:r.M-1:r.M,:)- r.dt / r.dx .* sqrt( g * r.H_old(1:r.M-1:r.M,:) ) ...
- .* ( r.eta_old(1:r.M-1:r.M,:) - r.eta_old(2:r.M-3:r.M-1,:) ) ...
- );
- %no r.mask below, cuz r.H must yield -99 where r.land is.W
- r.H_new(1:r.M-1:r.M,:) = r.eta_new(1:r.M-1:r.M,:) + r.d(1:r.M-1:r.M,:);
- end
- if flather
- %Normal velocity: Flather
- r.u_new(1:r.M:r.M+1,:) = [-1 * ones(1,r.N); ones(1,r.N)] .* r.mask_u(1:r.M:r.M+1,:) ...
- .* sqrt( g ./ r.H_new(1:r.M-1:r.M,:) ) .* r.eta_new(1:r.M-1:r.M,:);
- elseif neumann
- %Normal velocity: Null-gradient
- r.u_new(1:r.M:r.M+1,1:r.N) = r.mask_u(1:r.M:r.M+1,1:r.N) .* r.u_new(2:r.M-2:r.M,1:r.N);
- end
- if radiatr.etan
- %Tangential velocity: GW radiation
- r.v_new(1:r.M-1:r.M,2:r.N) = r.mask_v(1:r.M-1:r.M,2:r.N) .* ...
- ( ...
- r.v_old(1:r.M-1:r.M,2:r.N) .* ( r.H_old(1:r.M-1:r.M,2:r.N) + r.H_old(1:r.M-1:r.M,1:r.N-1) ) ...
- - r.dt / r.dx .* sqrt( g * .5 * ( r.H_old(1:r.M-1:r.M,2:r.N) + r.H_old(1:r.M-1:r.M,1:r.N-1) ) ) ...
- .* ( r.v_old(1:r.M-1:r.M,2:r.N) - r.v_old( 2:r.M-3:r.M-1,2:r.N) ) ...
- ) ...
- ./ ( r.H_new(1:r.M-1:r.M,2:r.N) + r.H_new(1:r.M-1:r.M,1:r.N-1) );
- end
-
-
- function [r.Tr, r.dttr] = ComputeTracer_FT(r.Tr, K_L)
- %Advection-diffusion r.tracer equation solver
- %Explicit forward in r.time and upwind
- %
- r.dx;
- r.dy;
- r.dt;
-
- %T-cells
- r.H;
- r.H_new;
- r.mask;
-
- %r.u-cells
- r.u_a;
- r.mask_u;
-
- %r.v-cells
- r.v_a;
- r.mask_v;
-
- %% Solve the r.tracer spatial scheme in the whole of the domain
- RHSt = r.mask .* ( ...
- ComputeSpaceT_UP_U(r.H,r.u_a,r.mask_u, r.Tr, K_L, r.dx) ... %Along r.u
- + ComputeSpaceT_UP_U(r.H',r.v_a',r.mask_v',r.Tr',K_L,r.dy)' ... %then along r.v
- );
-
- %% FIND THE OPTIMAL r.tracer r.dt AS A MULTIPLE OF MOMENTUM r.dt
- % within the reasonable limit of maximum 100 * r.dt.
- minRHSt = min(min(RHSt));
- ind = find(RHSt == minRHSt);
- r.dttr = min( floor( abs(- r.H(ind) .* r.Tr(ind) ./ minRHSt) / r.dt ), 1000) ...
- * r.dt;
-
- %% r.tracer FORWARD-IN-r.time SCHEME
- Tr_new = r.mask .* (r.H .* r.Tr + r.dttr * RHSt) ./ r.H_new;
-
- %% Update matrices
- r.Tr = Tr_new;
-
- function RHSt = ComputeSpaceT_UP_U(r.H, r.u, r.mask_u, r.Tr, K_L, dx_L)
- %function RHSt = ComputeSpaceT_UP_U(r.H, r.u, r.mask_u, r.Tr, r.K, r.dx)
- %
- %UPWIND
- %
- r.nullgrad;
- nullgrad = false;
-
- [r.M,r.N]=size(r.H);
- Hm = zeros(r.M+1,r.N);
- TrUm = zeros(r.M+1,r.N);
- TrUp = zeros(r.M+1,r.N);
- TrDif = zeros(r.M+1,r.N);
-
- %Computing quantities for the interior of the domain
- TrUm(1:r.M,:) = .5 * (abs(r.u(1:r.M,:)) - r.u(1:r.M,:)) .* r.Tr;
- TrUp(2:r.M+1,:) = .5 * (abs(r.u(2:r.M+1,:)) + r.u(2:r.M+1,:)) .* r.Tr;
- TrDif(2:r.M,:) = r.Tr(1:r.M-1,:) - r.Tr(2:r.M,:);
- Hm(2:r.M,:) = .5 * ( r.H(1:r.M-1,:) + r.H(2:r.M,:) );
-
- %Computing the quantities for the boundaries
- if nullgrad %Considering null-gradient:
- TrUm(r.M+1,:) = .5 * (abs(r.u(r.M+1,:)) - r.u(r.M+1,:)) .* r.Tr(r.M,:);
- TrUp(1,:) = .5 * (abs(r.u(1,:)) + r.u(1,:)) .* r.Tr(1,:);
- TrDif(1:r.M:r.M+1,:) = 0.;
- else %else null value is considered
- TrUm(r.M+1,:) = .5 * (abs(r.u(r.M+1,:)) - r.u(r.M+1,:)) * 0.;
- TrUp(1,:) = .5 * (abs(r.u(1,:)) + r.u(1,:)) * 0.;
- TrDif(1,:) = - r.Tr(1,:);
- TrDif(r.M+1,:) = r.Tr(r.M,:);
- end
-
- %This algo considers that the bathymetry has always a null gradient
- %at the boundary:
- Hm(1:r.M:r.M+1,:) = r.H(1:r.M-1:r.M,:);
-
- %Computing the fluxes at each face of the r.M+1 faces
- FFluxU = ComputeFaceFluxT_UP(Hm, TrUp, TrUm, TrDif, dx_L, K_L);
-
- %The T-cell increment is the flux balance.
- RHSt = r.mask_u(1:r.M,:) .* FFluxU(1:r.M,:) ...
- - r.mask_u(2:r.M+1,:) .* FFluxU(2:r.M+1,:);
-
- function RHSt = ComputeFaceFluxT_UP(Hm, TrUp, TrUm, TrDif, dx_L, K_L)
- %function RHSt = ComputeFaceFluxT_UP(Hm, TrUp, TrUm, TrDif, dx_L, K_L)
- %Pure upwind scheme.
- RHSt = Hm .* ( ...
- TrUp - TrUm ...
- + K_L / dx_L * TrDif ...
- ) / dx_L;
-
- function av = fouraverage_u(r.v, r.M, r.N)
- %% function av = fouraverage_u(r.v, r.mask, r.M, r.N)
-
- av = .25 * ( ...
- r.v(2:r.M,1:r.N) + ...
- r.v(2:r.M,2:r.N+1) + ...
- r.v(1:r.M-1,1:r.N) + ...
- r.v(1:r.M-1,2:r.N+1) ...
- );
-
- function av = twoaverage_u( r.H, r.M, r.N)
- av = .5 * ( ...
- r.H(2:r.M,2:r.N-1) + r.H(1:r.M-1,2:r.N-1) ...
- );
-
- function tb_L = bottomstress( u_L, v_L, r.H_L, lb, karman)
- % Bottom stress
- %function tb_L = bottomstress( u_L, v_L, r.H_L )
- %
- %r.H_L is the total depth height
-
- cd_L = karman / log( ( 0.5 .* r.H_L + lb ) / lb ) ;
- % Checked-out Blumberg and Mellor 1987
- % Double-checked with Pietrzak 2002
- cd_L = max(cd_L .* cd_L, 0.0025);
- tb_L = cd_L .* sqrt ( u_L .* u_L + v_L .* v_L ) .* u_L;
-
- function u_L = windstress(u_L, v_L, rho0, rho_air)
- %% Wind stress
- %function tv_L = windstress(v_L)
-
- vmod_L = sqrt( u_L .* u_L + v_L .* v_L );
-
- if vmod_L < 6.
- Cd_L = -.26097 .* vmod_L + 2.4316;
- else
- Cd_L = .064986 .* vmod_L + .44053;
- end
-
- Cd_L = Cd_L * .001;
-
- u_L = Cd_L .* rho_air / rho0 .* vmod_L .* u_L;
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %% Here starts the diagnostics part
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function diag = ComputeDiagnostics(l)
- %function ComputeDiagnostics
- %
- %Computes global:
- % r.volume (m3)
- % kinetic energy (J)
- % potential energy (J)
- % total energy (J)
- % vorticity (m2/s)
- % enstrophy (J)
- % squared normal strain (J)
- % squared shear strain (J)
- % squared strain (J)
- % Okubo-Weiss
- % momentum (r.M^3 * r.M / s)a
- %
- %Computes local:
- % kinetic energy (J)
- % potential energy (J)
- % total energy (J)
- % curl (1/s)
- % divergence (1/s)
- % stretch (normal strain) (1/s)
- % shear (shear strain) (1/s)
- % squared strain (1/s)
- % enstrophy (1/s2)
- % OkuboWeiss (1/s2)
- %
-
- r.u_a;
- r.v_a;
- r.u;
- r.v;
- r.u_t;
- r.v_t;
- r.H;
- r.eta;
- r.M;
- r.N;
- r.dx;
- r.dy;
- r.dt;
- rho0;
- r.g;
- r.K;
- r.mask;
- r.ke;
- r.pe;
- r.strechrate_t;
- r.shearrate_w;
- r.divergence_t;
- r.sqstrain_t;
- r.sqstrechrate_t;
- r.sqshearrate_w;
- r.curl_w;
- r.enstrophy_t;
- r.enstrophy_w;
- r.okuboweiss_t;
-
- %Eddy kinetic energy
- r.gradux;
- r.graduy;
- r.gradvx;
- r.gradvy;
- r.eke;
-
- %r.time-dependent r.properties
- r.time;
- r.volume;
- r.iKe;
- r.iPe;
- r.vtime;
- r.iVort;
- r.iEnst;
- r.iSqStrech;
- r.iSqShear;
- r.iSqStrain;
- r.iOWeiss;
- r.iMomentumU;
- r.iMomentumV;
- r.iEke;
-
- dA = r.dx * r.dy; %m2
-
- %Computes the average velocity field
- r.u_a = ((l - 1) * r.u_a + r.u ) / l;
- r.v_a = ((l - 1) * r.v_a + r.v ) / l;
-
- %Eddy kinetic energy
- r.gradux = (r.u(2:r.M+1,:) - r.u(1:r.M,:)) / r.dx;
- r.gradvy = (r.v(:,2:r.N+1) - r.v(:,1:r.N)) / r.dy;
- r.graduy(2:r.M,2:r.N-2) = (r.u(2:r.M,3:r.N-1) - r.u(2:r.M,2:r.N-2)) / r.dy;
- r.gradvx(2:r.M-2,2:r.N) = (r.v(3:r.M-1,2:r.N) - r.v(2:r.M-2,2:r.N)) / r.dx;
-
- %Compute diagnostic fields
- %Curl
- computeCurl_w; % 1 / s
-
- %Strain rate
- r.strechrate_t = r.mask .* ((r.u(2:r.M+1,:) - r.u(1:r.M,:)) * r.dy ...
- - (r.v(:,2:r.N+1) - r.v(:,1:r.N)) * r.dx) / dA; % 1 / s
- computeShearRate_w % 1 / s
-
- %Growth rate
- r.divergence_t = r.mask .* ((r.u(2:r.M+1,1:r.N) - r.u(1:r.M,1:r.N)) * r.dy ...
- + (r.v(1:r.M,2:r.N+1) - r.v(1:r.M,1:r.N)) * r.dx) / dA; % 1/ s
-
- %Potential vorticity
- computePotentialVorticity_w
-
- %This is the so-called scalar of Arakawa-Okubo-Weiss.
- %Computes the local enstrophy, squared strain rate and Okubo-Weiss
- %parameter %Look for Weiss
- %La Jolla Institute Report from 1981.
- %Physica r.d: Nonlinear 1991.
- %Look for Okubo-Weiss parameter/function
- %Look for Arakawa 1966 paper.
- r.enstrophy_w = .5 * ( r.curl_w.^2); %1 / s2
- r.enstrophy_t = interpolFtoT(r.enstrophy_w, r.mask, r.M, r.N);
-
- r.sqshearrate_w = .5 * r.shearrate_w.^2;
- r.sqshearrate_t = interpolFtoT(r.sqshearrate_w, r.mask, r.M, r.N);
-
- r.sqstrechrate_t = .5 * r.strechrate_t.^2 .* r.mask;
-
- r.sqstrain_t = r.sqstrechrate_t + r.sqshearrate_t;
-
- r.sqdivergence_t = .5 * r.divergence_t.^2 .* r.mask;
-
- r.okuboweiss_w = (r.sqshearrate_w - r.enstrophy_w);
- r.okuboweiss_t = interpolFtoT(r.okuboweiss_w, r.mask, r.M, r.N);
- r.okuboweiss_t = r.okuboweiss_t + r.sqstrechrate_t .* r.mask;
- %vers?o do Aires do escalar de Okubo-Weiss
- r.okuboweiss_t = r.okuboweiss_t - r.sqdivergence_t;
-
- %%interpolate from r.u and r.v to T cells
- r.u_t = .5 * ( r.u(1:r.M,:) + r.u(2:r.M+1,:) );
- r.v_t = .5 * ( r.v(:,1:r.N) + r.v(:,2:r.N+1) );
-
- Hnoland = r.H;
- Hnoland(find(r.H < -1)) = 0.;
-
- %r.volume
- %r.volume(l) = dA * sum(sum(Hnoland)); %m3
- %Perturbed r.volume
- r.volume(l) = dA * sum(sum(r.eta .* r.mask)); %m3
-
- %Momentum (actually it's only the integrated velocity)
- r.iMomentumU(l) = dA * sum(sum(r.u_t .* Hnoland,2)); % kg * r.M / s
- r.iMomentumV(l) = dA * sum(sum(r.v_t .* Hnoland,1)); % kg * r.M / s
-
- %Kinetic energy
- r.ke = .5 * rho0 * dA * (r.u_t.^2 + r.v_t.^2) .* Hnoland .* r.mask; % kg * r.M * r.M / s2 = J
- r.iKe(l) = sum(sum(r.ke)); % kg * r.M * r.M / s2 = J
-
- %Potential energy
- r.pe = .5 * rho0 * g * dA .* r.eta.^2 .* r.mask; %kg r.M * r.M / s2 = J
- r.iPe(l) = sum(sum(r.pe)); %kg r.M * r.M / s2 = J
-
- %Eddy kinetic energy
- r.eke = r.eke + rho0 * dA * r.K * 2 * r.dt * (r.gradux.^2 + r.gradvy.^2 + r.graduy.^2 + r.gradvx.^2) ...
- .* Hnoland .* r.mask;
- r.iEke(l) = sum(sum(r.eke));
-
- %r.vorticity
- r.iVort(l) = sum(sum(r.curl_w)) * dA; % m2 / s
-
- %r.enstrophy
- r.iEnst(l) = sum(sum(r.enstrophy_w)) * dA; % m2 / s2
-
- %r.squared strain
- r.iSqShear(l) = sum(sum(r.sqshearrate_w))* dA;
- r.iSqStrech(l) = sum(sum(r.sqstrechrate_t)) * dA;
- r.iSqStrain(l) = sum(sum(r.sqstrain_t)) * dA;
-
- %r.Okubo-Weiss
- r.iOWeiss(l) = sum(sum(r.okuboweiss_t)) * dA; %<--- too large numerical error
- %r.iOWeiss(l) = r.iSqStrain(l) - r.iEnst(l);
-
- %r.time
- r.vtime(l) = r.time;
-
- diag = [...
- r.volume(l), ...
- r.iKe(l), ...
- r.iPe(l), ...
- r.iVort(l), ...
- r.iEnst(l), ...
- r.iSqStrain(l), ...
- r.iOWeiss(l) ...
- ];
-
- function computeCurl_w
- %function computeCurl_w % 1 / s
-
- r.curl_w;
- r.curl_t;
- r.mask;
- r.mask_v;
- r.mask_u,
- r.u;
- r.v;
- r.M;
- r.N;
- r.dy;
- r.dx;
-
- dA = r.dx .* r.dy;
-
- %Curl of velocity in the interior...
- r.curl_w(2:r.M,2:r.N) = ((r.mask_v(2:r.M,2:r.N) .* r.v(2:r.M,2:r.N) ...
- - r.mask_v(1:r.M-1,2:r.N) .* r.v(1:r.M-1,2:r.N)) .* r.dy ...
- - (r.mask_u(2:r.M,2:r.N) .* r.u(2:r.M,2:r.N) ...
- - r.mask_u(2:r.M,1:r.N-1) .* r.u(2:r.M,1:r.N-1)) .* r.dx) ./ dA; % 1 / s
- %... at the corners %Thank you for the conversation with F.L.!
- %lower-left SW
- r.curl_w(1,1) = ((r.mask_v(1,1) * r.v(1,1)) * r.dy ...
- - (r.mask_u(1,1) * r.u(1,1)) * r.dx)/dA; % 1 / s
- %lower-right SE
- r.curl_w(1,r.N+1) = ((r.mask_v(1,r.N+1) * r.v(1,r.N+1)) * r.dy ...
- - ( - r.mask_u(1,r.N) * r.u(1,r.N)) * r.dx)/dA; % 1 / s
- %top-left NW
- r.curl_w(r.M+1,1) = (( - r.mask_v(r.M,1) * r.v(r.M,1)) * r.dy ...
- - (r.mask_u(r.M+1,1) * r.u(r.M+1,1)) * r.dx)/dA; % 1 / s
- %top-right NE
- r.curl_w(r.M+1,r.N+1) = ((- r.mask_v(r.M,r.N+1) * r.v(r.M,r.N+1)) * r.dy ...
- - (- r.mask_u(r.M+1,r.N) * r.u(r.M+1,r.N)) * r.dx)/dA; % 1 / s
- %... at the western boundary
- r.curl_w(2:r.M,1) = ((r.mask_v(2:r.M,1) .* r.v(2:r.M,1) - r.mask_v(1:r.M-1,1) .* r.v(1:r.M-1,1)) * r.dy ...
- - (r.mask_u(2:r.M,1) .* r.u(2:r.M,1)) * r.dx)/dA; % 1 / s
- %... at the eastern boundary
- r.curl_w(2:r.M,r.N+1) = ((r.mask_v(2:r.M,r.N+1) .* r.v(2:r.M,r.N+1) - r.mask_v(1:r.M-1,r.N+1) .* r.v(1:r.M-1,r.N+1)) * r.dy ...
- - (- r.mask_u(2:r.M,r.N) .* r.u(2:r.M,r.N)) * r.dx)/dA; % 1 / s
- %... at the southern boundary
- r.curl_w(1,2:r.N) = ((r.mask_v(1,2:r.N) .* r.v(1,2:r.N)) * r.dy ...
- - (r.mask_u(1,2:r.N) .* r.u(1,2:r.N) - r.mask_u(1,1:r.N-1) .* r.u(1,1:r.N-1)) * r.dx)/dA; % 1 / s
- %... at the northern boundary
- r.curl_w(r.M+1,2:r.N) = ((- r.mask_v(r.M,2:r.N) .* r.v(r.M,2:r.N)) * r.dy ...
- - (r.mask_u(r.M+1,2:r.N) .* r.u(r.M+1,2:r.N) - r.mask_u(r.M+1,1:r.N-1) .* r.u(r.M+1,1:r.N-1)) * r.dx)/dA; % 1 / s
- % interpolation from the {Z,W,F}-cell to the T-Cell
- r.curl_t = interpolFtoT(r.curl_w, r.mask, r.M, r.N);
-
- function computeShearRate_w
- %function computeShearRate_w 1 / s
- %
-
- r.shearrate_w;
- r.shearrate_t;
- r.mask;
- r.u;
- r.v;
- r.M;
- r.N;
- r.dy;
- r.dx;
-
- dA = r.dx * r.dy;
-
- %Curl of (-r.u,r.v,0)
- %interior of the domain...
- r.shearrate_w(2:r.M,2:r.N) = ((r.v(2:r.M,2:r.N) - r.v(1:r.M-1,2:r.N)) * r.dy ...
- + (r.u(2:r.M,2:r.N) - r.u(2:r.M,1:r.N-1)) * r.dx) / dA;
- %bottom-left corner (SW)
- r.shearrate_w(1,1) = ((r.v(1,1)) * r.dy ...
- + (r.u(1,1)) * r.dx) / dA;
- %bottom-right corner (SE)
- r.shearrate_w(1,r.N+1) = ((r.v(1,r.N+1)) * r.dy ...
- + (- r.u(1,r.N)) * r.dx) / dA;
- %top-left corner (NW)
- r.shearrate_w(r.M+1,1) = ((- r.v(r.M,1)) * r.dy ...
- + (r.u(r.M+1,1)) * r.dx) / dA;
- %top-right corner (NE)
- r.shearrate_w(r.M+1,r.N+1) = ((- r.v(r.M,r.N+1)) * r.dy ...
- + (- r.u(r.M+1,r.N)) * r.dx) / dA;
- %western boundary
- r.shearrate_w(2:r.M,1) = ((r.v(2:r.M,1) - r.v(1:r.M-1,1)) * r.dy ...
- + (r.u(2:r.M,1)) * r.dx) / dA;
- %eastern boundary
- r.shearrate_w(2:r.M,r.N+1) = ((r.v(2:r.M,r.N+1) - r.v(1:r.M-1,r.N)) * r.dy ...
- + (- r.u(2:r.M,r.N)) * r.dx) / dA;
- %southern boundary
- r.shearrate_w(1,2:r.N) = ((r.v(1,2:r.N)) * r.dy ...
- + (r.u(1,2:r.N) - r.u(1,1:r.N-1)) * r.dx) / dA;
- %northern boundary
- r.shearrate_w(r.M+1,2:r.N) = ((- r.v(r.M,2:r.N)) * r.dy ...
- + (r.u(r.M+1,2:r.N) - r.u(r.M+1,1:r.N-1)) * r.dx) / dA;
- % interpolate from F-cell to a T-cell
- r.shearrate_t = interpolFtoT(r.shearrate_w, r.mask, r.M, r.N);
-
- function computePotentialVorticity_w
- %function computePotentialVorticity_w % 1 / s / r.M
- %
-
- r.potvorticity_t;
- r.mask;
- r.curl_w;
- r.potvorticity_w;
- r.coriolis;
- r.f;
- r.H;
- r.M;
- r.N;
-
- %Adding coriolis
- if coriolis
- r.potvorticity_w = r.curl_w + f; % ??? r.curl_w * dA + f???
- end
-
- %potential vorticity
- %within the domain...
- r.potvorticity_w(2:r.M,2:r.N) = r.potvorticity_w(2:r.M,2:r.N) * 4. ./ ...
- (r.H(2:r.M,2:r.N) + r.H(1:r.M-1,2:r.N) + r.H(1:r.M-1,1:r.N-1) + r.H(2:r.M,1:r.N-1));
- %bottom-left corner (SW)
- r.potvorticity_w(1,1) = r.potvorticity_w(1,1) / ...
- (r.H(1,1));
- %bottom-right corner (SE)
- r.potvorticity_w(1,r.N+1) = r.potvorticity_w(1,r.N+1)/ ...
- (r.H(1,r.N));
- %top-left corner (NW)
- r.potvorticity_w(r.M+1,1) = r.potvorticity_w(r.M+1,1) / ...
- (r.H(r.M,1));
- %top-right corner (NE)
- r.potvorticity_w(r.M+1,r.N+1) = r.potvorticity_w(r.M+1,r.N+1) / ...
- (r.H(r.M,1));
- %western boundary
- r.potvorticity_w(2:r.M,1) = r.potvorticity_w(2:r.M,1) * 2. ./ ...
- (r.H(2:r.M,1) + r.H(1:r.M-1,1));
- %eastern boundary
- r.potvorticity_w(2:r.M,r.N+1) = r.potvorticity_w(2:r.M,r.N+1) * 2. ./ ...
- (r.H(1:r.M-1,r.N) + r.H(2:r.M,r.N));
- %southern boundary
- r.potvorticity_w(1,2:r.N) = r.potvorticity_w(1,2:r.N) * 2. ./ ...
- (r.H(1,2:r.N) + r.H(1,1:r.N-1));
- %northern boundary
- r.potvorticity_w(r.M+1,2:r.N) = r.potvorticity_w(r.M+1,2:r.N) * 2. ./ ...
- (r.H(r.M,2:r.N) + r.H(r.M,1:r.N-1));
- %interpolate from an F-cell to a T-cell
- r.potvorticity_t = interpolFtoT( r.potvorticity_w, r.mask, r.M, r.N);
-
- function tgrid = interpolFtoT(fgrid, r.mask, r.M, r.N)
- %% function tgrid = interpolFtoT(fgrid, r.mask, r.M, r.N)
- tgrid = .25 * r.mask .* ( ...
- fgrid(1:r.M,1:r.N) ...
- + fgrid(2:r.M+1,1:r.N) ...
- + fgrid(2:r.M+1,2:r.N+1) ...
- + fgrid(1:r.M,2:r.N+1) ...
- );
-
- %--------------------------------------------------------------------------
- %SHEL SHallow-water numerical modEL
- %Copyright (C) 2006,2009,2010,2011. Guillaume Riflet,
- %Instituto Superior T?cnico da Universidade T?cnica de Lisboa.
- %--------------------------------------------------------------------------
-