MITgcm_contrib/tides/create_tide_bc_4bdry.m

% function  mitgcm_timec = create_tide_bc_4bdry(nx, ny, cbdry_used,tide_num, time_nodal) ; 
%
% Input:  nx, the size of first dimension of MITgcm 
%         ny, the size of second dimension of MITgcm 
%     cbdry_used, cell array specify which boundary to add tidal bc 
%                 {'s', 'n', 'w', 'e'}, South, North, West, and East 
%             If there is no ocean points on one boundary, donot specify it. 
%             Otherwise it will generate an error in terpolation. 
%             For this case, use tide_bc_empty to generate a zero amplitude tide BC. 
%     time_nodal: time for nodal correction 
%     
%
% 
% Output: mitgcm_timec, notal correction time constant used in MITgcm 
%
% 
% The script is to create tidal current amplitude and phase  for MITgcm  
% 
% 1) Read boundary grid file bdry_grid_[swne].mat generated by create_bc_grid.m 
% 2) Use CAT2008b to generate tidal bc  with names like 
%              OB[NSWE]am.seaice_obcs, current amplitude 
%              OB[NSWE]ph.seaice_obcs, phase 
%
%  
% 
%         tide_num [nc < 10] ten tidal constituents can be included 
%                          m2  s2  n2  k2  k1  o1  p1  q1 mf mm  
%                          1   2   3   4   5   6   7   8  9  10 
%
% OB[N,S,E,W][am,ph]File :: Files with boundary conditions,
%                           the letter combinations mean:
%              N/S/E/W   :: northern/southern/eastern/western boundary
%              am/ph     :: tidal amplitude (m/s) / radian 
%
%
%
% XC WANG 7dec2012 
%
%addpath /data17/home/xiao/CAT/Tool/TMD2.03   

function  mitgcm_timec = create_tide_bc_4bdry(nx, ny, cbdry_used,tide_num, time_nodal) ;  

addpath /data17/home/xiao/CAT/Tool/TMD2.03    
addpath /data17/home/xiao/Mat_Tools/MITgcm_Tool

sec_in_day = 86400 ; 
c8 = 8;  c4 =4 ; 

%model_file =['Model_CATS2008a_opt'];  
model_file =['Model_CATS2008b_km'];   

% Dimension of Boundary files  
% NX, grid points in horizontal direction
% NH, grid points in vertical direction 
% Use the first records to create tidal BC 
% NX=640; NY=640; NH=70;  month_in = 1;  
NX =nx; NY =ny; 

% 
rad2deg = 180/pi ; 
deg2rad = pi/180 ;   

% Set the Boundary to add tides  
% AngleCS_bc, AngleSN_bc are the Cos and Sin and angle of U axis of MITgcm 
% from due east 
%ibc = 3 ; 
%c_bdry3 = {'s', 'n', 'w', 'e'} ;   
c_bdry3 = cbdry_used ;   

for ibc = 1:length(c_bdry3)     

c_bdry=c_bdry3{ibc} ; 
bdry_file=['bdry_grid_',c_bdry, '.mat'];  
load(bdry_file, 'lat_bc', 'lon_bc', 'depth_bc', 'AngleCS_bc', 'AngleSN_bc') ;     


% Make sure it is column vector 
[nx_t1 ny_t1]  = size(lat_bc) ; 
if(nx_t1 ~= 1) 
lat_bc = lat_bc'; 
lon_bc = lon_bc';  
depth_bc = depth_bc' ; 
AngleCS_bc = AngleCS_bc'  ; 
AngleSN_bc = AngleSN_bc'  ;  
end 


% Characters for Model Tidal BC file 
cob=upper(c_bdry) ;
fld={'am','ph'}     ; 


if any(strcmp(cob, {'N','S'})) 
 OBlength = NX ; 
 else 
 OBlength = NY ; 
end 


%  Specify what tidal components to include  
%tide_num = tide_all  ; 
%tide_num = 5 ; 

