shelfice_remeshing/input/rdmds_init.m

%This program creates open boundary prescription files for the PIG
%experiment based on the 10 yrs spin-up run using OBCS with 
%U,V = 0, Tref, Sref

%following variables 
%-----3D fields-----
% T Temperature (C)
% S Salinity (psu)
% U u-velocity (m/s)
% PH ocean pressure (or atm geopotential)


%Set the grid size;
nx = 3;    delx = 1;   X = nx*delx;
ny = 200;    dely = 1;   Y = ny*dely;
nz = 90;    delz = 10;  Z = nz*delz;

%x axis 
x = 1:delx:X;
%y axis
y = 1:dely:Y;
%z axis [m]
z = delz:delz:Z;


%Note: 
%Ensure that the volume flux is zero at the open boundary


%% ________________________________________________________________________
%  First try: linear velocity profile. 
%  From U = -0.025 ms-1 (outwards) to  0.025 ms-1  
%  Constant in y direction (meridionally)
%  Do not take exactly Umin and Umax from the OBC of the original
%  experiment: Take modified numbers to that Umin = -Umax and volume is conserved

% Initialize variable
V_OBC(1:nx,1:nz) = NaN;


v_sfc = 0.025;                         %U at sfc
v_bottom = -0.025;                       %U at bottom
del_v = (v_bottom - v_sfc)/((nz-1)*delz);   %delu/delz
    
for iz = 1:nz;
    V_OBC(1:nx,iz) = v_sfc + del_v*((iz-1)*delz);    
end


%  Linear T profile
%  From Tmin at sfc increasing to Tmax at bottom
%  Take Tmin and Tmax from the OBC of the original PIG experiment

% Initialize variable
T_OBC(1:nx,1:nz) = NaN;


T_sfc = -1.9;
T_bottom = 2;
del_T = (T_bottom - T_sfc)/(59*delz);

for iz = 1:nz;
    
    
    tref(iz) = T_sfc + del_T*((iz-30)*delz); 
    if iz<=30;
        tref(iz)=-1.9;
    end
    if iz>=90
        tref(iz) =2;
    end
end


for iz = 1:nz;
    T_OBC(1:nx,iz) = tref(iz);    
end


%  Linear S profile
%  From Smin at sfc increasing to Smax at bottom
%  Take Smin and Smax from the OBC of the original PIG experiment

% Initialize variable
S_OBC(1:nx,1:nz) = NaN;

S_sfc = 34.2;
S_bottom = 34.7;
del_S = (S_bottom - S_sfc)/(59*delz);

for iz = 1:nz;
    
    
    sref(iz) = S_sfc + del_S*((iz-30)*delz); 
    if iz<=30;
        sref(iz)=34.2;
    end
    if iz>=90
        sref(iz) =34.7;
    end
end

for iz = 1:nz;
    S_OBC(1:nx,iz) = sref(iz);    
end


%% Initial T & S conditions

% Take western open boundary conditions for T & S 
% and assume no change in x-direction

T_init = zeros(nx,ny,nz);
S_init = zeros(nx,ny,nz);

for iy = 1:ny;
    T_init(:,iy,:) = T_OBC;
    S_init(:,iy,:) = S_OBC;
end


%% Print OBCS files for T, S

% fid=fopen('uvel.obw','w','b');  fprintf(fid,'%10.4f',U_OBC);fclose(fid);
fid=fopen('vvel.obw','w','b');  fwrite(fid,V_OBC,'real*8'); fclose(fid);
% fid=fopen('uvel.obw','w','b');  fwrite(fid,U_OBC,'real*8'); fclose(fid);
fid=fopen('theta.obw','w','b');  fwrite(fid,T_OBC,'real*8'); fclose(fid);
fid=fopen('salt.obw','w','b');  fwrite(fid,S_OBC,'real*8'); fclose(fid);
% fid=fopen('theta.obw','w');  fprintf(fid,'%10.4f',T_OBC);fclose(fid);
% fid=fopen('salt.obw','w');  fprintf(fid,'%10.4f',S_OBC);fclose(fid);


%% Print init files for T, S

fid_T=fopen('theta.init','w','b');fwrite(fid,T_init,'real*8');fclose(fid);
fid_T=fopen('salt.init','w','b');fwrite(fid,S_init,'real*8');fclose(fid);


