% function  create_tide_bc_second(nx, ny, cbdry_used,tide_num, time_nodal) ; 
%
% Input:  nx, the size of first dimension of MITgcm domain  
%         ny, the size of second dimension of MITgcm domain 
%     cbdry_used, cell array specify which boundary to add tidal boundary condition,  
%                 {'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 interpolation. 
%             For this case, use tide_bc_empty to generate a zero amplitude tide BC. 
%     tide_num: tidal components to be included in tidal BC file 
%     time_nodal: time for nodal correction 
% 
% Output: Tidal boundary condition files will be generated with names like
%          OB[N,S,E,W][am,ph].seaice_obcs
%
% 
% Note: 
% 
% 1) The script will read boundary grid file bdry_grid_[swne].mat generated by create_bc_grid.m 
% 2) The script will use  CAT2008b to generate tidal bc files  with names like 
%              OB[NSWE]am.seaice_obcs, for tidal current amplitude   (m/s) 
%              OB[NSWE]ph.seaice_obcs, phase in seconds  (s)  
% 
% 3) The following 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 
%
% 4) The phase generated here inludes the Greenwich Phase of tidal constituent, 
%     nodal correction term, and a time constant in order to use nodal correction formula of CAT2008b  
%
%          phase_in_second = ( tidal_period*phase_in_radian/2*pi - mitgcm_timec) 
%
% 5) CAT2008b directory needs to be specified in test_tide_bc.m 
% as following,  addpath /data17/home/xiao/CAT/Tool/TMD2.03   
%
% XC WANG 25mar2013 
%         17dec2012  
%

function  create_tide_bc_second(nx, ny, cbdry_used,tide_num, time_nodal) ;  

%The path of CATS2008b needed to be specified either here or in test_tide_bc.m 
%addpath /data17/home/xiao/CAT/Tool/TMD2.03    

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

%name of tide file of CATS2008b 
model_file =['Model_CATS2008b_km'];   

% Dimension of Boundary files  
% NX, grid points in first  dimension of MITgcm 
% NY, grid points in second dimension of MITgcm 
% NH, grid points in vertical direction 
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 
%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 


%  Use tide_num to specify what tidal components to include  
[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 = datenum(1979, 1, 1); 
time0 = time_nodal  ; 
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, :));    
 period(k) = 2*pi/omega(k) ; 
end  

% Time constant for nodal correction 
mitgcm_timec = (time_nodal-time1992)*sec_in_day    ; 

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 ;  


% 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 MITgcm  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_rad(:, itide) = vtide_ph1(:)*deg2rad ;   
% pha_all(:, itide) = vtide_ph1(:) ;   
 else 
 amp_all(:, itide) = utide_am1(:) ; 
 pha_all_rad(:, itide) = utide_ph1(:)*deg2rad ;   
% pha_all(:, itide) = utide_ph1(:) ;    
end  

 pha_all_sec(:, itide) = period(itide) .* pha_all_rad(:, itide)/(2*pi) - mitgcm_timec ;

end     % end of itide 

% Conver Nan to Zero 
badp  = ~isfinite(amp_all) ; 
amp_all(badp) = 0 ; 
badp  = ~isfinite(pha_all_sec) ; 
pha_all_sec(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_sec) ;    
%end 

end  


% mitgcm_timec = (time_nodal-time1992)*sec_in_day ; 

return