[uamp upha depth_tmd conlist0] = tmd_extract_HC(model_file, lat_bc, lon_bc, 'U',  tide_num);  
[vamp vpha depth_tmd conlist0] = tmd_extract_HC(model_file, lat_bc, lon_bc, 'V',  tide_num);      

% Make sure  uamp etc are [length(tide_num), OBlength) ]  
% When length(tide_num) == 1, the amp is (OBlength, 1), not (1, OBlength) 
[nxsize nysize] = size(uamp); 
if( nxsize ~= length(tide_num) ) 
 uamp= uamp' ; 
 vamp= vamp' ; 
 upha= upha' ; 
 vpha= vpha' ; 
end 
 

conlist=conlist0(tide_num, :) ;  

% Nodal Correction Using CAT2008b formula 
time0 = time_nodal  ; 
% time0 = datenum(1979, 1, 1); 
%time0 = datenum(1978, 12, 15); 
time1992 = datenum(1992, 1,1) ;  
[pu pf] = nodal(time0-time1992+48622, conlist);  
[pu pf] = nodal(time0-time1992+48622, conlist);   

% Get Ph parameters 
for k=1:length(tide_num) 
 [ispec(k), amp(k), ph(k), omega(k), alpha(k), cNum] = constit(conlist(k, :));   
end 


amp_all = zeros(OBlength, length(tide_num)) ; 
pha_all = zeros(OBlength, length(tide_num)) ; 

uamp_all = zeros(OBlength, length(tide_num)) ; 
upha_all = zeros(OBlength, length(tide_num)) ; 
vamp_all = zeros(OBlength, length(tide_num)) ; 
vpha_all = zeros(OBlength, length(tide_num)) ; 

for itide =1:length(tide_num)  
%for itide =1:1 

for itemp = 1:OBlength
utide_am0(itemp)  = uamp(itide, itemp) ;  
utide_ph0(itemp)  = upha(itide, itemp) ;   

vtide_am0(itemp)  = vamp(itide, itemp) ;  
vtide_ph0(itemp)  = vpha(itide, itemp) ;     
end 

% Interpolate to all grid  
utide_am0 = intp_line_safe(utide_am0, lat_bc) ; 
utide_ph0 = intp_line_safe(utide_ph0, lat_bc)  ; 
vtide_am0 = intp_line_safe(vtide_am0, lat_bc) ; 
vtide_ph0 = intp_line_safe(vtide_ph0, lat_bc) ;   

% Conduct nodal correction based  CMT2008b formulae and scripts 
%
% Tidal prediction is 
% pf*amp*cos[omega(t-time0+time0-time1992) -(pha - ph - pu)] 
utide_ph0 = utide_ph0*deg2rad-ph(itide)-pu(itide) ; 
vtide_ph0 = vtide_ph0*deg2rad-ph(itide)-pu(itide) ; 
utide_am0 = utide_am0*pf(itide) ;
vtide_am0 = vtide_am0*pf(itide) ;

% Convert to current based on Depth from GCM 
utide_ph0 = utide_ph0*rad2deg ; 
vtide_ph0 = vtide_ph0*rad2deg ;  


%[nx_t1 ny_t1] = size(utide_am0); 
%[nx_t2 ny_t2] = size(depth_bc) ;  

% make sure is [1 OBlength]; 
%if(nx_t1 ~= 1) 
% utide_am0 = utide_am0' ;  
% vtide_am0 = vtide_am0' ; 
%end  

%if(nx_t2 ~= 1) 
% depth_bc = depth_bc'; 
%end  

% Depth can have zeros 
utide_am0 = utide_am0 ./ depth_bc ; 
vtide_am0 = vtide_am0 ./ depth_bc ; 
%utide_am0(:) = utide_am0(:) ./ depth_bc(:) ; 
%vtide_am0(:) = vtide_am0(:) ./ depth_bc(:) ; 