% for iz = 1:nz;
%     fprintf(fid_T,'%10.4f',T_init(:,:,iz));
%     fprintf(fid_S,'%10.4f',S_init(:,:,iz)); 
% end
% 
% fclose(fid_T);
% fclose(fid_S);
1	dgoldberg	1.1	%This program creates open boundary prescription files for the PIG
2			%experiment based on the 10 yrs spin-up run using OBCS with
3			%U,V = 0, Tref, Sref
4
5			%following variables
6			%-----3D fields-----
7			% T Temperature (C)
8			% S Salinity (psu)
9			% U u-velocity (m/s)
10			% PH ocean pressure (or atm geopotential)
11
12
13
14
15			%Set the grid size;
16			nx = 3; delx = 1; X = nx*delx;
17			ny = 200; dely = 1; Y = ny*dely;
18			nz = 90; delz = 10; Z = nz*delz;
19
20			%x axis
21			x = 1:delx:X;
22			%y axis
23			y = 1:dely:Y;
24			%z axis [m]
25			z = delz:delz:Z;
26
27
28			%Note:
29			%Ensure that the volume flux is zero at the open boundary
30
31
32			%% ________________________________________________________________________
33			% First try: linear velocity profile.
34			% From U = -0.025 ms-1 (outwards) to 0.025 ms-1
35			% Constant in y direction (meridionally)
36			% Do not take exactly Umin and Umax from the OBC of the original
37			% experiment: Take modified numbers to that Umin = -Umax and volume is conserved
38
39			% Initialize variable
40			V_OBC(1:nx,1:nz) = NaN;
41
42
43			v_sfc = 0.025; %U at sfc
44			v_bottom = -0.025; %U at bottom
45			del_v = (v_bottom - v_sfc)/((nz-1)*delz); %delu/delz
46
47			for iz = 1:nz;
48			V_OBC(1:nx,iz) = v_sfc + del_v((iz-1)delz);
49			end
50
51
52			% Linear T profile
53			% From Tmin at sfc increasing to Tmax at bottom
54			% Take Tmin and Tmax from the OBC of the original PIG experiment
55
56			% Initialize variable
57			T_OBC(1:nx,1:nz) = NaN;
58
59
60			T_sfc = -1.9;
61			T_bottom = 2;
62			del_T = (T_bottom - T_sfc)/(59*delz);
63
64			for iz = 1:nz;
65
66
67			tref(iz) = T_sfc + del_T((iz-30)delz);
68			if iz<=30;
69			tref(iz)=-1.9;
70			end
71			if iz>=90
72			tref(iz) =2;
73			end
74			end
75
76
77
78			for iz = 1:nz;
79			T_OBC(1:nx,iz) = tref(iz);
80			end
81
82
83			% Linear S profile
84			% From Smin at sfc increasing to Smax at bottom
85			% Take Smin and Smax from the OBC of the original PIG experiment
86
87			% Initialize variable
88			S_OBC(1:nx,1:nz) = NaN;
89
90			S_sfc = 34.2;
91			S_bottom = 34.7;
92			del_S = (S_bottom - S_sfc)/(59*delz);
93
94			for iz = 1:nz;
95
96
97			sref(iz) = S_sfc + del_S((iz-30)delz);
98			if iz<=30;
99			sref(iz)=34.2;
100			end
101			if iz>=90
102			sref(iz) =34.7;
103			end
104			end
105
106			for iz = 1:nz;
107			S_OBC(1:nx,iz) = sref(iz);
108			end
109
110
111
112			%% Initial T & S conditions
113
114			% Take western open boundary conditions for T & S
115			% and assume no change in x-direction
116
117			T_init = zeros(nx,ny,nz);
118			S_init = zeros(nx,ny,nz);
119
120			for iy = 1:ny;
121			T_init(:,iy,:) = T_OBC;
122			S_init(:,iy,:) = S_OBC;
123			end
124
125
126
127			%% Print OBCS files for T, S
128
129			% fid=fopen('uvel.obw','w','b'); fprintf(fid,'%10.4f',U_OBC);fclose(fid);
130			fid=fopen('vvel.obw','w','b'); fwrite(fid,V_OBC,'real*8'); fclose(fid);
131			% fid=fopen('uvel.obw','w','b'); fwrite(fid,U_OBC,'real*8'); fclose(fid);
132			fid=fopen('theta.obw','w','b'); fwrite(fid,T_OBC,'real*8'); fclose(fid);
133			fid=fopen('salt.obw','w','b'); fwrite(fid,S_OBC,'real*8'); fclose(fid);
134			% fid=fopen('theta.obw','w'); fprintf(fid,'%10.4f',T_OBC);fclose(fid);
135			% fid=fopen('salt.obw','w'); fprintf(fid,'%10.4f',S_OBC);fclose(fid);
136
137
138
139			%% Print init files for T, S
140
141			fid_T=fopen('theta.init','w','b');fwrite(fid,T_init,'real*8');fclose(fid);
142			fid_T=fopen('salt.init','w','b');fwrite(fid,S_init,'real*8');fclose(fid);
143
144
145			% for iz = 1:nz;
146			% fprintf(fid_T,'%10.4f',T_init(:,:,iz));
147			% fprintf(fid_S,'%10.4f',S_init(:,:,iz));
148			% end
149			%
150			% fclose(fid_T);
151			% fclose(fid_S);