| 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 |
% pf*amp*cos[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 |
Parent Directory
| 
Revision Log
| 
Revision Graph