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