% 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