% rotate to a new coordinate
[utide_am1 utide_ph1 vtide_am1 vtide_ph1] = tide_rot(utide_am0,utide_ph0,vtide_am0,vtide_ph0,AngleCS_bc, AngleSN_bc)  ;    

uamp_all(:, itide) = utide_am1(:) ; 
upha_all(:, itide) = utide_ph1(:) ; 
vamp_all(:, itide) = vtide_am1(:) ; 
vpha_all(:, itide) = vtide_ph1(:) ; 

if any(strcmp(cob, {'N', 'S'}))   
 amp_all(:, itide) = vtide_am1(:) ; 
 pha_all(:, itide) = vtide_ph1(:)*deg2rad ;   
% pha_all(:, itide) = vtide_ph1(:) ;   
 else 
 amp_all(:, itide) = utide_am1(:) ; 
 pha_all(:, itide) = utide_ph1(:)*deg2rad ;   
% pha_all(:, itide) = utide_ph1(:) ;   
end 

end    

% Conver Nan to Zero 
badp  = ~isfinite(amp_all) ; 
amp_all(badp) = 0 ; 
badp  = ~isfinite(pha_all) ; 
pha_all(badp) = 0 ; 


%if(1==0) 
% Write to file  
cfld = fld{1} ; 
%fnm=['OB' cob cfld '.seaice_obcs.k1']; 
fnm=['OB' cob cfld '.seaice_obcs']; 
writebin(fnm, amp_all) ;  


cfld = fld{2} ; 
%fnm=['OB' cob cfld '.seaice_obcs.k1']; 
fnm=['OB' cob cfld '.seaice_obcs']; 
writebin(fnm, pha_all) ;    
%end 

end  


mitgcm_timec = (time_nodal-time1992)*sec_in_day ; 

