Dear Community,
Thank you for your kind input.
I have worked more on my bathymetry, as I realize my rx1 number was way too high. Using GRID_SmoothPositive_ROMS_rx1.m from the LP-bathymetry package with rx1max=5 i obtained a rx0-value of 0.078 which is way lower than I need. In addition I have simplified the sigmalayer distribution to Vstretch=1, Vtrans=2, thetab=thetas=1.
The max currents are now smaller but still several cm/s. When i run without salinity differences the currents are <1e.5 cm/s everywhere so it is definetely something with how the stratification is handled.
When checking my initialization file I realize that the I did not succeed in creating perfectly horizontal isolines. It may be that these errors induce the currents? Below are some images from the xi=130 crossection. From the first figure it is clear that the initial isolines are not horizontal everywhere. After 24 hours the salinity is smoothed out and apparently diffused somewhat downward in the lower layers, and tilted somewhat near the coast. The resulting pressure gradients differ in the vertical and induce the currents as shown.
https://www.myroms.org/links/t_1692/xi130_saltstart.jpg
https://www.myroms.org/links/t_1692/xi130_salt24.jpg
https://www.myroms.org/links/t_1692/xi130_u24.jpg
Below I attach the matlab used to generate the initial file.
Is there a better way to do this? What I have is a vector mysalt with salt values at set z-levels, constant salinity below 90m depth. Then I run though all the grid points, extract the z-values of each sigma layer in that point (using scoord.m from the utility package), and insert the correct salt value from the salt vector entry corresponding to the nearmost z-value. mk_vanm_init() is the main routine, and getSalt() is where the actual salt values are generated.
Eivind
Code: Select all
function [salt]=getSalt(zvec,mysalt)
% This function extracts salt values from mysalt, placing
% them at sigmalayers given by zvec.
% Returns a vector of length ll=length(zvec)
%
% zvec: Vector of length ll, negative values where z(1) is the
% deepmost value.
% mysalt: kk*2 matrix, first column gives salt values at negative z-values
% given in last column. mysalt(1,:) is the surface value
ll = length(zvec);
kk = length(mysalt(:,2));
salt = zeros(1,ll);
dz = abs(mysalt(1,2)-mysalt(2,2)); %Constant
for l=1:ll
k=0;
while (abs(zvec(l)-mysalt(kk-k,2))>=dz)
k = k+1;
end
if (k==kk || k==kk-1) %Surface
salt(l) = mysalt(1,1);
else
dzz = abs(zvec(l)-mysalt(kk-k,2))/dz;
salt(l) = dzz*mysalt(kk-k,1)+(1-dzz)*mysalt(kk-k+1,1); %Linear interpolation
end
end
Code: Select all
function mk_vanm_init(mysalt)
fname = 'vanmini_hydro.nc';
h = nc_varget('vanm_grid.nc','h');
%Grid parametres, applies to Van Mijen only
dx = 200;
dy = 200;
Lx = dx*373;
Ly = dy*248;
xi_rho = 375 ;
xi_u = 374 ;
xi_v = 375 ;
eta_rho = 250 ;
eta_u = 250 ;
eta_v = 249 ;
s_rho = 32; % =N
s_w = 33; % = N+1
%Stretching parametres (must be similar to those in ocean_*.in):
Vtransf = 2;
Vstretch = 1;
theta_s = 1.0;
theta_b = 1.0;
Tcline = 25;
% Make x,y grid
xr = [-dx/2:dx:Lx+dx/2];
yr = [-dy/2:dy:Ly+dy/2]';
[xr,yr] = meshgrid(xr,yr);
% Get z, sigma and Cs-values
column = 0;
z_rho = zeros(eta_rho,xi_rho,s_rho);
z_w = zeros(eta_rho,xi_rho,s_w);
for index = 1:xi_rho
[z_rho(:,index,:),sig_rho,Cs_r] = scoord(h, xr, yr, Vtransf, Vstretch,...
theta_s, theta_b, Tcline, s_rho, 0, column, index, 0);
[z_w(:,index,:),sig_w,Cs_w] = scoord(h, xr, yr, Vtransf, Vstretch,...
theta_s, theta_b, Tcline, s_rho, 1, column, index, 0);
end
sig_rho = sig_rho';
sig_w = sig_w';
Cs_r = Cs_r';
Cs_w = Cs_w';
%Temp data (constant = 2C everywhere):
temp = 2*ones(1,s_rho,eta_rho,xi_rho);
%Salt data:
% This setup:
% I have a vector mysalt (257*2)
% mysalt(:,1): salinity values, constant below 90 m
% mysalt(:,2): z-value, dz=-0.625m
%
% Need salt value at each s_rho, dependent on posistion.
% z-values at each sigma layer l is given by
% z_rho(i,j,l)
% For each point i,j and layer l:
% - Go through mysalt, find where the z-value deviates by less than dz
% - Extract saltvalue
disp('Filling salt matrix..');
salt = zeros (1,s_rho,eta_rho,xi_rho);
landsalt = getSalt(squeeze(z_rho(1,1,:)),mysalt); % Standard vector over land
for j=1:xi_rho
for i=1:eta_rho
if h(i,j) == 3
salt(1,:,i,j) = landsalt;
else
salt(1,:,i,j) = getSalt(squeeze(z_rho(i,j,:)),mysalt);
end
end
end
% Write to NetCDF file *************************
nc_varput ...
etc