return 
1	dimitri	1.1	% function mitgcm_timec = create_tide_bc_4bdry(nx, ny, cbdry_used,tide_num, time_nodal) ;
2			%
3			% Input: nx, the size of first dimension of MITgcm
4			% ny, the size of second dimension of MITgcm
5			% cbdry_used, cell array specify which boundary to add tidal bc
6			% {'s', 'n', 'w', 'e'}, South, North, West, and East
7			% If there is no ocean points on one boundary, donot specify it.
8			% Otherwise it will generate an error in terpolation.
9			% For this case, use tide_bc_empty to generate a zero amplitude tide BC.
10			% time_nodal: time for nodal correction
11			%
12			%
13			%
14			% Output: mitgcm_timec, notal correction time constant used in MITgcm
15			%
16			%
17			% The script is to create tidal current amplitude and phase for MITgcm
18			%
19			% 1) Read boundary grid file bdry_grid_[swne].mat generated by create_bc_grid.m
20			% 2) Use CAT2008b to generate tidal bc with names like
21			% OB[NSWE]am.seaice_obcs, current amplitude
22			% OB[NSWE]ph.seaice_obcs, phase
23			%
24			%
25			%
26			% tide_num [nc < 10] ten tidal constituents can be included
27			% m2 s2 n2 k2 k1 o1 p1 q1 mf mm
28			% 1 2 3 4 5 6 7 8 9 10
29			%
30			% OB[N,S,E,W][am,ph]File :: Files with boundary conditions,
31			% the letter combinations mean:
32			% N/S/E/W :: northern/southern/eastern/western boundary
33			% am/ph :: tidal amplitude (m/s) / radian
34			%
35			%
36			%
37			% XC WANG 7dec2012
38			%
39			%addpath /data17/home/xiao/CAT/Tool/TMD2.03
40
41			function mitgcm_timec = create_tide_bc_4bdry(nx, ny, cbdry_used,tide_num, time_nodal) ;
42
43			addpath /data17/home/xiao/CAT/Tool/TMD2.03
44			addpath /data17/home/xiao/Mat_Tools/MITgcm_Tool
45
46			sec_in_day = 86400 ;
47			c8 = 8; c4 =4 ;
48
49			%model_file =['Model_CATS2008a_opt'];
50			model_file =['Model_CATS2008b_km'];
51
52			% Dimension of Boundary files
53			% NX, grid points in horizontal direction
54			% NH, grid points in vertical direction
55			% Use the first records to create tidal BC
56			% NX=640; NY=640; NH=70; month_in = 1;
57			NX =nx; NY =ny;
58
59			%
60			rad2deg = 180/pi ;
61			deg2rad = pi/180 ;
62
63			% Set the Boundary to add tides
64			% AngleCS_bc, AngleSN_bc are the Cos and Sin and angle of U axis of MITgcm
65			% from due east
66			%ibc = 3 ;
67			%c_bdry3 = {'s', 'n', 'w', 'e'} ;
68			c_bdry3 = cbdry_used ;
69
70			for ibc = 1:length(c_bdry3)
71
72			c_bdry=c_bdry3{ibc} ;
73			bdry_file=['bdry_grid_',c_bdry, '.mat'];
74			load(bdry_file, 'lat_bc', 'lon_bc', 'depth_bc', 'AngleCS_bc', 'AngleSN_bc') ;
75
76
77			% Make sure it is column vector
78			[nx_t1 ny_t1] = size(lat_bc) ;
79			if(nx_t1 ~= 1)
80			lat_bc = lat_bc';
81			lon_bc = lon_bc';
82			depth_bc = depth_bc' ;
83			AngleCS_bc = AngleCS_bc' ;
84			AngleSN_bc = AngleSN_bc' ;
85			end
86
87
88			% Characters for Model Tidal BC file
89			cob=upper(c_bdry) ;
90			fld={'am','ph'} ;
91
92
93			if any(strcmp(cob, {'N','S'}))
94			OBlength = NX ;
95			else
96			OBlength = NY ;
97			end
98
99
100			% Specify what tidal components to include
101			%tide_num = tide_all ;
102			%tide_num = 5 ;
103
104			[uamp upha depth_tmd conlist0] = tmd_extract_HC(model_file, lat_bc, lon_bc, 'U', tide_num);
105			[vamp vpha depth_tmd conlist0] = tmd_extract_HC(model_file, lat_bc, lon_bc, 'V', tide_num);
106
107			% Make sure uamp etc are [length(tide_num), OBlength) ]
108			% When length(tide_num) == 1, the amp is (OBlength, 1), not (1, OBlength)
109			[nxsize nysize] = size(uamp);
110			if( nxsize ~= length(tide_num) )
111			uamp= uamp' ;
112			vamp= vamp' ;
113			upha= upha' ;
114			vpha= vpha' ;
115			end
116
117
118			conlist=conlist0(tide_num, :) ;
119
120			% Nodal Correction Using CAT2008b formula
121			time0 = time_nodal ;
122			% time0 = datenum(1979, 1, 1);
123			%time0 = datenum(1978, 12, 15);
124			time1992 = datenum(1992, 1,1) ;
125			[pu pf] = nodal(time0-time1992+48622, conlist);
126			[pu pf] = nodal(time0-time1992+48622, conlist);
127
128			% Get Ph parameters
129			for k=1:length(tide_num)
130			[ispec(k), amp(k), ph(k), omega(k), alpha(k), cNum] = constit(conlist(k, :));
131			end
132
133
134			amp_all = zeros(OBlength, length(tide_num)) ;
135			pha_all = zeros(OBlength, length(tide_num)) ;
136
137			uamp_all = zeros(OBlength, length(tide_num)) ;
138			upha_all = zeros(OBlength, length(tide_num)) ;
139			vamp_all = zeros(OBlength, length(tide_num)) ;
140			vpha_all = zeros(OBlength, length(tide_num)) ;
141
142			for itide =1:length(tide_num)
143			%for itide =1:1
144
145			for itemp = 1:OBlength
146			utide_am0(itemp) = uamp(itide, itemp) ;
147			utide_ph0(itemp) = upha(itide, itemp) ;
148
149			vtide_am0(itemp) = vamp(itide, itemp) ;
150			vtide_ph0(itemp) = vpha(itide, itemp) ;
151			end
152
153			% Interpolate to all grid
154			utide_am0 = intp_line_safe(utide_am0, lat_bc) ;
155			utide_ph0 = intp_line_safe(utide_ph0, lat_bc) ;
156			vtide_am0 = intp_line_safe(vtide_am0, lat_bc) ;
157			vtide_ph0 = intp_line_safe(vtide_ph0, lat_bc) ;
158
159			% Conduct nodal correction based CMT2008b formulae and scripts
160			%
161			% Tidal prediction is
162			% pfampcos[omega(t-time0+time0-time1992) -(pha - ph - pu)]
163			utide_ph0 = utide_ph0*deg2rad-ph(itide)-pu(itide) ;
164			vtide_ph0 = vtide_ph0*deg2rad-ph(itide)-pu(itide) ;
165			utide_am0 = utide_am0*pf(itide) ;
166			vtide_am0 = vtide_am0*pf(itide) ;
167
168			% Convert to current based on Depth from GCM
169			utide_ph0 = utide_ph0*rad2deg ;
170			vtide_ph0 = vtide_ph0*rad2deg ;
171
172
173			%[nx_t1 ny_t1] = size(utide_am0);
174			%[nx_t2 ny_t2] = size(depth_bc) ;
175
176			% make sure is [1 OBlength];
177			%if(nx_t1 ~= 1)
178			% utide_am0 = utide_am0' ;
179			% vtide_am0 = vtide_am0' ;
180			%end
181
182			%if(nx_t2 ~= 1)
183			% depth_bc = depth_bc';
184			%end
185
186			% Depth can have zeros
187			utide_am0 = utide_am0 ./ depth_bc ;
188			vtide_am0 = vtide_am0 ./ depth_bc ;
189			%utide_am0(:) = utide_am0(:) ./ depth_bc(:) ;
190			%vtide_am0(:) = vtide_am0(:) ./ depth_bc(:) ;
191
192
193			% rotate to a new coordinate
194			[utide_am1 utide_ph1 vtide_am1 vtide_ph1] = tide_rot(utide_am0,utide_ph0,vtide_am0,vtide_ph0,AngleCS_bc, AngleSN_bc) ;
195
196			uamp_all(:, itide) = utide_am1(:) ;
197			upha_all(:, itide) = utide_ph1(:) ;
198			vamp_all(:, itide) = vtide_am1(:) ;
199			vpha_all(:, itide) = vtide_ph1(:) ;
200
201			if any(strcmp(cob, {'N', 'S'}))
202			amp_all(:, itide) = vtide_am1(:) ;
203			pha_all(:, itide) = vtide_ph1(:)*deg2rad ;
204			% pha_all(:, itide) = vtide_ph1(:) ;
205			else
206			amp_all(:, itide) = utide_am1(:) ;
207			pha_all(:, itide) = utide_ph1(:)*deg2rad ;
208			% pha_all(:, itide) = utide_ph1(:) ;
209			end
210
211			end
212
213			% Conver Nan to Zero
214			badp = ~isfinite(amp_all) ;
215			amp_all(badp) = 0 ;
216			badp = ~isfinite(pha_all) ;
217			pha_all(badp) = 0 ;
218
219
220			%if(1==0)
221			% Write to file
222			cfld = fld{1} ;
223			%fnm=['OB' cob cfld '.seaice_obcs.k1'];
224			fnm=['OB' cob cfld '.seaice_obcs'];
225			writebin(fnm, amp_all) ;
226
227
228			cfld = fld{2} ;
229			%fnm=['OB' cob cfld '.seaice_obcs.k1'];
230			fnm=['OB' cob cfld '.seaice_obcs'];
231			writebin(fnm, pha_all) ;
232			%end
233
234			end
235
236
237			mitgcm_timec = (time_nodal-time1992)*sec_in_day ;
238
239			return